htcp

test_htcp_load_from_file.c (6242B)

---

 /* Copyright (C) 2018, 2020-2023, 2025, 2025 |Méso|Star> (contact@meso-star.com)
  * Copyright (C) 2018 Centre National de la Recherche Scientifique
  * Copyright (C) 2018 Université Paul Sabatier
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #define _POSIX_C_SOURCE 200112L /* snprintf support */
 
 #include "htcp.h"
 #include "test_htcp_utils.h"
 
 #include <rsys/cstr.h>
 #include <rsys/double3.h>
 #include <rsys/math.h>
 #include <string.h>
 
 /* Potential temperature to temperature
  *  T = THT * (PABST/P0)^(R/(M_DA*Cp0)) */
 static FINLINE double
 THT_to_T(const double THT, const size_t i /*Id of THT*/, const void* ctx)
 {
   const double* PABST = ctx;
   const double P0 = 101325;
 #if 0
   const double R = 8.3144598; /* In kg.m^2.s^-2.mol^-1.K */
   const double M_DA = 28.9644E-3; /* In kg.mol^-1 */
   const double Cp0 = 1.005E+3; /* In J/kg/K */
 #endif
   const double exponent = 0.28562974993698609; /* R/M_DA*cp0 */
   return THT * pow(PABST[i]/P0, exponent);
 }
 
 static void
 check_descriptor
   (const struct htcp_desc* desc,
    const char* path, /* Path where the reference values are */
    const char* basename) /* Basename of the reference files */
 {
   char buf[128];
   double nx, ny, nz, ntimes;
   FILE* fp = NULL;
   size_t i;
 
   printf("Check the descriptor\n");
 
   CHK(desc && path && basename);
 
   i = (size_t)snprintf(buf, sizeof(buf), "%s/%s_desc", path, basename);
   CHK(i < sizeof(buf));
   CHK(fp = fopen(buf, "r"));
   CHK(fscanf(fp, "%lg %lg %lg %lg", &nx, &ny, &nz, &ntimes) == 4);
   CHK(fscanf(fp, "%*g") == EOF);
   CHK(fclose(fp) == 0);
 
   CHK(desc->spatial_definition[0] == nx);
   CHK(desc->spatial_definition[1] == ny);
   CHK(desc->spatial_definition[2] == nz);
   CHK(desc->time_definition == ntimes);
 }
 
 static void
 check_variable
   (const double* loaded_data,
    const size_t ndata,
    const char* var, /* Name of the variable to check */
    const char* path, /* Path where the reference values are */
    const char* basename, /* Basename of the reference files */
    double (*convert)(const double, const size_t, const void*), /* May be NULL */
    const void* ctx)
 {
   char buf[128];
   FILE* fp = NULL;
   size_t i;
 
   printf("Check the %s variable\n", var);
 
   CHK(loaded_data && ndata && var && path && basename);
 
   i = (size_t)snprintf(buf, sizeof(buf), "%s/%s_%s", path, basename, var);
   CHK(i < sizeof(buf));
   CHK(fp = fopen(buf, "r"));
   FOR_EACH(i, 0, ndata) {
     double val;
     CHK(fscanf(fp, "%lg", &val) == 1);
     if(convert) val = convert(val, i, ctx);
     CHK(eq_eps(val, loaded_data[i], loaded_data[i]*1.e-6));
   }
   CHK(fscanf(fp, "%*g") == EOF);
   CHK(fclose(fp) == 0);
 }
 
 static void
 check_misc(const struct htcp* htcp)
 {
   struct htcp_desc desc;
   double low[3], upp[3];
   size_t x, y, z;
 
   printf("Check miscellaneous\n");
 
   CHK(htcp_get_desc(htcp, &desc) == RES_OK);
   CHK(desc.irregular_z || desc.coord_z == NULL);
 
   FOR_EACH(x, 0, desc.spatial_definition[0]) {
     low[0] = (double)x * desc.vxsz_x;
     upp[0] = low[0] + desc.vxsz_x;
     FOR_EACH(y, 0, desc.spatial_definition[1]) {
       low[1] = (double)y * desc.vxsz_y;
       upp[1] = low[1] + desc.vxsz_y;
       FOR_EACH(z, 0, desc.spatial_definition[2]) {
         double vox_low[3];
         double vox_upp[3];
         if(!desc.irregular_z) {
           low[2] = (double)z * desc.vxsz_z[0];
           upp[2] = low[2] + desc.vxsz_z[0];
         } else {
           size_t i;
           low[2] = desc.lower[2];
           FOR_EACH(i, 0, z) low[2] += desc.vxsz_z[i];
           CHK(eq_eps(low[2], desc.coord_z[z], 1.e-6));
           upp[2] = low[2] + desc.vxsz_z[z];
         }
         htcp_desc_get_voxel_aabb(&desc, x, y, z, vox_low, vox_upp);
         CHK(d3_eq_eps(vox_low, low, 1.e-6));
         CHK(d3_eq_eps(vox_upp, upp, 1.e-6));
       }
     }
   }
 
   upp[0] = desc.lower[0] + desc.vxsz_x * (double)desc.spatial_definition[0];
   upp[1] = desc.lower[1] + desc.vxsz_y * (double)desc.spatial_definition[1];
   if(!desc.irregular_z) {
     upp[2] = desc.lower[2] + desc.vxsz_z[0] * (double)desc.spatial_definition[2];
   } else {
     upp[2] = desc.lower[2];
     FOR_EACH(z, 0, desc.spatial_definition[2]) {
       upp[2] += desc.vxsz_z[z];
     }
   }
   CHK(d3_eq_eps(upp, desc.upper, 1.e-6));
 }
 
 int
 main(int argc, char** argv)
 {
   struct htcp* htcp = NULL;
   struct htcp_desc desc = HTCP_DESC_NULL;
   double fp_to_meter = 1.0;
   char* filename = NULL;
   char* path = NULL;
   char* base = NULL;
   char* p = NULL;
   size_t n;
 
   if(argc < 3) {
     fprintf(stderr, "Usage: %s <htcp> <ref-data-path> [fp-to-meter]\n", argv[0]);
     return -1;
   }
   if(argc > 3) {
     CHK(cstr_to_double(argv[3], &fp_to_meter) == RES_OK);
   }
 
   filename = argv[1];
   path = argv[2];
 
   CHK(htcp_create(NULL, &mem_default_allocator, 1, &htcp) == RES_OK);
   CHK(htcp_load(htcp, filename) == RES_OK);
 
   CHK(htcp_get_desc(htcp, &desc) == RES_OK);
 
   /* Compute the basename of the ref file from the submitted htcp file */
   base = filename;
   p = strrchr(filename, '/');
   if(p) base = p+1;
   p = strrchr(filename, '.');
   if(p) *p = '\0';
 
   n = desc.spatial_definition[0]
     * desc.spatial_definition[1]
     * desc.spatial_definition[2]
     * desc.time_definition;
 
   printf("Irregular Z: %i\n", desc.irregular_z);
   check_descriptor(&desc, path, base);
   check_misc(htcp);
   check_variable(desc.RVT, n, "RVT", path, base, NULL, NULL);
   check_variable(desc.RCT, n, "RCT", path, base, NULL, NULL);
   check_variable(desc.PABST, n, "PABST", path, base, NULL, NULL);
   check_variable(desc.T, n, "THT", path, base, THT_to_T, desc.PABST);
 
   CHK(htcp_ref_put(htcp) == RES_OK);
   check_memory_allocator(&mem_default_allocator);
   CHK(mem_allocated_size() == 0);
   return 0;
 }