test_sdis_flux_with_h.c (9250B)
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 configuration is a rectangle whose the conductivity is lambda and its 21 * temperature is unknown. Its top and bottom boundaries rectangle are 22 * adiabatics while its left and right ones have a fixed fluxes of phi1 and 23 * phi2 respectively. The left boundary has also a convective exchange with the 24 * surrounding fluid whose temperature is Text. At stationnary, the 25 * temperature at a given position into the solid rectangle is: 26 * 27 * T(x) = phi2/lambda*x + (Text + (phi1 + phi2)/h) 28 * 29 * with h the convective coefficient on the left boundary 30 * 31 * 32 * Text ///// (0.2,0.5) 33 * +---------+ 34 * | | 35 * h _\ |--> <--| 36 * / / phi1-> <-phi2 37 * \__/ |--> <--| 38 * | | 39 * +---------+ 40 * (0,0) /////// 41 */ 42 43 #define LAMBDA 25.0 44 #define RHO 7500.0 45 #define CP 500.0 46 #define DELTA 0.01 47 #define Text 373.15 48 #define PHI1 1000.0 49 #define PHI2 5000.0 50 #define H 10 51 52 #define N 10000 /* #realisations */ 53 54 /******************************************************************************* 55 * Geometry 56 ******************************************************************************/ 57 static void 58 get_position(const size_t ivert, double pos[2], void* ctx) 59 { 60 square_get_position(ivert, pos, ctx); 61 pos[0] *= 0.2; 62 pos[1] *= 0.5; 63 } 64 65 /******************************************************************************* 66 * Media 67 ******************************************************************************/ 68 static double 69 solid_get_calorific_capacity 70 (const struct sdis_rwalk_vertex* vtx, 71 struct sdis_data* data) 72 { 73 (void)data, (void)vtx; 74 return CP; 75 } 76 77 static double 78 solid_get_thermal_conductivity 79 (const struct sdis_rwalk_vertex* vtx, 80 struct sdis_data* data) 81 { 82 (void)data, (void)vtx; 83 return LAMBDA; 84 } 85 86 static double 87 solid_get_volumic_mass 88 (const struct sdis_rwalk_vertex* vtx, 89 struct sdis_data* data) 90 { 91 (void)data, (void)vtx; 92 return RHO; 93 } 94 95 static double 96 solid_get_delta 97 (const struct sdis_rwalk_vertex* vtx, 98 struct sdis_data* data) 99 { 100 (void)data, (void)vtx; 101 return DELTA; 102 } 103 104 static double 105 solid_get_temperature 106 (const struct sdis_rwalk_vertex* vtx, 107 struct sdis_data* data) 108 { 109 (void)data, (void)vtx; 110 return SDIS_TEMPERATURE_NONE; 111 } 112 113 static double 114 fluid_get_temperature 115 (const struct sdis_rwalk_vertex* vtx, 116 struct sdis_data* data) 117 { 118 (void)data, (void)vtx; 119 return Text; 120 } 121 122 static struct sdis_medium* 123 create_solid(struct sdis_device* dev) 124 { 125 struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL; 126 struct sdis_medium* solid = NULL; 127 128 /* Create the solid_medium */ 129 shader.calorific_capacity = solid_get_calorific_capacity; 130 shader.thermal_conductivity = solid_get_thermal_conductivity; 131 shader.volumic_mass = solid_get_volumic_mass; 132 shader.delta = solid_get_delta; 133 shader.temperature = solid_get_temperature; 134 OK(sdis_solid_create(dev, &shader, NULL, &solid)); 135 return solid; 136 } 137 138 static struct sdis_medium* 139 create_fluid(struct sdis_device* dev) 140 { 141 struct sdis_fluid_shader shader = DUMMY_FLUID_SHADER; 142 struct sdis_medium* fluid = NULL; 143 144 /* Create the solid_medium */ 145 shader.temperature = fluid_get_temperature; 146 OK(sdis_fluid_create(dev, &shader, NULL, &fluid)); 147 return fluid; 148 } 149 150 /******************************************************************************* 151 * Interfaces 152 ******************************************************************************/ 153 struct interf { 154 double h; 155 double phi; 156 }; 157 158 static double 159 interface_get_convection_coef 160 (const struct sdis_interface_fragment* frag, 161 struct sdis_data* data) 162 { 163 const struct interf* interf = sdis_data_cget(data); 164 (void)frag; 165 return interf->h; 166 } 167 168 static double 169 interface_get_flux 170 (const struct sdis_interface_fragment* frag, 171 struct sdis_data* data) 172 { 173 const struct interf* interf = sdis_data_cget(data); 174 (void)frag; 175 return interf->phi; 176 } 177 178 static struct sdis_interface* 179 interface_create 180 (struct sdis_device* sdis, 181 struct sdis_medium* front, 182 struct sdis_medium* back, 183 const double h, 184 const double phi) 185 { 186 struct sdis_interface* interf = NULL; 187 struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL; 188 struct sdis_data* data = NULL; 189 struct interf* props = NULL; 190 191 shader.front.flux = interface_get_flux; 192 shader.convection_coef = interface_get_convection_coef; 193 shader.convection_coef_upper_bound = h; 194 195 OK(sdis_data_create(sdis, sizeof(struct interf), 16, NULL, &data)); 196 props = sdis_data_get(data); 197 props->h = h; 198 props->phi = phi; 199 OK(sdis_interface_create(sdis, front, back, &shader, data, &interf)); 200 OK(sdis_data_ref_put(data)); 201 return interf; 202 } 203 204 /******************************************************************************* 205 * Test 206 ******************************************************************************/ 207 int 208 main(int argc, char** argv) 209 { 210 struct sdis_device* dev = NULL; 211 struct sdis_estimator* estimator = NULL; 212 struct sdis_interface* interf_left = NULL; 213 struct sdis_interface* interf_right = NULL; 214 struct sdis_interface* interf_adiab = NULL; 215 struct sdis_interface* interfaces[4]; 216 struct sdis_medium* solid = NULL; 217 struct sdis_medium* fluid = NULL; 218 struct sdis_scene* scn = NULL; 219 struct sdis_mc mc = SDIS_MC_NULL; 220 struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; 221 struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; 222 struct sdis_solve_boundary_flux_args flux_args = 223 SDIS_SOLVE_BOUNDARY_FLUX_ARGS_DEFAULT; 224 struct sdis_solve_probe_boundary_flux_args probe_flux_args = 225 SDIS_SOLVE_PROBE_BOUNDARY_FLUX_ARGS_DEFAULT; 226 double Tref = 0; 227 size_t prim = 0; 228 (void)argc, (void)argv; 229 230 OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); 231 232 solid = create_solid(dev); 233 fluid = create_fluid(dev); 234 235 interf_left = interface_create(dev, solid, fluid, H, PHI1); 236 interf_right = interface_create(dev, solid, fluid, 0, PHI2); 237 interf_adiab = interface_create(dev, solid, fluid, 0, SDIS_FLUX_NONE); 238 interfaces[0] = interf_adiab; /* Bottom */ 239 interfaces[1] = interf_left; 240 interfaces[2] = interf_adiab; /* Top */ 241 interfaces[3] = interf_right; 242 243 scn_args.get_indices = square_get_indices; 244 scn_args.get_interface = square_get_interface; 245 scn_args.get_position = get_position; 246 scn_args.nprimitives = square_nsegments; 247 scn_args.nvertices = square_nvertices; 248 scn_args.context = interfaces; 249 OK(sdis_scene_2d_create(dev, &scn_args, &scn)); 250 251 prim = 1; /* Left */ 252 flux_args.nrealisations = N; 253 flux_args.nprimitives = 1; 254 flux_args.primitives = &prim; 255 OK(sdis_solve_boundary_flux(scn, &flux_args, &estimator)); 256 OK(sdis_estimator_get_convective_flux(estimator, &mc)); 257 printf("Convective flux (left side) ~ %g W/m² +/- %g\n", mc.E, mc.SE); 258 OK(sdis_estimator_get_imposed_flux(estimator, &mc)); 259 printf("Imposed flux (left side) ~ %g W/m² +/- %g\n", mc.E, mc.SE); 260 OK(sdis_estimator_get_total_flux(estimator, &mc)); 261 printf("Total flux (left side) ~ %g W/m² +/- %g\n", mc.E, mc.SE); 262 CHK(eq_eps(-PHI2, mc.E, 3*mc.SE)); 263 OK(sdis_estimator_ref_put(estimator)); 264 265 probe_flux_args.nrealisations = N; 266 probe_flux_args.iprim = 1; /* Left */ 267 probe_flux_args.uv[0] = 0.5; 268 OK(sdis_solve_probe_boundary_flux(scn, &probe_flux_args, &estimator)); 269 OK(sdis_estimator_get_convective_flux(estimator, &mc)); 270 printf("Convective flux (probe on left side) ~ %g W/m² +/- %g\n", mc.E, mc.SE); 271 OK(sdis_estimator_get_imposed_flux(estimator, &mc)); 272 printf("Imposed flux (probe on left side) ~ %g W/m² +/- %g\n", mc.E, mc.SE); 273 OK(sdis_estimator_get_total_flux(estimator, &mc)); 274 printf("Total flux (probe on left side) ~ %g W/m² +/- %g\n", mc.E, mc.SE); 275 CHK(eq_eps(-PHI2, mc.E, 3*mc.SE)); 276 OK(sdis_estimator_ref_put(estimator)); 277 278 probe_args.nrealisations = N; 279 probe_args.position[0] = 0.2; 280 probe_args.position[1] = 0.25; 281 OK(sdis_solve_probe(scn, &probe_args, &estimator)); 282 OK(sdis_estimator_get_temperature(estimator, &mc)); 283 OK(sdis_estimator_ref_put(estimator)); 284 285 Tref = PHI2/LAMBDA*probe_args.position[0] + (Text + (PHI1 + PHI2)/H); 286 printf("Temperature at %g = %g ~ %g +/- %g\n", 287 probe_args.position[0], Tref, mc.E, mc.SE); 288 CHK(eq_eps(mc.E, Tref, 3*mc.SE)); 289 290 OK(sdis_device_ref_put(dev)); 291 OK(sdis_interface_ref_put(interf_left)); 292 OK(sdis_interface_ref_put(interf_right)); 293 OK(sdis_interface_ref_put(interf_adiab)); 294 OK(sdis_medium_ref_put(solid)); 295 OK(sdis_medium_ref_put(fluid)); 296 OK(sdis_scene_ref_put(scn)); 297 CHK(mem_allocated_size() == 0); 298 return 0; 299 }