sdis_heat_path_boundary_Xd_solid_fluid_picard1.h (12452B)
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_radiative_env_c.h" 22 #include "sdis_scene_c.h" 23 24 #include <star/ssp.h> 25 26 #include "sdis_Xd_begin.h" 27 28 /******************************************************************************* 29 * Helper functions 30 ******************************************************************************/ 31 static INLINE res_T 32 XD(rwalk_get_Tref) 33 (const struct sdis_scene* scn, 34 const struct rwalk* rwalk, 35 const struct temperature* T, 36 double* out_Tref) 37 { 38 double Tref = SDIS_TEMPERATURE_NONE; 39 res_T res = RES_OK; 40 ASSERT(rwalk && T && out_Tref); 41 42 if(T->done) { 43 /* The path reaches a limit condition, i.e. it goes to the infinity and 44 * fetches the ambient radiative temperature. We do not use the limit 45 * conditions as the reference temperature to make the sampled paths 46 * independant of them. */ 47 struct sdis_radiative_ray ray = SDIS_RADIATIVE_RAY_NULL; 48 ray.dir[0] = rwalk->dir[0]; 49 ray.dir[1] = rwalk->dir[1]; 50 ray.dir[2] = rwalk->dir[2]; 51 ray.time = rwalk->vtx.time; 52 Tref = radiative_env_get_reference_temperature(scn->radenv, &ray); 53 } else { 54 struct sdis_interface_fragment frag; 55 struct sdis_interface* interf = NULL; 56 ASSERT(!SXD_HIT_NONE(&rwalk->XD(hit))); 57 58 /* Fetch the interface where the random walk ends */ 59 interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id); 60 ASSERT(rwalk->hit_side!=SDIS_FRONT || interf->medium_front->type==SDIS_FLUID); 61 ASSERT(rwalk->hit_side!=SDIS_BACK || interf->medium_back->type==SDIS_FLUID); 62 63 /* Fragment on the fluid side of the boundary onto which the rwalk ends */ 64 XD(setup_interface_fragment) 65 (&frag, &rwalk->vtx, &rwalk->XD(hit), rwalk->hit_side); 66 67 Tref = interface_side_get_reference_temperature(interf, &frag); 68 } 69 70 res = XD(check_Tref)(scn, rwalk->vtx.P, Tref, FUNC_NAME); 71 if(res != RES_OK) goto error; 72 73 exit: 74 *out_Tref = Tref; 75 return res; 76 error: 77 Tref = -1; 78 goto exit; 79 } 80 81 /******************************************************************************* 82 * Boundary path between a solid and a fluid 83 ******************************************************************************/ 84 res_T 85 XD(solid_fluid_boundary_picard1_path) 86 (struct sdis_scene* scn, 87 struct rwalk_context* ctx, 88 const struct sdis_interface_fragment* frag, 89 struct rwalk* rwalk, 90 struct ssp_rng* rng, 91 struct temperature* T) 92 { 93 /* Input argument used to handle the net flux */ 94 struct handle_net_flux_args handle_net_flux_args = HANDLE_NET_FLUX_ARGS_NULL; 95 96 /* Input argument used to handle the external net flux */ 97 struct handle_external_net_flux_args handle_external_net_flux_args = 98 HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL; 99 100 /* Input/output arguments of the function used to sample a reinjection */ 101 struct sample_reinjection_step_args samp_reinject_step_args = 102 SAMPLE_REINJECTION_STEP_ARGS_NULL; 103 struct reinjection_step reinject_step = REINJECTION_STEP_NULL; 104 105 /* Temperature and random walk state of the sampled radiative path */ 106 struct temperature T_s; 107 struct rwalk rwalk_s; 108 109 /* Fragment on the fluid side of the boundary */ 110 struct sdis_interface_fragment frag_fluid; 111 112 /* The enclosures split by the boundary */ 113 unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; 114 115 /* Data attached to the boundary */ 116 struct sdis_interface* interf = NULL; 117 struct sdis_medium* solid = NULL; 118 struct sdis_medium* fluid = NULL; 119 120 double h_cond; /* Conductive coefficient */ 121 double h_conv; /* Convective coefficient */ 122 double h_radi_hat; /* Radiative coefficient with That */ 123 double h_hat; /* Sum of h_<conv|cond|rad_hat> */ 124 double p_conv; /* Convective proba */ 125 double p_cond; /* Conductive proba */ 126 127 double epsilon; /* Interface emissivity */ 128 double Tref; /* Reference temperature */ 129 double Tref_s; /* Reference temperature of the sampled radiative path */ 130 double lambda; /* Solid conductivity */ 131 double delta_boundary; /* Orthogonal reinjection dst at the boundary */ 132 double delta; /* Orthogonal fitted reinjection dst at the boundary */ 133 double delta_m; /* delta in meter */ 134 135 double r; 136 struct sdis_heat_vertex hvtx = SDIS_HEAT_VERTEX_NULL; 137 enum sdis_side solid_side = SDIS_SIDE_NULL__; 138 enum sdis_side fluid_side = SDIS_SIDE_NULL__; 139 res_T res = RES_OK; 140 141 ASSERT(scn && rwalk && rng && T && ctx); 142 ASSERT(XD(check_rwalk_fragment_consistency)(rwalk, frag) == RES_OK); 143 144 /* Retrieve the solid and the fluid split by the boundary */ 145 interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id); 146 solid = interface_get_medium(interf, SDIS_FRONT); 147 fluid = interface_get_medium(interf, SDIS_BACK); 148 solid_side = SDIS_FRONT; 149 fluid_side = SDIS_BACK; 150 if(solid->type != SDIS_SOLID) { 151 SWAP(struct sdis_medium*, solid, fluid); 152 SWAP(enum sdis_side, solid_side, fluid_side); 153 ASSERT(fluid->type == SDIS_FLUID); 154 } 155 156 /* Get the enclosures split by the boundary */ 157 scene_get_enclosure_ids(scn, rwalk->XD(hit).prim.prim_id, enc_ids); 158 159 /* Setup a fragment for the fluid side */ 160 frag_fluid = *frag; 161 frag_fluid.side = fluid_side; 162 163 /* Fetch the solid properties */ 164 lambda = solid_get_thermal_conductivity(solid, &rwalk->vtx); 165 delta = solid_get_delta(solid, &rwalk->vtx); 166 167 /* Fetch the boundary emissivity */ 168 epsilon = interface_side_get_emissivity 169 (interf, SDIS_INTERN_SOURCE_ID, &frag_fluid); 170 171 if(epsilon <= 0) { 172 Tref = 0; 173 } else { 174 /* Check the Tref */ 175 Tref = interface_side_get_reference_temperature(interf, &frag_fluid); 176 res = XD(check_Tref)(scn, frag_fluid.P, Tref, FUNC_NAME); 177 if(res != RES_OK) goto error; 178 } 179 180 /* Note that the reinjection distance is *FIXED*. It MUST ensure that the 181 * orthogonal distance from the boundary to the reinjection point is at most 182 * equal to delta. */ 183 delta_boundary = sqrt(DIM) * delta; 184 185 /* Sample a reinjection step */ 186 samp_reinject_step_args.rng = rng; 187 samp_reinject_step_args.solid_enc_id = enc_ids[solid_side]; 188 samp_reinject_step_args.rwalk = rwalk; 189 samp_reinject_step_args.distance = delta_boundary; 190 samp_reinject_step_args.side = solid_side; 191 res = XD(sample_reinjection_step_solid_fluid) 192 (scn, &samp_reinject_step_args, &reinject_step); 193 if(res != RES_OK) goto error; 194 195 /* Define the orthogonal dst from the reinjection pos to the interface */ 196 delta = reinject_step.distance / sqrt(DIM); 197 delta_m = delta * scn->fp_to_meter; 198 199 /* Compute the convective, conductive and the upper bound radiative coef */ 200 h_conv = interface_get_convection_coef(interf, frag); 201 h_cond = lambda / delta_m; 202 if(epsilon <= 0) { 203 h_radi_hat = 0; /* No radiative transfert */ 204 } else { 205 res = scene_check_temperature_range(scn); 206 if(res != RES_OK) { res = RES_BAD_OP_IRRECOVERABLE; goto error; } 207 h_radi_hat = 4.0 * BOLTZMANN_CONSTANT * ctx->That3 * epsilon; 208 } 209 210 /* Compute a global upper bound coefficient */ 211 h_hat = h_conv + h_cond + h_radi_hat; 212 213 /* Compute the probas to switch in solid, fluid or radiative random walk */ 214 p_conv = h_conv / h_hat; 215 p_cond = h_cond / h_hat; 216 217 /* Handle the net flux if any */ 218 handle_net_flux_args.interf = interf; 219 handle_net_flux_args.frag = frag; 220 handle_net_flux_args.green_path = ctx->green_path; 221 handle_net_flux_args.picard_order = get_picard_order(ctx); 222 handle_net_flux_args.h_cond = h_cond; 223 handle_net_flux_args.h_conv = h_conv; 224 handle_net_flux_args.h_radi = h_radi_hat; 225 res = XD(handle_net_flux)(scn, &handle_net_flux_args, T); 226 if(res != RES_OK) goto error; 227 228 /* Handle the external net flux if any */ 229 handle_external_net_flux_args.interf = interf; 230 handle_external_net_flux_args.frag = frag; 231 handle_external_net_flux_args.XD(hit) = &rwalk->XD(hit); 232 handle_external_net_flux_args.green_path = ctx->green_path; 233 handle_external_net_flux_args.picard_order = get_picard_order(ctx); 234 handle_external_net_flux_args.h_cond = h_cond; 235 handle_external_net_flux_args.h_conv = h_conv; 236 handle_external_net_flux_args.h_radi = h_radi_hat; 237 res = XD(handle_external_net_flux)(scn, rng, &handle_external_net_flux_args, T); 238 if(res != RES_OK) goto error; 239 240 /* Fetch the last registered heat path vertex */ 241 if(ctx->heat_path) hvtx = *heat_path_get_last_vertex(ctx->heat_path); 242 243 /* Null collision */ 244 for(;;) { 245 double h_radi; /* Radiative coefficient */ 246 double p_radi; /* Radiative proba */ 247 248 /* Indices of the registered vertex of the sampled radiative path */ 249 size_t ihvtx_radi_begin = 0; 250 size_t ihvtx_radi_end = 0; 251 252 r = ssp_rng_canonical(rng); 253 254 /* Switch in convective path */ 255 if(r < p_conv) { 256 T->func = XD(convective_path); 257 rwalk->enc_id = enc_ids[fluid_side]; 258 rwalk->hit_side = fluid_side; 259 break; 260 } 261 262 /* Switch in conductive path */ 263 if(r < p_conv + p_cond) { 264 struct solid_reinjection_args solid_reinject_args = 265 SOLID_REINJECTION_ARGS_NULL; 266 267 /* Perform the reinjection into the solid */ 268 solid_reinject_args.reinjection = &reinject_step; 269 solid_reinject_args.rwalk_ctx = ctx; 270 solid_reinject_args.rwalk = rwalk; 271 solid_reinject_args.rng = rng; 272 solid_reinject_args.T = T; 273 solid_reinject_args.fp_to_meter = scn->fp_to_meter; 274 res = XD(solid_reinjection)(scn, enc_ids[solid_side], &solid_reinject_args); 275 if(res != RES_OK) goto error; 276 break; 277 } 278 279 /* From there, we know the path is either a radiative path or a 280 * null-collision */ 281 282 if(ctx->heat_path) { 283 /* Fetch the index of the first vertex of the radiative path that is 284 * going to be traced i.e. the last registered vertex */ 285 ihvtx_radi_begin = heat_path_get_vertices_count(ctx->heat_path) - 1; 286 } 287 288 /* Sample a radiative path and get the Tref at its end. */ 289 T_s = *T; 290 rwalk_s = *rwalk; 291 rwalk_s.enc_id = enc_ids[fluid_side]; 292 rwalk_s.hit_side = fluid_side; 293 res = XD(radiative_path)(scn, ctx, &rwalk_s, rng, &T_s); 294 if(res != RES_OK) goto error; 295 296 /* Get the Tref at the end of the candidate radiative path */ 297 res = XD(rwalk_get_Tref)(scn, &rwalk_s, &T_s, &Tref_s); 298 if(res != RES_OK) goto error; 299 300 /* The reference temperatures must be known, as this is a radiative path. 301 * If this is not the case, an error should be reported before this point. 302 * Hence these assertions to detect unexpected behavior */ 303 ASSERT(SDIS_TEMPERATURE_IS_KNOWN(Tref)); 304 ASSERT(SDIS_TEMPERATURE_IS_KNOWN(Tref_s)); 305 306 h_radi = BOLTZMANN_CONSTANT * epsilon * 307 ( Tref*Tref*Tref 308 + Tref*Tref * Tref_s 309 + Tref * Tref_s*Tref_s 310 + Tref_s*Tref_s*Tref_s); 311 312 p_radi = h_radi / h_hat; 313 if(r < p_conv + p_cond + p_radi) { /* Radiative path */ 314 *rwalk = rwalk_s; 315 *T = T_s; 316 break; 317 318 /* Null collision: the sampled path is rejected. */ 319 } else { 320 321 if(ctx->green_path) { 322 /* The limit condition of the green path could be set by the rejected 323 * sampled radiative path. Reset this limit condition. */ 324 green_path_reset_limit(ctx->green_path); 325 } 326 327 if(ctx->heat_path) { 328 /* Set the sampled radiative path as a branch of the current path */ 329 ihvtx_radi_end = heat_path_get_vertices_count(ctx->heat_path); 330 heat_path_increment_sub_path_branch_id 331 (ctx->heat_path, ihvtx_radi_begin, ihvtx_radi_end); 332 333 /* Add a break into the heat path geometry and restart it from the 334 * position of the input random walk. */ 335 res = heat_path_restart(ctx->heat_path, &hvtx); 336 if(res != RES_OK) goto error; 337 } 338 } 339 340 /* Null-collision, looping at the beginning */ 341 } 342 343 exit: 344 return res; 345 error: 346 goto exit; 347 } 348 349 #include "sdis_Xd_end.h"