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