sdis_heat_path_conductive_delta_sphere_Xd.h (16415B)
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_log.h" 17 #include "sdis_green.h" 18 #include "sdis_interface_c.h" 19 #include "sdis_medium_c.h" 20 #include "sdis_misc.h" 21 #include "sdis_scene_c.h" 22 23 #include "sdis_Xd_begin.h" 24 25 /******************************************************************************* 26 * Helper functions 27 ******************************************************************************/ 28 /* Sample the next direction to walk toward and compute the distance to travel. 29 * Return the sampled direction `dir0', the distance to travel along this 30 * direction, the hit `hit0' along `dir0' wrt to the returned distance, the 31 * direction `dir1' used to adjust the displacement distance, and the hit 32 * `hit1' along `dir1' used to adjust the displacement distance. */ 33 static float 34 XD(sample_next_step) 35 (struct sdis_scene* scn, 36 struct ssp_rng* rng, 37 const float pos[DIM], 38 const float delta_solid, 39 float dir0[DIM], /* Sampled direction */ 40 float dir1[DIM], /* Direction used to adjust delta */ 41 struct sXd(hit)* hit0, /* Hit along the sampled direction */ 42 struct sXd(hit)* hit1) /* Hit used to adjust delta */ 43 { 44 struct sXd(hit) hits[2]; 45 float dirs[2][DIM]; 46 float range[2]; 47 float delta; 48 ASSERT(scn && rng && pos && delta_solid>0 && dir0 && dir1 && hit0 && hit1); 49 50 *hit0 = SXD_HIT_NULL; 51 *hit1 = SXD_HIT_NULL; 52 53 #if DIM == 2 54 /* Sample a main direction around 2PI */ 55 ssp_ran_circle_uniform_float(rng, dirs[0], NULL); 56 #else 57 /* Sample a main direction around 4PI */ 58 ssp_ran_sphere_uniform_float(rng, dirs[0], NULL); 59 #endif 60 61 /* Negate the sampled dir */ 62 fX(minus)(dirs[1], dirs[0]); 63 64 /* Use the previously sampled direction to estimate the minimum distance from 65 * `pos' to the scene boundary */ 66 f2(range, FLT_MIN, delta_solid*RAY_RANGE_MAX_SCALE); 67 SXD(scene_view_trace_ray(scn->sXd(view), pos, dirs[0], range, NULL, &hits[0])); 68 SXD(scene_view_trace_ray(scn->sXd(view), pos, dirs[1], range, NULL, &hits[1])); 69 if(SXD_HIT_NONE(&hits[0]) && SXD_HIT_NONE(&hits[1])) { 70 delta = delta_solid; 71 } else { 72 delta = MMIN(hits[0].distance, hits[1].distance); 73 } 74 75 if(!SXD_HIT_NONE(&hits[0]) 76 && delta != hits[0].distance 77 && eq_eps(hits[0].distance, delta, delta_solid*0.1)) { 78 /* Set delta to the main hit distance if it is roughly equal to it in order 79 * to avoid numerical issues on moving along the main direction. */ 80 delta = hits[0].distance; 81 } 82 83 /* Setup outputs */ 84 if(delta <= delta_solid*0.1 && hits[1].distance == delta) { 85 /* Snap the random walk to the boundary if delta is too small */ 86 fX(set)(dir0, dirs[1]); 87 *hit0 = hits[1]; 88 fX(splat)(dir1, (float)INF); 89 *hit1 = SXD_HIT_NULL; 90 } else { 91 fX(set)(dir0, dirs[0]); 92 *hit0 = hits[0]; 93 if(delta == hits[0].distance) { 94 fX(set)(dir1, dirs[0]); 95 *hit1 = hits[0]; 96 } else if(delta == hits[1].distance) { 97 fX(set)(dir1, dirs[1]); 98 *hit1 = hits[1]; 99 } else { 100 fX(splat)(dir1, 0); 101 *hit1 = SXD_HIT_NULL; 102 } 103 } 104 105 return delta; 106 } 107 108 /* Sample the next direction to walk toward and compute the distance to travel. 109 * If the targeted position does not lie inside the current medium, reject it 110 * and sample a new next step. */ 111 static res_T 112 XD(sample_next_step_robust) 113 (struct sdis_scene* scn, 114 const unsigned current_enc_id, 115 struct ssp_rng* rng, 116 const double pos[DIM], 117 const float delta_solid, 118 float dir0[DIM], /* Sampled direction */ 119 float dir1[DIM], /* Direction used to adjust delta */ 120 struct sXd(hit)* hit0, /* Hit along the sampled direction */ 121 struct sXd(hit)* hit1, /* Hit used to adjust delta */ 122 float* out_delta) 123 { 124 unsigned enc_id = ENCLOSURE_ID_NULL; 125 float delta; 126 float org[DIM]; 127 const size_t MAX_ATTEMPTS = 100; 128 size_t iattempt = 0; 129 res_T res = RES_OK; 130 ASSERT(scn && rng && pos && delta_solid > 0); 131 ASSERT(current_enc_id != ENCLOSURE_ID_NULL); 132 ASSERT(dir0 && dir1 && hit0 && hit1 && out_delta); 133 134 fX_set_dX(org, pos); 135 do { 136 double pos_next[DIM]; 137 138 /* Compute the next step */ 139 delta = XD(sample_next_step) 140 (scn, rng, org, delta_solid, dir0, dir1, hit0, hit1); 141 142 /* Retrieve the medium of the next step */ 143 if(hit0->distance > delta) { 144 XD(move_pos)(dX(set)(pos_next, pos), dir0, delta); 145 res = scene_get_enclosure_id_in_closed_boundaries(scn, pos_next, &enc_id); 146 if(res == RES_BAD_OP) { enc_id = ENCLOSURE_ID_NULL; res = RES_OK; } 147 if(res != RES_OK) goto error; 148 } else { 149 unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; 150 scene_get_enclosure_ids(scn, hit0->prim.prim_id, enc_ids); 151 enc_id = fX(dot)(dir0, hit0->normal) < 0 ? enc_ids[0] : enc_ids[1]; 152 } 153 154 /* Check medium consistency */ 155 if(current_enc_id != enc_id) { 156 #if 0 157 log_warn(scn->dev, 158 "%s: inconsistent medium during the solid random walk -- " 159 "pos=("FORMAT_VECX")\n", FUNC_NAME, SPLITX(pos)); 160 #endif 161 } 162 } while(current_enc_id != enc_id && ++iattempt < MAX_ATTEMPTS); 163 164 /* Handle error */ 165 if(iattempt >= MAX_ATTEMPTS) { 166 log_warn(scn->dev, 167 "%s: could not find a next valid conductive -- pos=("FORMAT_VECX")\n", 168 FUNC_NAME, SPLITX(pos)); 169 res = RES_BAD_OP; 170 goto error; 171 } 172 173 *out_delta = delta; 174 175 exit: 176 return res; 177 error: 178 goto exit; 179 } 180 181 /******************************************************************************* 182 * Handle the volumic power at a given diffusive step 183 ******************************************************************************/ 184 struct XD(handle_volumic_power_args) { 185 /* Forward/backward direction of the sampled diffusive step */ 186 const float* dir0; 187 const float* dir1; 188 189 /* Forward/backward intersections along the sampled diffusive step */ 190 const struct sXd(hit)* hit0; 191 const struct sXd(hit)* hit1; 192 193 /* Physical properties */ 194 double power; /* Volumic power */ 195 double lambda; /* Conductivity */ 196 197 float delta_solid; /* Challenged length of a diffusive step */ 198 float delta; /* Current length of the current diffusive step */ 199 200 size_t picard_order; 201 }; 202 static const struct XD(handle_volumic_power_args) 203 XD(HANDLE_VOLUMIC_POWER_ARGS_NULL) = { 204 NULL, NULL, NULL, NULL, -1, -1, -1, -1, 0 205 }; 206 207 static INLINE int 208 XD(check_handle_volumic_power_args) 209 (const struct XD(handle_volumic_power_args)* args) 210 { 211 ASSERT(args); 212 return args 213 && args->dir0 214 && args->dir1 215 && args->hit0 216 && args->hit1 217 && args->lambda >= 0 218 && args->delta_solid > 0 219 && args->delta >= 0 220 && args->delta_solid >= 0 221 && args->picard_order > 0; 222 } 223 224 static res_T 225 XD(handle_volumic_power) 226 (const struct sdis_scene* scn, 227 const struct XD(handle_volumic_power_args)* args, 228 double* out_power_term, 229 struct temperature* T) 230 { 231 double power_term = 0; 232 res_T res = RES_OK; 233 ASSERT(scn && out_power_term && T && XD(check_handle_volumic_power_args)(args)); 234 235 /* No volumic power. Do nothing */ 236 if(args->power == SDIS_VOLUMIC_POWER_NONE) goto exit; 237 238 /* Check that picardN is not enabled when a volumic power is set since in 239 * this situation the upper bound of the Monte-Carlo weight required by 240 * picardN cannot be known */ 241 if(args->picard_order > 1) { 242 log_err(scn->dev, 243 "%s: invalid not null volumic power '%g' kg/m^3. Could not manage a " 244 "volumic power when the picard order is not equal to 1; Picard order is " 245 "currently set to %lu.\n", 246 FUNC_NAME, args->power, (unsigned long)args->picard_order); 247 res = RES_BAD_ARG; 248 goto error; 249 } 250 251 /* No forward/backward intersection along the sampled direction */ 252 if(SXD_HIT_NONE(args->hit0) && SXD_HIT_NONE(args->hit1)) { 253 const double delta_in_meter = args->delta * scn->fp_to_meter; 254 power_term = delta_in_meter * delta_in_meter / (2.0 * DIM * args->lambda); 255 T->value += args->power * power_term; 256 257 /* An intersection along this diffusive step is find. Use it to statically 258 * correct the power term currently registered */ 259 } else { 260 const double delta_s_adjusted = args->delta_solid * RAY_RANGE_MAX_SCALE; 261 const double delta_s_in_meter = args->delta_solid * scn->fp_to_meter; 262 double h; 263 double h_in_meter; 264 double cos_U_N; 265 float N[DIM] = {0}; 266 267 if(args->delta == args->hit0->distance) { 268 fX(normalize)(N, args->hit0->normal); 269 cos_U_N = fX(dot)(args->dir0, N); 270 271 } else { 272 ASSERT(args->delta == args->hit1->distance); 273 fX(normalize)(N, args->hit1->normal); 274 cos_U_N = fX(dot)(args->dir1, N); 275 } 276 277 h = args->delta * fabs(cos_U_N); 278 h_in_meter = h * scn->fp_to_meter; 279 280 /* The regular power term */ 281 power_term = h_in_meter * h_in_meter / (2.0*args->lambda); 282 283 /* Add the power corrective term. Be careful to use the adjusted 284 * delta_solid to correctly handle the RAY_RANGE_MAX_SCALE factor in the 285 * computation of the limit angle. But keep going with the unmodified 286 * delta_solid in the corrective term since it was the one that was 287 * "wrongly" used in the previous step and that must be corrected. */ 288 if(h == delta_s_adjusted) { 289 power_term += 290 -(delta_s_in_meter * delta_s_in_meter) / (2.0*DIM*args->lambda); 291 292 } else if(h < delta_s_adjusted) { 293 const double sin_a = h / delta_s_adjusted; 294 #if DIM==2 295 /* tmp = sin(2a) / (PI - 2*a) */ 296 const double tmp = sin_a * sqrt(1 - sin_a*sin_a) / acos(sin_a); 297 power_term += 298 -(delta_s_in_meter * delta_s_in_meter) / (4.0*args->lambda) * tmp; 299 #else 300 const double tmp = (sin_a*sin_a*sin_a - sin_a) / (1-sin_a); 301 power_term += 302 (delta_s_in_meter * delta_s_in_meter) / (6.0*args->lambda) * tmp; 303 #endif 304 } 305 T->value += args->power * power_term; 306 } 307 308 exit: 309 *out_power_term = power_term; 310 return res; 311 error: 312 goto exit; 313 } 314 315 /******************************************************************************* 316 * Local function 317 ******************************************************************************/ 318 res_T 319 XD(conductive_path_delta_sphere) 320 (struct sdis_scene* scn, 321 struct rwalk_context* ctx, 322 struct rwalk* rwalk, 323 struct ssp_rng* rng, 324 struct temperature* T) 325 { 326 /* Enclosure/medium in which the conductive path starts */ 327 struct sdis_medium* mdm = NULL; 328 unsigned enc_id = ENCLOSURE_ID_NULL; 329 330 /* Physical properties */ 331 struct solid_props props_ref = SOLID_PROPS_NULL; 332 double green_power_term = 0; 333 334 /* Miscellaneous */ 335 double position_start[DIM] = {0}; 336 size_t istep = 0; /* Help for debug */ 337 res_T res = RES_OK; 338 339 /* Check pre-conditions */ 340 ASSERT(scn && rwalk && rng && T); 341 342 (void)ctx, (void)istep; /* Avoid "unsued variable" warnings */ 343 344 res = scene_get_enclosure_id_in_closed_boundaries(scn, rwalk->vtx.P, &enc_id); 345 if(res != RES_OK) goto error; 346 res = scene_get_enclosure_medium(scn, scene_get_enclosure(scn, enc_id), &mdm); 347 if(res != RES_OK) goto error; 348 ASSERT(sdis_medium_get_type(mdm) == SDIS_SOLID); 349 350 /* Check the random walk consistency */ 351 if(enc_id != rwalk->enc_id) { 352 log_err(scn->dev, "%s: invalid solid random walk. " 353 "Unexpected enclosure -- pos=("FORMAT_VECX")\n", 354 FUNC_NAME, SPLITX(rwalk->vtx.P)); 355 res = RES_BAD_OP_IRRECOVERABLE; 356 goto error; 357 } 358 359 /* Save the submitted position */ 360 dX(set)(position_start, rwalk->vtx.P); 361 362 /* Retrieve the solid properties at the current position. Use them to verify 363 * that those that are supposed to be constant by the conductive random walk 364 * remain the same. Note that we take care of the same constraints on the 365 * solid reinjection since once reinjected, the position of the random walk 366 * is that at the beginning of the conductive random walk. Thus, after a 367 * reinjection, the next line retrieves the properties of the reinjection 368 * position. By comparing them to the properties along the random walk, we 369 * thus verify that the properties are constant throughout the random walk 370 * with respect to the properties of the reinjected position. */ 371 solid_get_properties(mdm, &rwalk->vtx, &props_ref); 372 373 do { /* Solid random walk */ 374 struct XD(handle_volumic_power_args) handle_volpow_args = 375 XD(HANDLE_VOLUMIC_POWER_ARGS_NULL); 376 struct sXd(hit) hit0, hit1; 377 struct solid_props props = SOLID_PROPS_NULL; 378 double power_term = 0; 379 double mu; 380 float delta; /* Random walk numerical parameter */ 381 double delta_m; 382 float dir0[DIM], dir1[DIM]; 383 float org[DIM]; 384 385 /* Fetch solid properties */ 386 res = solid_get_properties(mdm, &rwalk->vtx, &props); 387 if(res != RES_OK) goto error; 388 389 res = check_solid_constant_properties 390 (scn->dev, ctx->green_path != NULL, 0/*use WoS?*/, &props_ref, &props); 391 if(res != RES_OK) goto error; 392 393 /* Check the limit condition 394 * REVIEW Rfo: This can be a bug if the random walk comes from a boundary */ 395 if(SDIS_TEMPERATURE_IS_KNOWN(props.temperature)) { 396 T->value += props.temperature; 397 T->done = 1; 398 399 if(ctx->green_path) { 400 res = green_path_set_limit_vertex 401 (ctx->green_path, mdm, &rwalk->vtx, rwalk->elapsed_time); 402 if(res != RES_OK) goto error; 403 } 404 405 if(ctx->heat_path) { 406 heat_path_get_last_vertex(ctx->heat_path)->weight = T->value; 407 } 408 409 break; 410 } 411 412 fX_set_dX(org, rwalk->vtx.P); 413 414 /* Sample the direction to walk toward and compute the distance to travel */ 415 res = XD(sample_next_step_robust)(scn, enc_id, rng, rwalk->vtx.P, 416 (float)props.delta, dir0, dir1, &hit0, &hit1, &delta); 417 if(res != RES_OK) goto error; 418 419 /* Add the volumic power density to the measured temperature */ 420 handle_volpow_args.dir0 = dir0; 421 handle_volpow_args.dir1 = dir1; 422 handle_volpow_args.hit0 = &hit0; 423 handle_volpow_args.hit1 = &hit1; 424 handle_volpow_args.power = props.power; 425 handle_volpow_args.lambda = props.lambda; 426 handle_volpow_args.delta_solid = (float)props.delta; 427 handle_volpow_args.delta = delta; 428 handle_volpow_args.picard_order = get_picard_order(ctx); 429 res = XD(handle_volumic_power)(scn, &handle_volpow_args, &power_term, T); 430 if(res != RES_OK) goto error; 431 432 /* Register the power term for the green function. Delay its registration 433 * until the end of the conductive path, i.e. the path is valid */ 434 if(ctx->green_path && props.power != SDIS_VOLUMIC_POWER_NONE) { 435 green_power_term += power_term; 436 } 437 438 /* Rewind the time */ 439 delta_m = delta * scn->fp_to_meter; 440 mu = (2*DIM*props.lambda)/(props.rho*props.cp*delta_m*delta_m); 441 res = time_rewind(scn, mu, props.t0, rng, rwalk, ctx, T); 442 if(res != RES_OK) goto error; 443 if(T->done) break; /* Limit condition was reached */ 444 445 /* Define if the random walk hits something along dir0 */ 446 if(hit0.distance > delta) { 447 rwalk->XD(hit) = SXD_HIT_NULL; 448 rwalk->hit_side = SDIS_SIDE_NULL__; 449 } else { 450 rwalk->XD(hit) = hit0; 451 rwalk->hit_side = fX(dot)(hit0.normal, dir0) < 0 ? SDIS_FRONT : SDIS_BACK; 452 } 453 454 /* Update the random walk position */ 455 XD(move_pos)(rwalk->vtx.P, dir0, delta); 456 457 /* Register the new vertex against the heat path */ 458 res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value, 459 SDIS_HEAT_VERTEX_CONDUCTION, (int)ctx->nbranchings); 460 if(res != RES_OK) goto error; 461 462 ++istep; 463 464 /* Keep going while the solid random walk does not hit an interface */ 465 } while(SXD_HIT_NONE(&rwalk->XD(hit))); 466 467 /* Register the power term for the green function */ 468 if(ctx->green_path && props_ref.power != SDIS_VOLUMIC_POWER_NONE) { 469 res = green_path_add_power_term 470 (ctx->green_path, mdm, &rwalk->vtx, green_power_term); 471 if(res != RES_OK) goto error; 472 } 473 474 T->func = XD(boundary_path); 475 rwalk->enc_id = ENCLOSURE_ID_NULL; /* At the interface between 2 media */ 476 477 exit: 478 return res; 479 error: 480 goto exit; 481 } 482 483 #include "sdis_Xd_end.h"