sdis_heat_path_convective_Xd.h (10416B)
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_green.h" 18 #include "sdis_heat_path.h" 19 #include "sdis_medium_c.h" 20 #include "sdis_scene_c.h" 21 22 #include <star/ssp.h> 23 24 #include "sdis_Xd_begin.h" 25 26 /******************************************************************************* 27 * Non generic helper functions 28 ******************************************************************************/ 29 #ifndef SDIS_HEAT_PATH_CONVECTIVE_XD_H 30 #define SDIS_HEAT_PATH_CONVECTIVE_XD_H 31 32 static res_T 33 check_fluid_constant_properties 34 (struct sdis_device* dev, 35 const struct fluid_props* props_ref, 36 const struct fluid_props* props) 37 { 38 res_T res = RES_OK; 39 ASSERT(dev && props_ref && props); 40 41 if(props_ref->rho != props->rho) { 42 log_err(dev, 43 "%s: invalid volumic mass. One assumes a constant volumic mass for " 44 "the whole fluid.\n", FUNC_NAME); 45 res = RES_BAD_ARG; 46 goto error; 47 } 48 49 if(props_ref->cp != props->cp) { 50 log_err(dev, 51 "%s: invalid calorific capacity. One assumes a constant calorific " 52 "capacity for the whole fluid.\n", FUNC_NAME); 53 res = RES_BAD_ARG; 54 goto error; 55 } 56 57 exit: 58 return res; 59 error: 60 goto exit; 61 } 62 63 #endif /* SDIS_HEAT_PATH_CONVECTIVE_XD_H */ 64 65 /******************************************************************************* 66 * Helper functions 67 ******************************************************************************/ 68 static res_T 69 XD(handle_known_fluid_temperature) 70 (struct rwalk_context* ctx, 71 struct rwalk* rwalk, 72 struct sdis_medium* mdm, 73 struct temperature* T) 74 { 75 double temperature; 76 int known_temperature; 77 res_T res = RES_OK; 78 ASSERT(ctx && rwalk && T); 79 ASSERT(sdis_medium_get_type(mdm) == SDIS_FLUID); 80 81 temperature = fluid_get_temperature(mdm, &rwalk->vtx); 82 83 /* Check if the temperature is known */ 84 known_temperature = SDIS_TEMPERATURE_IS_KNOWN(temperature); 85 if(!known_temperature) goto exit; 86 87 T->value += temperature; 88 T->done = 1; 89 90 if(ctx->green_path) { 91 res = green_path_set_limit_vertex 92 (ctx->green_path, mdm, &rwalk->vtx, rwalk->elapsed_time); 93 if(res != RES_OK) goto error; 94 } 95 96 if(ctx->heat_path) { 97 heat_path_get_last_vertex(ctx->heat_path)->weight = T->value; 98 } 99 100 exit: 101 return res; 102 error: 103 goto exit; 104 } 105 106 static res_T 107 XD(handle_convective_path_startup) 108 (struct sdis_scene* scn, 109 struct rwalk* rwalk, 110 int* path_starts_in_fluid) 111 { 112 const float range[2] = {FLT_MIN, FLT_MAX}; 113 float dir[DIM] = {0}; 114 float org[DIM] = {0}; 115 res_T res = RES_OK; 116 ASSERT(scn && rwalk && path_starts_in_fluid); 117 118 *path_starts_in_fluid = SXD_HIT_NONE(&rwalk->XD(hit)); 119 if(*path_starts_in_fluid == 0) goto exit; /* Nothing to do */ 120 121 dir[DIM-1] = 1; 122 fX_set_dX(org, rwalk->vtx.P); 123 124 /* Init the path hit field required to define the current enclosure and 125 * fetch the interface data */ 126 SXD(scene_view_trace_ray(scn->sXd(view), org, dir, range, NULL, &rwalk->XD(hit))); 127 if(SXD_HIT_NONE(&rwalk->XD(hit))) { 128 log_err(scn->dev, 129 "%s: the position %g %g %g lies in the surrounding fluid whose " 130 "temperature must be known.\n", FUNC_NAME, SPLIT3(rwalk->vtx.P)); 131 res = RES_BAD_OP; 132 goto error; 133 } 134 135 rwalk->hit_side = fX(dot)(rwalk->XD(hit).normal, dir) < 0 136 ? SDIS_FRONT : SDIS_BACK; 137 138 exit: 139 return res; 140 error: 141 goto exit; 142 } 143 144 static res_T 145 XD(check_enclosure) 146 (struct sdis_scene* scn, 147 const struct rwalk* rwalk, 148 const struct enclosure* enc) 149 { 150 res_T res = RES_OK; 151 ASSERT(scn && rwalk && enc); 152 153 if(enc->medium_id == MEDIUM_ID_MULTI) { 154 /* The enclosures with multiple media are used to describe limit 155 * conditions and therefore they cannot be fetched */ 156 log_err(scn->dev, 157 "%s: enclosure with multiple media at ("FORMAT_VECX"). " 158 "The path should be reached a limit condition before.\n", 159 FUNC_NAME, SPLITX(rwalk->vtx.P)); 160 res = RES_BAD_ARG; 161 goto error; 162 } 163 164 exit: 165 return res; 166 error: 167 goto exit; 168 } 169 170 /******************************************************************************* 171 * Local functions 172 ******************************************************************************/ 173 res_T 174 XD(convective_path) 175 (struct sdis_scene* scn, 176 struct rwalk_context* ctx, 177 struct rwalk* rwalk, 178 struct ssp_rng* rng, 179 struct temperature* T) 180 { 181 /* Properties */ 182 struct fluid_props props_ref = FLUID_PROPS_NULL; 183 const struct sdis_interface* interf = NULL; 184 struct sdis_medium* mdm = NULL; 185 186 /* Enclosure */ 187 unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; 188 const struct enclosure* enc = NULL; 189 190 /* Miscellaneous */ 191 struct sXd(attrib) attr_P, attr_N; 192 struct sXd(hit)* rwalk_hit = NULL; 193 double r; 194 #if SDIS_XD_DIMENSION == 2 195 float st; 196 #else 197 float st[2]; 198 #endif 199 int path_starts_in_fluid; 200 res_T res = RES_OK; 201 202 ASSERT(scn && ctx && rwalk && rng && T); 203 (void)rng, (void)ctx; /* Avoid "unsued variable" warnings */ 204 205 rwalk_hit = &rwalk->XD(hit); 206 207 /* Get the enclosure medium */ 208 enc = scene_get_enclosure(scn, rwalk->enc_id); 209 if((res = XD(check_enclosure)(scn, rwalk, enc)) != RES_OK) goto error; 210 if((res = scene_get_enclosure_medium(scn, enc, &mdm)) != RES_OK) goto error; 211 ASSERT(sdis_medium_get_type(mdm) == SDIS_FLUID); 212 213 res = XD(handle_known_fluid_temperature)(ctx, rwalk, mdm, T); 214 if(res != RES_OK) goto error; 215 if(T->done) goto exit; /* The fluid temperature is known */ 216 217 /* Setup the missing random walk member variables when the convective path 218 * starts from the fluid */ 219 res = XD(handle_convective_path_startup)(scn, rwalk, &path_starts_in_fluid); 220 if(res != RES_OK) goto error; 221 222 /* Retrieve the fluid properties at the current position. Use them to verify 223 * that those that are supposed to be constant by the convective random walk 224 * remain the same. */ 225 res = fluid_get_properties(mdm, &rwalk->vtx, &props_ref); 226 if(res != RES_OK) goto error; 227 228 /* The hc upper bound can be 0 if h is uniformly 0. In that case the result 229 * is the initial condition. */ 230 if(enc->hc_upper_bound == 0) { 231 ASSERT(path_starts_in_fluid); /* Cannot be in the fluid without starting there. */ 232 rwalk->vtx.time = props_ref.t0; 233 res = XD(handle_known_fluid_temperature)(ctx, rwalk, mdm, T); 234 if(res != RES_OK) goto error; 235 if(T->done) { 236 goto exit; /* Stop the random walk */ 237 } else { 238 log_err(scn->dev, "%s: undefined initial condition.", FUNC_NAME); 239 res = RES_BAD_OP; 240 goto error; 241 } 242 } 243 244 /* Sample time until init condition is reached or a true convection occurs. */ 245 for(;;) { 246 struct sdis_interface_fragment frag; 247 struct sXd(primitive) prim; 248 struct fluid_props props = FLUID_PROPS_NULL; 249 double hc; 250 double mu; 251 252 /* Fetch fluid properties */ 253 res = fluid_get_properties(mdm, &rwalk->vtx, &props); 254 if(res != RES_OK) goto error; 255 256 res = check_fluid_constant_properties(scn->dev, &props_ref, &props); 257 if(res != RES_OK) goto error; 258 259 /* Sample the time using the upper bound. */ 260 mu = enc->hc_upper_bound / (props.rho * props.cp) * enc->S_over_V; 261 res = time_rewind(scn, mu, props.t0, rng, rwalk, ctx, T); 262 if(res != RES_OK) goto error; 263 if(T->done) break; /* Limit condition was reached */ 264 265 /* Uniformly sample the enclosure. */ 266 #if DIM == 2 267 SXD(scene_view_sample 268 (enc->sXd(view), 269 ssp_rng_canonical_float(rng), 270 ssp_rng_canonical_float(rng), 271 &prim, &rwalk_hit->u)); 272 st = rwalk_hit->u; 273 #else 274 SXD(scene_view_sample 275 (enc->sXd(view), 276 ssp_rng_canonical_float(rng), 277 ssp_rng_canonical_float(rng), 278 ssp_rng_canonical_float(rng), 279 &prim, rwalk_hit->uv)); 280 f2_set(st, rwalk_hit->uv); 281 #endif 282 /* Map the sampled primitive id from the enclosure space to the scene 283 * space. Note that the overall scene has only one shape. As a consequence 284 * neither the geom_id nor the inst_id needs to be updated */ 285 rwalk_hit->prim.prim_id = enclosure_local2global_prim_id(enc, prim.prim_id); 286 287 SXD(primitive_get_attrib(&rwalk_hit->prim, SXD_POSITION, st, &attr_P)); 288 SXD(primitive_get_attrib(&rwalk_hit->prim, SXD_GEOMETRY_NORMAL, st, &attr_N)); 289 dX_set_fX(rwalk->vtx.P, attr_P.value); 290 fX(set)(rwalk_hit->normal, attr_N.value); 291 292 /* Define the interface side */ 293 scene_get_enclosure_ids(scn, rwalk_hit->prim.prim_id, enc_ids); 294 if(rwalk->enc_id == enc_ids[SDIS_BACK]) { 295 rwalk->hit_side = SDIS_BACK; 296 } else if(rwalk->enc_id == enc_ids[SDIS_FRONT]) { 297 rwalk->hit_side = SDIS_FRONT; 298 } else { 299 FATAL("Unexpected fluid interface\n"); 300 } 301 302 /* Get the interface properties */ 303 interf = scene_get_interface(scn, rwalk_hit->prim.prim_id); 304 305 /* Register the new vertex against the heat path */ 306 res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value, 307 SDIS_HEAT_VERTEX_CONVECTION, (int)ctx->nbranchings); 308 if(res != RES_OK) goto error; 309 310 /* Setup the fragment of the sampled position into the enclosure. */ 311 XD(setup_interface_fragment)(&frag, &rwalk->vtx, rwalk_hit, rwalk->hit_side); 312 313 /* Fetch the convection coefficient of the sampled position */ 314 hc = interface_get_convection_coef(interf, &frag); 315 if(hc > enc->hc_upper_bound) { 316 log_err(scn->dev, 317 "%s: hc (%g) exceeds its provided upper bound (%g) at %g %g %g.\n", 318 FUNC_NAME, hc, enc->hc_upper_bound, SPLIT3(rwalk->vtx.P)); 319 res = RES_BAD_OP; 320 goto error; 321 } 322 323 r = ssp_rng_canonical_float(rng); 324 if(r < hc / enc->hc_upper_bound) { 325 /* True convection. Always true if hc == bound. */ 326 break; 327 } 328 } 329 330 rwalk_hit->distance = 0; 331 T->func = XD(boundary_path); 332 rwalk->enc_id = ENCLOSURE_ID_NULL; /* Interface between 2 enclosures */ 333 334 exit: 335 return res; 336 error: 337 goto exit; 338 } 339 340 #include "sdis_Xd_end.h"