Add convergence test for Simpson integration method

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Multipole/trunk@94 4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843

6 files changed, 120 insertions, 2 deletions

diff --git a/interface.ccl b/interface.ccl
index a8eda59..6c12ff7 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -13,3 +13,8 @@ CCTK_REAL harmonics type=GF timelevels=1
 {
   harmonic_re, harmonic_im
 } "Spherical harmonics"
+
+CCTK_REAL test_integration_convergence_orders type=SCALAR
+{
+  test_simpson_convergence_order
+} "Test integration convergence orders"
diff --git a/schedule.ccl b/schedule.ccl
index 9c94f21..d460e0f 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,5 +1,7 @@
 #schedule.ccl for thorn Multipole
 
+STORAGE: test_integration_convergence_orders
+
 if (enable_test)
 {
   STORAGE: harmonics[1]
@@ -24,3 +26,10 @@ schedule Multipole_ParamCheck at CCTK_PARAMCHECK
   LANG: C
   OPTIONS: GLOBAL
 } "Check Multipole parameters"
+
+# # Tests
+
+schedule Multipole_TestSimpson at CCTK_PARAMCHECK
+{
+  LANG: C
+} "Test simpson integration method"
diff --git a/src/make.code.defn b/src/make.code.defn
index cccefbd..4d49656 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -1 +1 @@
-SRCS=interpolate.cc multipole.cc utils.cc sphericalharmonic.cc integrate.cc
+SRCS=interpolate.cc multipole.cc utils.cc sphericalharmonic.cc integrate.cc tests.cc
diff --git a/src/tests.cc b/src/tests.cc
new file mode 100644
index 0000000..b63e18b
--- /dev/null
+++ b/src/tests.cc
@@ -0,0 +1,97 @@
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "utils.hh"
+#include "integrate.hh"
+
+static CCTK_REAL test_simpson_integral(int n)
+{
+  const int nx = n;
+  const int ny = n;
+  const int array_size=(nx+1)*(ny+1);
+
+  CCTK_REAL *f = new CCTK_REAL[array_size];
+
+  const CCTK_REAL dx = 1./nx;
+  const CCTK_REAL dy = 1./ny;
+
+  for (int ix = 0; ix <= nx; ix++)
+  {
+    for (int iy = 0; iy <= ny; iy++)
+    {
+      const int i = Multipole_Index(ix, iy, nx);
+
+      const CCTK_REAL x = ix*dx;
+      const CCTK_REAL y = iy*dy;
+      const CCTK_REAL PI = acos(-1.0);
+
+      f[i] = x*pow(y,2)*pow(cos(2*PI*y),2)*pow(sin(2*PI*x),2);
+    }
+  }
+
+  const CCTK_REAL result = Simpson2DIntegral(f, nx, ny, dx, dy);
+  delete [] f;
+  return result;
+}
+
+
+void Multipole_TestSimpson(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+
+  const int n1 = 100;
+  const int n2 = 200;
+  const CCTK_REAL PI = acos(-1.0);
+  const CCTK_REAL result1 = test_simpson_integral(100);
+  const CCTK_REAL result2 = test_simpson_integral(200);
+  const CCTK_REAL exact = 1./24 + 1./(64 * pow(PI,2));
+  const CCTK_REAL error1 = fabs(result1 - exact);
+  const CCTK_REAL error2 = fabs(result2 - exact);
+  *test_simpson_convergence_order = log10(error1/error2) / log10((CCTK_REAL) n2/n1);
+}
+
+// void Multipole_TestIntegrate(CCTK_ARGUMENTS)
+// {
+//   const int n = 100;
+
+//   const int nth = n;
+//   const int nph = n;
+//   const int array_size=(nth+1)*(nph+1);
+
+//   CCTK_REAL *f = new CCTK_REAL[array_size];
+
+//   const CCTK_REAL dth = 1/nx;
+//   const CCTK_REAL dph = 1/ny;
+
+//   for (int ix = 0; ix <= nx; ix++)
+//   {
+//     for (int iy = 0; iy <= ny; iy++)
+//     {
+//       const int i = Multipole_Index(ix, iy, nx);
+
+//       const CCTK_REAL x = ix*dx;
+//       const CCTK_REAL y = iy*dy;
+
+//       f[i] = sin(2*PI*x)*cos(2*PI*y);
+//     }
+//   }
+
+
+
+
+//   CCTK_REAL result = Multipole_Integrate(int array_size, int nthetap,
+//                                          CCTK_REAL const array1r[], CCTK_REAL const array1i[],
+//                                          CCTK_REAL const array2r[], CCTK_REAL const array2i[],
+//                                          CCTK_REAL const th[], CCTK_REAL const ph[], 
+//                                          CCTK_REAL *outre, CCTK_REAL *outim)
+
+
+
+//   printf("Integration result: %.19g\n", result);
+// }
diff --git a/test/test_simpson.par b/test/test_simpson.par
index 78daaf4..772d2a9 100644
--- a/test/test_simpson.par
+++ b/test/test_simpson.par
@@ -1,5 +1,5 @@
 
-ActiveThorns = "CoordBase SymBase Boundary CartGrid3d IOUtil Carpet CarpetLib  CarpetInterp AEILocalInterp InitBase Multipole LoopControl"
+ActiveThorns = "CoordBase SymBase Boundary CartGrid3d IOUtil Carpet CarpetLib  CarpetInterp AEILocalInterp InitBase Multipole LoopControl CarpetIOASCII"
 
 #############################################################
 # Grid
@@ -87,3 +87,6 @@ IO::out_fileinfo                     = "none"
 # CarpetIOASCII::out1d_x = yes
 # CarpetIOASCII::out1d_y = yes
 # CarpetIOASCII::out1d_z = yes
+
+CarpetIOASCII::out0d_vars = "Multipole::test_simpson_convergence_order"
+CarpetIOASCII::out0d_every = 1
diff --git a/test/test_simpson/test_simpson_convergence_order..asc b/test/test_simpson/test_simpson_convergence_order..asc
new file mode 100644
index 0000000..446999a
--- /dev/null
+++ b/test/test_simpson/test_simpson_convergence_order..asc
@@ -0,0 +1,4 @@
+# 0D ASCII output created by CarpetIOASCII
+#
+0	0	0 0 0	0 0 0	0	0 0 0	4.00339319061966
+