rnatm

commit d027c2ef0473a7a7dbd7e5a3fc8d918ea395b66a
parent 4c45085fcbd988b9d4abe19dee7e8f22997e099d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri,  3 Mar 2023 09:28:15 +0100

Fix the rnatm_sammple_component function

The cumulative radiative coefficients per component were calculated from
the components that had a tetrahedron at the position queried whereas
during sampling we assumed that it was constructed from all components.
This could only work for configurations where all components shared the
same volumetric mesh. This commit solves this problem by calculating the
cumulative for all components, assuming a zero radiative coefficient
when a component has no tetrahedron at the given position.

Diffstat:
M
src/rnatm_properties.c
 | 
60
+++++++++++++++++++++++++++++++++---------------------------

1 file changed, 33 insertions(+), 27 deletions(-)
diff --git a/src/rnatm_properties.c b/src/rnatm_properties.c
@@ -994,19 +994,18 @@ compute_unnormalized_cumulative_radcoef
    const size_t iband,
    const size_t iquad,
    float cumulative[RNATM_MAX_COMPONENTS_COUNT],
-   size_t* cumulative_sz,
    /* For debug */
    const double k_min,
    const double k_max)
 {
   struct rnatm_cell_get_radcoef_args cell_args = RNATM_CELL_GET_RADCOEF_ARGS_NULL;
   size_t cpnt = RNATM_GAS;
-  size_t naerosols = 0;
   size_t icumul = 0;
+  size_t naerosols = 0;
   float k = 0;
   res_T res = RES_OK;
   ASSERT(atm && cells && (unsigned)radcoef < RNATM_RADCOEFS_COUNT__);
-  ASSERT(cumulative && cumulative_sz);
+  ASSERT(cumulative);
   (void)k_min;
 
   naerosols = darray_aerosol_size_get(&atm->aerosols);
@@ -1019,26 +1018,26 @@ compute_unnormalized_cumulative_radcoef
   cell_args.k_max = k_max; /* For Debug */
 
   do {
-    double per_cell_k;
 
     cell_args.cell = cells[cpnt+1];
 
-    /* This component does not exist here */
-    if(SUVM_PRIMITIVE_NONE(&cell_args.cell.prim)) continue;
+    /* Test if the component exists here */
+    if(!SUVM_PRIMITIVE_NONE(&cell_args.cell.prim)) {
+      double per_cell_k;
 
-    /* Add the component's contribution to the radiative coefficient */
-    res = rnatm_cell_get_radcoef(atm, &cell_args, &per_cell_k);
-    if(res != RES_OK) goto error;
+      /* Add the component's contribution to the radiative coefficient */
+      res = rnatm_cell_get_radcoef(atm, &cell_args, &per_cell_k);
+      if(res != RES_OK) goto error;
+      k += (float)per_cell_k;
+    }
 
-    k += (float)per_cell_k;
+    /* Update the cumulative */
     cumulative[icumul] = k;
     ++icumul;
-
   } while(++cpnt < naerosols);
-  *cumulative_sz = icumul;
 
-  ASSERT(!icumul || (float)k_min <= cumulative[icumul-1]);
-  ASSERT(!icumul || (float)k_max >= cumulative[icumul-1]);
+  ASSERT(cumulative[icumul-1] == 0 || (float)k_min <= cumulative[icumul-1]);
+  ASSERT(cumulative[icumul-1] == 0 || (float)k_max >= cumulative[icumul-1]);
 
 exit:
   return res;
@@ -1056,7 +1055,7 @@ rnatm_get_radcoef
    double* out_k)
 {
   float cumul[RNATM_MAX_COMPONENTS_COUNT];
-  size_t cumul_sz;
+  size_t ncpnts;
   double k = 0;
   res_T res = RES_OK;
 
@@ -1064,18 +1063,20 @@ rnatm_get_radcoef
   res = check_rnatm_get_radcoef_args(atm, args);
   if(res != RES_OK) goto error;
 
+  ncpnts = 1/*gas*/ + darray_aerosol_size_get(&atm->aerosols);
+  ASSERT(ncpnts <= RNATM_MAX_COMPONENTS_COUNT);
+
   /* Calculate the cumulative (unnormalized) of radiative coefficients. Its last
    * entry is the sum of the radiative coefficients of each component, which is
    * the atmospheric radiative coefficient to be returned. */
   res = compute_unnormalized_cumulative_radcoef(atm, args->radcoef,
-    args->cells, args->iband, args->iquad, cumul, &cumul_sz, args->k_min,
-    args->k_max);
+    args->cells, args->iband, args->iquad, cumul, args->k_min, args->k_max);
   if(res != RES_OK) goto error;
 
-  if(cumul_sz == 0) {
+  if(cumul[ncpnts-1] == 0) {
     k = 0; /* No atmospheric data */
   } else {
-    k = cumul[cumul_sz-1];
+    k = cumul[ncpnts-1];
     ASSERT(args->k_min <= k && k <= args->k_max);
   }
 
@@ -1094,7 +1095,7 @@ rnatm_sample_component
 {
   float cumul[RNATM_MAX_COMPONENTS_COUNT];
   float norm;
-  size_t cumul_sz;
+  size_t ncpnts;
   size_t i;
   res_T res = RES_OK;
 
@@ -1102,28 +1103,33 @@ rnatm_sample_component
   res = check_rnatm_sample_component_args(atm, args);
   if(res != RES_OK) goto error;
 
+  ncpnts = 1/*gas*/ + darray_aerosol_size_get(&atm->aerosols);
+  ASSERT(ncpnts <= RNATM_MAX_COMPONENTS_COUNT);
+
   /* Discard the calculation of the cumulative if there is only gas */
-  if(!darray_aerosol_size_get(&atm->aerosols)) {
+  if(ncpnts == 1) {
+    ASSERT(!SUVM_PRIMITIVE_NONE(&args->cells[0].prim));
     *cpnt = RNATM_GAS;
     goto exit;
   }
 
   res = compute_unnormalized_cumulative_radcoef(atm, args->radcoef, args->cells,
-    args->iband, args->iquad, cumul, &cumul_sz, -DBL_MAX, DBL_MAX);
+    args->iband, args->iquad, cumul, -DBL_MAX, DBL_MAX);
   if(res != RES_OK) goto error;
-  ASSERT(cumul_sz >= 1);
+  ASSERT(cumul[ncpnts-1] > 0);
 
   /* Normalize the cumulative */
-  norm = cumul[cumul_sz-1];
-  FOR_EACH(i, 0, cumul_sz) cumul[i] /= norm;
-  cumul[cumul_sz-1] = 1.f; /* Handle precision issues */
+  norm = cumul[ncpnts-1];
+  FOR_EACH(i, 0, ncpnts) cumul[i] /= norm;
+  cumul[ncpnts-1] = 1.f; /* Handle precision issues */
 
   /* Use a simple linear search to sample the component since there are
    * usually very few aerosols to consider */
-  FOR_EACH(i, 0, cumul_sz) if(args->r < cumul[i]) break;
+  FOR_EACH(i, 0, ncpnts) if(args->r < cumul[i]) break;
 
   *cpnt = i-1;
   ASSERT(*cpnt == RNATM_GAS || *cpnt < darray_aerosol_size_get(&atm->aerosols));
+  ASSERT(!SUVM_PRIMITIVE_NONE(&args->cells[i].prim));
 
 exit:
   return res;