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