Add tests for symmetry of integration

It is desirable that integrands with exact pi-antisymmetry integrate to zero; we add measures of this here.  This is the case for all but the midpoint method.

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

3 files changed, 60 insertions, 0 deletions

diff --git a/interface.ccl b/interface.ccl
index 2e893ab..26d7ef7 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -20,3 +20,11 @@ CCTK_REAL test_integration_convergence_orders type=SCALAR
   test_trapezoidal_convergence_order,
   test_simpson_convergence_order
 } "Test integration convergence orders"
+
+CCTK_REAL test_integration_symmetries type=SCALAR
+{
+  test_midpoint_pi_symmetry,
+  test_trapezoidal_pi_symmetry,
+  test_simpson_pi_symmetry,
+  test_driscollhealy_pi_symmetry
+} "Test integration symmetries"
diff --git a/schedule.ccl b/schedule.ccl
index 328e1de..3b46cae 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -4,6 +4,7 @@ if (enable_test)
 {
   STORAGE: harmonics[1]
   STORAGE: test_integration_convergence_orders
+  STORAGE: test_integration_symmetries
 }
 
 schedule Multipole_Calc at CCTK_ANALYSIS after (calc_np,PsiKadelia,Accelerator_CopyBack)
@@ -34,4 +35,9 @@ if (enable_test)
   {
     LANG: C
   } "Test convergence of integration"
+
+  schedule Multipole_TestIntegrationSymmetry at CCTK_PARAMCHECK
+  {
+    LANG: C
+  } "Test symmetry of integration"
 }
diff --git a/src/tests.cc b/src/tests.cc
index 42e79b3..de711eb 100644
--- a/src/tests.cc
+++ b/src/tests.cc
@@ -40,6 +40,38 @@ static CCTK_REAL test_integral(int n, CCTK_REAL (*integration_fn) (const CCTK_RE
   return result;
 }
 
+
+static CCTK_REAL test_pi_symmetric_sphere_integral(CCTK_REAL (*integration_fn) (const CCTK_REAL *, int, int, CCTK_REAL, CCTK_REAL))
+{
+  const int n = 100;
+  const int nth = n;
+  const int nph = n;
+  const int array_size=(nth+1)*(nph+1);
+  const CCTK_REAL PI = acos(-1.0);
+
+  CCTK_REAL *f = new CCTK_REAL[array_size];
+
+  const CCTK_REAL dth = PI/nth;
+  const CCTK_REAL dph = 2*PI/nph;
+
+  for (int ith = 0; ith <= nth; ith++)
+  {
+    for (int iph = 0; iph <= nph; iph++)
+    {
+      const int i = Multipole_Index(ith, iph, nth);
+
+      const CCTK_REAL th = ith*dth;
+      const CCTK_REAL ph = iph*dph;
+
+      f[i] = -(cos(ph)*sqrt(5/PI)*pow(cos(th/2.),3)*sin(th/2.));
+    }
+  }
+
+  const CCTK_REAL result = integration_fn(f, nth, nph, dth, dph);
+  delete [] f;
+  return result;
+}
+
 CCTK_REAL integration_convergence_order(CCTK_REAL (*integration_fn) (const CCTK_REAL *, int, int, CCTK_REAL, CCTK_REAL))
 {
   const int n1 = 100;
@@ -61,6 +93,20 @@ void Multipole_TestIntegrationConvergence(CCTK_ARGUMENTS)
   *test_midpoint_convergence_order = integration_convergence_order(&Midpoint2DIntegral);
 }
 
+void Multipole_TestIntegrationSymmetry(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+
+  *test_simpson_pi_symmetry = test_pi_symmetric_sphere_integral(&Simpson2DIntegral);
+  *test_midpoint_pi_symmetry = test_pi_symmetric_sphere_integral(&Midpoint2DIntegral);
+  *test_trapezoidal_pi_symmetry = test_pi_symmetric_sphere_integral(&Trapezoidal2DIntegral);
+  *test_driscollhealy_pi_symmetry = test_pi_symmetric_sphere_integral(&DriscollHealy2DIntegral);
+  printf("Pi symmetry Simpson integral: %.19g\n", *test_simpson_pi_symmetry);
+  printf("Pi symmetry midpoint integral: %.19g\n", *test_midpoint_pi_symmetry);
+  printf("Pi symmetry trapezoidal integral: %.19g\n", *test_trapezoidal_pi_symmetry);
+  printf("Pi symmetry Driscoll and Healy integral: %.19g\n", *test_driscollhealy_pi_symmetry);
+}
+
 // void Multipole_TestIntegrate(CCTK_ARGUMENTS)
 // {
 //   const int n = 100;