test_sdis_unsteady_analytic_profile_2d.c (13300B)
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/mem_allocator.h> 20 #include <rsys/stretchy_array.h> 21 22 /* 23 * The system is an unsteady-state temperature profile, meaning that at any 24 * point, at any: time, we can analytically calculate the temperature. We 25 * immerse in this temperature field a supershape representing a solid in which 26 * we want to evaluate the temperature by Monte Carlo at a given position and 27 * observation time. On the Monte Carlo side, the temperature of the supershape 28 * is unknown. Only the boundary temperature is fixed at the temperature of the 29 * unsteady trilinear profile mentioned above. Monte Carlo would then have to 30 * find the temperature defined by the unsteady profile. 31 * 32 * T(y) /\ <-- T(x,y,t) 33 * | ___/ \___ 34 * | \ . T=? / 35 * o--- T(x) /_ __ _\ 36 * \/ \/ 37 * 38 */ 39 40 #define LAMBDA 0.1 /* [W/(m.K)] */ 41 #define RHO 25.0 /* [kg/m^3] */ 42 #define CP 2.0 /* [J/K/kg)] */ 43 #define DELTA 1.0/20.0 /* [m/fp_to_meter] */ 44 45 #define NREALISATIONS 100000 46 47 struct super_shape { 48 double* positions; 49 size_t* indices; 50 }; 51 #define SUPER_SHAPE_NULL__ {NULL, NULL} 52 static const struct super_shape SUPER_SHAPE_NULL = SUPER_SHAPE_NULL__; 53 54 /******************************************************************************* 55 * Helper functions 56 ******************************************************************************/ 57 /* Unsteady-state temperature profile. Its corresponding power term is : 58 * P(x, y) = A*[12*x^2*y^3 + 5*x^4*y] */ 59 static double temperature(const double pos[2], const double time) 60 { 61 const double kx = PI/4; 62 const double ky = PI/4; 63 const double alpha = LAMBDA / (RHO*CP); /* Diffusivity */ 64 65 const double A = 0; /* No termal source */ 66 const double B1 = 0; 67 const double B2 = 1000; 68 69 double x, y, t; 70 double a, b, c; 71 double temp; 72 ASSERT(pos); 73 74 x = pos[0]; 75 y = pos[1]; 76 t = time; 77 78 a = B1*(x*x*x - 3*x*y*y); 79 b = B2*sin(kx*x)*sin(ky*y)*exp(-alpha*(kx*kx + ky*ky)*t); 80 c = A * x*x*x*x + y*y*y; 81 82 temp = (a + b - c) / LAMBDA; 83 return temp; 84 } 85 86 static void 87 dump_paths 88 (FILE* fp, 89 struct sdis_scene* scn, 90 const enum sdis_diffusion_algorithm diff_algo, 91 const double pos[2], 92 const double time, 93 const size_t npaths) 94 { 95 struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 96 struct sdis_estimator* estimator = NULL; 97 98 args.nrealisations = npaths; 99 args.position[0] = pos[0]; 100 args.position[1] = pos[1]; 101 args.time_range[0] = time; 102 args.time_range[1] = time; 103 args.diff_algo = diff_algo; 104 args.register_paths = SDIS_HEAT_PATH_ALL; 105 OK(sdis_solve_probe(scn, &args, &estimator)); 106 107 dump_heat_paths(fp, estimator); 108 109 OK(sdis_estimator_ref_put(estimator)); 110 } 111 112 /******************************************************************************* 113 * Solid, i.e. medium of the super shape 114 ******************************************************************************/ 115 #define SOLID_PROP(Prop, Val) \ 116 static double \ 117 solid_get_##Prop \ 118 (const struct sdis_rwalk_vertex* vtx, \ 119 struct sdis_data* data) \ 120 { \ 121 (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \ 122 return Val; \ 123 } 124 SOLID_PROP(calorific_capacity, CP) /* [J/K/kg] */ 125 SOLID_PROP(thermal_conductivity, LAMBDA) /* [W/m/K] */ 126 SOLID_PROP(volumic_mass, RHO) /* [kg/m^3] */ 127 SOLID_PROP(delta, DELTA) /* [m/fp_to_meter] */ 128 SOLID_PROP(temperature, SDIS_TEMPERATURE_NONE/*<=> unknown*/) /* [K] */ 129 #undef SOLID_PROP 130 131 static struct sdis_medium* 132 create_solid(struct sdis_device* sdis) 133 { 134 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 135 struct sdis_medium* solid = NULL; 136 137 shader.calorific_capacity = solid_get_calorific_capacity; 138 shader.thermal_conductivity = solid_get_thermal_conductivity; 139 shader.volumic_mass = solid_get_volumic_mass; 140 shader.delta = solid_get_delta; 141 shader.temperature = solid_get_temperature; 142 shader.t0 = -INF; 143 OK(sdis_solid_create(sdis, &shader, NULL, &solid)); 144 return solid; 145 } 146 147 /******************************************************************************* 148 * Dummy environment, i.e. environment surrounding the super shape. It is 149 * defined only for Stardis compliance: in Stardis, an interface must divide 2 150 * media. 151 ******************************************************************************/ 152 static struct sdis_medium* 153 create_dummy(struct sdis_device* sdis) 154 { 155 struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL; 156 struct sdis_medium* dummy = NULL; 157 158 shader.calorific_capacity = dummy_medium_getter; 159 shader.volumic_mass = dummy_medium_getter; 160 shader.temperature = dummy_medium_getter; 161 OK(sdis_fluid_create(sdis, &shader, NULL, &dummy)); 162 return dummy; 163 } 164 165 /******************************************************************************* 166 * Interface: its temperature is fixed to the unsteady-state profile 167 ******************************************************************************/ 168 static double 169 interface_get_temperature 170 (const struct sdis_interface_fragment* frag, 171 struct sdis_data* data) 172 { 173 (void)data; 174 return temperature(frag->P, frag->time); 175 } 176 177 static struct sdis_interface* 178 create_interface 179 (struct sdis_device* sdis, 180 struct sdis_medium* front, 181 struct sdis_medium* back) 182 { 183 struct sdis_interface* interf = NULL; 184 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 185 186 shader.front.temperature = interface_get_temperature; 187 shader.back.temperature = interface_get_temperature; 188 OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf)); 189 return interf; 190 } 191 192 /******************************************************************************* 193 * Super shape 194 ******************************************************************************/ 195 static struct super_shape 196 create_super_shape(void) 197 { 198 struct super_shape sshape = SUPER_SHAPE_NULL; 199 200 const unsigned nslices = 128; 201 const double a = 1.0; 202 const double b = 1.0; 203 const double n1 = 1.0; 204 const double n2 = 1.0; 205 const double n3 = 1.0; 206 const double m = 6.0; 207 size_t i = 0; 208 209 FOR_EACH(i, 0, nslices) { 210 const double theta = (double)i * (2.0*PI / (double)nslices); 211 const double tmp0 = pow(fabs(1.0/a * cos(m/4.0*theta)), n2); 212 const double tmp1 = pow(fabs(1.0/b * sin(m/4.0*theta)), n3); 213 const double tmp2 = pow(tmp0 + tmp1, 1.0/n1); 214 const double r = 1.0 / tmp2; 215 const double x = cos(theta) * r; 216 const double y = sin(theta) * r; 217 const size_t j = (i + 1) % nslices; 218 219 sa_push(sshape.positions, x); 220 sa_push(sshape.positions, y); 221 sa_push(sshape.indices, i); 222 sa_push(sshape.indices, j); 223 } 224 225 return sshape; 226 } 227 228 static void 229 release_super_shape(struct super_shape* sshape) 230 { 231 CHK(sshape); 232 sa_release(sshape->positions); 233 sa_release(sshape->indices); 234 } 235 236 static INLINE size_t 237 super_shape_nsegments(const struct super_shape* sshape) 238 { 239 CHK(sshape); 240 return sa_size(sshape->indices) / 2/*#indices per segment*/; 241 } 242 243 static INLINE size_t 244 super_shape_nvertices(const struct super_shape* sshape) 245 { 246 CHK(sshape); 247 return sa_size(sshape->positions) / 2/*#coords per vertex*/; 248 } 249 250 /******************************************************************************* 251 * Scene, i.e. the system to simulate 252 ******************************************************************************/ 253 struct scene_context { 254 const struct super_shape* sshape; 255 struct sdis_interface* interf; 256 }; 257 static const struct scene_context SCENE_CONTEXT_NULL = {NULL, NULL}; 258 259 static void 260 scene_get_indices(const size_t iseg, size_t ids[2], void* ctx) 261 { 262 struct scene_context* context = ctx; 263 CHK(ids && context); 264 /* Flip the indices to ensure that the normal points into the super shape */ 265 ids[0] = (unsigned)context->sshape->indices[iseg*2+1]; 266 ids[1] = (unsigned)context->sshape->indices[iseg*2+0]; 267 } 268 269 static void 270 scene_get_interface(const size_t iseg, struct sdis_interface** interf, void* ctx) 271 { 272 struct scene_context* context = ctx; 273 (void)iseg; 274 CHK(interf && context); 275 *interf = context->interf; 276 } 277 278 static void 279 scene_get_position(const size_t ivert, double pos[2], void* ctx) 280 { 281 struct scene_context* context = ctx; 282 CHK(pos && context); 283 pos[0] = context->sshape->positions[ivert*2+0]; 284 pos[1] = context->sshape->positions[ivert*2+1]; 285 } 286 287 static struct sdis_scene* 288 create_scene 289 (struct sdis_device* sdis, 290 const struct super_shape* sshape, 291 struct sdis_interface* interf) 292 { 293 struct sdis_scene* scn = NULL; 294 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 295 struct scene_context context = SCENE_CONTEXT_NULL; 296 297 context.interf = interf; 298 context.sshape = sshape; 299 300 scn_args.get_indices = scene_get_indices; 301 scn_args.get_interface = scene_get_interface; 302 scn_args.get_position = scene_get_position; 303 scn_args.nprimitives = super_shape_nsegments(sshape); 304 scn_args.nvertices = super_shape_nvertices(sshape); 305 scn_args.context = &context; 306 OK(sdis_scene_2d_create(sdis, &scn_args, &scn)); 307 return scn; 308 } 309 310 /******************************************************************************* 311 * Validations 312 ******************************************************************************/ 313 static void 314 check_probe 315 (struct sdis_scene* scn, 316 const enum sdis_diffusion_algorithm diff_algo, 317 const double pos[2], 318 const double time, /* [s] */ 319 const int green) 320 { 321 struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 322 struct sdis_mc T = SDIS_MC_NULL; 323 struct sdis_estimator* estimator = NULL; 324 double ref = 0; 325 326 args.nrealisations = NREALISATIONS; 327 args.position[0] = pos[0]; 328 args.position[1] = pos[1]; 329 args.time_range[0] = time; 330 args.time_range[1] = time; 331 args.diff_algo = diff_algo; 332 333 if(!green) { 334 OK(sdis_solve_probe(scn, &args, &estimator)); 335 } else { 336 struct sdis_green_function* greenfn = NULL; 337 338 OK(sdis_solve_probe_green_function(scn, &args, &greenfn)); 339 OK(sdis_green_function_solve(greenfn, &estimator)); 340 OK(sdis_green_function_ref_put(greenfn)); 341 } 342 OK(sdis_estimator_get_temperature(estimator, &T)); 343 344 ref = temperature(pos, time); 345 printf("T(%g, %g, %g) = %g ~ %g +/- %g\n", 346 SPLIT2(pos), time, ref, T.E, T.SE); 347 CHK(eq_eps(ref, T.E, 3*T.SE)); 348 349 OK(sdis_estimator_ref_put(estimator)); 350 } 351 352 /******************************************************************************* 353 * The test 354 ******************************************************************************/ 355 int 356 main(int argc, char** argv) 357 { 358 /* Stardis */ 359 struct sdis_device* sdis = NULL; 360 struct sdis_interface* interf = NULL; 361 struct sdis_medium* solid = NULL; 362 struct sdis_medium* dummy = NULL; /* Medium surrounding the solid */ 363 struct sdis_scene* scn = NULL; 364 365 /* Miscellaneous */ 366 FILE* fp = NULL; 367 struct super_shape sshape = SUPER_SHAPE_NULL; 368 const double pos[3] = {0.2,0.3}; /* [m/fp_to_meter] */ 369 const double time = 5; /* [s] */ 370 (void)argc, (void)argv; 371 372 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &sdis)); 373 374 sshape = create_super_shape(); 375 376 /* Save the super shape geometry for debug and visualisation */ 377 CHK(fp = fopen("super_shape_2d.obj", "w")); 378 dump_segments(fp, sshape.positions, super_shape_nvertices(&sshape), 379 sshape.indices, super_shape_nsegments(&sshape)); 380 CHK(fclose(fp) == 0); 381 382 solid = create_solid(sdis); 383 dummy = create_dummy(sdis); 384 interf = create_interface(sdis, solid, dummy); 385 scn = create_scene(sdis, &sshape, interf); 386 387 check_probe(scn, SDIS_DIFFUSION_DELTA_SPHERE, pos, time, 0/*green*/); 388 check_probe(scn, SDIS_DIFFUSION_WOS, pos, time, 0/*green*/); 389 check_probe(scn, SDIS_DIFFUSION_WOS, pos, time, 1/*green*/); 390 391 /* Write 10 heat paths sampled by the delta sphere algorithm */ 392 CHK(fp = fopen("paths_delta_sphere_2d.vtk", "w")); 393 dump_paths(fp, scn, SDIS_DIFFUSION_DELTA_SPHERE, pos, time, 10); 394 CHK(fclose(fp) == 0); 395 396 /* Write 10 heat paths sampled by the WoS algorithm */ 397 CHK(fp = fopen("paths_wos_2d.vtk", "w")); 398 dump_paths(fp, scn, SDIS_DIFFUSION_WOS, pos, time, 10); 399 CHK(fclose(fp) == 0); 400 401 release_super_shape(&sshape); 402 OK(sdis_device_ref_put(sdis)); 403 OK(sdis_interface_ref_put(interf)); 404 OK(sdis_medium_ref_put(solid)); 405 OK(sdis_medium_ref_put(dummy)); 406 OK(sdis_scene_ref_put(scn)); 407 408 CHK(mem_allocated_size() == 0); 409 return 0; 410 }