sdis_heat_path_boundary_Xd_handle_external_net_flux.h (14170B)
1 /* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #include "sdis_brdf.h" 17 #include "sdis_heat_path_boundary_c.h" 18 #include "sdis_interface_c.h" 19 #include "sdis_log.h" 20 #include "sdis_scene_c.h" 21 #include "sdis_source_c.h" 22 23 #include <rsys/cstr.h> /* res_to_cstr */ 24 25 #include "sdis_Xd_begin.h" 26 27 /******************************************************************************* 28 * Non generic data types and function 29 ******************************************************************************/ 30 #ifndef SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H 31 #define SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H 32 33 /* Incident diffuse flux is made up of two components. One corresponds to the 34 * diffuse flux due to the reflection of the source on surfaces. The other is 35 * the diffuse flux due to the source's radiation scattering at least once in 36 * the environment. */ 37 struct incident_diffuse_flux { 38 double reflected; /* [W/m^2] */ 39 double scattered; /* [W/m^2] */ 40 double dir[3]; /* Direction along wich the scattered part was retrieved */ 41 }; 42 #define INCIDENT_DIFFUSE_FLUX_NULL__ {0, 0, {0,0,0}} 43 static const struct incident_diffuse_flux INCIDENT_DIFFUSE_FLUX_NULL = 44 INCIDENT_DIFFUSE_FLUX_NULL__; 45 46 #endif /* SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H */ 47 48 /******************************************************************************* 49 * Generic helper functions 50 ******************************************************************************/ 51 static INLINE res_T 52 XD(check_handle_external_net_flux_args) 53 (const struct sdis_scene* scn, 54 const char* func_name, 55 const struct handle_external_net_flux_args* args) 56 { 57 int net_flux = 0; 58 res_T res = RES_OK; 59 60 /* Handle bugs */ 61 ASSERT(scn && func_name && args); 62 ASSERT(args->interf && args->frag); 63 ASSERT(!SXD_HIT_NONE(args->XD(hit))); 64 ASSERT(args->h_cond >= 0 && args->h_conv >= 0 && args->h_radi >= 0); 65 ASSERT(args->h_cond + args->h_conv + args->h_radi > 0); 66 67 net_flux = interface_side_is_external_flux_handled(args->interf, args->frag); 68 net_flux = net_flux && (scn->source != NULL); 69 70 if(net_flux && args->picard_order != 1) { 71 log_err(scn->dev, 72 "%s: Impossible to process external fluxes when Picard order is not " 73 "equal to 1; Picard order is currently set to %lu.\n", 74 func_name, (unsigned long)args->picard_order); 75 res = RES_BAD_ARG; 76 return res; 77 } 78 79 if(sdis_medium_get_type(args->interf->medium_back) 80 == sdis_medium_get_type(args->interf->medium_front)) { 81 log_err(scn->dev, 82 "%s: external fluxes can only be processed on fluid/solid interfaces.\n", 83 func_name); 84 res = RES_BAD_ARG; 85 return res; 86 } 87 88 return RES_OK; 89 } 90 91 static INLINE double /* [W/m^2/sr] */ 92 XD(direct_contribution) 93 (struct sdis_scene* scn, 94 struct source_sample* sample, 95 const double pos[DIM], 96 const unsigned enc_id, /* Current enclosure */ 97 const struct sXd(hit)* hit_from) 98 { 99 struct sXd(hit) hit = SXD_HIT_NULL; 100 ASSERT(scn && sample && pos && hit_from); 101 102 /* Is the source hidden */ 103 XD(trace_ray)(scn, pos, sample->dir, sample->dst, enc_id, hit_from, &hit); 104 if(!SXD_HIT_NONE(&hit)) return 0; /* [W/m^2/sr] */ 105 106 /* Note that the value returned is not the source's actual radiance, but the 107 * radiance relative to the source's power. Care must therefore be taken to 108 * multiply it by the power of the source to obtain its real contribution. 109 * This trick makes it possible to manage the external flux in the green 110 * function. */ 111 return sample->radiance_term; /* [W/m^2/sr] */ 112 } 113 114 static res_T 115 XD(compute_incident_diffuse_flux) 116 (struct sdis_scene* scn, 117 struct ssp_rng* rng, 118 const struct source_props* props, 119 const double in_pos[DIM], /* position */ 120 const double in_N[DIM], /* Surface normal. (Away from the surface) */ 121 const double time, 122 const unsigned enc_id, /* Current enclosure */ 123 const struct sXd(hit)* in_hit, /* Current intersection */ 124 struct incident_diffuse_flux* diffuse_flux) /* [W/m^2] */ 125 { 126 struct sXd(hit) hit = SXD_HIT_NULL; 127 double pos[3] = {0}; /* In 3D for ray tracing ray to the source */ 128 double dir[3] = {0}; /* Incident direction (toward the surface). Always 3D.*/ 129 double N[3] = {0}; /* Surface normal. Always 3D */ 130 size_t nbounces = 0; /* For debug */ 131 res_T res = RES_OK; 132 ASSERT(props && in_pos && in_N && in_hit && diffuse_flux); 133 134 /* Local copy of input argument */ 135 dX(set)(pos, in_pos); 136 dX(set)(N, in_N); 137 hit = *in_hit; 138 139 /* Sample a diffusive direction in 3D */ 140 ssp_ran_hemisphere_cos(rng, N, dir, NULL); 141 142 *diffuse_flux = INCIDENT_DIFFUSE_FLUX_NULL; 143 144 for(;;) { 145 /* External sources */ 146 struct source_sample samp = SOURCE_SAMPLE_NULL; 147 148 /* Interface */ 149 struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; 150 struct sdis_interface* interf = NULL; 151 152 /* BRDF */ 153 struct brdf brdf = BRDF_NULL; 154 struct brdf_sample bounce = BRDF_SAMPLE_NULL; 155 struct brdf_setup_args brdf_setup_args = BRDF_SETUP_ARGS_NULL; 156 157 /* Miscellaneous */ 158 double L = 0; /* incident flux to bounce position */ 159 double wi[3] = {0}; /* Incident direction (outward the surface). Always 3D */ 160 161 d3_minus(wi, dir); /* Always in 3D */ 162 163 res = XD(find_next_fragment) 164 (scn, pos, dir, &hit, time, enc_id, &hit, &interf, &frag); 165 if(res != RES_OK) goto error; 166 167 if(SXD_HIT_NONE(&hit)) { 168 /* No surface. Handle the radiance emitted by the source and scattered at 169 * least once in the environment. Note that the value returned is not the 170 * actual scattered component of the incident diffuse flux: it relates 171 * to the radiance of the source scattered along the input dir at the 172 * given instant. It must therefore be multiplied by this radiance to 173 * obtain its real contribution. This trick makes it possible to manage 174 * the external flux in the green function. */ 175 diffuse_flux->scattered = PI; 176 diffuse_flux->dir[0] = dir[0]; 177 diffuse_flux->dir[1] = dir[1]; 178 diffuse_flux->dir[2] = dir[2]; 179 break; 180 } 181 182 d3_set(pos, frag.P); 183 184 /* Retrieve BRDF at current interface position */ 185 brdf_setup_args.interf = interf; 186 brdf_setup_args.frag = &frag; 187 brdf_setup_args.source_id = sdis_source_get_id(scn->source); 188 res = brdf_setup(scn->dev, &brdf_setup_args, &brdf); 189 if(res != RES_OK) goto error; 190 191 /* Check if path is absorbed */ 192 if(ssp_rng_canonical(rng) < brdf.emissivity) break; 193 194 /* Sample rebound direction */ 195 switch(frag.side) { 196 case SDIS_FRONT: dX(set)(N, frag.Ng); break; 197 case SDIS_BACK: dX(minus)(N, frag.Ng); break; 198 default: FATAL("Unreachable code\n"); 199 } 200 brdf_sample(&brdf, rng, wi, N, &bounce); 201 d3_set(dir, bounce.dir); /* Always in 3D */ 202 203 /* Calculate the direct contribution if the rebound is specular */ 204 if(bounce.cpnt == BRDF_SPECULAR) { 205 res = source_trace_to(scn->source, props, pos, bounce.dir, &samp); 206 if(res != RES_OK) goto error; 207 208 if(!SOURCE_SAMPLE_NONE(&samp)) { 209 double Ld = XD(direct_contribution)(scn, &samp, pos, enc_id, &hit); 210 L = Ld; /* [W/m^2/sr] */ 211 } 212 213 /* Calculate the direct contribution of the rebound is diffuse */ 214 } else { 215 double cos_theta = 0; 216 ASSERT(bounce.cpnt == BRDF_DIFFUSE); 217 218 /* Sample an external source to handle its direct contribution at the 219 * bounce position */ 220 res = source_sample(scn->source, props, rng, pos, &samp); 221 CHK(res == RES_OK); 222 cos_theta = d3_dot(samp.dir, N); 223 224 /* The source is behind the surface */ 225 if(cos_theta <= 0) { 226 L = 0; /* [W/m^2/sr] */ 227 228 /* The source is above the surface */ 229 } else { 230 double Ld = XD(direct_contribution)(scn, &samp, pos, enc_id, &hit); 231 L = Ld * cos_theta / (PI * samp.pdf); /* [W/m^2/sr] */ 232 } 233 } 234 diffuse_flux->reflected += L; /* [W/m^2/sr] */ 235 ++nbounces; 236 } 237 diffuse_flux->reflected *= PI; /* [W/m^2] */ 238 239 exit: 240 return res; 241 error: 242 goto exit; 243 } 244 245 /******************************************************************************* 246 * Local functions 247 ******************************************************************************/ 248 res_T 249 XD(handle_external_net_flux) 250 (struct sdis_scene* scn, 251 struct ssp_rng* rng, 252 const struct handle_external_net_flux_args* args, 253 struct temperature* T) 254 { 255 /* Terms to be registered in the green function */ 256 struct sdis_green_external_flux_terms green = 257 SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL; 258 259 /* External source */ 260 struct source_props src_props = SOURCE_PROPS_NULL; 261 struct source_sample src_sample = SOURCE_SAMPLE_NULL; 262 263 /* External flux */ 264 struct incident_diffuse_flux incident_flux_diffuse = INCIDENT_DIFFUSE_FLUX_NULL; 265 double incident_flux = 0; /* [W/m^2] */ 266 double incident_flux_direct = 0; /* [W/m^2] */ 267 double net_flux = 0; /* [W/m^2] */ 268 double net_flux_sc = 0; /* [W/m^2] */ 269 270 /* Sampled path */ 271 double N[3] = {0}; /* Normal. Always in 3D */ 272 273 /* On the fluid side */ 274 struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; 275 276 /* Miscellaneous */ 277 unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; 278 double sum_h = 0; 279 double emissivity = 0; /* Emissivity */ 280 double Ld = 0; /* Incident radiance [W/m^2/sr] */ 281 double cos_theta = 0; 282 unsigned src_id = 0; 283 int handle_flux = 0; 284 res_T res = RES_OK; 285 ASSERT(scn && args && T); 286 287 res = XD(check_handle_external_net_flux_args)(scn, FUNC_NAME, args); 288 if(res != RES_OK) goto error; 289 290 /* Setup the interface fragment on fluid side */ 291 frag = *args->frag; 292 if(sdis_medium_get_type(args->interf->medium_front) == SDIS_FLUID) { 293 frag.side = SDIS_FRONT; 294 } else { 295 ASSERT(sdis_medium_get_type(args->interf->medium_back) == SDIS_FLUID); 296 frag.side = SDIS_BACK; 297 } 298 299 /* Retrieve the enclosures */ 300 scene_get_enclosure_ids(scn, args->XD(hit)->prim.prim_id, enc_ids); 301 302 /* No external sources <=> no external fluxes. Nothing to do */ 303 handle_flux = interface_side_is_external_flux_handled(args->interf, &frag); 304 handle_flux = handle_flux && (scn->source != NULL); 305 if(!handle_flux) goto exit; 306 307 /* Emissivity is null <=> external flux is null. Nothing to do */ 308 src_id = sdis_source_get_id(scn->source); 309 emissivity = interface_side_get_emissivity(args->interf, src_id, &frag); 310 res = interface_side_check_emissivity(scn->dev, emissivity, frag.P, frag.time); 311 if(res != RES_OK) goto error; 312 313 if(emissivity == 0) goto exit; 314 315 res = source_get_props(scn->source, frag.time, &src_props); 316 if(res != RES_OK) goto error; 317 318 /* Sample a direction toward the source to add its direct contribution */ 319 res = source_sample(scn->source, &src_props, rng, frag.P, &src_sample); 320 if(res != RES_OK) goto error; 321 322 /* Setup the normal to ensure that it points toward the fluid medium */ 323 dX(set)(N, frag.Ng); 324 if(frag.side == SDIS_BACK) dX(minus)(N, N); 325 326 /* Calculate the incident direct flux if the external source is above the 327 * interface side */ 328 cos_theta = d3_dot(N, src_sample.dir); 329 if(cos_theta > 0) { 330 Ld = XD(direct_contribution) 331 (scn, &src_sample, frag.P, enc_ids[frag.side], args->XD(hit)); 332 incident_flux_direct = cos_theta * Ld / src_sample.pdf; /* [W/m^2] */ 333 } 334 335 /* Calculate the incident diffuse flux [W/m^2] */ 336 res = XD(compute_incident_diffuse_flux)(scn, rng, &src_props, frag.P, N, 337 frag.time, enc_ids[frag.side], args->XD(hit), &incident_flux_diffuse); 338 if(res != RES_OK) goto error; 339 340 /* Calculate the incident flux without the part scattered by the environment. 341 * The latter depends on the source's diffuse radiance, not on its power. On 342 * the other hand, both the direct incident flux and the diffuse incident flux 343 * reflected by surfaces are linear with respect to the source power. This 344 * term can therefore be recorded in the green function in relation to this 345 * power, whereas the incident diffused flux coming from the scattered source 346 * radiance depends on the diffuse radiance of the source */ 347 incident_flux = /* [W/m^2] */ 348 incident_flux_direct + incident_flux_diffuse.reflected; 349 350 /* Calculate the net flux [W/m^2] */ 351 net_flux = incident_flux * emissivity; /* [W/m^2] */ 352 353 /* Calculate the net flux from the radiance source scattered at least once by 354 * the medium */ 355 net_flux_sc = incident_flux_diffuse.scattered * emissivity; /* [W/m^2] */ 356 357 /* Until now, net flux has been calculated on the basis of source power and 358 * source diffuse radiance. What is actually calculated are the external flux 359 * terms of the green function. These must be multiplied by the source power 360 * and the source diffuse radiance, then added together to give the actual 361 * external flux */ 362 sum_h = (args->h_radi + args->h_conv + args->h_cond); 363 green.term_wrt_power = net_flux / sum_h; /* [K/W] */ 364 green.term_wrt_diffuse_radiance = net_flux_sc / sum_h; /* [K/W/m^2/sr] */ 365 green.time = frag.time; /* [s] */ 366 green.dir[0] = incident_flux_diffuse.dir[0]; 367 green.dir[1] = incident_flux_diffuse.dir[1]; 368 green.dir[2] = incident_flux_diffuse.dir[2]; 369 370 T->value += green.term_wrt_power * src_props.power; 371 if(green.term_wrt_diffuse_radiance) { 372 T->value += 373 green.term_wrt_diffuse_radiance 374 * source_get_diffuse_radiance(scn->source, green.time, green.dir); 375 } 376 377 /* Register the external net flux terms */ 378 if(args->green_path) { 379 res = green_path_add_external_flux_terms(args->green_path, &green); 380 if(res != RES_OK) goto error; 381 } 382 383 exit: 384 return res; 385 error: 386 goto exit; 387 } 388 389 #include "sdis_Xd_end.h"