star-ck

commit 15190372e57c4c9f09693939094ec7560d101814
parent b8e1c8fcca2db7d211ff03acf46d5959d016c2c3
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed,  4 Nov 2020 16:19:31 +0100

Test the atrck_band_sample_quad_pt function

Diffstat:
M
src/atrck.c
 | 
1
+
M
src/atrck.h
 | 
1
+
M
src/test_atrck_load.c
 | 
67
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---

3 files changed, 66 insertions(+), 3 deletions(-)
diff --git a/src/atrck.c b/src/atrck.c
@@ -483,6 +483,7 @@ atrck_band_get_quad_pt
   atrck_quad_pt->ka_list = quad_pt->ka_list;
   atrck_quad_pt->abscissa = quad_pt->abscissa;
   atrck_quad_pt->weight = quad_pt->weight;
+  atrck_quad_pt->id = iquad_pt;
 
 exit:
   return res;
diff --git a/src/atrck.h b/src/atrck.h
@@ -52,6 +52,7 @@ struct atrck_quad_pt {
   float* ka_list; /* Per node ka */
   double abscissa; /* m^-1 */
   double weight;
+  size_t id;
 };
 static const struct atrck_quad_pt ATRCK_QUAD_PT_NULL;
 
diff --git a/src/test_atrck_load.c b/src/test_atrck_load.c
@@ -15,7 +15,17 @@
 
 #include "atrck.h"
 
+#include <rsys/math.h>
 #include <rsys/mem_allocator.h>
+#include <math.h>
+
+static INLINE float
+rand_canonic(void)
+{
+  int r;
+  while((r = rand()) == RAND_MAX);
+  return (float)r / (float)RAND_MAX;
+}
 
 /*******************************************************************************
  * Helper functions
@@ -26,9 +36,12 @@ check_atrck_load
    const size_t nbands,
    const size_t nnodes)
 {
+  #define NMOMENTS 4
   struct atrck_band band = ATRCK_BAND_NULL;
   struct atrck_quad_pt qpt = ATRCK_QUAD_PT_NULL;
   size_t iband;
+  const size_t nsamps = 1000000;
+  size_t isamp;
 
   CHK(atrck);
   CHK(nbands);
@@ -43,15 +56,26 @@ check_atrck_load
   CHK(atrck_get_band(atrck, 0, &band) == RES_OK);
 
   CHK(band.quad_pts_count == 1);
-  CHK(atrck_band_get_quad_pt( NULL, 0, &qpt) == RES_BAD_ARG);
-  CHK(atrck_band_get_quad_pt( &band, 1, &qpt) == RES_BAD_ARG);
-  CHK(atrck_band_get_quad_pt( &band, 0, NULL) == RES_BAD_ARG);
+  CHK(atrck_band_get_quad_pt(NULL, 0, &qpt) == RES_BAD_ARG);
+  CHK(atrck_band_get_quad_pt(&band, 1, &qpt) == RES_BAD_ARG);
+  CHK(atrck_band_get_quad_pt(&band, 0, NULL) == RES_BAD_ARG);
+
+  CHK(atrck_band_sample_quad_pt(NULL, 0, &qpt) == RES_BAD_ARG);
+  CHK(atrck_band_sample_quad_pt(&band, 1, &qpt) == RES_BAD_ARG);
+  CHK(atrck_band_sample_quad_pt(&band, 0, NULL) == RES_BAD_ARG);
 
   FOR_EACH(iband, 0, nbands) {
     const double low = (double)iband;
     const double upp = (double)(iband+1);
     const size_t nqpts = iband + 1;
     size_t iqpt;
+    double sum[NMOMENTS] = {0};
+    double sum2[NMOMENTS] = {0};
+    double E[NMOMENTS] = {0};
+    double V[NMOMENTS] = {0};
+    double SE[NMOMENTS] = {0};
+    double ref;
+    int imoment;
 
     CHK(atrck_get_band(atrck, iband, &band) == RES_OK);
     CHK(band.lower == low);
@@ -67,12 +91,49 @@ check_atrck_load
       CHK(atrck_band_get_quad_pt(&band, iqpt, &qpt) == RES_OK);
       CHK(qpt.abscissa == abscissa);
       CHK(qpt.weight == weight);
+      CHK(qpt.id == iqpt);
 
       FOR_EACH(inode, 0, nnodes) {
         const float ka = (float)(iband*1000 + iqpt*100 + inode);
         CHK(qpt.ka_list[inode] == ka);
       }
     }
+
+    /* Check the sampling of the quadrature points */
+    FOR_EACH(isamp, 0, nsamps) {
+      const double r = rand_canonic();
+      double w[NMOMENTS];
+
+      CHK(atrck_band_sample_quad_pt(&band, r, &qpt) == RES_OK);
+      FOR_EACH(imoment, 0, NMOMENTS) {
+        if(imoment == 0) {
+          w[imoment] = (double)(qpt.id + 1);
+        } else {
+          w[imoment] = (double)(qpt.id + 1) * w[imoment-1];;
+        }
+        sum[imoment] += w[imoment];
+        sum2[imoment] += w[imoment]*w[imoment];
+      }
+    }
+
+    FOR_EACH(imoment, 0, NMOMENTS) {
+      E[imoment] = sum[imoment] / (double)nsamps;
+      V[imoment] = sum2[imoment] / (double)nsamps - E[imoment]*E[imoment];
+      SE[imoment] = sqrt(V[imoment]/(double)nsamps);
+    }
+
+    ref = (double)(nqpts+1)/2.0;
+    CHK(eq_eps(ref, E[0], 3*SE[0]));
+
+    ref = (double)(2*nqpts*nqpts + 3*nqpts + 1) / 6.0;
+    CHK(eq_eps(ref, E[1], 3*SE[1]));
+
+    ref = ((double)nqpts * pow((double)nqpts + 1.0, 2.0)) / 4.0;
+    CHK(eq_eps(ref, E[2], 3*SE[2]));
+
+    ref = (double)1.0/30.0 * (double)
+      (6*nqpts*nqpts*nqpts*nqpts + 15*nqpts*nqpts*nqpts + 10*nqpts*nqpts - 1);
+    CHK(eq_eps(ref, E[3], 3*SE[3]));
   }
 }