test_sdis_compute_power.c (11462B)
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 <star/s3dut.h> 21 22 #define N 100000ul /* #realisations */ 23 #define POWER0 10 24 #define POWER1 5 25 26 static INLINE void 27 check_intersection 28 (const double val0, 29 const double eps0, 30 const double val1, 31 const double eps1) 32 { 33 double interval0[2], interval1[2]; 34 double intersection[2]; 35 interval0[0] = val0 - eps0; 36 interval0[1] = val0 + eps0; 37 interval1[0] = val1 - eps1; 38 interval1[1] = val1 + eps1; 39 intersection[0] = MMAX(interval0[0], interval1[0]); 40 intersection[1] = MMIN(interval0[1], interval1[1]); 41 CHK(intersection[0] <= intersection[1]); 42 } 43 44 /******************************************************************************* 45 * Geometry 46 ******************************************************************************/ 47 struct context { 48 struct s3dut_mesh_data msh0; 49 struct s3dut_mesh_data msh1; 50 struct sdis_interface* interf0; 51 struct sdis_interface* interf1; 52 }; 53 54 static void 55 get_indices(const size_t itri, size_t ids[3], void* context) 56 { 57 const struct context* ctx = context; 58 /* Note that we swap the indices to ensure that the triangle normals point 59 * inward the super shape */ 60 if(itri < ctx->msh0.nprimitives) { 61 ids[0] = ctx->msh0.indices[itri*3+0]; 62 ids[2] = ctx->msh0.indices[itri*3+1]; 63 ids[1] = ctx->msh0.indices[itri*3+2]; 64 } else { 65 const size_t itri2 = itri - ctx->msh0.nprimitives; 66 ids[0] = ctx->msh1.indices[itri2*3+0] + ctx->msh0.nvertices; 67 ids[2] = ctx->msh1.indices[itri2*3+1] + ctx->msh0.nvertices; 68 ids[1] = ctx->msh1.indices[itri2*3+2] + ctx->msh0.nvertices; 69 } 70 } 71 72 static void 73 get_position(const size_t ivert, double pos[3], void* context) 74 { 75 const struct context* ctx = context; 76 if(ivert < ctx->msh0.nvertices) { 77 pos[0] = ctx->msh0.positions[ivert*3+0] - 2.0; 78 pos[1] = ctx->msh0.positions[ivert*3+1]; 79 pos[2] = ctx->msh0.positions[ivert*3+2]; 80 } else { 81 const size_t ivert2 = ivert - ctx->msh0.nvertices; 82 pos[0] = ctx->msh1.positions[ivert2*3+0] + 2.0; 83 pos[1] = ctx->msh1.positions[ivert2*3+1]; 84 pos[2] = ctx->msh1.positions[ivert2*3+2]; 85 } 86 } 87 88 static void 89 get_interface(const size_t itri, struct sdis_interface** bound, void* context) 90 { 91 const struct context* ctx = context; 92 *bound = itri < ctx->msh0.nprimitives ? ctx->interf0 : ctx->interf1; 93 } 94 95 /******************************************************************************* 96 * Interface 97 ******************************************************************************/ 98 static double 99 interface_get_convection_coef 100 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 101 { 102 CHK(frag != NULL); (void)data; 103 return 1.0; 104 } 105 106 /******************************************************************************* 107 * Fluid medium 108 ******************************************************************************/ 109 static double 110 fluid_get_temperature 111 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 112 { 113 CHK(vtx == NULL); (void)data; 114 return 300; 115 } 116 117 /******************************************************************************* 118 * Solid medium 119 ******************************************************************************/ 120 static double 121 solid_get_calorific_capacity 122 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 123 { 124 CHK(vtx != NULL); (void)data; 125 return 1.0; 126 } 127 128 static double 129 solid_get_thermal_conductivity 130 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 131 { 132 CHK(vtx != NULL); (void)data; 133 return 1.0; 134 } 135 136 static double 137 solid_get_volumic_mass 138 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 139 { 140 CHK(vtx != NULL); (void)data; 141 return 1.0; 142 } 143 144 static double 145 solid_get_delta 146 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 147 { 148 CHK(vtx != NULL); (void)data; 149 return 1.0 / 20.0; 150 } 151 152 static double 153 solid_get_volumic_power 154 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 155 { 156 CHK(vtx != NULL); (void)data; 157 return vtx->P[0] < 0 ? POWER0 : POWER1; 158 } 159 160 /******************************************************************************* 161 * Test 162 ******************************************************************************/ 163 int 164 main(int argc, char** argv) 165 { 166 struct context ctx; 167 struct s3dut_mesh* sphere = NULL; 168 struct s3dut_mesh* cylinder = NULL; 169 struct sdis_data* data = NULL; 170 struct sdis_device* dev = NULL; 171 struct sdis_estimator* estimator = NULL; 172 struct sdis_medium* fluid = NULL; 173 struct sdis_medium* solid0 = NULL; 174 struct sdis_medium* solid1 = NULL; 175 struct sdis_interface* interf0 = NULL; 176 struct sdis_interface* interf1 = NULL; 177 struct sdis_scene* scn = NULL; 178 struct sdis_mc mpow = SDIS_MC_NULL; 179 struct sdis_mc time = SDIS_MC_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 interf_shader = SDIS_INTERFACE_SHADER_NULL; 184 struct sdis_compute_power_args args = SDIS_COMPUTE_POWER_ARGS_DEFAULT; 185 size_t nverts = 0; 186 size_t ntris = 0; 187 double ref = 0; 188 int is_master_process; 189 (void)argc, (void) argv; 190 191 create_default_device(&argc, &argv, &is_master_process, &dev); 192 193 /* Setup the interface shader */ 194 interf_shader.convection_coef = interface_get_convection_coef; 195 196 /* Setup the fluid shader */ 197 fluid_shader.temperature = fluid_get_temperature; 198 199 /* Setup the solid shader */ 200 solid_shader.calorific_capacity = solid_get_calorific_capacity; 201 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 202 solid_shader.volumic_mass = solid_get_volumic_mass; 203 solid_shader.delta = solid_get_delta; 204 solid_shader.volumic_power = solid_get_volumic_power; 205 206 /* Create the fluid */ 207 OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); 208 209 /* Create the solids */ 210 OK(sdis_solid_create(dev, &solid_shader, data, &solid0)); 211 OK(sdis_solid_create(dev, &solid_shader, data, &solid1)); 212 213 /* Create the interface0 */ 214 OK(sdis_interface_create(dev, solid0, fluid, &interf_shader, NULL, &interf0)); 215 OK(sdis_interface_create(dev, solid1, fluid, &interf_shader, NULL, &interf1)); 216 ctx.interf0 = interf0; 217 ctx.interf1 = interf1; 218 219 /* Create the geometry */ 220 OK(s3dut_create_sphere(NULL, 1, 512, 256, &sphere)); 221 OK(s3dut_create_cylinder(NULL, 1, 10, 512, 8, &cylinder)); 222 OK(s3dut_mesh_get_data(sphere, &ctx.msh0)); 223 OK(s3dut_mesh_get_data(cylinder, &ctx.msh1)); 224 225 /* Create the scene */ 226 ntris = ctx.msh0.nprimitives + ctx.msh1.nprimitives; 227 nverts = ctx.msh0.nvertices + ctx.msh1.nvertices; 228 scn_args.get_indices = get_indices; 229 scn_args.get_interface = get_interface; 230 scn_args.get_position = get_position; 231 scn_args.nprimitives = ntris; 232 scn_args.nvertices = nverts; 233 scn_args.context = &ctx; 234 OK(sdis_scene_create(dev, &scn_args, &scn)); 235 236 /* Test sdis_compute_power function */ 237 args.nrealisations = N; 238 args.medium = solid0; 239 args.time_range[0] = INF; 240 args.time_range[1] = INF; 241 BA(sdis_compute_power(NULL, &args, &estimator)); 242 BA(sdis_compute_power(scn, NULL, &estimator)); 243 BA(sdis_compute_power(scn, &args, NULL)); 244 args.nrealisations = 0; 245 BA(sdis_compute_power(scn, &args, &estimator)); 246 args.nrealisations = N; 247 args.medium = NULL; 248 BA(sdis_compute_power(scn, &args, &estimator)); 249 args.medium = solid0; 250 args.time_range[0] = args.time_range[1] = -1; 251 BA(sdis_compute_power(scn, &args, &estimator)); 252 args.time_range[0] = 1; 253 BA(sdis_compute_power(scn, &args, &estimator)); 254 args.time_range[1] = 0; 255 BA(sdis_compute_power(scn, &args, &estimator)); 256 args.time_range[0] = args.time_range[1] = INF; 257 OK(sdis_compute_power(scn, &args, &estimator)); 258 259 if(!is_master_process) { 260 CHK(estimator == NULL); 261 } else { 262 BA(sdis_estimator_get_power(NULL, &mpow)); 263 BA(sdis_estimator_get_power(estimator, NULL)); 264 OK(sdis_estimator_get_power(estimator, &mpow)); 265 OK(sdis_estimator_get_realisation_time(estimator, &time)); 266 267 /* Check results for solid 0 */ 268 ref = 4.0/3.0 * PI * POWER0; 269 printf("Mean power of the solid0 = %g W ~ %g W +/- %g\n", 270 ref, mpow.E, mpow.SE); 271 check_intersection(ref, 1.e-3*ref, mpow.E, 3*mpow.SE); 272 OK(sdis_estimator_ref_put(estimator)); 273 } 274 275 args.medium = solid1; 276 OK(sdis_compute_power(scn, &args, &estimator)); 277 278 if(is_master_process) { 279 /* Check results for solid 1 */ 280 OK(sdis_estimator_get_power(estimator, &mpow)); 281 ref = PI * 10 * POWER1; 282 printf("Mean power of the solid1 = %g W ~ %g W +/- %g\n", 283 ref, mpow.E, mpow.SE); 284 check_intersection(ref, 1.e-3*ref, mpow.E, 3*mpow.SE); 285 OK(sdis_estimator_ref_put(estimator)); 286 } 287 288 args.time_range[0] = 0; 289 args.time_range[1] = 10; 290 OK(sdis_compute_power(scn, &args, &estimator)); 291 292 if(is_master_process) { 293 /* Check for a not null time range */ 294 OK(sdis_estimator_get_power(estimator, &mpow)); 295 ref = PI * 10 * POWER1; 296 printf("Mean power of the solid1 in [0, 10] s = %g W ~ %g W +/- %g\n", 297 ref, mpow.E, mpow.SE); 298 check_intersection(ref, 1.e-3*ref, mpow.E, 3*mpow.SE); 299 OK(sdis_estimator_ref_put(estimator)); 300 } 301 302 /* Reset the scene with only one solid medium */ 303 OK(sdis_scene_ref_put(scn)); 304 ctx.interf0 = interf0; 305 ctx.interf1 = interf0; 306 OK(sdis_scene_create(dev, &scn_args, &scn)); 307 308 /* Check invalid medium */ 309 args.time_range[0] = args.time_range[1] = 1; 310 args.medium = solid1; 311 BA(sdis_compute_power(scn, &args, &estimator)); 312 313 /* Check non constant volumic power */ 314 args.medium = solid0; 315 OK(sdis_compute_power(scn, &args, &estimator)); 316 if(is_master_process) { 317 OK(sdis_estimator_get_power(estimator, &mpow)); 318 ref = 4.0/3.0*PI*POWER0 + PI*10*POWER1; 319 printf("Mean power of the sphere+cylinder = %g W ~ %g W +/- %g\n", 320 ref, mpow.E, mpow.SE); 321 check_intersection(ref, 1e-2*ref, mpow.E, 3*mpow.SE); 322 OK(sdis_estimator_ref_put(estimator)); 323 } 324 325 #if 0 326 { 327 double* vertices = NULL; 328 size_t* indices = NULL; 329 size_t i; 330 CHK(vertices = MEM_CALLOC(&allocator, nverts*3, sizeof(*vertices))); 331 CHK(indices = MEM_CALLOC(&allocator, ntris*3, sizeof(*indices))); 332 FOR_EACH(i, 0, ntris) get_indices(i, indices + i*3, &ctx); 333 FOR_EACH(i, 0, nverts) get_position(i, vertices + i*3, &ctx); 334 dump_mesh(stdout, vertices, nverts, indices, ntris); 335 MEM_RM(&allocator, vertices); 336 MEM_RM(&allocator, indices); 337 } 338 #endif 339 340 /* Clean up memory */ 341 OK(sdis_medium_ref_put(fluid)); 342 OK(sdis_medium_ref_put(solid0)); 343 OK(sdis_medium_ref_put(solid1)); 344 OK(sdis_interface_ref_put(interf0)); 345 OK(sdis_interface_ref_put(interf1)); 346 OK(sdis_scene_ref_put(scn)); 347 OK(s3dut_mesh_ref_put(sphere)); 348 OK(s3dut_mesh_ref_put(cylinder)); 349 350 free_default_device(dev); 351 352 CHK(mem_allocated_size() == 0); 353 return 0; 354 }