sdis_realisation_Xd.h (15329B)
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_device_c.h" 17 #include "sdis_heat_path.h" 18 #include "sdis_heat_path_boundary_c.h" 19 #include "sdis_interface_c.h" 20 #include "sdis_log.h" 21 #include "sdis_medium_c.h" 22 #include "sdis_misc.h" 23 #include "sdis_scene_c.h" 24 25 #include <rsys/stretchy_array.h> 26 #include <star/ssp.h> 27 28 #include "sdis_Xd_begin.h" 29 30 /******************************************************************************* 31 * Non generic helper functions 32 ******************************************************************************/ 33 #ifndef SDIS_REALISATION_XD_H 34 #define SDIS_REALISATION_XD_H 35 36 static INLINE res_T 37 check_probe_realisation_args(const struct probe_realisation_args* args) 38 { 39 return args 40 && args->rng 41 && args->enc_id != ENCLOSURE_ID_NULL 42 && args->time >= 0 43 && args->picard_order > 0 44 && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__ 45 ? RES_OK : RES_BAD_ARG; 46 } 47 48 static INLINE res_T 49 check_boundary_realisation_args(const struct boundary_realisation_args* args) 50 { 51 return args 52 && args->rng 53 && args->uv[0] >= 0 54 && args->uv[0] <= 1 55 && args->uv[1] >= 0 56 && args->uv[1] <= 1 57 && args->time >= 0 58 && args->picard_order > 0 59 && (args->side == SDIS_FRONT || args->side == SDIS_BACK) 60 && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__ 61 ? RES_OK : RES_BAD_ARG; 62 } 63 64 static INLINE res_T 65 check_boundary_flux_realisation_args 66 (const struct boundary_flux_realisation_args* args) 67 { 68 return args 69 && args->rng 70 && args->uv[0] >= 0 71 && args->uv[0] <= 1 72 && args->uv[1] >= 0 73 && args->uv[1] <= 1 74 && args->time >= 0 75 && args->picard_order > 0 76 && (args->solid_side == SDIS_FRONT || args->solid_side == SDIS_BACK) 77 && (unsigned)args->diff_algo < SDIS_DIFFUSION_ALGORITHMS_COUNT__ 78 ? RES_OK : RES_BAD_ARG; 79 } 80 #endif /* SDIS_REALISATION_XD_H */ 81 82 /******************************************************************************* 83 * Local functions 84 ******************************************************************************/ 85 res_T 86 XD(sample_coupled_path) 87 (struct sdis_scene* scn, 88 struct rwalk_context* ctx, 89 struct rwalk* rwalk, 90 struct ssp_rng* rng, 91 struct temperature* T) 92 { 93 #ifndef NDEBUG 94 /* Stack that saves the state of each recursion steps. */ 95 struct entry { 96 struct temperature temperature; 97 struct rwalk rwalk; 98 }* stack = NULL; 99 size_t istack = 0; 100 #endif 101 struct sdis_heat_vertex* heat_vtx = NULL; 102 /* Maximum accepted #failures before stopping the realisation */ 103 const size_t MAX_FAILS = 1; 104 res_T res = RES_OK; 105 ASSERT(scn && ctx && rwalk && rng && T); 106 107 ctx->nbranchings += 1; 108 CHK(ctx->nbranchings <= ctx->max_branchings); 109 110 if(ctx->heat_path && T->func == XD(boundary_path)) { 111 heat_vtx = heat_path_get_last_vertex(ctx->heat_path); 112 } 113 114 while(!T->done) { 115 /* Save the current random walk state */ 116 const struct rwalk rwalk_bkp = *rwalk; 117 const struct temperature T_bkp = *T; 118 size_t nfails = 0; /* #failures */ 119 120 #ifndef NDEBUG 121 struct entry e; 122 e.temperature = *T; 123 e.rwalk = *rwalk; 124 sa_push(stack, e); 125 ++istack; 126 #endif 127 128 /* Reject the step if a BAD_OP occurs and retry up to MAX_FAILS times */ 129 do { 130 res = T->func(scn, ctx, rwalk, rng, T); 131 if(res == RES_BAD_OP) { *rwalk = rwalk_bkp; *T = T_bkp; } 132 } while(res == RES_BAD_OP && ++nfails < MAX_FAILS); 133 if(res != RES_OK) { 134 log_err(scn->dev, "%s: reject path (realisation: %lu; branch: %lu)\n", 135 FUNC_NAME, 136 (unsigned long)ctx->irealisation, 137 (unsigned long)ctx->nbranchings); 138 goto error; 139 } 140 141 /* Update the type of the first vertex of the random walks that begin on a 142 * boundary. Indeed, one knows the "right" type of the first vertex only 143 * after the boundary_path execution that defines the sub path to resolve 144 * from the submitted boundary position. Note that if the boundary 145 * temperature is known, the type is let as it. */ 146 if(heat_vtx && !T->done && T->func != XD(boundary_path)) { 147 heat_vtx = heat_path_get_last_vertex(ctx->heat_path); 148 if(T->func == XD(conductive_path)) { 149 heat_vtx->type = SDIS_HEAT_VERTEX_CONDUCTION; 150 } else if(T->func == XD(convective_path)) { 151 heat_vtx->type = SDIS_HEAT_VERTEX_CONVECTION; 152 } else if(T->func == XD(radiative_path)) { 153 heat_vtx->type = SDIS_HEAT_VERTEX_RADIATIVE; 154 } else { 155 FATAL("Unreachable code.\n"); 156 } 157 heat_vtx = NULL; /* Notify that the first vertex is finalized */ 158 } 159 } 160 161 162 exit: 163 #ifndef NDEBUG 164 sa_release(stack); 165 #endif 166 ctx->nbranchings -= 1; 167 return res == RES_BAD_OP_IRRECOVERABLE ? RES_BAD_OP : res; 168 error: 169 goto exit; 170 } 171 172 res_T 173 XD(probe_realisation) 174 (struct sdis_scene* scn, 175 struct probe_realisation_args* args, 176 double* weight) 177 { 178 /* Starting enclosure/medium */ 179 const struct enclosure* enc = NULL; 180 struct sdis_medium* mdm = NULL; 181 182 /* Random walk */ 183 struct rwalk_context ctx = RWALK_CONTEXT_NULL; 184 struct rwalk rwalk = RWALK_NULL; 185 struct temperature T = TEMPERATURE_NULL; 186 187 /* Miscellaneous */ 188 enum sdis_heat_vertex_type type; 189 double t0; 190 double (*get_initial_temperature) 191 (const struct sdis_medium* mdm, 192 const struct sdis_rwalk_vertex* vtx); 193 res_T res = RES_OK; 194 ASSERT(scn && weight && check_probe_realisation_args(args) == RES_OK); 195 196 /* Get the enclosure medium */ 197 enc = scene_get_enclosure(scn, args->enc_id); 198 res = scene_get_enclosure_medium(scn, enc, &mdm); 199 if(res != RES_OK) goto error; 200 201 switch(sdis_medium_get_type(mdm)) { 202 case SDIS_FLUID: 203 T.func = XD(convective_path); 204 get_initial_temperature = fluid_get_temperature; 205 t0 = fluid_get_t0(mdm); 206 break; 207 case SDIS_SOLID: 208 T.func = XD(conductive_path); 209 get_initial_temperature = solid_get_temperature; 210 t0 = solid_get_t0(mdm); 211 break; 212 default: FATAL("Unreachable code\n"); break; 213 } 214 215 dX(set)(rwalk.vtx.P, args->position); 216 rwalk.vtx.time = args->time; 217 218 /* Register the starting position against the heat path */ 219 type = sdis_medium_get_type(mdm) == SDIS_SOLID 220 ? SDIS_HEAT_VERTEX_CONDUCTION 221 : SDIS_HEAT_VERTEX_CONVECTION; 222 res = register_heat_vertex(args->heat_path, &rwalk.vtx, 0, type, 0); 223 if(res != RES_OK) goto error; 224 225 if(t0 >= rwalk.vtx.time) { 226 double tmp; 227 /* Check the initial condition. */ 228 rwalk.vtx.time = t0; 229 tmp = get_initial_temperature(mdm, &rwalk.vtx); 230 if(SDIS_TEMPERATURE_IS_KNOWN(tmp)) { 231 *weight = tmp; 232 goto exit; 233 } 234 /* The initial condition should have been reached */ 235 log_err(scn->dev, 236 "%s: undefined initial condition. " 237 "The time is %g but the temperature remains unknown.\n", 238 FUNC_NAME, t0); 239 res = RES_BAD_OP; 240 goto error; 241 } 242 243 rwalk.XD(hit) = SXD_HIT_NULL; 244 rwalk.enc_id = args->enc_id; 245 246 ctx.green_path = args->green_path; 247 ctx.heat_path = args->heat_path; 248 ctx.Tmin = scn->tmin; 249 ctx.Tmin2 = ctx.Tmin * ctx.Tmin; 250 ctx.Tmin3 = ctx.Tmin * ctx.Tmin2; 251 ctx.That = scn->tmax; 252 ctx.That2 = ctx.That * ctx.That; 253 ctx.That3 = ctx.That * ctx.That2; 254 ctx.max_branchings = args->picard_order - 1; 255 ctx.irealisation = args->irealisation; 256 ctx.diff_algo = args->diff_algo; 257 258 res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T); 259 if(res != RES_OK) goto error; 260 261 ASSERT(SDIS_TEMPERATURE_IS_KNOWN(T.value)); 262 *weight = T.value; 263 264 exit: 265 return res; 266 error: 267 goto exit; 268 } 269 270 res_T 271 XD(boundary_realisation) 272 (struct sdis_scene* scn, 273 struct boundary_realisation_args* args, 274 double* weight) 275 { 276 struct rwalk_context ctx = RWALK_CONTEXT_NULL; 277 struct rwalk rwalk = RWALK_NULL; 278 struct temperature T = TEMPERATURE_NULL; 279 struct sXd(attrib) attr; 280 #if SDIS_XD_DIMENSION == 2 281 float st; 282 #else 283 float st[2]; 284 #endif 285 res_T res = RES_OK; 286 ASSERT(scn && weight && check_boundary_realisation_args(args) == RES_OK); 287 288 T.func = XD(boundary_path); 289 rwalk.hit_side = args->side; 290 rwalk.XD(hit).distance = 0; 291 rwalk.vtx.time = args->time; 292 rwalk.enc_id = ENCLOSURE_ID_NULL; /* At an interface between 2 enclosures */ 293 294 #if SDIS_XD_DIMENSION == 2 295 st = (float)args->uv[0]; 296 #else 297 f2_set_d2(st, args->uv); 298 #endif 299 300 /* Fetch the primitive */ 301 SXD(scene_view_get_primitive 302 (scn->sXd(view), (unsigned int)args->iprim, &rwalk.XD(hit).prim)); 303 304 /* Retrieve the world space position of the probe onto the primitive */ 305 SXD(primitive_get_attrib(&rwalk.XD(hit).prim, SXD_POSITION, st, &attr)); 306 dX_set_fX(rwalk.vtx.P, attr.value); 307 308 /* Retrieve the primitive normal */ 309 SXD(primitive_get_attrib(&rwalk.XD(hit).prim, SXD_GEOMETRY_NORMAL, st, &attr)); 310 fX(set)(rwalk.XD(hit).normal, attr.value); 311 312 #if SDIS_XD_DIMENSION==2 313 rwalk.XD(hit).u = st; 314 #else 315 f2_set(rwalk.XD(hit).uv, st); 316 #endif 317 318 res = register_heat_vertex(args->heat_path, &rwalk.vtx, 0/*weight*/, 319 SDIS_HEAT_VERTEX_CONDUCTION, 0/*Branch id*/); 320 if(res != RES_OK) goto error; 321 322 ctx.green_path = args->green_path; 323 ctx.heat_path = args->heat_path; 324 ctx.Tmin = scn->tmin; 325 ctx.Tmin2 = ctx.Tmin * ctx.Tmin; 326 ctx.Tmin3 = ctx.Tmin * ctx.Tmin2; 327 ctx.That = scn->tmax; 328 ctx.That2 = ctx.That * ctx.That; 329 ctx.That3 = ctx.That * ctx.That2; 330 ctx.max_branchings = args->picard_order - 1; 331 ctx.irealisation = args->irealisation; 332 ctx.diff_algo = args->diff_algo; 333 334 res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T); 335 if(res != RES_OK) goto error; 336 337 *weight = T.value; 338 339 exit: 340 return res; 341 error: 342 goto exit; 343 } 344 345 res_T 346 XD(boundary_flux_realisation) 347 (struct sdis_scene* scn, 348 struct boundary_flux_realisation_args* args, 349 struct bound_flux_result* result) 350 { 351 /* Random walk */ 352 struct rwalk_context ctx = RWALK_CONTEXT_NULL; 353 struct rwalk rwalk = RWALK_NULL; 354 struct temperature T = TEMPERATURE_NULL; 355 356 /* Boundary */ 357 struct sXd(attrib) attr; 358 struct sXd(primitive) prim; 359 #if SDIS_XD_DIMENSION == 2 360 float st; 361 #else 362 float st[2]; 363 #endif 364 double P[SDIS_XD_DIMENSION]; 365 float N[SDIS_XD_DIMENSION]; 366 unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; 367 enum sdis_side fluid_side; 368 369 /* Miscellaneous */ 370 double Tmin, Tmin2, Tmin3; 371 double That, That2, That3; 372 res_T res = RES_OK; 373 char compute_radiative; 374 char compute_convective; 375 376 ASSERT(scn && result && check_boundary_flux_realisation_args(args) == RES_OK); 377 378 #if SDIS_XD_DIMENSION == 2 379 #define SET_PARAM(Dest, Src) (Dest).u = (Src); 380 st = (float)args->uv[0]; 381 #else 382 #define SET_PARAM(Dest, Src) f2_set((Dest).uv, (Src)); 383 f2_set_d2(st, args->uv); 384 #endif 385 386 Tmin = scn->tmin; 387 Tmin2 = Tmin * Tmin; 388 Tmin3 = Tmin * Tmin2; 389 That = scn->tmax; 390 That2 = That * That; 391 That3 = That * That2; 392 393 fluid_side = (args->solid_side/*solid*/==SDIS_FRONT) ? SDIS_BACK : SDIS_FRONT; 394 395 compute_radiative = (args->flux_mask & FLUX_FLAG_RADIATIVE) != 0; 396 compute_convective = (args->flux_mask & FLUX_FLAG_CONVECTIVE) != 0; 397 398 /* Fetch the primitive */ 399 SXD(scene_view_get_primitive(scn->sXd(view), (unsigned int)args->iprim, &prim)); 400 401 /* Retrieve the world space position of the probe onto the primitive */ 402 SXD(primitive_get_attrib(&prim, SXD_POSITION, st, &attr)); 403 dX_set_fX(P, attr.value); 404 405 /* Retrieve the primitive normal */ 406 SXD(primitive_get_attrib(&prim, SXD_GEOMETRY_NORMAL, st, &attr)); 407 fX(set)(N, attr.value); 408 409 #define RESET_WALK(Side, EncId) { \ 410 rwalk = RWALK_NULL; \ 411 rwalk.hit_side = (Side); \ 412 rwalk.XD(hit).distance = 0; \ 413 rwalk.vtx.time = args->time; \ 414 rwalk.enc_id = (EncId); \ 415 rwalk.XD(hit).prim = prim; \ 416 SET_PARAM(rwalk.XD(hit), st); \ 417 ctx.Tmin = Tmin; \ 418 ctx.Tmin3 = Tmin3; \ 419 ctx.That = That; \ 420 ctx.That2 = That2; \ 421 ctx.That3 = That3; \ 422 ctx.max_branchings = args->picard_order - 1; \ 423 ctx.irealisation = args->irealisation; \ 424 ctx.diff_algo = args->diff_algo; \ 425 dX(set)(rwalk.vtx.P, P); \ 426 fX(set)(rwalk.XD(hit).normal, N); \ 427 T = TEMPERATURE_NULL; \ 428 } (void)0 429 430 /* Compute boundary temperature */ 431 RESET_WALK(args->solid_side, ENCLOSURE_ID_NULL); 432 T.func = XD(boundary_path); 433 res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T); 434 if(res != RES_OK) return res; 435 result->Tboundary = T.value; 436 437 /* Get the enclosures */ 438 scene_get_enclosure_ids(scn, (unsigned)args->iprim, enc_ids); 439 440 /* Compute radiative temperature */ 441 if(compute_radiative) { 442 RESET_WALK(fluid_side, enc_ids[fluid_side]); 443 T.func = XD(radiative_path); 444 res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T); 445 if(res != RES_OK) return res; 446 ASSERT(SDIS_TEMPERATURE_IS_KNOWN(T.value)); 447 result->Tradiative = T.value; 448 } 449 450 /* Compute fluid temperature */ 451 if(compute_convective) { 452 RESET_WALK(fluid_side, enc_ids[fluid_side]); 453 454 /* Check whether the temperature of the fluid is known by querying it from 455 * its boundary. This makes it possible to handle situations where fluids 456 * are used as Robin's boundary condition. In this case, a geometric 457 * enclosure may have several fluids, each defining the temperature of a 458 * boundary condition. Sampling of convective paths in such an enclosure is 459 * forbidden since this enclosure does not exist for this path space: it is 460 * beyond its boundary */ 461 res = XD(query_medium_temperature_from_boundary)(scn, &ctx, &rwalk, &T); 462 if(res != RES_OK) return res; 463 464 /* Robin's boundary condition */ 465 if(T.done) { 466 result->Tfluid = T.value; 467 468 /* Sample a convective path */ 469 } else { 470 T.func = XD(convective_path); 471 res = XD(sample_coupled_path)(scn, &ctx, &rwalk, args->rng, &T); 472 if(res != RES_OK) return res; 473 result->Tfluid = T.value; 474 } 475 } 476 477 #undef SET_PARAM 478 #undef RESET_WALK 479 480 return RES_OK; 481 } 482 483 #include "sdis_Xd_end.h"