test_sdis_solve_probe_2d.c (7414B)
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 with unknown temperature. The 23 * surrounding fluid has a fixed constant temperature. 24 * 25 * (1,1) 26 * +-------+ _\ 27 * | | / / 28 * | | \__/ 300K 29 * | | 30 * +-------+ 31 * (0,0) 32 */ 33 34 /******************************************************************************* 35 * Geometry 36 ******************************************************************************/ 37 struct context { 38 const double* positions; 39 const size_t* indices; 40 struct sdis_interface* interf; 41 }; 42 43 static void 44 get_indices(const size_t iseg, size_t ids[2], void* context) 45 { 46 struct context* ctx = context; 47 ids[0] = ctx->indices[iseg*2+0]; 48 ids[1] = ctx->indices[iseg*2+1]; 49 } 50 51 static void 52 get_position(const size_t ivert, double pos[2], void* context) 53 { 54 struct context* ctx = context; 55 pos[0] = ctx->positions[ivert*2+0]; 56 pos[1] = ctx->positions[ivert*2+1]; 57 } 58 59 static void 60 get_interface(const size_t iseg, struct sdis_interface** bound, void* context) 61 { 62 struct context* ctx = context; 63 (void)iseg; 64 *bound = ctx->interf; 65 } 66 67 /******************************************************************************* 68 * Media & interface 69 ******************************************************************************/ 70 struct fluid { 71 double temperature; 72 }; 73 74 static double 75 fluid_get_temperature 76 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 77 { 78 CHK(data != NULL && vtx != NULL); 79 return ((const struct fluid*)sdis_data_cget(data))->temperature; 80 } 81 82 static double 83 solid_get_calorific_capacity 84 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 85 { 86 (void)vtx, (void)data; 87 return 1.0; 88 } 89 90 static double 91 solid_get_thermal_conductivity 92 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 93 { 94 (void)vtx, (void)data; 95 return 0.1; 96 } 97 98 static double 99 solid_get_volumic_mass 100 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 101 { 102 (void)vtx, (void)data; 103 return 1.0; 104 } 105 106 static double 107 solid_get_delta 108 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 109 { 110 (void)vtx, (void)data; 111 return 1.0/20.0; 112 } 113 114 static double 115 solid_get_temperature 116 (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) 117 { 118 (void)vtx, (void)data; 119 return SDIS_TEMPERATURE_NONE; 120 } 121 122 static double 123 interface_get_convection_coef 124 (const struct sdis_interface_fragment* frag, struct sdis_data* data) 125 { 126 (void)frag, (void)data; 127 return 0.5; 128 } 129 130 /******************************************************************************* 131 * Main test 132 ******************************************************************************/ 133 int 134 main(int argc, char** argv) 135 { 136 struct sdis_mc T = SDIS_MC_NULL; 137 struct sdis_mc time = SDIS_MC_NULL; 138 struct sdis_device* dev = NULL; 139 struct sdis_medium* solid = NULL; 140 struct sdis_medium* fluid = NULL; 141 struct sdis_interface* interf = NULL; 142 struct sdis_scene* scn = NULL; 143 struct sdis_data* data = NULL; 144 struct sdis_estimator* estimator = NULL; 145 struct sdis_estimator* estimator2 = NULL; 146 struct sdis_green_function* green = NULL; 147 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 148 struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; 149 struct sdis_solid_shader solid_shader = DUMMY_SOLID_SHADER; 150 struct sdis_interface_shader interface_shader = DUMMY_INTERFACE_SHADER; 151 struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 152 struct context ctx; 153 struct fluid* fluid_param; 154 double ref; 155 const size_t N = 1000; 156 size_t nreals; 157 size_t nfails; 158 (void)argc, (void)argv; 159 160 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 161 162 /* Create the fluid medium */ 163 OK(sdis_data_create 164 (dev, sizeof(struct fluid), ALIGNOF(struct fluid), NULL, &data)); 165 fluid_param = sdis_data_get(data); 166 fluid_param->temperature = 300; 167 fluid_shader.temperature = fluid_get_temperature; 168 OK(sdis_fluid_create(dev, &fluid_shader, data, &fluid)); 169 OK(sdis_data_ref_put(data)); 170 171 /* Create the solid medium */ 172 solid_shader.calorific_capacity = solid_get_calorific_capacity; 173 solid_shader.thermal_conductivity = solid_get_thermal_conductivity; 174 solid_shader.volumic_mass = solid_get_volumic_mass; 175 solid_shader.delta = solid_get_delta; 176 solid_shader.temperature = solid_get_temperature; 177 OK(sdis_solid_create(dev, &solid_shader, NULL, &solid)); 178 179 /* Create the solid/fluid interface */ 180 interface_shader.convection_coef = interface_get_convection_coef; 181 interface_shader.front = SDIS_INTERFACE_SIDE_SHADER_NULL; 182 interface_shader.back = SDIS_INTERFACE_SIDE_SHADER_NULL; 183 OK(sdis_interface_create 184 (dev, solid, fluid, &interface_shader, NULL, &interf)); 185 186 /* Release the media */ 187 OK(sdis_medium_ref_put(solid)); 188 OK(sdis_medium_ref_put(fluid)); 189 190 /* Create the scene */ 191 ctx.positions = square_vertices; 192 ctx.indices = square_indices; 193 ctx.interf = interf; 194 scn_args.get_indices = get_indices; 195 scn_args.get_interface = get_interface; 196 scn_args.get_position = get_position; 197 scn_args.nprimitives = square_nsegments; 198 scn_args.nvertices = square_nvertices; 199 scn_args.context = &ctx; 200 OK(sdis_scene_2d_create(dev, &scn_args, &scn)); 201 202 OK(sdis_interface_ref_put(interf)); 203 204 /* Test the solver */ 205 solve_args.nrealisations = N; 206 solve_args.position[0] = 0.5; 207 solve_args.position[1] = 0.5; 208 solve_args.time_range[0] = INF; 209 solve_args.time_range[1] = INF; 210 211 OK(sdis_solve_probe(scn, &solve_args, &estimator)); 212 OK(sdis_estimator_get_realisation_count(estimator, &nreals)); 213 OK(sdis_estimator_get_failure_count(estimator, &nfails)); 214 OK(sdis_estimator_get_temperature(estimator, &T)); 215 OK(sdis_estimator_get_realisation_time(estimator, &time)); 216 217 ref = 300; 218 printf("Temperature at (%g, %g) = %g ~ %g +/- %g\n", 219 SPLIT2(solve_args.position), ref, T.E, T.SE); 220 printf("Time per realisation (in usec) = %g +/- %g\n", time.E, time.SE); 221 printf("#failures = %lu/%lu\n", (unsigned long)nfails, (unsigned long)N); 222 223 CHK(nfails + nreals == N); 224 CHK(nfails < N/1000); 225 CHK(eq_eps(T.E, ref, T.SE)); 226 227 OK(sdis_solve_probe_green_function(scn, &solve_args, &green)); 228 OK(sdis_green_function_solve(green, &estimator2)); 229 check_green_function(green); 230 check_estimator_eq(estimator, estimator2); 231 232 OK(sdis_estimator_ref_put(estimator)); 233 OK(sdis_estimator_ref_put(estimator2)); 234 OK(sdis_green_function_ref_put(green)); 235 236 /* The external fluid cannot have an unknown temperature */ 237 fluid_param->temperature = SDIS_TEMPERATURE_NONE; 238 239 BA(sdis_solve_probe(scn, &solve_args, &estimator)); 240 241 OK(sdis_scene_ref_put(scn)); 242 OK(sdis_device_ref_put(dev)); 243 244 CHK(mem_allocated_size() == 0); 245 return 0; 246 }