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