test_sdis_unsteady.c (9075B)
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 /* 20 * The scene is a solid cube whose temperature is fixed on its faces. The test 21 * calculates its temperature at a given position at different observation times 22 * and validates the result against the calculated values by analytically 23 * solving the green functions. 24 * 25 * T3 (0.1,0.1,0.1) 26 * +-------+ 27 * /' T4 /| 28 * +-------+ | T1 29 * T0 | +.....|.+ 30 * |, T5 |/ 31 * +-------+ 32 * (0,0,0) T2 33 */ 34 35 #define T_INIT 280 /* [K] */ 36 #define T0 310 /* [K] */ 37 #define T1 320 /* [K] */ 38 #define T2 330 /* [K] */ 39 #define T3 310 /* [K] */ 40 #define T4 320 /* [K] */ 41 #define T5 300 /* [K] */ 42 43 #define NREALISATIONS 10000 44 #define FP_TO_METER 0.1 45 46 struct reference { 47 double pos[3]; /* [m/FP_TO_METER] */ 48 double time; /* [s] */ 49 double temp; /* [K] */ 50 }; 51 52 static const struct reference references[] = { 53 {{0.3, 0.4, 0.6}, 1000.0000, 281.33455593977152}, 54 {{0.3, 0.4, 0.6}, 2000.0000, 286.90151817350699}, 55 {{0.3, 0.4, 0.6}, 3000.0000, 292.84330866161531}, 56 {{0.3, 0.4, 0.6}, 4000.0000, 297.81444160746452}, 57 {{0.3, 0.4, 0.6}, 5000.0000, 301.70787295764546}, 58 {{0.3, 0.4, 0.6}, 10000.000, 310.78920179442139}, 59 {{0.3, 0.4, 0.6}, 20000.000, 313.37629443163121}, 60 {{0.3, 0.4, 0.6}, 30000.000, 313.51064004438581}, 61 {{0.3, 0.4, 0.6}, 1000000.0, 313.51797642855502} 62 }; 63 static const size_t nreferences = sizeof(references)/sizeof(*references); 64 65 /******************************************************************************* 66 * Solid, i.e. medium of the cube 67 ******************************************************************************/ 68 #define SOLID_PROP(Prop, Val) \ 69 static double \ 70 solid_get_##Prop \ 71 (const struct sdis_rwalk_vertex* vtx, \ 72 struct sdis_data* data) \ 73 { \ 74 (void)vtx, (void)data; /* Avoid the "unused variable" warning */ \ 75 return Val; \ 76 } 77 SOLID_PROP(calorific_capacity, 2000.0) /* [J/K/kg] */ 78 SOLID_PROP(thermal_conductivity, 0.5) /* [W/m/K] */ 79 SOLID_PROP(volumic_mass, 2500.0) /* [kg/m^3] */ 80 SOLID_PROP(delta, 1.0/60.0) 81 #undef SOLID_PROP 82 83 static double /* [K] */ 84 solid_get_temperature 85 (const struct sdis_rwalk_vertex* vtx, 86 struct sdis_data* data) 87 { 88 (void)data; 89 ASSERT(vtx); 90 if(vtx->time <= 0) return T_INIT; /* Initial temperature [K] */ 91 return SDIS_TEMPERATURE_NONE; 92 } 93 94 static struct sdis_medium* 95 create_solid(struct sdis_device* sdis) 96 { 97 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 98 struct sdis_medium* solid = NULL; 99 100 shader.calorific_capacity = solid_get_calorific_capacity; 101 shader.thermal_conductivity = solid_get_thermal_conductivity; 102 shader.volumic_mass = solid_get_volumic_mass; 103 shader.delta = solid_get_delta; 104 shader.temperature = solid_get_temperature; 105 OK(sdis_solid_create(sdis, &shader, NULL, &solid)); 106 return solid; 107 } 108 109 /******************************************************************************* 110 * Dummy environment, i.e. environment surrounding the cube. It is defined only 111 * for Stardis compliance: in Stardis, an interface must divide 2 media. 112 ******************************************************************************/ 113 static struct sdis_medium* 114 create_dummy(struct sdis_device* sdis) 115 { 116 struct sdis_fluid_shader shader = SDIS_FLUID_SHADER_NULL; 117 struct sdis_medium* dummy = NULL; 118 119 shader.calorific_capacity = dummy_medium_getter; 120 shader.volumic_mass = dummy_medium_getter; 121 shader.temperature = dummy_medium_getter; 122 OK(sdis_fluid_create(sdis, &shader, NULL, &dummy)); 123 return dummy; 124 } 125 126 /******************************************************************************* 127 * Interface: its temperature is known 128 ******************************************************************************/ 129 static double /* [K] */ 130 interface_get_temperature 131 (const struct sdis_interface_fragment* frag, 132 struct sdis_data* data) 133 { 134 (void)data; 135 ASSERT(frag); 136 137 if(frag->Ng[0] == 1) return T0; 138 else if(frag->Ng[0] == -1) return T1; 139 else if(frag->Ng[1] == 1) return T2; 140 else if(frag->Ng[1] == -1) return T3; 141 else if(frag->Ng[2] == 1) return T4; 142 else if(frag->Ng[2] == -1) return T5; 143 else FATAL("Unreachable code\n"); 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.temperature = interface_get_temperature; 156 shader.back.temperature = interface_get_temperature; 157 OK(sdis_interface_create(sdis, front, back, &shader, NULL, &interf)); 158 return interf; 159 } 160 161 /******************************************************************************* 162 * The scene 163 ******************************************************************************/ 164 static void 165 get_interface(const size_t itri, struct sdis_interface** interf, void* ctx) 166 { 167 (void)itri; 168 ASSERT(interf && ctx); 169 *interf = ctx; 170 } 171 172 static struct sdis_scene* 173 create_scene 174 (struct sdis_device* sdis, 175 struct sdis_interface* interf) 176 { 177 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 178 struct sdis_scene* scn = NULL; 179 180 scn_args.get_indices = box_get_indices; 181 scn_args.get_interface = get_interface; 182 scn_args.get_position = box_get_position; 183 scn_args.nprimitives = box_ntriangles; 184 scn_args.nvertices = box_nvertices; 185 scn_args.context = interf; 186 scn_args.fp_to_meter = FP_TO_METER; 187 OK(sdis_scene_create(sdis, &scn_args, &scn)); 188 189 return scn; 190 } 191 192 /******************************************************************************* 193 * Validations 194 ******************************************************************************/ 195 static void 196 check_probe 197 (struct sdis_scene* scn, 198 const struct reference* ref, 199 const enum sdis_diffusion_algorithm diff_algo) 200 { 201 struct sdis_solve_probe_args args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 202 struct sdis_mc T = SDIS_MC_NULL; 203 struct sdis_estimator* estimator = NULL; 204 CHK(ref); 205 206 args.position[0] = ref->pos[0]; 207 args.position[1] = ref->pos[1]; 208 args.position[2] = ref->pos[2]; 209 args.time_range[0] = 210 args.time_range[1] = ref->time; /* [s] */ 211 args.diff_algo = diff_algo; 212 args.nrealisations = NREALISATIONS; 213 OK(sdis_solve_probe(scn, &args, &estimator)); 214 OK(sdis_estimator_get_temperature(estimator, &T)); 215 216 printf("T(%g, %g, %g) at %g s = %g ~ %g +/- %g\n", 217 ref->pos[0]*FP_TO_METER, ref->pos[1]*FP_TO_METER, ref->pos[2]*FP_TO_METER, 218 ref->time, ref->temp, T.E, T.SE); 219 CHK(eq_eps(ref->temp, T.E, 3*T.SE)); 220 221 OK(sdis_estimator_ref_put(estimator)); 222 } 223 224 static void 225 check 226 (struct sdis_scene* scn, 227 const enum sdis_diffusion_algorithm diff_algo) 228 { 229 const char* str_algo = NULL; 230 size_t i = 0; 231 232 switch(diff_algo) { 233 case SDIS_DIFFUSION_WOS: str_algo = "Walk on Sphere"; break; 234 case SDIS_DIFFUSION_DELTA_SPHERE: str_algo = "Delta sphere"; break; 235 default: FATAL("Unreachable code.\n"); break; 236 } 237 238 printf("\n\x1b[1m\x1b[35m%s\x1b[0m\n", str_algo); 239 FOR_EACH(i, 0, nreferences) { 240 check_probe(scn, references + i, diff_algo); 241 } 242 } 243 244 /******************************************************************************* 245 * The test 246 ******************************************************************************/ 247 int 248 main(int argc, char** argv) 249 { 250 struct sdis_device* sdis = NULL; 251 struct sdis_interface* interf = NULL; 252 struct sdis_medium* solid = NULL; 253 struct sdis_medium* dummy = NULL; 254 struct sdis_scene* scene = NULL; 255 (void)argc, (void)argv; 256 257 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &sdis)); 258 259 solid = create_solid(sdis); 260 dummy = create_dummy(sdis); 261 interf = create_interface(sdis, solid, dummy); 262 scene = create_scene(sdis, interf); 263 264 check(scene, SDIS_DIFFUSION_DELTA_SPHERE); 265 check(scene, SDIS_DIFFUSION_WOS); 266 267 OK(sdis_interface_ref_put(interf)); 268 OK(sdis_medium_ref_put(solid)); 269 OK(sdis_medium_ref_put(dummy)); 270 OK(sdis_scene_ref_put(scene)); 271 OK(sdis_device_ref_put(sdis)); 272 273 CHK(mem_allocated_size() == 0); 274 return 0; 275 }