Implement Driscoll&Healy integration

Implements a more accurate integration over the sphere using an
algorithm by Driscoll & Healy. This algorithm uses Gaussian
integration weights, leading (almost) to exponential convergence.


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

19 files changed, 588 insertions, 6 deletions

diff --git a/param.ccl b/param.ccl
index 9af3650..b8ed861 100644
--- a/param.ccl
+++ b/param.ccl
@@ -29,8 +29,9 @@ CCTK_STRING coord_system "What is the coord system?"
 
 KEYWORD integration_method "How to do surface integrals" STEERABLE=always
 {
-  "midpoint" :: "Midpoint rule (1st order)"
-  "Simpson"  :: "Simpson's rule (4th order)"
+  "midpoint"      :: "Midpoint rule (1st order)"
+  "Simpson"       :: "Simpson's rule (4th order) [requires even ntheta and nphi]"
+  "DriscollHealy" :: "Driscoll & Healy (exponentially convergent) [requires even ntheta]"
 } "midpoint"
 
 CCTK_INT out_every "How often to output" \
diff --git a/src/integrate.cc b/src/integrate.cc
index 0f1c392..1a9f44b 100644
--- a/src/integrate.cc
+++ b/src/integrate.cc
@@ -107,6 +107,52 @@ CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, CCTK_REAL hx, CC
   return ((double) 1 / (double) 9) * hx * hy * integrand_sum;
 }
 
+// See: J.R. Driscoll and D.M. Healy Jr., Computing Fourier transforms
+// and convolutions on the 2-sphere. Advances in Applied Mathematics,
+// 15(2):202–250, 1994.
+CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *const f,
+                                  int const nx, int const ny,
+                                  CCTK_REAL const hx, CCTK_REAL const hy)
+{
+  assert(f);
+  assert(nx >= 0);
+  assert(ny >= 0);
+  assert(nx % 2 == 0);
+  
+  CCTK_REAL integrand_sum = 0.0;
+  
+  // Skip the poles (ix=0 and ix=nx), since the weight there is zero
+  // anyway
+#pragma omp parallel for reduction(+: integrand_sum)
+  for (int ix = 1; ix < nx; ++ ix)
+  {
+    
+    // These weights lead to an almost spectral convergence
+    CCTK_REAL const theta = M_PI * ix / nx;
+    CCTK_REAL weight = 0.0;
+    for (int l = 0; l < nx/2; ++ l)
+    {
+      weight += sin((2*l+1)*theta)/(2*l+1);
+    }
+    weight *= 4.0 / M_PI;
+    
+    CCTK_REAL local_sum = 0.0;
+    // Skip the last point (iy=ny), since we assume periodicity and
+    // therefor it has the same value as the first point. We don't use
+    // weights in this direction, which leads to spectral convergence.
+    // (Yay periodicity!)
+    for (int iy = 0; iy < ny; ++ iy)
+    {
+      local_sum += f[idx(ix,iy)];
+    }
+    
+    integrand_sum += weight * local_sum;
+    
+  }
+  
+  return hx * hy * integrand_sum;
+}
+
 // 1D integrals
 
 static CCTK_REAL Simpson1DIntegral(CCTK_REAL const *f, int n, CCTK_REAL h)
@@ -114,7 +160,8 @@ static CCTK_REAL Simpson1DIntegral(CCTK_REAL const *f, int n, CCTK_REAL h)
   CCTK_REAL integrand_sum = 0;
   int i = 0;
 
-  assert(n > 0); assert(f);
+  assert(f);
+  assert(n > 0);
   assert(n % 2 == 0);
 
   int p = n / 2;
diff --git a/src/integrate.hh b/src/integrate.hh
index c45917e..5bcd43f 100644
--- a/src/integrate.hh
+++ b/src/integrate.hh
@@ -6,4 +6,7 @@
 CCTK_REAL Simpson2DIntegral(CCTK_REAL const *f, int nx, int ny, 
                             CCTK_REAL hx, CCTK_REAL hy);
 
+CCTK_REAL DriscollHealy2DIntegral(CCTK_REAL const *f, int nx, int ny,
+                                  CCTK_REAL hx, CCTK_REAL hy);
+
 #endif
diff --git a/src/utils.cc b/src/utils.cc
index 74d9eef..cb25bdf 100644
--- a/src/utils.cc
+++ b/src/utils.cc
@@ -346,7 +346,8 @@ void Multipole_Integrate(int array_size, int nthetap,
       *outim = tempi;
     } 
   }
-  else if (CCTK_Equals(integration_method, "simpson"))
+  else if (CCTK_Equals(integration_method, "Simpson") ||
+           CCTK_Equals(integration_method, "DriscollHealy"))
   {
     static CCTK_REAL *fr = 0;
     static CCTK_REAL *fi = 0;
@@ -368,7 +369,31 @@ void Multipole_Integrate(int array_size, int nthetap,
                array1i[i] * array2r[i] ) * sin(th[i]);
     }
     
-    *outre = Simpson2DIntegral(fr, ntheta, nphi, dth, dph);
-    *outim = Simpson2DIntegral(fi, ntheta, nphi, dth, dph);
+    if (CCTK_Equals(integration_method, "Simpson"))
+    {
+      if (nphi % 2 != 0 || ntheta % 2 != 0)
+      {
+        CCTK_WARN (CCTK_WARN_ABORT, "The Simpson integration method requires even ntheta and even nphi");
+      }
+      *outre = Simpson2DIntegral(fr, ntheta, nphi, dth, dph);
+      *outim = Simpson2DIntegral(fi, ntheta, nphi, dth, dph);
+    }
+    else if (CCTK_Equals(integration_method, "DriscollHealy"))
+    {
+      if (ntheta % 2 != 0)
+      {
+        CCTK_WARN (CCTK_WARN_ABORT, "The Driscoll&Healy integration method requires even ntheta");
+      }
+      *outre = DriscollHealy2DIntegral(fr, ntheta, nphi, dth, dph);
+      *outim = DriscollHealy2DIntegral(fi, ntheta, nphi, dth, dph);
+    }
+    else
+    {
+      CCTK_WARN(CCTK_WARN_ABORT, "internal error");
+    }
+  }
+  else
+  {
+    CCTK_WARN(CCTK_WARN_ABORT, "internal error");
   }
 } 
diff --git a/test/test_driscollhealy.par b/test/test_driscollhealy.par
new file mode 100644
index 0000000..c1fee6f
--- /dev/null
+++ b/test/test_driscollhealy.par
@@ -0,0 +1,89 @@
+
+ActiveThorns = "CoordBase SymBase Boundary CartGrid3d IOUtil Carpet CarpetLib  CarpetInterp AEILocalInterp InitBase Multipole LoopControl"
+
+#############################################################
+# Grid
+#############################################################
+
+CartGrid3D::type                        = "coordbase"
+CartGrid3D::domain                      = "full"
+CartGrid3D::avoid_origin                = "no"
+
+CoordBase::domainsize                   = minmax
+CoordBase::xmin                         = -10
+CoordBase::ymin                         = -10
+CoordBase::zmin                         = -10
+CoordBase::xmax                         = 10
+CoordBase::ymax                         = 10
+CoordBase::zmax                         = 10
+CoordBase::dx                           = 0.2
+CoordBase::dy                           = 0.2
+CoordBase::dz                           = 0.2
+CoordBase::boundary_size_x_lower        = 2
+CoordBase::boundary_size_y_lower        = 2
+CoordBase::boundary_size_z_lower        = 2
+CoordBase::boundary_shiftout_x_lower    = 0
+CoordBase::boundary_shiftout_y_lower    = 0
+CoordBase::boundary_shiftout_z_lower    = 0
+CoordBase::boundary_size_x_upper        = 2
+CoordBase::boundary_size_y_upper        = 2
+CoordBase::boundary_size_z_upper        = 2
+CoordBase::boundary_shiftout_x_upper    = 0
+CoordBase::boundary_shiftout_y_upper    = 0
+CoordBase::boundary_shiftout_z_upper    = 0
+
+#############################################################
+# Carpet
+#############################################################
+
+Carpet::ghost_size                      = 3
+Carpet::domain_from_coordbase           = "yes"
+Carpet::poison_new_timelevels           = "yes"
+Carpet::check_for_poison                = "no"
+CarpetLib::poison_value                    = 113
+Carpet::init_fill_timelevels            = "yes"
+
+#############################################################
+# CarpetLib
+#############################################################
+
+CarpetLib::poison_new_memory            = "yes"
+
+#############################################################
+# Cactus
+#############################################################
+
+Cactus::terminate                     = "iteration"
+Cactus::cctk_itlast                   = 0
+
+#############################################################
+# Multipole
+#############################################################
+
+Multipole::nradii = 1
+Multipole::radius[0] = 8.0
+Multipole::variables = "Multipole::harmonic_re{sw=-2 cmplx='Multipole::harmonic_im' name='harmonic'}"
+Multipole::integration_method = "DriscollHealy"
+
+Multipole::enable_test = "yes"
+Multipole::test_l = 2
+Multipole::test_m = 2
+Multipole::out_1d_every = 1
+
+#############################################################
+# Output
+#############################################################
+
+IO::out_dir                          = $parfile
+IO::out_fileinfo                     = "none"
+
+# Enabling 1D output for the test grid functions would be helpful for
+# localising any failures but it makes the tests dependent on the
+# number of processors, as CarpetIOASCII's output is dependent on
+# this.
+
+# CarpetIOASCII::out1d_vars = "Multipole::harmonics"
+# CarpetIOASCII::out1d_every = 1
+# CarpetIOASCII::out1d_x = yes
+# CarpetIOASCII::out1d_y = yes
+# CarpetIOASCII::out1d_z = yes
diff --git a/test/test_driscollhealy/mp_harmonic_im_r8.00.ph.asc b/test/test_driscollhealy/mp_harmonic_im_r8.00.ph.asc
new file mode 100644
index 0000000..b689c89
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_im_r8.00.ph.asc
@@ -0,0 +1,104 @@
+"Time = 0
+0.000000 -6.907573430162997911e-17
+0.062832 0.0597709107978601703
+0.125664 0.1185991090836183148
+0.188496 0.1755571071177308706
+0.251327 0.22974631087048486
+0.314159 0.2803123461157627783
+0.376991 0.326457761226630172
+0.439823 0.3674546220110177508
+0.502655 0.4026566798676182168
+0.565487 0.4315084643619017113
+0.628319 0.4535551023743125687
+0.691150 0.4684489721639752968
+0.753982 0.4759551277630860722
+0.816814 0.4759551277630863497
+0.879646 0.4684489721639752968
+0.942478 0.4535551023743126797
+1.005310 0.4315084643619018778
+1.068142 0.4026566798676179393
+1.130973 0.3674546220110178063
+1.193805 0.3264577612266306161
+1.256637 0.2803123461157628893
+1.319469 0.2297463108704850265
+1.382301 0.1755571071177310649
+1.445133 0.1185991090836186479
+1.507964 0.05977091079786012173
+1.570796 -1.401964684891840811e-16
+1.633628 -0.05977091079786012867
+1.696460 -0.1185991090836181483
+1.759292 -0.1755571071177305376
+1.822124 -0.229746310870484749
+1.884956 -0.2803123461157626672
+1.947787 -0.326457761226630172
+2.010619 -0.3674546220110179173
+2.073451 -0.4026566798676175507
+2.136283 -0.4315084643619017668
+2.199115 -0.4535551023743126242
+2.261947 -0.4684489721639758519
+2.324779 -0.4759551277630864607
+2.387610 -0.4759551277630861277
+2.450442 -0.4684489721639749082
+2.513274 -0.4535551023743128463
+2.576106 -0.4315084643619019333
+2.638938 -0.4026566798676183834
+2.701770 -0.3674546220110177508
+2.764602 -0.3264577612266314488
+2.827433 -0.2803123461157625007
+2.890265 -0.2297463108704854706
+2.953097 -0.1755571071177309539
+3.015929 -0.1185991090836186895
+3.078761 -0.05977091079786037847
+3.141593 1.425743635182407722e-16
+3.204425 0.05977091079786019112
+3.267256 0.1185991090836186063
+3.330088 0.1755571071177305098
+3.392920 0.2297463108704850543
+3.455752 0.2803123461157625562
+3.518584 0.3264577612266305606
+3.581416 0.3674546220110170847
+3.644247 0.4026566798676178283
+3.707079 0.4315084643619019888
+3.769911 0.4535551023743124577
+3.832743 0.4684489721639751858
+3.895575 0.4759551277630860167
+3.958407 0.4759551277630860722
+4.021239 0.4684489721639746862
+4.084070 0.4535551023743122911
+4.146902 0.4315084643619018778
+4.209734 0.4026566798676178838
+4.272566 0.3674546220110170291
+4.335398 0.3264577612266304496
+4.398230 0.2803123461157624452
+4.461062 0.2297463108704845824
+4.523893 0.1755571071177310372
+4.586725 0.1185991090836189532
+4.649557 0.05977091079785964989
+4.712389 2.43532537710590483e-16
+4.775221 -0.05977091079785951111
+4.838053 -0.1185991090836182593
+4.900885 -0.1755571071177308429
+4.963716 -0.229746310870485082
+5.026548 -0.2803123461157622787
+5.089380 -0.3264577612266307272
+5.152212 -0.3674546220110169736
+5.215044 -0.4026566798676180503
+5.277876 -0.4315084643619020444
+5.340708 -0.4535551023743122911
+5.403539 -0.4684489721639749082
+5.466371 -0.4759551277630857391
+5.529203 -0.4759551277630865163
+5.592035 -0.4684489721639766291
+5.654867 -0.4535551023743119581
+5.717699 -0.4315084643619015448
+5.780530 -0.4026566798676182724
+5.843362 -0.3674546220110180839
+5.906194 -0.3264577612266302276
+5.969026 -0.2803123461157628338
+6.031858 -0.2297463108704853041
+6.094690 -0.1755571071177306763
+6.157522 -0.1185991090836190087
+6.220353 -0.059770910797860205
+6.283185 -4.802722211955578559e-16
+
+
diff --git a/test/test_driscollhealy/mp_harmonic_im_r8.00.th.asc b/test/test_driscollhealy/mp_harmonic_im_r8.00.th.asc
new file mode 100644
index 0000000..41a6a88
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_im_r8.00.th.asc
@@ -0,0 +1,55 @@
+"Time = 0
+0.000000 0.6307831305050407567
+0.061600 0.007024422862686013555
+0.123200 -8.818650445389855207e-16
+0.184800 2.386973079634007145e-16
+0.246399 3.090954546896031658e-16
+0.307999 1.88689857105083995e-16
+0.369599 1.787601751083378799e-16
+0.431199 6.42809675982812758e-17
+0.492799 -5.223354783714025765e-17
+0.554399 1.052874518799022369e-16
+0.615999 1.117372074496677921e-16
+0.677598 8.829044623742264276e-17
+0.739198 -6.907573430162997911e-17
+0.800798 -3.90179774949959511e-17
+0.862398 6.159252871967265588e-17
+0.923998 5.627127698640937574e-17
+0.985598 4.652774069628360505e-17
+1.047198 7.285297644863500048e-17
+1.108797 -6.542430094850337162e-18
+1.170397 6.465233305208766899e-17
+1.231997 1.779484214503371193e-17
+1.293597 -7.602576948937617888e-18
+1.355197 1.390238123941617177e-17
+1.416797 2.070976566529800156e-17
+1.478397 1.584136196322225341e-17
+1.539996 2.282432166526202462e-17
+1.601596 -3.480861916600505987e-17
+1.663196 4.246513194032316906e-18
+1.724796 8.325741691803758443e-18
+1.786396 8.248138457401257405e-18
+1.847996 1.165433598435465015e-17
+1.909596 6.66306970460519502e-18
+1.971195 5.865616208071797884e-18
+2.032795 4.518590889726565158e-18
+2.094395 5.639724593917154933e-19
+2.155995 8.462818399384888541e-19
+2.217595 7.220434734079161545e-18
+2.279195 -1.940108903329609695e-18
+2.340795 1.091953982035055254e-18
+2.402394 9.65339816474644947e-19
+2.463994 5.865854807296045742e-19
+2.525594 1.106952583858313206e-18
+2.587194 1.616933492592317859e-18
+2.648794 -8.596749700189567001e-19
+2.710394 3.992357306555509448e-19
+2.771994 1.582759388194540812e-19
+2.833593 2.217055110251093269e-20
+2.895193 2.827594302876126491e-20
+2.956793 2.496238261006691607e-20
+3.018393 7.633200804732358595e-21
+3.079993 1.074058560061705883e-21
+3.141593 -7.22781544709519842e-37
+
+
diff --git a/test/test_driscollhealy/mp_harmonic_l0_m0_r8.00.asc b/test/test_driscollhealy/mp_harmonic_l0_m0_r8.00.asc
new file mode 100644
index 0000000..23eaf79
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_l0_m0_r8.00.asc
@@ -0,0 +1 @@
+0.000000 0 0
diff --git a/test/test_driscollhealy/mp_harmonic_l1_m-1_r8.00.asc b/test/test_driscollhealy/mp_harmonic_l1_m-1_r8.00.asc
new file mode 100644
index 0000000..23eaf79
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_l1_m-1_r8.00.asc
@@ -0,0 +1 @@
+0.000000 0 0
diff --git a/test/test_driscollhealy/mp_harmonic_l1_m0_r8.00.asc b/test/test_driscollhealy/mp_harmonic_l1_m0_r8.00.asc
new file mode 100644
index 0000000..23eaf79
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_l1_m0_r8.00.asc
@@ -0,0 +1 @@
+0.000000 0 0
diff --git a/test/test_driscollhealy/mp_harmonic_l1_m1_r8.00.asc b/test/test_driscollhealy/mp_harmonic_l1_m1_r8.00.asc
new file mode 100644
index 0000000..23eaf79
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_l1_m1_r8.00.asc
@@ -0,0 +1 @@
+0.000000 0 0
diff --git a/test/test_driscollhealy/mp_harmonic_l2_m-1_r8.00.asc b/test/test_driscollhealy/mp_harmonic_l2_m-1_r8.00.asc
new file mode 100644
index 0000000..eab32f1
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_l2_m-1_r8.00.asc
@@ -0,0 +1 @@
+0.000000 -8.332830142340763872e-17 1.022973937129801201e-16
diff --git a/test/test_driscollhealy/mp_harmonic_l2_m-2_r8.00.asc b/test/test_driscollhealy/mp_harmonic_l2_m-2_r8.00.asc
new file mode 100644
index 0000000..2fa5550
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_l2_m-2_r8.00.asc
@@ -0,0 +1 @@
+0.000000 2.138002145330596485e-10 2.446208802313476189e-17
diff --git a/test/test_driscollhealy/mp_harmonic_l2_m0_r8.00.asc b/test/test_driscollhealy/mp_harmonic_l2_m0_r8.00.asc
new file mode 100644
index 0000000..08d5282
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_l2_m0_r8.00.asc
@@ -0,0 +1 @@
+0.000000 -1.931250071812291289e-16 2.833087430018189041e-07
diff --git a/test/test_driscollhealy/mp_harmonic_l2_m1_r8.00.asc b/test/test_driscollhealy/mp_harmonic_l2_m1_r8.00.asc
new file mode 100644
index 0000000..5b6f793
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_l2_m1_r8.00.asc
@@ -0,0 +1 @@
+0.000000 -8.269870904790199682e-17 -5.485377474579588873e-17
diff --git a/test/test_driscollhealy/mp_harmonic_l2_m2_r8.00.asc b/test/test_driscollhealy/mp_harmonic_l2_m2_r8.00.asc
new file mode 100644
index 0000000..cfe52d1
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_l2_m2_r8.00.asc
@@ -0,0 +1 @@
+0.000000 3.141510003733607892 -2.084879755699846063e-17
diff --git a/test/test_driscollhealy/mp_harmonic_re_r8.00.ph.asc b/test/test_driscollhealy/mp_harmonic_re_r8.00.ph.asc
new file mode 100644
index 0000000..219f02c
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_re_r8.00.ph.asc
@@ -0,0 +1,104 @@
+"Time = 0
+0.000000 0.4768960190949341205
+0.062832 0.4731354400496181056
+0.125664 0.4619133783331734744
+0.188496 0.443406739337559197
+0.251327 0.417907108617889822
+0.314159 0.3858170683438172932
+0.376991 0.3476422227979080182
+0.439823 0.3039851008795564358
+0.502655 0.255533712494027232
+0.565487 0.2030525812388497831
+0.628319 0.1473690189180452692
+0.691150 0.08936138578181727909
+0.753982 0.02994455985679102608
+0.816814 -0.02994455985679060975
+0.879646 -0.08936138578181719583
+0.942478 -0.1473690189180451304
+1.005310 -0.2030525812388500884
+1.068142 -0.255533712494027232
+1.130973 -0.3039851008795565468
+1.193805 -0.3476422227979084623
+1.256637 -0.3858170683438171267
+1.319469 -0.4179071086178900996
+1.382301 -0.443406739337558975
+1.445133 -0.461913378333173974
+1.507964 -0.4731354400496178836
+1.570796 -0.4768960190949332878
+1.633628 -0.4731354400496186607
+1.696460 -0.461913378333174196
+1.759292 -0.4434067393375594746
+1.822124 -0.4179071086178899885
+1.884956 -0.3858170683438170712
+1.947787 -0.3476422227979082957
+2.010619 -0.3039851008795566023
+2.073451 -0.255533712494027232
+2.136283 -0.2030525812388493112
+2.199115 -0.1473690189180452137
+2.261947 -0.08936138578181727909
+2.324779 -0.02994455985679057852
+2.387610 0.02994455985679041893
+2.450442 0.08936138578181734848
+2.513274 0.1473690189180449917
+2.576106 0.2030525812388492279
+2.638938 0.2555337124940273985
+2.701770 0.3039851008795561027
+2.764602 0.3476422227979088508
+2.827433 0.3858170683438172932
+2.890265 0.4179071086178900996
+2.953097 0.443406739337559308
+3.015929 0.4619133783331736964
+3.078761 0.4731354400496182722
+3.141593 0.4768960190949315114
+3.204425 0.4731354400496181611
+3.267256 0.4619133783331736409
+3.330088 0.4434067393375588639
+3.392920 0.4179071086178901551
+3.455752 0.3858170683438176263
+3.518584 0.3476422227979086288
+3.581416 0.3039851008795569354
+3.644247 0.2555337124940272875
+3.707079 0.203052581238849561
+3.769911 0.1473690189180455468
+3.832743 0.08936138578181757053
+3.895575 0.02994455985679117527
+3.958407 -0.02994455985679039117
+4.021239 -0.08936138578181766767
+4.084070 -0.1473690189180449084
+4.146902 -0.2030525812388496165
+4.209734 -0.2555337124940270654
+4.272566 -0.3039851008795565468
+4.335398 -0.3476422227979082402
+4.398230 -0.3858170683438171822
+4.461062 -0.4179071086178904326
+4.523893 -0.443406739337559086
+4.586725 -0.4619133783331737519
+4.649557 -0.4731354400496182722
+4.712389 -0.4768960190949321221
+4.775221 -0.4731354400496184387
+4.838053 -0.4619133783331735299
+4.900885 -0.4434067393375591415
+4.963716 -0.4179071086178902106
+5.026548 -0.3858170683438174597
+5.089380 -0.3476422227979087953
+5.152212 -0.3039851008795573795
+5.215044 -0.2555337124940266214
+5.277876 -0.2030525812388495333
+5.340708 -0.1473690189180454357
+5.403539 -0.08936138578181807013
+5.466371 -0.02994455985679158813
+5.529203 0.02994455985679015872
+5.592035 0.08936138578181786196
+5.654867 0.1473690189180450194
+5.717699 0.2030525812388488394
+5.780530 0.2555337124940263438
+5.843362 0.3039851008795563803
+5.906194 0.3476422227979081847
+5.969026 0.3858170683438172932
+6.031858 0.417907108617889711
+6.094690 0.4434067393375583643
+6.157522 0.4619133783331743626
+6.220353 0.4731354400496186052
+6.283185 0.476896019094934509
+
+
diff --git a/test/test_driscollhealy/mp_harmonic_re_r8.00.th.asc b/test/test_driscollhealy/mp_harmonic_re_r8.00.th.asc
new file mode 100644
index 0000000..36c31b5
--- /dev/null
+++ b/test/test_driscollhealy/mp_harmonic_re_r8.00.th.asc
@@ -0,0 +1,55 @@
+"Time = 0
+0.000000 -4.885011147332195149e-16
+0.061600 0.6225628802689283647
+0.123200 0.626011174765530698
+0.184800 0.6200885652671891402
+0.246399 0.6118754000725465936
+0.307999 0.601448949309237535
+0.369599 0.588906848052695775
+0.431199 0.5743658428948120731
+0.492799 0.5579603039292125866
+0.554399 0.5398405168177836666
+0.615999 0.5201707733415493751
+0.677598 0.4991272890280229468
+0.739198 0.4768960190949341205
+0.800798 0.4536703256043045362
+0.862398 0.4296485914087362912
+0.923998 0.4050318059467630238
+0.985598 0.3800211327412000295
+1.047198 0.3548155100560987929
+1.108797 0.3296093187267329361
+1.170397 0.3045901516673371057
+1.231997 0.2799367090384122236
+1.293597 0.2558168423155763316
+1.355197 0.2323857922780937779
+1.416797 0.2097846387827586034
+1.478397 0.1881389819792761153
+1.539996 0.1675578340283152279
+1.601596 0.1481328320161413303
+1.663196 0.1299376592113383189
+1.724796 0.1130277740794584929
+1.786396 0.09744041746750978228
+1.847996 0.0831948475902995338
+1.909596 0.07029291427378291102
+1.971195 0.05871980508769748808
+2.032795 0.0484450861834871474
+2.094395 0.03942394482988391402
+2.155995 0.03159862529246817292
+2.217595 0.02490004252143376146
+2.279195 0.01924954067796907226
+2.340795 0.01456077249599862025
+2.402394 0.01074165773964111219
+2.463994 0.007696404127104090713
+2.525594 0.005327525361079997121
+2.587194 0.003537890131800171169
+2.648794 0.002232672092008527748
+2.710394 0.00132127104105327531
+2.771994 0.0007190932719948640991
+2.833593 0.000349210492446978024
+2.895193 0.0001438537891083224067
+2.956793 4.571824601967693921e-05
+3.018393 9.059240056068421681e-06
+3.079993 5.663901966311766602e-07
+3.141593 -1.077203552463963929e-22
+
+
diff --git a/test/test_driscollhealy/test_driscollhealy.par b/test/test_driscollhealy/test_driscollhealy.par
new file mode 100644
index 0000000..4edac91
--- /dev/null
+++ b/test/test_driscollhealy/test_driscollhealy.par
@@ -0,0 +1,90 @@
+
+ActiveThorns = "CoordBase SymBase Boundary CartGrid3d IOUtil Carpet CarpetLib  CarpetInterp AEILocalInterp InitBase Multipole LoopControl"
+
+#############################################################
+# Grid
+#############################################################
+
+CartGrid3D::type                        = "coordbase"
+CartGrid3D::domain                      = "full"
+CartGrid3D::avoid_origin                = "no"
+
+CoordBase::domainsize                   = minmax
+CoordBase::xmin                         = -10
+CoordBase::ymin                         = -10
+CoordBase::zmin                         = -10
+CoordBase::xmax                         = 10
+CoordBase::ymax                         = 10
+CoordBase::zmax                         = 10
+CoordBase::dx                           = 0.2
+CoordBase::dy                           = 0.2
+CoordBase::dz                           = 0.2
+CoordBase::boundary_size_x_lower        = 2
+CoordBase::boundary_size_y_lower        = 2
+CoordBase::boundary_size_z_lower        = 2
+CoordBase::boundary_shiftout_x_lower    = 0
+CoordBase::boundary_shiftout_y_lower    = 0
+CoordBase::boundary_shiftout_z_lower    = 0
+CoordBase::boundary_size_x_upper        = 2
+CoordBase::boundary_size_y_upper        = 2
+CoordBase::boundary_size_z_upper        = 2
+CoordBase::boundary_shiftout_x_upper    = 0
+CoordBase::boundary_shiftout_y_upper    = 0
+CoordBase::boundary_shiftout_z_upper    = 0
+
+#############################################################
+# Carpet
+#############################################################
+
+Carpet::ghost_size                      = 3
+Carpet::domain_from_coordbase           = "yes"
+Carpet::poison_new_timelevels           = "yes"
+Carpet::check_for_poison                = "no"
+CarpetLib::poison_value                    = 113
+Carpet::init_fill_timelevels            = "yes"
+
+#############################################################
+# CarpetLib
+#############################################################
+
+CarpetLib::poison_new_memory            = "yes"
+
+#############################################################
+# Cactus
+#############################################################
+
+Cactus::terminate                     = "iteration"
+Cactus::cctk_itlast                   = 0
+
+#############################################################
+# Multipole
+#############################################################
+
+Multipole::nradii = 1
+Multipole::ntheta = 51
+Multipole::radius[0] = 8.0
+Multipole::variables = "Multipole::harmonic_re{sw=-2 cmplx='Multipole::harmonic_im' name='harmonic'}"
+Multipole::integration_method = "DriscollHealy"
+
+Multipole::enable_test = "yes"
+Multipole::test_l = 2
+Multipole::test_m = 2
+Multipole::out_1d_every = 1
+
+#############################################################
+# Output
+#############################################################
+
+IO::out_dir                          = $parfile
+IO::out_fileinfo                     = "none"
+
+# Enabling 1D output for the test grid functions would be helpful for
+# localising any failures but it makes the tests dependent on the
+# number of processors, as CarpetIOASCII's output is dependent on
+# this.
+
+# CarpetIOASCII::out1d_vars = "Multipole::harmonics"
+# CarpetIOASCII::out1d_every = 1
+# CarpetIOASCII::out1d_x = yes
+# CarpetIOASCII::out1d_y = yes
+# CarpetIOASCII::out1d_z = yes