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