test_sgas.c (12529B)
1 /* Copyright (C) 2025, 2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2025, 2026 Université de Lorraine 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #define _POSIX_C_SOURCE 200809L /* snprintf */ 18 19 #include "sgas.h" 20 21 #include <rsys/mem_allocator.h> 22 23 #include <stdio.h> 24 25 /* Indentifier of thermodynamic properties */ 26 enum { PRESSURE, TEMPERATURE, XH2O, XCO2, XCO }; 27 28 /* Define the value of a property according to its type of property, the 29 * tetrahedron to which it is attached and the time to which it is defined */ 30 #define PROP_VAL(Prop, Tetra, Time) ((double)(Prop+Tetra*100+Time*10)) 31 32 /* Vertices of the volumetric mesh of a unit cube with a lower corner in 33 * (0.0.0), and therefore its upper corner in (1.1.1). The vertices are those 34 * meshes of a cube to which is added an additional vertex in its center. So 35 * the base of the tetrahedra are the original triangles and have a shared 36 * vertex in the center of the cube. Adding a vertex to the centre of the cybe 37 * is not necessary but simplifies the volumetric mesh of the cube and the 38 * indentification of the tetraheron in which a position lies. */ 39 static const double box_vertices[9/*#vertices*/*3/*#coords per vertex*/] = { 40 0.0, 0.0, 0.0, 41 1.0, 0.0, 0.0, 42 0.0, 1.0, 0.0, 43 1.0, 1.0, 0.0, 44 0.0, 0.0, 1.0, 45 1.0, 0.0, 1.0, 46 0.0, 1.0, 1.0, 47 1.0, 1.0, 1.0, 48 0.5, 0.5, 0.5 49 }; 50 static const uint64_t box_nverts = 51 (uint64_t)(sizeof(box_vertices) / (sizeof(double)*3)); 52 53 /* The following array lists the indices toward the 3D vertices of each 54 * boundary triangle. The index 8 references the vertex at the center of the 55 * box. 56 * ,2---,3 ,2----3 57 * ,' | ,'/| ,'/| \ | Y 58 * 6----7' / | 6' / | \ | | 59 * |', | / ,1 | / ,0---,1 o--X 60 * | ',|/,' |/,' | ,' / 61 * 4----5' 4----5' Z */ 62 static const uint64_t box_indices[12/*#tetras*/*4/*#indices per tetra*/] = { 63 0, 1, 2, 8, 1, 3, 2, 8, /* -Z */ 64 0, 2, 4, 8, 2, 6, 4, 8, /* -X */ 65 4, 6, 5, 8, 6, 7, 5, 8, /* +Z */ 66 3, 5, 7, 8, 5, 3, 1, 8, /* +X */ 67 2, 7, 6, 8, 7, 2, 3, 8, /* +Y */ 68 0, 5, 1, 8, 5, 0, 4, 8 /* -Y */ 69 }; 70 static const uint64_t box_ntetras = 71 (uint64_t)(sizeof(box_indices) / (sizeof(uint64_t)*4)); 72 73 /******************************************************************************* 74 * Helper functions 75 ******************************************************************************/ 76 static void 77 write_smsh 78 (FILE* fp, 79 const double* verts, 80 const uint64_t nverts, 81 const uint64_t* ids, 82 const uint64_t ntetras) 83 { 84 const uint64_t pagesize = 8192; 85 const uint32_t dimnode = 3; 86 const uint32_t dimcell = 4; 87 88 89 #define WRITE(Item, Count) \ 90 CHK(fwrite(Item, sizeof(*Item), (Count), fp) == (Count)) 91 92 #define PADDING { \ 93 char byte = 0; \ 94 CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); \ 95 CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); /* EOF indicator */ \ 96 } (void)0 97 98 /* Header */ 99 WRITE(&pagesize, 1); 100 WRITE(&nverts, 1); 101 WRITE(&ntetras, 1); 102 WRITE(&dimnode, 1); 103 WRITE(&dimcell, 1); 104 PADDING; 105 106 /* Nodes */ 107 WRITE(verts, nverts*dimnode); 108 PADDING; 109 110 /* Cells */ 111 WRITE(ids, ntetras*dimcell); 112 PADDING; 113 114 CHK(fflush(fp) == 0); 115 116 #undef WRITE 117 #undef PADDING 118 } 119 120 static void 121 write_atrtp 122 (FILE* fp, 123 const uint64_t ntetras, 124 const uint64_t time_index) 125 { 126 const uint64_t pagesize = 8192; 127 size_t i = 9; 128 129 #define WRITE(Item, Count) \ 130 CHK(fwrite(Item, sizeof(*Item), (Count), fp) == (Count)) 131 132 #define PADDING { \ 133 char byte = 0; \ 134 CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); \ 135 CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); /* EOF indicator */ \ 136 } (void)0 137 138 /* Header */ 139 WRITE(&pagesize, 1); 140 WRITE(&ntetras, 1); 141 PADDING; 142 143 FOR_EACH(i, 0, ntetras) { 144 const double pressure = PROP_VAL(PRESSURE, i, time_index); 145 const double temperature = PROP_VAL(TEMPERATURE, i, time_index); 146 const double xH2O = PROP_VAL(XH2O, i, time_index); 147 const double xCO2 = PROP_VAL(XCO2, i, time_index); 148 const double xCO = PROP_VAL(XCO, i, time_index); 149 150 /* Value of properties not exposed by the sgas library but required by the 151 * atrtp format. Its only requirement is that it be valid with respect to 152 * the properties for which it is used. */ 153 const double dummy = 0; 154 155 WRITE(&pressure, 1); 156 WRITE(&temperature, 1); 157 WRITE(&xH2O, 1); 158 WRITE(&xCO2, 1); 159 WRITE(&xCO, 1); 160 WRITE(&dummy, 1); /* vf-soot (not used) */ 161 WRITE(&dummy, 1); /* np-soot (not used) */ 162 WRITE(&dummy, 1); /* dp-soot (not used) */ 163 } 164 PADDING; 165 } 166 167 static void 168 check_props 169 (const struct sgas_props* props, 170 const uint64_t tetra, 171 const uint64_t time) 172 { 173 CHK(props); 174 CHK(props->pressure == PROP_VAL(PRESSURE, tetra, time)); 175 CHK(props->temperature == PROP_VAL(TEMPERATURE, tetra, time)); 176 CHK(props->xH2O == PROP_VAL(XH2O, tetra, time)); 177 CHK(props->xCO2 == PROP_VAL(XCO2, tetra, time)); 178 CHK(props->xCO == PROP_VAL(XCO, tetra, time)); 179 } 180 181 static void 182 test_gas_creation(void) 183 { 184 struct sgas_create_args args = SGAS_CREATE_ARGS_DEFAULT; 185 struct sgas* gas = NULL; 186 187 const char* name_mesh = "mesh.smsh"; 188 const char* name_props = "props.atrtp"; 189 const char* name_props_list = "props_list.txt"; 190 FILE* fp_mesh = NULL; 191 FILE* fp_props = NULL; 192 FILE* fp_props_list = NULL; 193 194 /* Create the gas mesh */ 195 CHK(fp_mesh = fopen(name_mesh, "w+")); 196 write_smsh(fp_mesh, box_vertices, box_nverts, box_indices, box_ntetras); 197 198 /* Create one set of thermodynamic properties */ 199 CHK(fp_props = fopen(name_props, "w")); 200 write_atrtp(fp_props, box_ntetras, 0/*time step*/); 201 CHK(fclose(fp_props) == 0); 202 203 /* List the thermodynamic properties, actually only one */ 204 CHK(fp_props_list = fopen(name_props_list, "w+")); 205 CHK(fprintf(fp_props_list, "#time thermo_properties\n")); 206 CHK(fprintf(fp_props_list, "0 %s\n", name_props)); 207 CHK(fflush(fp_props_list) == 0); 208 209 args.mesh.name = name_mesh; 210 args.therm_props_list.name = name_props_list; 211 args.verbose = 1; 212 213 /* Check the gas creation API */ 214 CHK(sgas_create(NULL, &gas) == RES_BAD_ARG); 215 CHK(sgas_create(&args, NULL) == RES_BAD_ARG); 216 CHK(sgas_create(&args, &gas) == RES_OK); 217 218 /* Check the reference counting API */ 219 CHK(sgas_ref_get(NULL) == RES_BAD_ARG); 220 CHK(sgas_ref_get(gas) == RES_OK); 221 CHK(sgas_ref_put(NULL) == RES_BAD_ARG); 222 CHK(sgas_ref_put(gas) == RES_OK); 223 CHK(sgas_ref_put(gas) == RES_OK); 224 225 /* The mesh is missing. */ 226 args.mesh.name = NULL; 227 CHK(sgas_create(&args, &gas) == RES_BAD_ARG); 228 229 /* The list of thermodynamic properties is missing */ 230 args.mesh.name = name_mesh; 231 args.therm_props_list.name = NULL; 232 CHK(sgas_create(&args, &gas) == RES_BAD_ARG); 233 234 /* Invalid stream for the list of thermodynamic properties */ 235 args.mesh.name = name_mesh; 236 args.therm_props_list.name = NULL; 237 args.therm_props_list.stream = fp_props_list; 238 CHK(sgas_create(&args, &gas) == RES_BAD_ARG); 239 240 /* The stream for the list of properties becomes valid */ 241 rewind(fp_props_list); 242 CHK(sgas_create(&args, &gas) == RES_OK); 243 CHK(sgas_ref_put(gas) == RES_OK); 244 245 /* Invalid stream for the gas mesh */ 246 rewind(fp_props_list); 247 args.mesh.stream = fp_mesh; 248 CHK(sgas_create(&args, &gas) != RES_OK); 249 250 /* The stream for the mesh becomes valid */ 251 rewind(fp_mesh); 252 CHK(sgas_create(&args, &gas) == RES_OK); 253 CHK(sgas_ref_put(gas) == RES_OK); 254 255 CHK(fclose(fp_mesh) == 0); 256 CHK(fclose(fp_props_list) == 0); 257 } 258 259 static void 260 test_gas_query(void) 261 { 262 struct sgas_create_args args = SGAS_CREATE_ARGS_DEFAULT; 263 struct sgas_props props = SGAS_PROPS_NULL; 264 struct sgas* gas = NULL; 265 266 const char* name_mesh = "mesh.smsh"; 267 const char* name_props_list = "props_list.txt"; 268 FILE* fp_mesh = NULL; 269 FILE* fp_props_list = NULL; 270 271 double pos[3] = {0}; 272 const int ntime_steps = 4; 273 int i = 0; 274 275 /* Create the gas mesh */ 276 CHK(fp_mesh = fopen(name_mesh, "w")); 277 write_smsh(fp_mesh, box_vertices, box_nverts, box_indices, box_ntetras); 278 CHK(fclose(fp_mesh) == 0); 279 280 CHK(fp_props_list = fopen(name_props_list, "w")); 281 282 /* Write the time varying thermodynamic properties 283 * and the file that lists them */ 284 FOR_EACH(i, 0, ntime_steps) { 285 char name_props[32] = {0}; 286 FILE* fp_props = NULL; 287 288 CHK((size_t)snprintf(name_props, sizeof(name_props), "props%d.atrtp", i) 289 < sizeof(name_props)-1/*NULL byte*/); 290 291 CHK(fp_props = fopen(name_props, "w")); 292 write_atrtp(fp_props, box_ntetras, (uint64_t)i/*time step*/); 293 CHK(fclose(fp_props) == 0); 294 295 CHK(fprintf(fp_props_list, "0.%de1 %s\n", i, name_props) != 0); 296 } 297 CHK(fclose(fp_props_list) == 0); 298 299 /* Create the gas from the previous data */ 300 args.mesh.name = name_mesh; 301 args.therm_props_list.name = name_props_list; 302 args.verbose = 1; 303 CHK(sgas_create(&args, &gas) == RES_OK); 304 305 CHK(sgas_get_props(NULL, pos, 0, &props) == RES_BAD_ARG); 306 CHK(sgas_get_props(gas, NULL, 0, &props) == RES_BAD_ARG); 307 CHK(sgas_get_props(gas, pos, 0, NULL) == RES_BAD_ARG); 308 309 /* The position is outside the mesh 310 * and therefore all returned properties are null */ 311 pos[0] = pos[1] = pos[2] = -1; 312 CHK(sgas_get_props(gas, pos, 0, &props) == RES_OK); 313 CHK(props.pressure == SGAS_PROPS_NULL.pressure); 314 CHK(props.temperature == SGAS_PROPS_NULL.temperature); 315 CHK(props.xH2O == SGAS_PROPS_NULL.xH2O); 316 CHK(props.xCO2 == SGAS_PROPS_NULL.xCO2); 317 CHK(props.xCO == SGAS_PROPS_NULL.xCO); 318 319 /* Query the first tetrahedron at the first time */ 320 pos[0] = 0.25; 321 pos[1] = 0.50; 322 pos[2] = 0.25; 323 CHK(sgas_get_props(gas, pos, 0, &props) == RES_OK); 324 check_props(&props, 0/*tetra*/, 0/*time*/); 325 326 /* Query the second tetrahedron, always at the first time */ 327 pos[0] = 0.75; 328 pos[1] = 0.50; 329 pos[2] = 0.25; 330 CHK(sgas_get_props(gas, pos, 0, &props) == RES_OK); 331 check_props(&props, 1/*tetra*/, 0/*time*/); 332 333 /* Do not move but query a time beyond the limit, i.e. time is truncated at 334 * the last defined time */ 335 CHK(sgas_get_props(gas, pos, ntime_steps, &props) == RES_OK); 336 check_props(&props, 1/*tetra*/, (uint64_t)(ntime_steps-1)/*time*/); 337 338 /* Query a time near the third time step */ 339 CHK(sgas_get_props(gas, pos, 1.51, &props) == RES_OK); 340 check_props(&props, 1/*tetra*/, 2/*time*/); 341 CHK(sgas_get_props(gas, pos, 2.49, &props) == RES_OK); 342 check_props(&props, 1/*tetra*/, 2/*time*/); 343 344 /* Query a position in the 9^th tetrahedron near the second time step. 345 * 346 * Note that the position was calculated from the vertices of the tetrahedron 347 * using the ssp_ran_tetrahedron_uniform function of the Star-SamPling 348 * library, via an adhoc executable. Do not use Star-SamPling directly here 349 * avoid making it a project dependency. */ 350 pos[0] = 0.23; 351 pos[1] = 0.98; 352 pos[2] = 0.52; 353 CHK(sgas_get_props(gas, pos, 0.7, &props) == RES_OK); 354 check_props(&props, 8/*tetra*/, 1/*time*/); 355 356 /* Query properties in the same position but exactly between 2 steps of time. 357 * In this case, none of the two times is better than the other. The queried 358 * time step depends only on a choice of arbitrary implementation. Below is 359 * therefore verified that the properties returned correspond to those of the 360 * tetrahedron for one of the 2 valid time steps */ 361 CHK(sgas_get_props(gas, pos, 1.5, &props) == RES_OK); 362 if(props.pressure == PROP_VAL(PRESSURE, 8/*tetra*/, 1/*time*/)) { 363 check_props(&props, 8/*tetra*/, 1/*time*/); 364 } else if(props.pressure == PROP_VAL(PRESSURE, 8/*tetra*/, 2/*time*/)) { 365 check_props(&props, 8/*tetra*/, 2/*time*/); 366 } else { 367 CHK(0); /* Unexpected properties values */ 368 } 369 370 CHK(sgas_ref_put(gas) == RES_OK); 371 } 372 373 /******************************************************************************* 374 * The test 375 ******************************************************************************/ 376 int 377 main(void) 378 { 379 test_gas_creation(); 380 test_gas_query(); 381 CHK(mem_allocated_size() == 0); 382 return 0; 383 }