test_sdis_convection_non_uniform.c (13189B)
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.h" 17 #include "test_sdis_utils.h" 18 19 #include <rsys/double3.h> 20 #include <rsys/math.h> 21 22 /* 23 * The scene is composed of an unit fluid cube/square whose temperature is 24 * unknown. The convection coefficient with the surrounding solid is H 25 * everywhere the temperature of the -/+X, -/+Y and -/+Z faces are fixed to T0 26 * and T1, T2, T3, T4 and T5, respectively. This test computes the temperature 27 * of the fluid Tf at an observation time t. This temperature is equal to: 28 * 29 * Tf(t) = Tf(0) * e^(-nu*t) + Tinf*(1-e^(-nu*t)) 30 * 31 * nu = (Sum_{i=0..5}(H*Si)) / (RHO*CP*V) 32 * Tinf = (Sum_{i=0..5}(H*Si*Ti) / (Sum_{i=0..5}(H*Si)); 33 * 34 * with Si surface of the faces (i.e. one), V the volume of the cube (i.e. 35 * one), RHO the volumic mass of the fluid and CP its calorific capacity. 36 * 37 * 3D 2D 38 * 39 * (1,1,1) (1,1) 40 * +---------+ +-----T3----+ 41 * /' T3 /|T4 | | 42 * +---------+ | | H _\ | 43 * | ' H _\ |T1 T0 / / T1 44 * |T0 / / | | | \__/ | 45 * | +..\__/.|.+ | | 46 * T5|, T2 |/ +-----------+ 47 * +---------+ (0,0) 48 * (0,0,0) 49 */ 50 51 #define N 100000 /* #realisations */ 52 53 #define Tf_0 280.0 54 55 #define T0 300.0 56 #define T1 310.0 57 #define T2 320.0 58 #define T3 330.0 59 #define T4 340.0 60 #define T5 350.0 61 62 #define HC0 100.0 63 #define HC1 30.0 64 #define HC2 3020.0 65 #define HC3 7300.0 66 #define HC4 3400.0 67 #define HC5 50.0 68 69 #define RHO 25.0 70 #define CP 2.0 71 72 /******************************************************************************* 73 * Media 74 ******************************************************************************/ 75 static double 76 fluid_get_temperature 77 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* is_stationary) 78 { 79 CHK(vtx != NULL); 80 CHK(is_stationary != NULL); 81 if(*((int*)sdis_data_cget(is_stationary))) { 82 return SDIS_TEMPERATURE_NONE; 83 } else { 84 return vtx->time <= 0 ? Tf_0 : SDIS_TEMPERATURE_NONE; 85 } 86 } 87 88 static double 89 fluid_get_volumic_mass 90 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 91 { 92 (void)data; 93 CHK(vtx != NULL); 94 return RHO; 95 } 96 97 static double 98 fluid_get_calorific_capacity 99 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 100 { 101 (void)data; 102 CHK(vtx != NULL); 103 return CP; 104 } 105 106 /******************************************************************************* 107 * Interface 108 ******************************************************************************/ 109 struct interf { 110 double temperature; 111 double hc; 112 }; 113 114 static double 115 interface_get_temperature 116 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 117 { 118 const struct interf* interf = sdis_data_cget(data); 119 CHK(frag && data); 120 return interf->temperature; 121 } 122 123 static double 124 interface_get_convection_coef 125 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 126 { 127 const struct interf* interf = sdis_data_cget(data); 128 CHK(frag && data); 129 return interf->hc; 130 } 131 132 static double 133 interface_get_emissivity 134 (const struct sdis_interface_fragment* frag, 135 const unsigned source_id, 136 struct sdis_data* data) 137 { 138 (void)source_id; 139 CHK(frag && data); 140 return 0; 141 } 142 143 static double 144 interface_get_specular_fraction 145 (const struct sdis_interface_fragment* frag, 146 const unsigned source_id, 147 struct sdis_data* data) 148 { 149 (void)source_id; 150 CHK(frag && data); 151 return 0; 152 } 153 154 static struct sdis_interface* 155 create_interface 156 (struct sdis_device* dev, 157 struct sdis_medium* front, 158 struct sdis_medium* back, 159 const struct sdis_interface_shader* interf_shader, 160 const double temperature, 161 const double hc) 162 { 163 struct sdis_data* data; 164 struct sdis_interface* interf; 165 struct interf* interf_props; 166 167 OK(sdis_data_create(dev, sizeof(struct interf), ALIGNOF(struct interf), 168 NULL, &data)); 169 interf_props = sdis_data_get(data); 170 interf_props->temperature = temperature; 171 interf_props->hc = hc; 172 OK(sdis_interface_create(dev, front, back, interf_shader, data, &interf)); 173 OK(sdis_data_ref_put(data)); 174 return interf; 175 } 176 177 /******************************************************************************* 178 * Test 179 ******************************************************************************/ 180 int 181 main(int argc, char** argv) 182 { 183 struct sdis_mc T = SDIS_MC_NULL; 184 struct sdis_mc mc_time = SDIS_MC_NULL; 185 struct sdis_device* dev = NULL; 186 struct sdis_medium* fluid = NULL; 187 struct sdis_medium* solid = NULL; 188 struct sdis_data* is_stationary = NULL; 189 struct sdis_interface* interf_T0 = NULL; 190 struct sdis_interface* interf_T1 = NULL; 191 struct sdis_interface* interf_T2 = NULL; 192 struct sdis_interface* interf_T3 = NULL; 193 struct sdis_interface* interf_T4 = NULL; 194 struct sdis_interface* interf_T5 = NULL; 195 struct sdis_scene* box_scn = NULL; 196 struct sdis_scene* square_scn = NULL; 197 struct sdis_estimator* estimator = NULL; 198 struct sdis_estimator* estimator2 = NULL; 199 struct sdis_green_function* green = NULL; 200 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 201 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 202 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 203 struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER; 204 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 205 struct sdis_interface* box_interfaces[12/*#triangles*/]; 206 struct sdis_interface* square_interfaces[4/*#segments*/]; 207 double ref; 208 double Tinf; 209 double nu; 210 size_t nreals; 211 size_t nfails; 212 int i; 213 (void)argc, (void)argv; 214 215 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 216 217 OK(sdis_data_create(dev, sizeof(int), ALIGNOF(int), NULL, &is_stationary)); 218 *((int*)sdis_data_get(is_stationary)) = 0; 219 220 /* Create the fluid medium */ 221 fluid_shader.temperature = fluid_get_temperature; 222 fluid_shader.calorific_capacity = fluid_get_calorific_capacity; 223 fluid_shader.volumic_mass = fluid_get_volumic_mass; 224 OK(sdis_fluid_create(dev, &fluid_shader, is_stationary, &fluid)); 225 226 /* Create the solid_medium */ 227 OK(sdis_solid_create(dev, &solid_shader, NULL, &solid)); 228 229 /* Setup the interface shader */ 230 interf_shader.convection_coef = interface_get_convection_coef; 231 interf_shader.front.temperature = interface_get_temperature; 232 interf_shader.front.emissivity = interface_get_emissivity; 233 interf_shader.front.specular_fraction = interface_get_specular_fraction; 234 235 /* Create the interfaces */ 236 interf_shader.convection_coef_upper_bound = HC0; 237 interf_T0 = create_interface(dev, fluid, solid, &interf_shader, T0, HC0); 238 interf_shader.convection_coef_upper_bound = HC1; 239 interf_T1 = create_interface(dev, fluid, solid, &interf_shader, T1, HC1); 240 interf_shader.convection_coef_upper_bound = HC2; 241 interf_T2 = create_interface(dev, fluid, solid, &interf_shader, T2, HC2); 242 interf_shader.convection_coef_upper_bound = HC3; 243 interf_T3 = create_interface(dev, fluid, solid, &interf_shader, T3, HC3); 244 interf_shader.convection_coef_upper_bound = HC4; 245 interf_T4 = create_interface(dev, fluid, solid, &interf_shader, T4, HC4); 246 interf_shader.convection_coef_upper_bound = HC5; 247 interf_T5 = create_interface(dev, fluid, solid, &interf_shader, T5, HC5); 248 249 /* Release the media */ 250 OK(sdis_medium_ref_put(solid)); 251 OK(sdis_medium_ref_put(fluid)); 252 253 /* Map the interfaces to their box triangles */ 254 box_interfaces[0] = box_interfaces[1] = interf_T5; /* Front */ 255 box_interfaces[2] = box_interfaces[3] = interf_T0; /* Left */ 256 box_interfaces[4] = box_interfaces[5] = interf_T4; /* Back */ 257 box_interfaces[6] = box_interfaces[7] = interf_T1; /* Right */ 258 box_interfaces[8] = box_interfaces[9] = interf_T3; /* Top */ 259 box_interfaces[10]= box_interfaces[11]= interf_T2; /* Bottom */ 260 261 /* Map the interfaces to their square segments */ 262 square_interfaces[0] = interf_T2; /* Bottom */ 263 square_interfaces[1] = interf_T0; /* Left */ 264 square_interfaces[2] = interf_T3; /* Top */ 265 square_interfaces[3] = interf_T1; /* Right */ 266 267 /* Create the box scene */ 268 scn_args.get_indices = box_get_indices; 269 scn_args.get_interface = box_get_interface; 270 scn_args.get_position = box_get_position; 271 scn_args.nprimitives = box_ntriangles; 272 scn_args.nvertices = box_nvertices; 273 scn_args.context = box_interfaces; 274 OK(sdis_scene_create(dev, &scn_args, &box_scn)); 275 276 /* Create the square scene */ 277 scn_args.get_indices = square_get_indices; 278 scn_args.get_interface = square_get_interface; 279 scn_args.get_position = square_get_position; 280 scn_args.nprimitives = square_nsegments; 281 scn_args.nvertices = square_nvertices; 282 scn_args.context = square_interfaces; 283 OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); 284 285 /* Release the interfaces */ 286 OK(sdis_interface_ref_put(interf_T0)); 287 OK(sdis_interface_ref_put(interf_T1)); 288 OK(sdis_interface_ref_put(interf_T2)); 289 OK(sdis_interface_ref_put(interf_T3)); 290 OK(sdis_interface_ref_put(interf_T4)); 291 OK(sdis_interface_ref_put(interf_T5)); 292 293 solve_args.nrealisations = N; 294 solve_args.position[0] = 0.25; 295 solve_args.position[1] = 0.25; 296 solve_args.position[2] = 0.25; 297 298 /* Test in 3D for various time values. */ 299 nu = (HC0 + HC1 + HC2 + HC3 + HC4 + HC5) / (RHO * CP); 300 Tinf = (HC0 * T0 + HC1 * T1 + HC2 * T2 + HC3 * T3 + HC4 * T4 + HC5 * T5) 301 / (HC0 + HC1 + HC2 + HC3 + HC4 + HC5); 302 printf(">>> Temperature of the box at (%g %g %g)\n\n", 303 SPLIT3(solve_args.position)); 304 FOR_EACH(i, 0, 5) { 305 double time = i ? (double)i / nu : INF; 306 solve_args.time_range[0] = time; 307 solve_args.time_range[1] = time; 308 309 ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); 310 311 *((int*)sdis_data_get(is_stationary)) = IS_INF(time); 312 313 /* Solve in 3D */ 314 OK(sdis_solve_probe(box_scn, &solve_args, &estimator)); 315 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 316 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 317 CHK(nfails + nreals == N); 318 OK(sdis_estimator_get_temperature(estimator, &T)); 319 OK(sdis_estimator_get_realisation_time(estimator, &mc_time)); 320 printf("Temperature at %g = %g ~ %g +/- %g\n", time, ref, T.E, T.SE); 321 printf("Time per realisation (in usec) = %g +/- %g\n", mc_time.E, mc_time.SE); 322 if(nfails) 323 printf("#failures = %lu/%lu\n", (unsigned long)nfails,(unsigned long)N); 324 CHK(eq_eps(T.E, ref, T.SE * 3)); 325 326 OK(sdis_solve_probe_green_function(box_scn, &solve_args, &green)); 327 OK(sdis_green_function_solve(green, &estimator2)); 328 check_green_function(green); 329 check_estimator_eq(estimator, estimator2); 330 check_green_serialization(green, box_scn); 331 OK(sdis_estimator_ref_put(estimator2)); 332 OK(sdis_green_function_ref_put(green)); 333 334 OK(sdis_estimator_ref_put(estimator)); 335 printf("\n"); 336 } 337 338 /* Test in 2D for various time values. */ 339 nu = (HC0 + HC1 + HC2 + HC3) / (RHO * CP); 340 Tinf = (HC0 * T0 + HC1 * T1 + HC2 * T2 + HC3 * T3) / (HC0 + HC1 + HC2 + HC3); 341 printf(">>> Temperature of the square at (%g %g)\n\n", 342 SPLIT2(solve_args.position)); 343 FOR_EACH(i, 0, 5) { 344 double time = i ? (double)i / nu : INF; 345 346 solve_args.time_range[0] = time; 347 solve_args.time_range[1] = time; 348 349 ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); 350 351 *((int*)sdis_data_get(is_stationary)) = IS_INF(time); 352 353 OK(sdis_solve_probe(square_scn, &solve_args, &estimator)); 354 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 355 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 356 CHK(nfails + nreals == N); 357 OK(sdis_estimator_get_temperature(estimator, &T)); 358 OK(sdis_estimator_get_realisation_time(estimator, &mc_time)); 359 printf("Temperature at %g = %g ~ %g +/- %g\n", time, ref, T.E, T.SE); 360 printf("Time per realisation (in usec) = %g +/- %g\n", mc_time.E, mc_time.SE); 361 if(nfails) 362 printf("#failures = %lu/%lu\n", (unsigned long)nfails,(unsigned long)N); 363 CHK(eq_eps(T.E, ref, T.SE * 3)); 364 365 OK(sdis_solve_probe_green_function(square_scn, &solve_args, &green)); 366 OK(sdis_green_function_solve(green, &estimator2)); 367 check_green_function(green); 368 check_estimator_eq(estimator, estimator2); 369 check_green_serialization(green, square_scn); 370 OK(sdis_estimator_ref_put(estimator2)); 371 OK(sdis_green_function_ref_put(green)); 372 373 OK(sdis_estimator_ref_put(estimator)); 374 printf("\n"); 375 } 376 377 OK(sdis_scene_ref_put(box_scn)); 378 OK(sdis_scene_ref_put(square_scn)); 379 OK(sdis_device_ref_put(dev)); 380 OK(sdis_data_ref_put(is_stationary)); 381 382 CHK(mem_allocated_size() == 0); 383 return 0; 384 } 385