atrstm_radcoefs.c (10666B)
1 /* Copyright (C) 2022, 2023 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2020, 2021 Centre National de la Recherche Scientifique 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #include "atrstm_c.h" 18 #include "atrstm_log.h" 19 #include "atrstm_radcoefs.h" 20 21 #ifdef ATRSTM_USE_SIMD 22 #include "atrstm_radcoefs_simd4.h" 23 #endif 24 25 #include <astoria/atrri.h> 26 #include <astoria/atrtp.h> 27 28 #include <star/suvm.h> 29 30 #include <rsys/cstr.h> 31 32 #include <string.h> 33 34 /******************************************************************************* 35 * Exported functions 36 ******************************************************************************/ 37 res_T 38 atrstm_fetch_radcoefs 39 (const struct atrstm* atrstm, 40 const struct atrstm_fetch_radcoefs_args* args, 41 atrstm_radcoefs_T radcoefs) 42 { 43 res_T res = RES_OK; 44 45 #ifndef ATRSTM_USE_SIMD 46 { 47 ASSERT(atrstm->use_simd == 0); 48 res = fetch_radcoefs(atrstm, args, radcoefs); 49 if(res != RES_OK) goto error; 50 } 51 #else 52 { 53 if(!atrstm->use_simd) { 54 res = fetch_radcoefs(atrstm, args, radcoefs); 55 if(res != RES_OK) goto error; 56 } else { 57 res = fetch_radcoefs_simd4(atrstm, args, radcoefs); 58 if(res != RES_OK) goto error; 59 } 60 } 61 #endif /* ATRSTM_USE_SIMD */ 62 63 exit: 64 return res; 65 error: 66 goto exit; 67 } 68 69 /******************************************************************************* 70 * Local functions functions 71 ******************************************************************************/ 72 int 73 check_fetch_radcoefs_args 74 (const struct atrstm* atrstm, 75 const struct atrstm_fetch_radcoefs_args* args, 76 const char* func_name) 77 { 78 int i; 79 ASSERT(atrstm && args && func_name); 80 81 if(!args 82 || SUVM_PRIMITIVE_NONE(&args->prim) 83 || args->wavelength < atrstm->wlen_range[0] 84 || args->wavelength > atrstm->wlen_range[1] 85 || !(args->radcoefs_mask & ATRSTM_RADCOEFS_MASK_ALL) 86 || !(args->components_mask & ATRSTM_CPNTS_MASK_ALL)) 87 return 0; 88 89 if(atrstm->spectral_type != ATRSTM_SPECTRAL_SW) { 90 log_err(atrstm, "%s: only shortwave is currently supported.\n", 91 func_name); 92 return 0; 93 } 94 95 FOR_EACH(i, 0, ATRSTM_RADCOEFS_COUNT__) { 96 if(args->radcoefs_mask & BIT(i)) { 97 if(args->k_min[i] > args->k_max[i]) return 0; 98 } 99 } 100 101 return 1; 102 } 103 104 res_T 105 fetch_radcoefs 106 (const struct atrstm* atrstm, 107 const struct atrstm_fetch_radcoefs_args* args, 108 atrstm_radcoefs_T radcoefs) 109 { 110 struct radcoefs prim_radcoefs[4]; 111 res_T res = RES_OK; 112 113 if(!atrstm 114 || !check_fetch_radcoefs_args(atrstm, args, FUNC_NAME) 115 || !radcoefs) { 116 res = RES_BAD_ARG; 117 goto error; 118 } 119 memset(radcoefs, 0, sizeof(double[ATRSTM_RADCOEFS_COUNT__])); 120 121 /* Compute the radiative properties of the primitive nodes */ 122 res = primitive_compute_radcoefs(atrstm, &atrstm->refract_id, &args->prim, 123 args->radcoefs_mask, prim_radcoefs); 124 if(res != RES_OK) { 125 log_err(atrstm, 126 "Error computing the radiative properties for the primitive %lu -- %s.\n", 127 (unsigned long)args->prim.iprim, 128 res_to_cstr(res)); 129 goto error; 130 } 131 132 /* Linearly interpolate the node radiative properties */ 133 if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) { 134 const double prim_ka_min = MMIN 135 (MMIN(prim_radcoefs[0].ka, prim_radcoefs[1].ka), 136 MMIN(prim_radcoefs[2].ka, prim_radcoefs[3].ka)); 137 const double prim_ka_max = MMAX 138 (MMAX(prim_radcoefs[0].ka, prim_radcoefs[1].ka), 139 MMAX(prim_radcoefs[2].ka, prim_radcoefs[3].ka)); 140 radcoefs[ATRSTM_RADCOEF_Ka] = 141 prim_radcoefs[0].ka * args->bcoords[0] 142 + prim_radcoefs[1].ka * args->bcoords[1] 143 + prim_radcoefs[2].ka * args->bcoords[2] 144 + prim_radcoefs[3].ka * args->bcoords[3]; 145 radcoefs[ATRSTM_RADCOEF_Ka] = /* Handle numerical uncertainty */ 146 CLAMP(radcoefs[ATRSTM_RADCOEF_Ka], prim_ka_min, prim_ka_max); 147 ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] >= args->k_min[ATRSTM_RADCOEF_Ka]); 148 ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] <= args->k_max[ATRSTM_RADCOEF_Ka]); 149 } 150 if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) { 151 const double prim_ks_min = MMIN 152 (MMIN(prim_radcoefs[0].ks, prim_radcoefs[1].ks), 153 MMIN(prim_radcoefs[2].ks, prim_radcoefs[3].ks)); 154 const double prim_ks_max = MMAX 155 (MMAX(prim_radcoefs[0].ks, prim_radcoefs[1].ks), 156 MMAX(prim_radcoefs[2].ks, prim_radcoefs[3].ks)); 157 radcoefs[ATRSTM_RADCOEF_Ks] = 158 prim_radcoefs[0].ks * args->bcoords[0] 159 + prim_radcoefs[1].ks * args->bcoords[1] 160 + prim_radcoefs[2].ks * args->bcoords[2] 161 + prim_radcoefs[3].ks * args->bcoords[3]; 162 radcoefs[ATRSTM_RADCOEF_Ks] = /* Handle numerical uncertainty */ 163 CLAMP(radcoefs[ATRSTM_RADCOEF_Ks], prim_ks_min, prim_ks_max); 164 ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] >= args->k_min[ATRSTM_RADCOEF_Ks]); 165 ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] <= args->k_max[ATRSTM_RADCOEF_Ks]); 166 } 167 if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Kext) { 168 const double prim_kext_min = MMIN 169 (MMIN(prim_radcoefs[0].kext, prim_radcoefs[1].kext), 170 MMIN(prim_radcoefs[2].kext, prim_radcoefs[3].kext)); 171 const double prim_kext_max = MMAX 172 (MMAX(prim_radcoefs[0].kext, prim_radcoefs[1].kext), 173 MMAX(prim_radcoefs[2].kext, prim_radcoefs[3].kext)); 174 radcoefs[ATRSTM_RADCOEF_Kext] = 175 prim_radcoefs[0].kext * args->bcoords[0] 176 + prim_radcoefs[1].kext * args->bcoords[1] 177 + prim_radcoefs[2].kext * args->bcoords[2] 178 + prim_radcoefs[3].kext * args->bcoords[3]; 179 radcoefs[ATRSTM_RADCOEF_Kext] = /* Handle numerical uncertainty */ 180 CLAMP(radcoefs[ATRSTM_RADCOEF_Kext], prim_kext_min, prim_kext_max); 181 ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] >= args->k_min[ATRSTM_RADCOEF_Kext]); 182 ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] <= args->k_max[ATRSTM_RADCOEF_Kext]); 183 } 184 185 exit: 186 return res; 187 error: 188 goto exit; 189 } 190 191 res_T 192 primitive_compute_radcoefs 193 (const struct atrstm* atrstm, 194 const struct atrri_refractive_index* refract_id, 195 const struct suvm_primitive* prim, 196 const int radcoefs_mask, /* Combination of atrstm_radcoef_flag */ 197 struct radcoefs radcoefs[]) 198 { 199 struct radcoefs_compute_args args = RADCOEFS_COMPUTE_ARGS_NULL; 200 struct atrtp_desc desc = ATRTP_DESC_NULL; 201 size_t inode; 202 res_T res = RES_OK; 203 ASSERT(atrstm && prim && prim->nvertices == 4 && refract_id && radcoefs); 204 205 res = atrtp_get_desc(atrstm->atrtp, &desc); 206 if(res != RES_OK) goto error; 207 208 /* Setup the constant input arguments of the optical properties computation 209 * routine */ 210 args.lambda = refract_id->wavelength; 211 args.n = refract_id->n; 212 args.kappa = refract_id->kappa; 213 args.fractal_prefactor = atrstm->fractal_prefactor; 214 args.fractal_dimension = atrstm->fractal_dimension; 215 args.radcoefs_mask = radcoefs_mask; 216 217 FOR_EACH(inode, 0, 4) { 218 const double* node; 219 220 /* Fetch the thermodynamic properties of the node */ 221 node = atrtp_desc_get_node_properties(&desc, prim->indices[inode]); 222 223 /* Setup the per node input arguments of the optical properties computation 224 * routine */ 225 args.soot_volumic_fraction = 226 node[ATRTP_SOOT_VOLFRAC]; 227 args.soot_primary_particles_count = 228 node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT]; 229 args.soot_primary_particles_diameter = 230 node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER]; 231 232 /* Compute the node optical properties */ 233 radcoefs_compute(&radcoefs[inode], &args); 234 } 235 236 exit: 237 return res; 238 error: 239 goto exit; 240 } 241 242 res_T 243 dump_radcoefs 244 (const struct atrstm* atrstm, 245 const struct atrri_refractive_index* refract_id, 246 FILE* stream) 247 { 248 struct radcoefs_compute_args args = RADCOEFS_COMPUTE_ARGS_NULL; 249 struct atrtp_desc desc = ATRTP_DESC_NULL; 250 struct radcoefs props_min; 251 struct radcoefs props_max; 252 size_t inode; 253 254 res_T res = RES_OK; 255 ASSERT(atrstm && stream); 256 257 res = atrtp_get_desc(atrstm->atrtp, &desc); 258 if(res != RES_OK) goto error; 259 260 /* Setup the constant input params of the optical properties computation */ 261 args.lambda = refract_id->wavelength; 262 args.n = refract_id->n; 263 args.kappa = refract_id->kappa; 264 args.fractal_prefactor = atrstm->fractal_prefactor; 265 args.fractal_dimension = atrstm->fractal_dimension; 266 args.radcoefs_mask = ATRSTM_RADCOEFS_MASK_ALL; 267 268 props_min.ka = DBL_MAX; 269 props_max.ka =-DBL_MAX; 270 props_min.ks = DBL_MAX; 271 props_max.ks =-DBL_MAX; 272 props_min.kext = DBL_MAX; 273 props_max.kext =-DBL_MAX; 274 275 #define FPRINTF(Text, Args) { \ 276 const int n = fprintf(stream, Text COMMA_##Args LIST_##Args); \ 277 if(n < 0) { \ 278 log_err(atrstm, "Error writing optical properties.\n"); \ 279 res = RES_IO_ERR; \ 280 goto error; \ 281 } \ 282 } (void)0 283 284 FPRINTF("# Ka Ks Kext\n", ARG0()); 285 286 FOR_EACH(inode, 0, desc.nnodes) { 287 struct radcoefs props; 288 const double* node; 289 290 /* Fetch the thermodynamic properties of the node */ 291 node = atrtp_desc_get_node_properties(&desc, inode); 292 293 /* Setup the per node input args of the optical properties computation */ 294 args.soot_volumic_fraction = 295 node[ATRTP_SOOT_VOLFRAC]; 296 args.soot_primary_particles_count = 297 node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT]; 298 args.soot_primary_particles_diameter = 299 node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER]; 300 301 /* Compute the node optical properties */ 302 radcoefs_compute(&props, &args); 303 304 FPRINTF("%g %g %g\n", ARG3(props.ka, props.ks, props.kext)); 305 306 props_min.ka = MMIN(props_min.ka, props.ka); 307 props_max.ka = MMAX(props_max.ka, props.ka); 308 props_min.ks = MMIN(props_min.ks, props.ks); 309 props_max.ks = MMAX(props_max.ks, props.ks); 310 props_min.kext = MMIN(props_min.kext, props.kext); 311 props_max.kext = MMAX(props_max.kext, props.kext); 312 313 } 314 315 FPRINTF("# Bounds = [%g %g %g], [%g, %g, %g]\n", ARG6 316 (props_min.ka, props_min.ks, props_min.kext, 317 props_max.ka, props_max.ks, props_max.kext)); 318 319 #undef FPRINTF 320 321 exit: 322 return res; 323 error: 324 goto exit; 325 }