star-wf

commit e324a89de5800d68a10b84a5146c5dc059de7fc0
parent 7207a3746c7a45734ab515b474aa5bd9076af351
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Sun, 10 Mar 2024 11:51:30 +0100

Test tabulation accuracy

Check that with the default tabulation arguments, the relative error of
inversion is less than 1e-4 or 1e-6, whether using linear or quadratic
interpolation, respectively.

Diffstat:
M
Makefile
 | 
2
+-
M
src/test_swf_H3d.c
 | 
44
++++++++++++++++++++++++++++++++++++--------

2 files changed, 37 insertions(+), 9 deletions(-)
diff --git a/Makefile b/Makefile
@@ -138,4 +138,4 @@ $(TEST_OBJ): config.mk swf-local.pc
 	$(CC) $(CFLAGS_EXE) $(SWF_CFLAGS) $(RSYS_CFLAGS) -c $(@:.o=.c) -o $@
 
 test_swf_H3d: config.mk swf-local.pc $(LIBNAME)
-	$(CC) $(CFLAGS_EXE) -o $@ src/$@.o $(LDFLAGS_EXE) $(SWF_LIBS) $(RSYS_LIBS)
+	$(CC) $(CFLAGS_EXE) -o $@ src/$@.o $(LDFLAGS_EXE) $(SWF_LIBS) $(RSYS_LIBS) -lm
diff --git a/src/test_swf_H3d.c b/src/test_swf_H3d.c
@@ -19,10 +19,18 @@
 #include <rsys/math.h>
 #include <rsys/mem_allocator.h>
 
+#include <stdlib.h>
+
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
-static void
+static double
+rand_canonic(void)
+{
+  return (double)rand() / ((double)RAND_MAX+1);
+}
+
+static INLINE void
 check_tabulation_creation(void)
 {
   struct swf_H_tabulate_args args = SWF_H_TABULATE_ARGS_DEFAULT;
@@ -61,16 +69,35 @@ check_tabulation_creation(void)
 static void
 check_tabulation_inversion(void)
 {
+  const double delta_Hx = 1.0e-6;
+
+  struct swf_H_tabulate_args args = SWF_H_TABULATE_ARGS_DEFAULT;
   struct swf_tabulation* tab = NULL;
-  double x, y, ref, err;
+  double errl = 0;
+  double errq = 0;
+  size_t n = 0;
+  size_t i = 0;
+
+  CHK(swf_H3d_tabulate(&args, &tab) == RES_OK);
 
-  CHK(swf_H3d_tabulate(&SWF_H_TABULATE_ARGS_DEFAULT, &tab) == RES_OK);
+  n = (size_t)ceil(1.0/delta_Hx);
 
-  y = 1e-12;
-  x = swf_tabulation_inverse(tab, y);
-  ref = swf_H3d_eval(x);
-  err = fabs((ref - y) / ref);
-  CHK(err < 1.e-4);
+  FOR_EACH(i, 0, n) {
+    const double Hx = (double)i*delta_Hx + rand_canonic()*delta_Hx;
+    const double x_ref = swf_H3d_inverse(Hx);
+    const double xl = swf_tabulation_inverse(tab, SWF_LINEAR, Hx);
+    const double xq = swf_tabulation_inverse(tab, SWF_QUADRATIC, Hx);
+
+    if(!x_ref) {
+      CHK(x_ref == xl);
+    } else {
+      errl = fabs((x_ref - xl) / x_ref);
+      errq = fabs((x_ref - xq) / x_ref);
+      /*printf("%e %e %e %e %e\n", x_ref, xq, xl, errq, errl);*/
+      CHK(errl < 1.e-4);
+      CHK(errq < 1.e-6);
+    }
+  }
 
   CHK(swf_tabulation_ref_put(tab) == RES_OK);
 }
@@ -82,6 +109,7 @@ int
 main(int argc, char** argv)
 {
   (void)argc, (void)argv;
+
   check_tabulation_creation();
   check_tabulation_inversion();
   CHK(mem_allocated_size() == 0);