Added private STRING parameter Extract::interpolation_operator defaulting to

"Lagrange polynomial interpolation".

See LocalInterp for other interpolation_operator values.

At the moment LocalInterp will use
  "Lagrange polynomial interpolation (tensor product)".

Changed the testsuites to use
  "Lagrange polynomial interpolation (maximum degree)"
to get the same behaviour as before and the testsuites to pass again.

They didn't pass before, because the interpolator gives slightly different
answers for the "tensor product" version. The Moncrief Q's should be
zero for Schwarzschild. The result is for example 10^-6 rather than
10^-5 which was reported as a fail.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Extract/trunk@77 5301f0c2-dbc4-4cee-b2f5-8d7afba4d129

8 files changed, 35 insertions, 10 deletions

diff --git a/param.ccl b/param.ccl
index 13b1826..1173c70 100644
--- a/param.ccl
+++ b/param.ccl
@@ -22,7 +22,6 @@ KEYWORD timecoord "Which time coordinate to use"
   "both" :: ""
 } "both"
 
-
 BOOLEAN all_modes "Extract: all l,m modes up to l"
 {
 } "yes"
@@ -145,6 +144,11 @@ REAL Cauchy_dr "Gridspacing for Cauchy data extraction"
 *:* :: ""
 } 0.2
 
+STRING interpolation_operator "Interpolation operator to use (check LocalInterp)"
+{
+  .+ :: ""
+} "Lagrange polynomial interpolation"
+
 INT interpolation_order "Order for interpolation" STEERABLE = ALWAYS
 {
   1:4 :: "Choose between first and forth order interpolation"
diff --git a/src/D3_extract.F b/src/D3_extract.F
index 9f6fcad..47a7fb3 100644
--- a/src/D3_extract.F
+++ b/src/D3_extract.F
@@ -4,7 +4,7 @@
 c     ==================================================================
 
       SUBROUTINE D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,igrid,
-     &     origin,myproc,interpolation_order,Nt,Np,all_modes,
+     &     origin,myproc,interpolation_operator,interpolation_order,Nt,Np,all_modes,
      &     l,m,x,y,z,Dx,Dy,Dz,Psi_power,Psi,
      &     g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
      &     eta,ADMmass,momentum,spin,mass,rsch,Qodd,Qeven,Extract_temp3d,dtaudt)
@@ -85,6 +85,8 @@ c     Input variables
      &     x,y,z
       CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) ::
      &     Extract_temp3d
+      CCTK_STRING,INTENT(IN) ::
+     &     interpolation_operator
 
 c     Output variables
       CCTK_REAL,INTENT(OUT) :: 
@@ -162,7 +164,8 @@ c
 c     ------------------------------------------------------------------
 
       CALL D3_to_D2(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
-     &     Psi_power,origin,myproc,interpolation_order,Dx,Dy,Dz,Psi,
+     &     Psi_power,origin,myproc,interpolation_operator,interpolation_order,
+     &     Dx,Dy,Dz,Psi,
      &     g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
      &     x,y,z,eta,Nt,Np,theta,phi,Psis,g00s,gxxs,gxys,gxzs,
      &     gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,dgzzs,
diff --git a/src/D3_extract_int.F b/src/D3_extract_int.F
index 07d7dab..ce14b37 100644
--- a/src/D3_extract_int.F
+++ b/src/D3_extract_int.F
@@ -8,7 +8,7 @@ c     ------------------------------------------------------------------
       INTERFACE
 
       SUBROUTINE D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,igrid,
-     &        origin,myproc,interpolation_order,
+     &        origin,myproc,interpolation_operator,interpolation_order,
      &        Nt,Np,all_modes,l,m,x,y,z,Dx,Dy,Dz,Psi_power,Psi,
      &        g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
      &        eta,ADMmass,momentum,spin,mass,rsch,Qodd,Qeven,
@@ -36,6 +36,9 @@ c     ------------------------------------------------------------------
       CCTK_REAL,INTENT(OUT) :: 
      &     ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt,
      &     momentum(3),spin(3)
+      CCTK_STRING,INTENT(IN) ::
+     &     interpolation_operator
+
 
       END SUBROUTINE
          
diff --git a/src/D3_to_D2.F b/src/D3_to_D2.F
index 7bc7586..da2294e 100644
--- a/src/D3_to_D2.F
+++ b/src/D3_to_D2.F
@@ -4,7 +4,8 @@
 c     ========================================================================
 
       SUBROUTINE D3_to_D2(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
-     &     Psi_power,origin,myproc,interpolation_order,Dx,Dy,Dz,Psi,
+     &     Psi_power,origin,myproc,interpolation_operator,interpolation_order,
+     &     Dx,Dy,Dz,Psi,
      &     g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
      &     x,y,z,eta,Nt,Np,theta,phi,
      &     Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,
@@ -46,6 +47,8 @@ c     Input variables
      &     Extract_temp3d
       LOGICAL ::
      &     do_ADMmass(2)
+      CCTK_STRING,INTENT(IN) :: 
+     &     interpolation_operator
 
 c     Output variables
 
@@ -77,7 +80,9 @@ c     Local variables, here only
       CCTK_POINTER, dimension(8) :: out_arrays
       CCTK_INT, dimension(8)     :: out_array_type_codes
 
-      character(30) options_string
+      character(128) options_string
+      character(128) operator
+      CCTK_INT nchars
 
 c     ------------------------------------------------------------------------
 
@@ -120,7 +125,9 @@ c     --------------------------------------------------------------------
       param_table_handle = -1
       coord_system_handle = -1
 
-      call CCTK_InterpHandle (interp_handle,"Lagrange polynomial interpolation")
+
+      call CCTK_FortranString(nchars,interpolation_operator,operator)
+      call CCTK_InterpHandle (interp_handle,operator)
       if (interp_handle .lt. 0) then
         call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
       endif
diff --git a/src/D3_to_D2_int.F b/src/D3_to_D2_int.F
index 8ca896c..db013a2 100644
--- a/src/D3_to_D2_int.F
+++ b/src/D3_to_D2_int.F
@@ -8,7 +8,8 @@ c     ------------------------------------------------------------------
       INTERFACE
 
       SUBROUTINE D3_to_D2(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
-     &        Psi_power,origin,myproc,interpolation_order,Dx,Dy,Dz,Psi,
+     &        Psi_power,origin,myproc,interpolation_operator,interpolation_order,
+     &        Dx,Dy,Dz,Psi,
      &        g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
      &        x,y,z,eta,Nt,Np,theta,phi,Psis,g00s,gxxs,gxys,
      &        gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,
@@ -38,6 +39,9 @@ c     ------------------------------------------------------------------
      &     spin_int1,spin_int2,spin_int3
       CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) ::
      &     Extract_temp3d
+      CCTK_STRING,INTENT(IN) ::
+     &     interpolation_operator
+
       END SUBROUTINE
 
       END INTERFACE
diff --git a/src/Extract.F b/src/Extract.F
index 1ea5e86..ace99fb 100644
--- a/src/Extract.F
+++ b/src/Extract.F
@@ -378,7 +378,8 @@ c           Do extraction at each radius
             extract_at_each_radius: DO WHILE (radius < r2)
 
                CALL D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
-     &           igrid,orig,myproc,interpolation_order,Nt,Np,all_modes,lmode,
+     &           igrid,orig,myproc,interpolation_operator,interpolation_order,
+     &           Nt,Np,all_modes,lmode,
      &           mmode,x_1d,y_1d,z_1d,Dx,Dy,Dz,Psi_power,Psi,g00,
      &           gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,kyy,kyz,kzz,
      &           radius,ADMmass,momentum,spin,mass,rsch,
@@ -549,7 +550,8 @@ c     Cannot use the conformal equation for ADM mass now
             END IF
 
             CALL D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
-     &           igrid,orig,myproc,interpolation_order,Nt,Np,all_modes,
+     &           igrid,orig,myproc,interpolation_operator,interpolation_order,
+     &           Nt,Np,all_modes,
      &           lmode,mmode,x_1d,y_1d,z_1d,
      &           Dx,Dy,Dz,Psi_power,Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,
      &           kyy,kyz,kzz,radius,ADMmass,momentum,spin,mass,rsch,
diff --git a/test/extract_physical.par b/test/extract_physical.par
index b75eade..c381b9c 100644
--- a/test/extract_physical.par
+++ b/test/extract_physical.par
@@ -30,6 +30,7 @@ ioascii::out1d_vars = "admbase::alp"
 ### Extraction parameters
 
 extract::interpolation_order = 1
+extract::interpolation_operator="Lagrange polynomial interpolation (maximum degree)"
 extract::num_detectors = 2
 extract::itout         = 1
 
diff --git a/test/test_extract.par b/test/test_extract.par
index db9de22..e09551d 100644
--- a/test/test_extract.par
+++ b/test/test_extract.par
@@ -30,6 +30,7 @@ ioascii::out1d_vars = "admbase::alp"
 ### Extraction parameters
 
 extract::interpolation_order = 1
+extract::interpolation_operator="Lagrange polynomial interpolation (maximum degree)"
 extract::num_detectors = 2
 extract::itout         = 1