test_sdis_convection.c (12641B)
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 H 10.0 63 #define RHO 25.0 64 #define CP 2.0 65 66 /******************************************************************************* 67 * Media 68 ******************************************************************************/ 69 static double 70 fluid_get_temperature 71 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* is_stationary) 72 { 73 CHK(vtx != NULL); 74 if(*((int*)sdis_data_cget(is_stationary))) { 75 return SDIS_TEMPERATURE_NONE; 76 } else { 77 return vtx->time <= 0 ? Tf_0 : SDIS_TEMPERATURE_NONE; 78 } 79 } 80 81 static double 82 fluid_get_volumic_mass 83 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* is_stationary) 84 { 85 (void)is_stationary; 86 CHK(vtx != NULL); 87 return RHO; 88 } 89 90 static double 91 fluid_get_calorific_capacity 92 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* is_stationary) 93 { 94 (void)is_stationary; 95 CHK(vtx != NULL); 96 return CP; 97 } 98 99 /******************************************************************************* 100 * Interface 101 ******************************************************************************/ 102 struct interf { 103 double temperature; 104 }; 105 106 static double 107 interface_get_temperature 108 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 109 { 110 const struct interf* interf = sdis_data_cget(data); 111 CHK(frag && data); 112 return interf->temperature; 113 } 114 115 static double 116 interface_get_convection_coef 117 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 118 { 119 CHK(frag && data); 120 return H; 121 } 122 123 static double 124 interface_get_emissivity 125 (const struct sdis_interface_fragment* frag, 126 const unsigned source_id, 127 struct sdis_data* data) 128 { 129 (void)source_id; 130 CHK(frag && data); 131 return 0; 132 } 133 134 static double 135 interface_get_specular_fraction 136 (const struct sdis_interface_fragment* frag, 137 const unsigned source_id, 138 struct sdis_data* data) 139 { 140 (void)source_id; 141 CHK(frag && data); 142 return 0; 143 } 144 145 static struct sdis_interface* 146 create_interface 147 (struct sdis_device* dev, 148 struct sdis_medium* front, 149 struct sdis_medium* back, 150 const struct sdis_interface_shader* interf_shader, 151 const double temperature) 152 { 153 struct sdis_data* data; 154 struct sdis_interface* interf; 155 struct interf* interf_props; 156 157 OK(sdis_data_create 158 (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data)); 159 interf_props = sdis_data_get(data); 160 interf_props->temperature = temperature; 161 OK(sdis_interface_create 162 (dev, front, back, interf_shader, data, &interf)); 163 OK(sdis_data_ref_put(data)); 164 return interf; 165 } 166 167 /******************************************************************************* 168 * Test 169 ******************************************************************************/ 170 int 171 main(int argc, char** argv) 172 { 173 struct sdis_mc T = SDIS_MC_NULL; 174 struct sdis_mc mc_time = SDIS_MC_NULL; 175 struct sdis_device* dev = NULL; 176 struct sdis_medium* fluid = NULL; 177 struct sdis_medium* solid = NULL; 178 struct sdis_data* is_stationary = NULL; 179 struct sdis_interface* interf_T0 = NULL; 180 struct sdis_interface* interf_T1 = NULL; 181 struct sdis_interface* interf_T2 = NULL; 182 struct sdis_interface* interf_T3 = NULL; 183 struct sdis_interface* interf_T4 = NULL; 184 struct sdis_interface* interf_T5 = NULL; 185 struct sdis_scene* box_scn = NULL; 186 struct sdis_scene* square_scn = NULL; 187 struct sdis_estimator* estimator = NULL; 188 struct sdis_estimator* estimator2 = NULL; 189 struct sdis_green_function* green = NULL; 190 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 191 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 192 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 193 struct sdis_interface_shader interf_shader = DUMMY_INTERFACE_SHADER; 194 struct sdis_interface* box_interfaces[12/*#triangles*/]; 195 struct sdis_interface* square_interfaces[4/*#segments*/]; 196 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 197 double ref; 198 double Tinf; 199 double nu; 200 size_t nreals; 201 size_t nfails; 202 int i; 203 (void)argc, (void)argv; 204 205 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 206 207 /* Create the fluid medium */ 208 OK(sdis_data_create(dev, sizeof(int), ALIGNOF(int), NULL, &is_stationary)); 209 *((int*)sdis_data_get(is_stationary)) = 0; 210 fluid_shader.temperature = fluid_get_temperature; 211 fluid_shader.calorific_capacity = fluid_get_calorific_capacity; 212 fluid_shader.volumic_mass = fluid_get_volumic_mass; 213 OK(sdis_fluid_create(dev, &fluid_shader, is_stationary, &fluid)); 214 215 /* Create the solid_medium */ 216 OK(sdis_solid_create(dev, &solid_shader, NULL, &solid)); 217 218 /* Setup the interface shader */ 219 interf_shader.convection_coef = interface_get_convection_coef; 220 interf_shader.front.temperature = interface_get_temperature; 221 interf_shader.front.emissivity = interface_get_emissivity; 222 interf_shader.front.specular_fraction = interface_get_specular_fraction; 223 interf_shader.convection_coef_upper_bound = H; 224 225 /* Create the interfaces */ 226 interf_T0 = create_interface(dev, fluid, solid, &interf_shader, T0); 227 interf_T1 = create_interface(dev, fluid, solid, &interf_shader, T1); 228 interf_T2 = create_interface(dev, fluid, solid, &interf_shader, T2); 229 interf_T3 = create_interface(dev, fluid, solid, &interf_shader, T3); 230 interf_T4 = create_interface(dev, fluid, solid, &interf_shader, T4); 231 interf_T5 = create_interface(dev, fluid, solid, &interf_shader, T5); 232 233 /* Release the media */ 234 OK(sdis_medium_ref_put(solid)); 235 OK(sdis_medium_ref_put(fluid)); 236 237 /* Map the interfaces to their box triangles */ 238 box_interfaces[0] = box_interfaces[1] = interf_T5; /* Front */ 239 box_interfaces[2] = box_interfaces[3] = interf_T0; /* Left */ 240 box_interfaces[4] = box_interfaces[5] = interf_T4; /* Back */ 241 box_interfaces[6] = box_interfaces[7] = interf_T1; /* Right */ 242 box_interfaces[8] = box_interfaces[9] = interf_T3; /* Top */ 243 box_interfaces[10]= box_interfaces[11]= interf_T2; /* Bottom */ 244 245 /* Map the interfaces to their square segments */ 246 square_interfaces[0] = interf_T2; /* Bottom */ 247 square_interfaces[1] = interf_T0; /* Left */ 248 square_interfaces[2] = interf_T3; /* Top */ 249 square_interfaces[3] = interf_T1; /* Right */ 250 251 /* Create the box scene */ 252 scn_args.get_indices = box_get_indices; 253 scn_args.get_interface = box_get_interface; 254 scn_args.get_position = box_get_position; 255 scn_args.nprimitives = box_ntriangles; 256 scn_args.nvertices = box_nvertices; 257 scn_args.context = box_interfaces; 258 OK(sdis_scene_create(dev, &scn_args, &box_scn)); 259 260 /* Create the square scene */ 261 scn_args.get_indices = square_get_indices; 262 scn_args.get_interface = square_get_interface; 263 scn_args.get_position = square_get_position; 264 scn_args.nprimitives = square_nsegments; 265 scn_args.nvertices = square_nvertices; 266 scn_args.context = square_interfaces; 267 OK(sdis_scene_2d_create(dev, &scn_args, &square_scn)); 268 269 /* Release the interfaces */ 270 OK(sdis_interface_ref_put(interf_T0)); 271 OK(sdis_interface_ref_put(interf_T1)); 272 OK(sdis_interface_ref_put(interf_T2)); 273 OK(sdis_interface_ref_put(interf_T3)); 274 OK(sdis_interface_ref_put(interf_T4)); 275 OK(sdis_interface_ref_put(interf_T5)); 276 277 solve_args.nrealisations = N; 278 solve_args.position[0] = 0.25; 279 solve_args.position[1] = 0.25; 280 solve_args.position[2] = 0.25; 281 282 /* Test in 3D for various time values. */ 283 nu = (6 * H) / (RHO*CP); 284 Tinf = (H*(T0 + T1 + T2 + T3 + T4 + T5)) / (6 * H); 285 printf(">>> Temperature of the box at (%g %g %g)\n\n", 286 SPLIT3(solve_args.position)); 287 FOR_EACH(i, 0, 5) { 288 double time = i ? (double)i / nu : INF; 289 ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); 290 291 solve_args.time_range[0] = time; 292 solve_args.time_range[1] = time; 293 294 /* Setup stationary state */ 295 *((int*)sdis_data_get(is_stationary)) = IS_INF(time); 296 297 /* Solve in 3D */ 298 OK(sdis_solve_probe(box_scn, &solve_args, &estimator)); 299 OK(sdis_estimator_get_temperature(estimator, &T)); 300 OK(sdis_estimator_get_realisation_time(estimator, &mc_time)); 301 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 302 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 303 CHK(nfails + nreals == N); 304 printf("Temperature at %g = %g ~ %g +/- %g\n", time, ref, T.E, T.SE); 305 printf("Time per realisation (in usec) = %g +/- %g\n", mc_time.E, mc_time.SE); 306 if(nfails) 307 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 308 CHK(eq_eps(T.E, ref, T.SE * 3)); 309 310 OK(sdis_solve_probe_green_function(box_scn, &solve_args, &green)); 311 OK(sdis_green_function_solve(green, &estimator2)); 312 check_green_function(green); 313 check_estimator_eq(estimator, estimator2); 314 check_green_serialization(green, box_scn); 315 OK(sdis_estimator_ref_put(estimator2)); 316 OK(sdis_green_function_ref_put(green)); 317 318 OK(sdis_estimator_ref_put(estimator)); 319 printf("\n"); 320 } 321 322 /* Test in 2D for various time values. */ 323 nu = (4 * H) / (RHO*CP); 324 Tinf = (H * (T0 + T1 + T2 + T3)) / (4 * H); 325 printf(">>> Temperature of the square at (%g %g)\n\n", 326 SPLIT2(solve_args.position)); 327 FOR_EACH(i, 0, 5) { 328 double time = i ? (double)i / nu : INF; 329 ref = Tf_0 * exp(-nu * time) + Tinf * (1 - exp(-nu * time)); 330 331 solve_args.time_range[0] = time; 332 solve_args.time_range[1] = time; 333 334 /* Setup stationnary state */ 335 *((int*)sdis_data_get(is_stationary)) = IS_INF(time); 336 337 /* Solve in 2D */ 338 OK(sdis_solve_probe(square_scn, &solve_args, &estimator)); 339 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 340 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 341 CHK(nfails + nreals == N); 342 OK(sdis_estimator_get_temperature(estimator, &T)); 343 OK(sdis_estimator_get_realisation_time(estimator, &mc_time)); 344 printf("Temperature at %g = %g ~ %g +/- %g\n", time, ref, T.E, T.SE); 345 printf("Time per realisation (in usec) = %g +/- %g\n", mc_time.E, mc_time.SE); 346 if(nfails) 347 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 348 CHK(eq_eps(T.E, ref, T.SE * 3)); 349 350 OK(sdis_solve_probe_green_function(square_scn, &solve_args, &green)); 351 OK(sdis_green_function_solve(green, &estimator2)); 352 check_green_function(green); 353 check_estimator_eq(estimator, estimator2); 354 check_green_serialization(green, square_scn); 355 OK(sdis_estimator_ref_put(estimator2)); 356 OK(sdis_green_function_ref_put(green)); 357 358 OK(sdis_estimator_ref_put(estimator)); 359 printf("\n"); 360 } 361 362 OK(sdis_scene_ref_put(box_scn)); 363 OK(sdis_scene_ref_put(square_scn)); 364 OK(sdis_device_ref_put(dev)); 365 OK(sdis_data_ref_put(is_stationary)); 366 367 CHK(mem_allocated_size() == 0); 368 return 0; 369 } 370