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