test_sdis_solve_probe3.c (11230B)
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/stretchy_array.h> 20 #include <rsys/math.h> 21 22 #include <star/s3dut.h> 23 24 /* 25 * The scene is composed of a solid cube whose temperature is unknown. The 26 * convection coefficient with the surrounding fluid is null. The temperature 27 * is fixed at the front and back face. At the center of the cube there is a 28 * solid sphere whose physical properties are the same of the solid cube; i.e. 29 * the sphere influences the random walks but not the result. 30 * 31 * (1,1,1) 32 * +----------------+ 33 * /' # # /| 34 * +----*--------*--+ | 35 * | ' # # | |350K 36 * | ' # # | | 37 * 300K| ' # # | | 38 * | +.....#..#.....|.+ 39 * |/ |/ 40 * +----------------+ 41 * (0,0,0) 42 */ 43 44 /******************************************************************************* 45 * Geometry 46 ******************************************************************************/ 47 struct context { 48 double* positions; 49 size_t* indices; 50 struct sdis_interface* solid_fluid_Tnone; 51 struct sdis_interface* solid_fluid_T300; 52 struct sdis_interface* solid_fluid_T350; 53 struct sdis_interface* solid_solid; 54 }; 55 static const struct context CONTEXT_NULL = {NULL, NULL, NULL, NULL, NULL, NULL}; 56 57 static void 58 get_indices(const size_t itri, size_t ids[3], void* context) 59 { 60 struct context* ctx = context; 61 ids[0] = ctx->indices[itri*3+0]; 62 ids[1] = ctx->indices[itri*3+1]; 63 ids[2] = ctx->indices[itri*3+2]; 64 } 65 66 static void 67 get_position(const size_t ivert, double pos[3], void* context) 68 { 69 struct context* ctx = context; 70 pos[0] = ctx->positions[ivert*3+0]; 71 pos[1] = ctx->positions[ivert*3+1]; 72 pos[2] = ctx->positions[ivert*3+2]; 73 } 74 75 static void 76 get_interface(const size_t itri, struct sdis_interface** bound, void* context) 77 { 78 struct context* ctx = context; 79 CHK(bound != NULL && context != NULL); 80 81 if(itri == 0 || itri == 1) { /* Box front face */ 82 *bound = ctx->solid_fluid_T300; 83 } else if(itri == 4 || itri == 5) { /* Box back face */ 84 *bound = ctx->solid_fluid_T350; 85 } else if(itri < box_ntriangles) { /* Box remaining faces */ 86 *bound = ctx->solid_fluid_Tnone; 87 } else { /* Faces of the internal geometry */ 88 *bound = ctx->solid_solid; 89 } 90 } 91 92 /******************************************************************************* 93 * Medium data 94 ******************************************************************************/ 95 static double 96 temperature_unknown(const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 97 { 98 (void)data; 99 CHK(vtx != NULL); 100 return SDIS_TEMPERATURE_NONE; 101 } 102 103 static double 104 solid_get_calorific_capacity 105 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 106 { 107 (void)data; 108 CHK(vtx != NULL && data == NULL); 109 return 2.0; 110 } 111 112 static double 113 solid_get_thermal_conductivity 114 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 115 { 116 (void)data; 117 CHK(vtx != NULL); 118 return 50.0; 119 } 120 121 static double 122 solid_get_volumic_mass 123 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 124 { 125 (void)data; 126 CHK(vtx != NULL); 127 return 25.0; 128 } 129 130 static double 131 solid_get_delta 132 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 133 { 134 (void)data; 135 CHK(vtx != NULL); 136 return 1.0/20.0; 137 } 138 139 /******************************************************************************* 140 * Interface 141 ******************************************************************************/ 142 struct interf { 143 double temperature; 144 }; 145 146 static double 147 null_interface_value 148 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 149 { 150 CHK(frag != NULL); 151 (void)data; 152 return 0; 153 } 154 155 static double 156 interface_get_temperature 157 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 158 { 159 CHK(data != NULL && frag != NULL); 160 return ((const struct interf*)sdis_data_cget(data))->temperature; 161 } 162 163 /******************************************************************************* 164 * Test 165 ******************************************************************************/ 166 int 167 main(int argc, char** argv) 168 { 169 struct sdis_mc T = SDIS_MC_NULL; 170 struct sdis_mc time = SDIS_MC_NULL; 171 struct sdis_device* dev = NULL; 172 struct sdis_data* data = NULL; 173 struct sdis_estimator* estimator = NULL; 174 struct sdis_estimator* estimator2 = NULL; 175 struct sdis_medium* solid = NULL; 176 struct sdis_medium* fluid = NULL; 177 struct sdis_interface* Tnone = NULL; 178 struct sdis_interface* T300 = NULL; 179 struct sdis_interface* T350 = NULL; 180 struct sdis_interface* solid_solid = NULL; 181 struct sdis_scene* scn = NULL; 182 struct sdis_green_function* green = NULL; 183 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 184 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 185 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 186 struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; 187 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 188 struct s3dut_mesh* msh = NULL; 189 struct s3dut_mesh_data msh_data; 190 struct context ctx = CONTEXT_NULL; 191 struct interf* interface_param = NULL; 192 double ref; 193 const size_t N = 10000; 194 size_t ntris; 195 size_t nverts; 196 size_t nreals; 197 size_t nfails; 198 size_t i; 199 (void)argc, (void)argv; 200 201 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 202 203 /* Create the fluid medium */ 204 fluid_shader.temperature = temperature_unknown; 205 OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); 206 207 /* Create the solid medium */ 208 solid_shader.calorific_capacity = solid_get_calorific_capacity; 209 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 210 solid_shader.volumic_mass = solid_get_volumic_mass; 211 solid_shader.delta = solid_get_delta; 212 solid_shader.temperature = temperature_unknown; 213 OK(sdis_solid_create(dev, &solid_shader, NULL, &solid)); 214 215 /* Create the fluid/solid interface with no limit conidition */ 216 interface_shader.convection_coef = null_interface_value; 217 interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; 218 interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; 219 OK(sdis_interface_create 220 (dev, solid, fluid, &interface_shader, NULL, &Tnone)); 221 222 /* Create the fluid/solid interface with a fixed temperature of 300K */ 223 OK(sdis_data_create(dev, sizeof(struct interf), 224 ALIGNOF(struct interf), NULL, &data)); 225 interface_param = sdis_data_get(data); 226 interface_param->temperature = 300; 227 interface_shader.front.temperature = interface_get_temperature; 228 OK(sdis_interface_create 229 (dev, solid, fluid, &interface_shader, data, &T300)); 230 OK(sdis_data_ref_put(data)); 231 232 /* Create the fluid/solid interface with a fixed temperature of 350K */ 233 OK(sdis_data_create(dev, sizeof(struct interf), 234 ALIGNOF(struct interf), NULL, &data)); 235 interface_param = sdis_data_get(data); 236 interface_param->temperature = 350; 237 OK(sdis_interface_create 238 (dev, solid, fluid, &interface_shader, data, &T350)); 239 OK(sdis_data_ref_put(data)); 240 241 /* Create the solid/solid interface */ 242 interface_shader = SDIS_INTERFACE_SHADER_NULL; 243 OK(sdis_interface_create 244 (dev, solid, solid, &interface_shader, NULL, &solid_solid)); 245 246 /* Release the media */ 247 OK(sdis_medium_ref_put(solid)); 248 OK(sdis_medium_ref_put(fluid)); 249 250 /* Register the box geometry */ 251 FOR_EACH(i, 0, box_nvertices) { 252 sa_push(ctx.positions, box_vertices[i*3+0]); 253 sa_push(ctx.positions, box_vertices[i*3+1]); 254 sa_push(ctx.positions, box_vertices[i*3+2]); 255 } 256 FOR_EACH(i, 0, box_ntriangles) { 257 sa_push(ctx.indices, box_indices[i*3+0]); 258 sa_push(ctx.indices, box_indices[i*3+1]); 259 sa_push(ctx.indices, box_indices[i*3+2]); 260 } 261 262 /* Setup a sphere at the center of the box */ 263 OK(s3dut_create_sphere(NULL, 0.25, 64, 32, &msh)); 264 OK(s3dut_mesh_get_data(msh, &msh_data)); 265 FOR_EACH(i, 0, msh_data.nvertices) { 266 sa_push(ctx.positions, msh_data.positions[i*3+0] + 0.5); 267 sa_push(ctx.positions, msh_data.positions[i*3+1] + 0.5); 268 sa_push(ctx.positions, msh_data.positions[i*3+2] + 0.5); 269 } 270 FOR_EACH(i, 0, msh_data.nprimitives) { 271 sa_push(ctx.indices, msh_data.indices[i*3+0] + box_nvertices); 272 sa_push(ctx.indices, msh_data.indices[i*3+1] + box_nvertices); 273 sa_push(ctx.indices, msh_data.indices[i*3+2] + box_nvertices); 274 } 275 OK(s3dut_mesh_ref_put(msh)); 276 277 /* Create the scene */ 278 ctx.solid_fluid_Tnone = Tnone; 279 ctx.solid_fluid_T300 = T300; 280 ctx.solid_fluid_T350 = T350; 281 ctx.solid_solid = solid_solid; 282 nverts = sa_size(ctx.positions) / 3; 283 ntris = sa_size(ctx.indices) / 3; 284 scn_args.get_indices = get_indices; 285 scn_args.get_interface = get_interface; 286 scn_args.get_position = get_position; 287 scn_args.nprimitives = ntris; 288 scn_args.nvertices = nverts; 289 scn_args.context = &ctx; 290 OK(sdis_scene_create(dev, &scn_args, &scn)); 291 292 /* Release the scene data */ 293 OK(sdis_interface_ref_put(Tnone)); 294 OK(sdis_interface_ref_put(T300)); 295 OK(sdis_interface_ref_put(T350)); 296 OK(sdis_interface_ref_put(solid_solid)); 297 sa_release(ctx.positions); 298 sa_release(ctx.indices); 299 300 /* Launch the solver */ 301 solve_args.nrealisations = N; 302 solve_args.position[0] = 0.5; 303 solve_args.position[1] = 0.5; 304 solve_args.position[2] = 0.5; 305 solve_args.time_range[0] = INF; 306 solve_args.time_range[1] = INF; 307 308 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 309 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 310 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 311 OK(sdis_estimator_get_temperature(estimator, &T)); 312 OK(sdis_estimator_get_realisation_time(estimator, &time)); 313 314 /* Print the estimation results */ 315 ref = 350 * solve_args.position[2] + (1-solve_args.position[2]) * 300; 316 printf("Temperature at (%g, %g, %g) = %g ~ %g +/- %g\n", 317 SPLIT3(solve_args.position), ref, T.E, T.SE); 318 printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); 319 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 320 321 /* Check the results */ 322 CHK(nfails + nreals == N); 323 CHK(nfails < N/1000); 324 CHK(eq_eps(T.E, ref, 3*T.SE)); 325 326 /* Check green function */ 327 OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); 328 OK(sdis_green_function_solve(green, &estimator2)); 329 check_green_function(green); 330 check_estimator_eq(estimator, estimator2); 331 check_green_serialization(green, scn); 332 333 /* Release data */ 334 OK(sdis_estimator_ref_put(estimator)); 335 OK(sdis_estimator_ref_put(estimator2)); 336 OK(sdis_green_function_ref_put(green)); 337 OK(sdis_scene_ref_put(scn)); 338 OK(sdis_device_ref_put(dev)); 339 340 CHK(mem_allocated_size() == 0); 341 return 0; 342 343 }