atrstm

commit f8454574d3116230c7b82ff64ce78fe2f12363a2
parent 319cd70f013dfc74bbfeb1c1ec2f5d30a17b1e58
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 26 May 2021 12:49:06 +0200

Handle numerical uncertainty at radcoef interpolation

Diffstat:
M
src/atrstm_radcoefs.c
 | 
35
++++++++++++++++++++++++-----------
M
src/atrstm_radcoefs_simd4.c
 | 
12
++++++++++++

2 files changed, 36 insertions(+), 11 deletions(-)
diff --git a/src/atrstm_radcoefs.c b/src/atrstm_radcoefs.c
@@ -55,17 +55,6 @@ atrstm_fetch_radcoefs
     } else {
       res = fetch_radcoefs_simd4(atrstm, args, radcoefs);
       if(res != RES_OK) goto error;
-      #ifndef NDEBUG
-      {
-        atrstm_radcoefs_T radcoefs_scalar;
-        int i;
-
-        ASSERT(fetch_radcoefs(atrstm, args, radcoefs_scalar) == RES_OK);
-        FOR_EACH(i, 0, ATRSTM_RADCOEFS_COUNT__) {
-          ASSERT(eq_eps(radcoefs[i], radcoefs_scalar[i], radcoefs[i]*1e-4));
-        }
-      }
-      #endif /* !NDEBUG */
     }
   }
   #endif /* ATRSTM_USE_SIMD */
@@ -141,29 +130,53 @@ fetch_radcoefs
 
   /* Linearly interpolate the node radiative properties */
   if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) {
+    const double prim_ka_min = MMIN
+      (MMIN(prim_radcoefs[0].ka, prim_radcoefs[1].ka),
+       MMIN(prim_radcoefs[2].ka, prim_radcoefs[3].ka));
+    const double prim_ka_max = MMAX
+      (MMAX(prim_radcoefs[0].ka, prim_radcoefs[1].ka),
+       MMAX(prim_radcoefs[2].ka, prim_radcoefs[3].ka));
     radcoefs[ATRSTM_RADCOEF_Ka] =
       prim_radcoefs[0].ka * args->bcoords[0]
     + prim_radcoefs[1].ka * args->bcoords[1]
     + prim_radcoefs[2].ka * args->bcoords[2]
     + prim_radcoefs[3].ka * args->bcoords[3];
+    radcoefs[ATRSTM_RADCOEF_Ka] = /* Handle numerical uncertainty */
+      CLAMP(radcoefs[ATRSTM_RADCOEF_Ka], prim_ka_min, prim_ka_max);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] >= args->k_min[ATRSTM_RADCOEF_Ka]);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] <= args->k_max[ATRSTM_RADCOEF_Ka]);
   }
   if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) {
+    const double prim_ks_min = MMIN
+      (MMIN(prim_radcoefs[0].ks, prim_radcoefs[1].ks),
+       MMIN(prim_radcoefs[2].ks, prim_radcoefs[3].ks));
+    const double prim_ks_max = MMAX
+      (MMAX(prim_radcoefs[0].ks, prim_radcoefs[1].ks),
+       MMAX(prim_radcoefs[2].ks, prim_radcoefs[3].ks));
     radcoefs[ATRSTM_RADCOEF_Ks] =
       prim_radcoefs[0].ks * args->bcoords[0]
     + prim_radcoefs[1].ks * args->bcoords[1]
     + prim_radcoefs[2].ks * args->bcoords[2]
     + prim_radcoefs[3].ks * args->bcoords[3];
+    radcoefs[ATRSTM_RADCOEF_Ks] = /* Handle numerical uncertainty */
+      CLAMP(radcoefs[ATRSTM_RADCOEF_Ks], prim_ks_min, prim_ks_max);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] >= args->k_min[ATRSTM_RADCOEF_Ks]);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] <= args->k_max[ATRSTM_RADCOEF_Ks]);
   }
   if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Kext) {
+    const double prim_kext_min = MMIN
+      (MMIN(prim_radcoefs[0].kext, prim_radcoefs[1].kext),
+       MMIN(prim_radcoefs[2].kext, prim_radcoefs[3].kext));
+    const double prim_kext_max = MMAX
+      (MMAX(prim_radcoefs[0].kext, prim_radcoefs[1].kext),
+       MMAX(prim_radcoefs[2].kext, prim_radcoefs[3].kext));
     radcoefs[ATRSTM_RADCOEF_Kext] =
       prim_radcoefs[0].kext * args->bcoords[0]
     + prim_radcoefs[1].kext * args->bcoords[1]
     + prim_radcoefs[2].kext * args->bcoords[2]
     + prim_radcoefs[3].kext * args->bcoords[3];
+    radcoefs[ATRSTM_RADCOEF_Kext] = /* Handle numerical uncertainty */
+      CLAMP(radcoefs[ATRSTM_RADCOEF_Kext], prim_kext_min, prim_kext_max);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] >= args->k_min[ATRSTM_RADCOEF_Kext]);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] <= args->k_max[ATRSTM_RADCOEF_Kext]);
   }
diff --git a/src/atrstm_radcoefs_simd4.c b/src/atrstm_radcoefs_simd4.c
@@ -63,17 +63,29 @@ fetch_radcoefs_simd4
 
   /* Linearly interpolate the node radiative properties */
   if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) {
+    const float prim_ka_min = v4f_x(v4f_reduce_min(prim_radcoefs.ka));
+    const float prim_ka_max = v4f_x(v4f_reduce_max(prim_radcoefs.ka));
     radcoefs[ATRSTM_RADCOEF_Ka] = v4f_x(v4f_dot(prim_radcoefs.ka, bcoords));
+    radcoefs[ATRSTM_RADCOEF_Ka] = /* Handle numerical uncertainty */
+      CLAMP(radcoefs[ATRSTM_RADCOEF_Ka], prim_ka_min, prim_ka_max);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] >= args->k_min[ATRSTM_RADCOEF_Ka]);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] <= args->k_max[ATRSTM_RADCOEF_Ka]);
   }
   if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) {
+    const float prim_ks_min = v4f_x(v4f_reduce_min(prim_radcoefs.ks));
+    const float prim_ks_max = v4f_x(v4f_reduce_max(prim_radcoefs.ks));
     radcoefs[ATRSTM_RADCOEF_Ks] = v4f_x(v4f_dot(prim_radcoefs.ks, bcoords));
+    radcoefs[ATRSTM_RADCOEF_Ks] = /* Handle numerical uncertainty */
+      CLAMP(radcoefs[ATRSTM_RADCOEF_Ks], prim_ks_min, prim_ks_max);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] >= args->k_min[ATRSTM_RADCOEF_Ks]);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] <= args->k_max[ATRSTM_RADCOEF_Ks]);
   }
   if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Kext) {
+    const float prim_kext_min = v4f_x(v4f_reduce_min(prim_radcoefs.kext));
+    const float prim_kext_max = v4f_x(v4f_reduce_max(prim_radcoefs.kext));
     radcoefs[ATRSTM_RADCOEF_Kext] = v4f_x(v4f_dot(prim_radcoefs.kext, bcoords));
+    radcoefs[ATRSTM_RADCOEF_Kext] = /* Handle numerical uncertainty */
+      CLAMP(radcoefs[ATRSTM_RADCOEF_Kext], prim_kext_min, prim_kext_max);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] >= args->k_min[ATRSTM_RADCOEF_Kext]);
     ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] <= args->k_max[ATRSTM_RADCOEF_Kext]);
   }