sdis_heat_path_boundary_Xd_solid_fluid_picardN.h (15359B)
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_green.h" 17 #include "sdis_heat_path_boundary_c.h" 18 #include "sdis_interface_c.h" 19 #include "sdis_medium_c.h" 20 #include "sdis_misc.h" 21 #include "sdis_realisation.h" 22 #include "sdis_scene_c.h" 23 24 #include <star/ssp.h> 25 26 #include "sdis_Xd_begin.h" 27 28 /******************************************************************************* 29 * Non generic helper functions 30 ******************************************************************************/ 31 #ifndef SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARD_N_H 32 #define SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARD_N_H 33 34 static INLINE res_T 35 check_net_flux 36 (struct sdis_scene* scn, 37 const struct sdis_interface* interf, 38 const struct sdis_interface_fragment* frag, 39 const size_t picard_order) 40 { 41 double phi; 42 res_T res = RES_OK; 43 ASSERT(scn && interf && frag && picard_order > 1); 44 45 phi = interface_side_get_flux(interf, frag); 46 if(phi != SDIS_FLUX_NONE && phi != 0) { 47 log_err(scn->dev, 48 "%s: invalid flux '%g' W/m^2. Could not manage a flux != 0 when the " 49 "picard order is not equal to 1; Picard order is currently set to %lu.\n", 50 FUNC_NAME, phi, (unsigned long)picard_order); 51 res = RES_BAD_ARG; 52 goto error; 53 } 54 55 exit: 56 return res; 57 error: 58 goto exit; 59 } 60 61 #endif /* SDIS_HEAT_PATH_BOUNDARY_XD_SOLID_FLUID_PICARD_N_H */ 62 63 /******************************************************************************* 64 * Generic helper functions 65 ******************************************************************************/ 66 static INLINE res_T 67 XD(sample_path) 68 (struct sdis_scene* scn, 69 const struct rwalk* rwalk_from, 70 struct rwalk_context* ctx, 71 struct ssp_rng* rng, 72 struct temperature* T) 73 { 74 struct rwalk rwalk = RWALK_NULL; 75 res_T res = RES_OK; 76 ASSERT(rwalk_from && rng && T); 77 78 /* Clean-up the output variable */ 79 *T = TEMPERATURE_NULL; 80 T->func = XD(boundary_path); 81 82 /* Init the random walk */ 83 rwalk.vtx = rwalk_from->vtx; 84 rwalk.enc_id = rwalk_from->enc_id; 85 rwalk.XD(hit) = rwalk_from->XD(hit); 86 rwalk.hit_side = rwalk_from->hit_side; 87 88 /* Start the registration of a new heat path */ 89 if(ctx->heat_path) { 90 struct sdis_heat_vertex heat_vtx = SDIS_HEAT_VERTEX_NULL; 91 92 heat_vtx.P[0] = rwalk.vtx.P[0]; 93 heat_vtx.P[1] = rwalk.vtx.P[1]; 94 heat_vtx.P[2] = rwalk.vtx.P[2]; 95 heat_vtx.time = rwalk.vtx.time; 96 heat_vtx.weight = 0; 97 heat_vtx.type = SDIS_HEAT_VERTEX_RADIATIVE; 98 heat_vtx.branch_id = (int)ctx->nbranchings + 1; 99 100 res = heat_path_restart(ctx->heat_path, &heat_vtx); 101 if(res != RES_OK) goto error; 102 } 103 104 /* Sample the path */ 105 res = XD(sample_coupled_path)(scn, ctx, &rwalk, rng, T); 106 if(res != RES_OK) goto error; 107 108 /* Check the returned temperature */ 109 ASSERT(T->done); 110 if(T->value < scn->tmin || scn->tmax < T->value) { 111 log_err(scn->dev, "%s: invalid temperature range `[%g, %g]K` regarding the " 112 "retrieved temperature %gK.\n", FUNC_NAME, scn->tmin, scn->tmax, T->value); 113 res = RES_BAD_OP_IRRECOVERABLE; 114 goto error; 115 } 116 117 exit: 118 return res; 119 error: 120 goto exit; 121 } 122 123 /******************************************************************************* 124 * Boundary path between a solid and a fluid 125 ******************************************************************************/ 126 res_T 127 XD(solid_fluid_boundary_picardN_path) 128 (struct sdis_scene* scn, 129 struct rwalk_context* ctx, 130 const struct sdis_interface_fragment* frag, 131 struct rwalk* rwalk, 132 struct ssp_rng* rng, 133 struct temperature* T) 134 { 135 /* Input/output arguments of the function used to sample a reinjection */ 136 struct sample_reinjection_step_args samp_reinject_step_args = 137 SAMPLE_REINJECTION_STEP_ARGS_NULL; 138 struct reinjection_step reinject_step = REINJECTION_STEP_NULL; 139 140 /* Fragment on the fluid side of the boundary */ 141 struct sdis_interface_fragment frag_fluid; 142 143 /* Vertex of the heat path */ 144 struct sdis_heat_vertex hvtx = SDIS_HEAT_VERTEX_NULL; 145 struct sdis_heat_vertex hvtx_s = SDIS_HEAT_VERTEX_NULL; 146 147 /* The enclosures split by the boundary */ 148 unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; 149 150 /* Data attached to the boundary */ 151 struct sdis_interface* interf = NULL; 152 struct sdis_medium* solid = NULL; 153 struct sdis_medium* fluid = NULL; 154 155 double h_cond; /* Conductive coefficient */ 156 double h_conv; /* Convective coefficient */ 157 double h_radi_hat; /* Radiative coefficient with That */ 158 double h_hat; /* Sum of h_<conv|cond|rad_hat> */ 159 double p_conv; /* Convective proba */ 160 double p_cond; /* Conductive proba */ 161 162 /* Min/Max Temperatures */ 163 double Tmin, Tmin2, Tmin3; 164 double That, That2, That3; 165 166 double epsilon; /* Interface emissivity */ 167 double lambda; /* Solid conductivity */ 168 double delta_boundary; /* Orthogonal reinjection dst at the boundary */ 169 double delta; /* Orthogonal fitted reinjection dst at the boundary */ 170 171 double r; 172 enum sdis_side solid_side = SDIS_SIDE_NULL__; 173 enum sdis_side fluid_side = SDIS_SIDE_NULL__; 174 res_T res = RES_OK; 175 176 ASSERT(scn && rwalk && rng && T && ctx); 177 ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag) == RES_OK); 178 179 /* Fetch the Min/max temperature */ 180 Tmin = ctx->Tmin; 181 Tmin2 = ctx->Tmin2; 182 Tmin3 = ctx->Tmin3; 183 That = ctx->That; 184 That2 = ctx->That2; 185 That3 = ctx->That3; 186 187 /* Retrieve the solid and the fluid split by the boundary */ 188 interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id); 189 solid = interface_get_medium(interf, SDIS_FRONT); 190 fluid = interface_get_medium(interf, SDIS_BACK); 191 solid_side = SDIS_FRONT; 192 fluid_side = SDIS_BACK; 193 if(solid->type != SDIS_SOLID) { 194 SWAP(struct sdis_medium*, solid, fluid); 195 SWAP(enum sdis_side, solid_side, fluid_side); 196 ASSERT(fluid->type == SDIS_FLUID); 197 } 198 199 /* Get the enclosures split by the boundary */ 200 scene_get_enclosure_ids(scn, rwalk->XD(hit).prim.prim_id, enc_ids); 201 202 /* Check that no net flux is set for this interface since the provided 203 * picardN algorithm does not handle it */ 204 res = check_net_flux(scn, interf, frag, get_picard_order(ctx)); 205 if(res != RES_OK) goto error; 206 207 /* Setup a fragment for the fluid side */ 208 frag_fluid = *frag; 209 frag_fluid.side = fluid_side; 210 211 /* Fetch the solid properties */ 212 lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx); 213 delta = solid_get_delta(solid, &rwalk->vtx); 214 215 /* Fetch the boundary emissivity */ 216 epsilon = interface_side_get_emissivity 217 (interf, SDIS_INTERN_SOURCE_ID, &frag_fluid); 218 219 /* Note that the reinjection distance is *FIXED*. It MUST ensure that the 220 * orthogonal distance from the boundary to the reinjection point is at most 221 * equal to delta. */ 222 delta_boundary = sqrt(DIM) * delta; 223 224 /* Sample a reinjection step */ 225 samp_reinject_step_args.rng = rng; 226 samp_reinject_step_args.solid_enc_id = enc_ids[solid_side]; 227 samp_reinject_step_args.rwalk = rwalk; 228 samp_reinject_step_args.distance = delta_boundary; 229 samp_reinject_step_args.side = solid_side; 230 res = XD(sample_reinjection_step_solid_fluid) 231 (scn, &samp_reinject_step_args, &reinject_step); 232 if(res != RES_OK) goto error; 233 234 /* Define the orthogonal dst from the reinjection pos to the interface */ 235 delta = reinject_step.distance / sqrt(DIM); 236 237 /* Compute the convective, conductive and the upper bound radiative coef */ 238 h_conv = interface_get_convection_coef(interf, frag); 239 h_cond = lambda / (delta * scn->fp_to_meter); 240 h_radi_hat = epsilon > 0 ? 4.0 * BOLTZMANN_CONSTANT * That3 * epsilon : 0; 241 242 if(epsilon <= 0) { 243 h_radi_hat = 0; /* No radiative transfert */ 244 } else { 245 res = scene_check_temperature_range(scn); 246 if(res != RES_OK) { res = RES_BAD_OP_IRRECOVERABLE; goto error; } 247 h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * That3 * epsilon; 248 } 249 250 /* Compute a global upper bound coefficient */ 251 h_hat = h_conv + h_cond + h_radi_hat; 252 253 /* Compute the probas to switch in convection or conduction */ 254 p_conv = h_conv / h_hat; 255 p_cond = h_cond / h_hat; 256 257 /* Fetch the last registered heat path vertex */ 258 if(ctx->heat_path) hvtx = *heat_path_get_last_vertex(ctx->heat_path); 259 260 /* Null collision main loop */ 261 for(;;) { 262 /* Temperature and random walk state of the sampled radiative path */ 263 struct temperature T_s; 264 struct rwalk rwalk_s; 265 266 double h_radi, h_radi_min, h_radi_max; /* Radiative coefficients */ 267 double p_radi, p_radi_min, p_radi_max; /* Radiative probas */ 268 double T0, T1, T2, T3, T4, T5; /* Computed temperatures */ 269 270 /* Indices of the registered vertex of the sampled radiative path */ 271 size_t ihvtx_radi_begin = 0; 272 size_t ihvtx_radi_end = 0; 273 274 r = ssp_rng_canonical(rng); 275 276 /* Switch in convective path */ 277 if(r < p_conv) { 278 T->func = XD(convective_path); 279 rwalk->enc_id = enc_ids[fluid_side]; 280 rwalk->hit_side = fluid_side; 281 break; 282 } 283 284 /* Switch in conductive path */ 285 if(r < p_conv + p_cond) { 286 struct solid_reinjection_args solid_reinject_args = 287 SOLID_REINJECTION_ARGS_NULL; 288 289 /* Perform the reinjection into the solid */ 290 solid_reinject_args.reinjection = &reinject_step; 291 solid_reinject_args.rwalk_ctx = ctx; 292 solid_reinject_args.rwalk = rwalk; 293 solid_reinject_args.rng = rng; 294 solid_reinject_args.T = T; 295 solid_reinject_args.fp_to_meter = scn->fp_to_meter; 296 res = XD(solid_reinjection)(scn, enc_ids[solid_side], &solid_reinject_args); 297 if(res != RES_OK) goto error; 298 break; 299 } 300 301 if(ctx->heat_path) { 302 /* Fetch the index of the first vertex of the radiative path that is 303 * going to be traced i.e. the last registered vertex */ 304 ihvtx_radi_begin = heat_path_get_vertices_count(ctx->heat_path) - 1; 305 } 306 307 /* Sample a radiative path */ 308 T_s = *T; 309 rwalk_s = *rwalk; 310 rwalk_s.enc_id = enc_ids[fluid_side]; 311 rwalk_s.hit_side = fluid_side; 312 res = XD(radiative_path)(scn, ctx, &rwalk_s, rng, &T_s); 313 if(res != RES_OK) goto error; 314 315 if(ctx->heat_path) { 316 /* Fetch the index after the last registered vertex of the sampled 317 * radiative path */ 318 ihvtx_radi_end = heat_path_get_vertices_count(ctx->heat_path); 319 } 320 321 /* Fetch the last registered heat path vertex of the radiative path */ 322 if(ctx->heat_path) hvtx_s = *heat_path_get_last_vertex(ctx->heat_path); 323 324 h_radi_min = 4.0 * BOLTZMANN_CONSTANT * Tmin3 * epsilon; 325 p_radi_min = h_radi_min / h_hat; 326 327 /* Switch in radiative path */ 328 if(r < p_conv + p_cond + p_radi_min) { *rwalk = rwalk_s; *T = T_s; break; } 329 330 /* Define some helper macros */ 331 #define SWITCH_IN_RADIATIVE { \ 332 *rwalk = rwalk_s; *T = T_s; \ 333 res = heat_path_restart(ctx->heat_path, &hvtx_s); \ 334 if(res != RES_OK) goto error; \ 335 } (void)0 336 337 #define NULL_COLLISION { \ 338 res = heat_path_restart(ctx->heat_path, &hvtx); \ 339 if(res != RES_OK) goto error; \ 340 if(ctx->heat_path) { \ 341 heat_path_increment_sub_path_branch_id \ 342 (ctx->heat_path, ihvtx_radi_begin, ihvtx_radi_end); \ 343 } \ 344 } (void)0 345 346 #define COMPUTE_TEMPERATURE(Result, RWalk, Temp) { \ 347 struct temperature T_p; \ 348 if((Temp)->done) { /* Ambient radiative temperature */ \ 349 ASSERT(SXD_HIT_NONE(&(RWalk)->XD(hit))); \ 350 T_p = *(Temp); \ 351 } else { \ 352 res = XD(sample_path)(scn, RWalk, ctx, rng, &T_p); \ 353 if(res != RES_OK) goto error; \ 354 } \ 355 Result = T_p.value; \ 356 } (void)0 357 358 #define CHECK_PMIN_PMAX { \ 359 p_radi_min = h_radi_min*epsilon / h_hat; \ 360 p_radi_max = h_radi_max*epsilon / h_hat; \ 361 if(r < p_conv + p_cond + p_radi_min) { SWITCH_IN_RADIATIVE; break; } \ 362 if(r > p_conv + p_cond + p_radi_max) { NULL_COLLISION; continue; } \ 363 } (void)0 364 365 /* Sample a 1st heat path at the end of the radiative path */ 366 COMPUTE_TEMPERATURE(T0, &rwalk_s, &T_s); 367 h_radi_min = BOLTZMANN_CONSTANT*(Tmin3 + 3*Tmin2*T0); 368 h_radi_max = BOLTZMANN_CONSTANT*(That3 + 3*That2*T0); 369 CHECK_PMIN_PMAX; 370 371 /* Sample a 2nd heat path at the end of the radiative path */ 372 COMPUTE_TEMPERATURE(T1, &rwalk_s, &T_s); 373 h_radi_min = BOLTZMANN_CONSTANT*(Tmin3 + Tmin2*T0 + 2*Tmin*T0*T1); 374 h_radi_max = BOLTZMANN_CONSTANT*(That3 + That2*T0 + 2*That*T0*T1); 375 CHECK_PMIN_PMAX; 376 377 /* Sample a 3rd heat path at the end of the radiative path */ 378 COMPUTE_TEMPERATURE(T2, &rwalk_s, &T_s); 379 h_radi_min = BOLTZMANN_CONSTANT*(Tmin3 + Tmin2*T0 + Tmin*T0*T1 + T0*T1*T2); 380 h_radi_max = BOLTZMANN_CONSTANT*(That3 + That2*T0 + That*T0*T1 + T0*T1*T2); 381 CHECK_PMIN_PMAX; 382 383 /* Sample a 1st heat path at the current position onto the interface */ 384 COMPUTE_TEMPERATURE(T3, rwalk, T); 385 h_radi_min = BOLTZMANN_CONSTANT*(Tmin2*T3 + Tmin*T0*T3 + T0*T1*T3 + T0*T1*T2); 386 h_radi_max = BOLTZMANN_CONSTANT*(That2*T3 + That*T0*T3 + T0*T1*T3 + T0*T1*T2); 387 CHECK_PMIN_PMAX; 388 389 /* Sample a 2nd heat path at the current position onto the interface */ 390 COMPUTE_TEMPERATURE(T4, rwalk, T); 391 h_radi_min = BOLTZMANN_CONSTANT*(Tmin*T3*T4 + T0*T3*T4 + T0*T1*T3 + T0*T1*T2); 392 h_radi_max = BOLTZMANN_CONSTANT*(That*T3*T4 + T0*T3*T4 + T0*T1*T3 + T0*T1*T2); 393 CHECK_PMIN_PMAX; 394 395 /* Sample a 3rd heat path at the current position onto the interface */ 396 COMPUTE_TEMPERATURE(T5, rwalk, T); 397 h_radi = BOLTZMANN_CONSTANT*(T3*T4*T5 + T0*T3*T4 + T0*T1*T3 + T0*T1*T2); 398 p_radi = h_radi * epsilon / h_hat; 399 400 /* Switch in radiative path */ 401 if(r < p_cond + p_conv + p_radi) { SWITCH_IN_RADIATIVE; break; } 402 403 /* Null-collision, looping at the beginning */ 404 NULL_COLLISION; 405 406 #undef SWITCH_IN_RADIATIVE 407 #undef NULL_COLLISION 408 #undef COMPUTE_TEMPERATURE 409 #undef CHECK_PMIN_PMAX 410 } 411 412 exit: 413 return res; 414 error: 415 goto exit; 416 } 417 418 #include "sdis_Xd_end.h" 419