test_sdis_external_flux_with_diffuse_radiance.c (14686B)
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 "test_sdis_utils.h" 17 #include "sdis.h" 18 19 #include <rsys/rsys.h> 20 #include <star/s3dut.h> 21 22 /* 23 * The system is a sphere with an emissivity of 1. It is illuminated by an 24 * external spherical source that is sufficiently far away from the sphere 25 * radius to assume a constant flux density received per m^2 of surface 26 * perpendicular to the direction towards the source center. 27 * 28 * In such situation, the value of the surface temperature is equal to : 29 * 30 * Ts = [q/(4*sigma)+Tamb^4]^0.25 31 * q = P/(4PI*d^2) 32 * 33 * with Tamb the temperature of the environment, sigma the Boltzmann constant, 34 * P the source power and d the distance of the source from the center of 35 * the sphere. 36 * 37 * If one adds a diffuse radiance Ldiff, the surface temperature becomes : 38 * 39 * Ts = [(q/4+Ldiff*PI)/sigma+Tamb^4]^0.25 40 * 41 * The test checks that we retrieved these analatycal results by Monte Carlo 42 * 43 * 44 * R = 0.01 m External source 45 * __ d = 10 m ### 46 * T = ? -> / .\ ..............................##### r = 0.3 m 47 * \__/ ##### 48 * E=1 ### 49 * P = 3^7 W 50 */ 51 52 /* The source */ 53 #define SOURCE_POWER 1.0e5 /* [W] */ 54 #define SOURCE_RADIUS 0.3 /* [m/fp_to_meter] */ 55 #define SOURCE_DISTANCE 10 /* [m/fp_to_meter] */ 56 57 /* Miscellaneous */ 58 #define SPHERE_RADIUS 0.01 /* [m/fp_to_meter] */ 59 #define SPHERE_T_REF 320 /* [K] */ 60 #define T_RAD 300.0 /* [K] */ 61 #define T_REF 300.0 /* [K] */ 62 #define T_RAD4 (T_RAD*T_RAD*T_RAD*T_RAD) 63 64 #define NREALISATIONS 10000 65 66 /******************************************************************************* 67 * Geometry 68 ******************************************************************************/ 69 static struct s3dut_mesh* 70 create_sphere(void) 71 { 72 struct s3dut_mesh* mesh = NULL; 73 OK(s3dut_create_sphere(NULL, SPHERE_RADIUS, 128, 64, &mesh)); 74 return mesh; 75 } 76 77 /******************************************************************************* 78 * Media 79 ******************************************************************************/ 80 #define MEDIUM_PROP(Type, Prop, Val) \ 81 static double \ 82 Type##_get_##Prop \ 83 (const struct sdis_rwalk_vertex* vtx, \ 84 struct sdis_data* data) \ 85 { \ 86 (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \ 87 return Val; \ 88 } 89 MEDIUM_PROP(medium, volumic_mass, 1) /* [kj/m^3] */ 90 MEDIUM_PROP(medium, calorific_capacity, 1) /* [J/K/Kg] */ 91 MEDIUM_PROP(medium, temperature, SDIS_TEMPERATURE_NONE) /* [K] */ 92 MEDIUM_PROP(solid, thermal_conductivity, 1) /* [W/m/K] */ 93 MEDIUM_PROP(solid, delta, (2*SPHERE_RADIUS)/10.0) /* [m] */ 94 #undef MEDIUM_PROP 95 96 static struct sdis_medium* 97 create_solid(struct sdis_device* sdis) 98 { 99 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 100 struct sdis_medium* solid = NULL; 101 102 shader.calorific_capacity = medium_get_calorific_capacity; 103 shader.thermal_conductivity = solid_get_thermal_conductivity; 104 shader.volumic_mass = medium_get_volumic_mass; 105 shader.delta = solid_get_delta; 106 shader.temperature = medium_get_temperature; 107 OK(sdis_solid_create(sdis, &shader, NULL, &solid)); 108 return solid; 109 } 110 111 static struct sdis_medium* 112 create_fluid(struct sdis_device* sdis) 113 { 114 struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL; 115 struct sdis_medium* fluid = NULL; 116 117 shader.calorific_capacity = medium_get_calorific_capacity; 118 shader.volumic_mass = medium_get_volumic_mass; 119 shader.temperature = medium_get_temperature; 120 OK(sdis_fluid_create(sdis, &shader, NULL, &fluid)); 121 return fluid; 122 } 123 124 /******************************************************************************* 125 * Interface 126 ******************************************************************************/ 127 static double 128 interface_get_emissivity 129 (const struct sdis_interface_fragment* frag, 130 const unsigned source_id, 131 struct sdis_data* data) 132 { 133 (void)frag, (void)source_id, (void)data; 134 return 1; 135 } 136 137 static double 138 interface_get_reference_temperature 139 (const struct sdis_interface_fragment* frag, 140 struct sdis_data* data) 141 { 142 (void)frag, (void)data; 143 return SPHERE_T_REF; /* [K] */ 144 } 145 146 static struct sdis_interface* 147 create_interface 148 (struct sdis_device* sdis, 149 struct sdis_medium* front, 150 struct sdis_medium* back) 151 { 152 struct sdis_interface* interf = NULL; 153 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 154 155 shader.front.emissivity = interface_get_emissivity; 156 shader.front.reference_temperature = interface_get_reference_temperature; 157 OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf)); 158 return interf; 159 } 160 161 /******************************************************************************* 162 * Radiative environment 163 ******************************************************************************/ 164 #define RADENV_PROP(Prop, Val) \ 165 static double \ 166 radenv_get_##Prop \ 167 (const struct sdis_radiative_ray* ray, \ 168 struct sdis_data* data) \ 169 { \ 170 (void)ray, (void)data; \ 171 return (Val); /* [K] */ \ 172 } 173 RADENV_PROP(temperature, T_RAD) 174 RADENV_PROP(reference_temperature, T_REF) 175 #undef RADENV_PROP 176 177 static struct sdis_radiative_env* 178 create_radenv(struct sdis_device* sdis) 179 { 180 struct sdis_radiative_env_shader shader = SDIS_RADIATIVE_ENV_SHADER_NULL; 181 struct sdis_radiative_env* radenv = NULL; 182 183 shader.temperature = radenv_get_temperature; 184 shader.reference_temperature = radenv_get_reference_temperature; 185 OK(sdis_radiative_env_create(sdis, &shader, NULL, &radenv)); 186 return radenv; 187 } 188 189 /******************************************************************************* 190 * External source 191 ******************************************************************************/ 192 static void 193 source_get_position 194 (const double time, 195 double pos[3], 196 struct sdis_data* data) 197 { 198 (void)time, (void)data; /* Avoid the "unusued variable" warning */ 199 pos[0] = SOURCE_DISTANCE; 200 pos[1] = 0; 201 pos[2] = 0; 202 } 203 204 static double 205 source_get_power(const double time, struct sdis_data* data) 206 { 207 (void)time, (void)data; /* Avoid the "unused variable" warning */ 208 return SOURCE_POWER; /* [W] */ 209 } 210 211 static double 212 source_get_diffuse_radiance 213 (const double time, 214 const double dir[3], 215 struct sdis_data* data) 216 { 217 const double* Ldiff = sdis_data_cget(data);; 218 (void)time, (void)dir; /* Avoid the "unused variable" warning */ 219 return *Ldiff; /* [W/m^2/sr] */ 220 } 221 222 static struct sdis_source* 223 create_source(struct sdis_device* sdis, double** diffuse_radiance) 224 { 225 struct sdis_spherical_source_shader shader = SDIS_SPHERICAL_SOURCE_SHADER_NULL; 226 struct sdis_data* data = NULL; 227 struct sdis_source* src = NULL; 228 229 OK(sdis_data_create(sdis, sizeof(double), ALIGNOF(double), NULL, &data)); 230 *diffuse_radiance = sdis_data_get(data); 231 **diffuse_radiance = 0; 232 233 shader.position = source_get_position; 234 shader.power = source_get_power; 235 shader.diffuse_radiance = source_get_diffuse_radiance; 236 shader.radius = SOURCE_RADIUS; 237 OK(sdis_spherical_source_create(sdis, &shader, data, &src)); 238 OK(sdis_data_ref_put(data)); 239 return src; 240 } 241 242 /******************************************************************************* 243 * The scene 244 ******************************************************************************/ 245 struct scene_context { 246 const struct s3dut_mesh_data* mesh; 247 struct sdis_interface* interf; 248 }; 249 static const struct scene_context SCENE_CONTEXT_NULL = {NULL, NULL}; 250 251 static void 252 scene_get_indices(const size_t itri, size_t ids[3], void* ctx) 253 { 254 struct scene_context* context = ctx; 255 CHK(ids && context && itri < context->mesh->nprimitives); 256 ids[0] = (unsigned)context->mesh->indices[itri*3+0]; 257 ids[1] = (unsigned)context->mesh->indices[itri*3+1]; 258 ids[2] = (unsigned)context->mesh->indices[itri*3+2]; 259 } 260 261 static void 262 scene_get_interface(const size_t itri, struct sdis_interface** interf, void* ctx) 263 { 264 struct scene_context* context = ctx; 265 CHK(interf && context && itri < context->mesh->nprimitives); 266 *interf = context->interf; 267 } 268 269 static void 270 scene_get_position(const size_t ivert, double pos[3], void* ctx) 271 { 272 struct scene_context* context = ctx; 273 CHK(pos && context && ivert < context->mesh->nvertices); 274 pos[0] = context->mesh->positions[ivert*3+0]; 275 pos[1] = context->mesh->positions[ivert*3+1]; 276 pos[2] = context->mesh->positions[ivert*3+2]; 277 } 278 279 static struct sdis_scene* 280 create_scene 281 (struct sdis_device* sdis, 282 const struct s3dut_mesh_data* mesh, 283 struct sdis_interface* interf, 284 struct sdis_source* source, 285 struct sdis_radiative_env* radenv) 286 { 287 struct sdis_scene* scn = NULL; 288 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 289 struct scene_context context = SCENE_CONTEXT_NULL; 290 291 context.mesh = mesh; 292 context.interf = interf; 293 294 scn_args.get_indices = scene_get_indices; 295 scn_args.get_interface = scene_get_interface; 296 scn_args.get_position = scene_get_position; 297 scn_args.nprimitives = mesh->nprimitives; 298 scn_args.nvertices = mesh->nvertices; 299 scn_args.t_range[0] = MMIN(T_REF, SPHERE_T_REF); 300 scn_args.t_range[1] = MMAX(T_REF, SPHERE_T_REF); 301 scn_args.source = source; 302 scn_args.radenv = radenv; 303 scn_args.context = &context; 304 OK(sdis_scene_create(sdis, &scn_args, &scn)); 305 return scn; 306 } 307 308 /******************************************************************************* 309 * Validations 310 ******************************************************************************/ 311 static double 312 analytic_temperature(const double Ldiff) 313 { 314 const double q = SOURCE_POWER/(4*PI*SOURCE_DISTANCE*SOURCE_DISTANCE); 315 const double Ts = pow((q/4 + Ldiff*PI)/BOLTZMANN_CONSTANT + T_RAD4, 0.25); 316 return Ts; 317 } 318 319 static void 320 check 321 (struct sdis_scene* scn, 322 const struct s3dut_mesh_data* mesh, 323 const double Ldiff) 324 { 325 struct sdis_solve_boundary_args args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT; 326 struct sdis_mc T = SDIS_MC_NULL; 327 struct sdis_estimator* estimator = NULL; 328 enum sdis_side* sides = NULL; 329 double ref = 0; 330 size_t i = 0; 331 332 sides = mem_alloc(mesh->nprimitives*sizeof(*sides)); 333 FOR_EACH(i, 0, mesh->nprimitives) sides[i] = SDIS_FRONT; 334 335 args.primitives = mesh->indices; 336 args.sides = sides; 337 args.nprimitives = mesh->nprimitives; 338 args.nrealisations = NREALISATIONS; 339 OK(sdis_solve_boundary(scn, &args, &estimator)); 340 341 OK(sdis_estimator_get_temperature(estimator, &T)); 342 ref = analytic_temperature(Ldiff); 343 printf("Ts = %g ~ %g +/- %g\n", ref, T.E, T.SE); 344 OK(sdis_estimator_ref_put(estimator)); 345 346 CHK(eq_eps(T.E, ref, T.SE * 3)); 347 348 mem_rm(sides); 349 } 350 351 static void 352 solve_green(struct sdis_green_function* green, const double Ldiff) 353 { 354 struct sdis_mc T = SDIS_MC_NULL; 355 struct sdis_estimator* estimator = NULL; 356 double ref = 0; 357 358 OK(sdis_green_function_solve(green, &estimator)); 359 OK(sdis_estimator_get_temperature(estimator, &T)); 360 361 ref = analytic_temperature(Ldiff); 362 printf("Ts = %g ~ %g +/- %g\n", ref, T.E, T.SE); 363 CHK(eq_eps(T.E, ref, T.SE * 3)); 364 365 OK(sdis_estimator_ref_put(estimator)); 366 } 367 368 static void 369 check_green 370 (struct sdis_scene* scn, 371 const struct s3dut_mesh_data* mesh, 372 double* Ldiff) 373 { 374 struct sdis_solve_boundary_args args = SDIS_SOLVE_BOUNDARY_ARGS_DEFAULT; 375 struct sdis_green_function* green = NULL; 376 enum sdis_side* sides = NULL; 377 size_t i = 0; 378 379 CHK(Ldiff); 380 381 sides = mem_alloc(mesh->nprimitives*sizeof(*sides)); 382 FOR_EACH(i, 0, mesh->nprimitives) sides[i] = SDIS_FRONT; 383 384 *Ldiff = 0; 385 386 args.primitives = mesh->indices; 387 args.sides = sides; 388 args.nprimitives = mesh->nprimitives; 389 args.nrealisations = NREALISATIONS; 390 OK(sdis_solve_boundary_green_function(scn, &args, &green)); 391 392 solve_green(green, (*Ldiff = 50)); 393 solve_green(green, (*Ldiff = 45)); 394 solve_green(green, (*Ldiff = 40)); 395 396 OK(sdis_green_function_ref_put(green)); 397 398 mem_rm(sides); 399 } 400 401 /******************************************************************************* 402 * The test 403 ******************************************************************************/ 404 int 405 main(int argc, char** argv) 406 { 407 /* Stardis */ 408 struct sdis_device* dev = NULL; 409 struct sdis_interface* interf = NULL; 410 struct sdis_medium* fluid = NULL; 411 struct sdis_medium* solid = NULL; 412 struct sdis_radiative_env* radenv = NULL; 413 struct sdis_scene* scene = NULL; 414 struct sdis_source* source = NULL; 415 416 /* Miscellaneous */ 417 struct s3dut_mesh_data mesh; 418 struct s3dut_mesh* sphere = NULL; 419 double* Ldiff = NULL; 420 421 (void)argc, (void)argv; 422 423 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 424 425 sphere = create_sphere(); 426 OK(s3dut_mesh_get_data(sphere, &mesh)); 427 428 fluid = create_fluid(dev); 429 solid = create_solid(dev); 430 interf = create_interface(dev, fluid, solid); 431 radenv = create_radenv(dev); 432 source = create_source(dev, &Ldiff); 433 434 scene = create_scene(dev, &mesh, interf, source, radenv); 435 436 check(scene, &mesh, (*Ldiff = 50)); 437 check_green(scene, &mesh, Ldiff); 438 439 OK(s3dut_mesh_ref_put(sphere)); 440 OK(sdis_device_ref_put(dev)); 441 OK(sdis_interface_ref_put(interf)); 442 OK(sdis_medium_ref_put(fluid)); 443 OK(sdis_medium_ref_put(solid)); 444 OK(sdis_radiative_env_ref_put(radenv)); 445 OK(sdis_scene_ref_put(scene)); 446 OK(sdis_source_ref_put(source)); 447 CHK(mem_allocated_size() == 0); 448 return 0; 449 }