[Commits] [svn:einsteintoolkit] Multipole/trunk/ (Rev. 102)

[ian.hinder at aei.mpg.de]

User: hinder
Date: 2013/10/09 12:57 PM

Modified:
 /trunk/
  interface.ccl, schedule.ccl
 /trunk/src/
  tests.cc

Log:
 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.

File Changes:

Directory: /trunk/src/
======================

File [modified]: tests.cc
Delta lines: +46 -0
===================================================================
--- trunk/src/tests.cc	2013-10-09 17:57:06 UTC (rev 101)
+++ trunk/src/tests.cc	2013-10-09 17:57:15 UTC (rev 102)
@@ -40,6 +40,38 @@
   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 @@
   *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;

Directory: /trunk/
==================

File [modified]: interface.ccl
Delta lines: +8 -0
===================================================================
--- trunk/interface.ccl	2013-10-09 17:57:06 UTC (rev 101)
+++ trunk/interface.ccl	2013-10-09 17:57:15 UTC (rev 102)
@@ -20,3 +20,11 @@
   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"

File [modified]: schedule.ccl
Delta lines: +6 -0
===================================================================
--- trunk/schedule.ccl	2013-10-09 17:57:06 UTC (rev 101)
+++ trunk/schedule.ccl	2013-10-09 17:57:15 UTC (rev 102)
@@ -4,6 +4,7 @@
 {
   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 @@
   {
     LANG: C
   } "Test convergence of integration"
+
+  schedule Multipole_TestIntegrationSymmetry at CCTK_PARAMCHECK
+  {
+    LANG: C
+  } "Test symmetry of integration"
 }