Merge /Users/eschnett/Cbeta/carpet

14 files changed, 266 insertions, 85 deletions

diff --git a/Carpet/CarpetLib/interface.ccl b/Carpet/CarpetLib/interface.ccl
index 44b7729d5..6e7462678 100644
--- a/Carpet/CarpetLib/interface.ccl
+++ b/Carpet/CarpetLib/interface.ccl
@@ -4,6 +4,7 @@ IMPLEMENTS: CarpetLib
 
 includes header: defs.hh in defs.hh
 includes header: dist.hh in dist.hh
+includes header: typeprops.hh in typeprops.hh
 
 includes header: bbox.hh in bbox.hh
 includes header: bboxset.hh in bboxset.hh
diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl
index c1d3754ec..45724f25d 100644
--- a/Carpet/CarpetLib/param.ccl
+++ b/Carpet/CarpetLib/param.ccl
@@ -85,6 +85,14 @@ STRING memstat_file "File name in which memstat output is collected (because std
 
 
 
+# Experimental recomposing parameters
+
+BOOLEAN combine_recompose "Recompose all grid functions of one refinement levels at once" STEERABLE=always
+{
+} "no"
+
+
+
 # Experimental communication parameters
 
 BOOLEAN interleave_communications "Try to interleave communications with each other; each processor begins to communicate with its 'right neighbour' in rank, instead of with the root processor" STEERABLE=always
diff --git a/Carpet/CarpetLib/src/commstate.cc b/Carpet/CarpetLib/src/commstate.cc
index b6ed4095a..7a00157c2 100644
--- a/Carpet/CarpetLib/src/commstate.cc
+++ b/Carpet/CarpetLib/src/commstate.cc
@@ -1,6 +1,7 @@
 #include <cassert>
 #include <cstdlib>
 #include <cstring>
+#include <iostream>
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
@@ -115,7 +116,7 @@ void comm_state::step ()
           if (commstate_verbose) {
             CCTK_INFO ("Finished MPI_Irecv");
           }
-          timer.stop (0);
+          timer.stop (procbuf.recvbufsize * datatypesize);
           num_posted_recvs++;
         }
       }
diff --git a/Carpet/CarpetLib/src/commstate.hh b/Carpet/CarpetLib/src/commstate.hh
index 8cc8d40b4..c01f732da 100644
--- a/Carpet/CarpetLib/src/commstate.hh
+++ b/Carpet/CarpetLib/src/commstate.hh
@@ -1,6 +1,7 @@
 #ifndef COMMSTATE_HH
 #define COMMSTATE_HH
 
+#include <cstdlib>
 #include <queue>
 #include <vector>
 
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index 1152c154a..d1d8b73f0 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -62,24 +62,32 @@ template<typename T, int D> class bbox;
 template<typename T, int D> class bboxset;
 template<typename T, int D, typename P> class fulltree;
 
-struct pseudoregion_t;
-struct region_t;
-
-typedef vect<bool,dim>   bvect;
-typedef vect<int,dim>    ivect;
-typedef bbox<int,dim>    ibbox;
-typedef bboxset<int,dim> ibset;
-typedef fulltree<int,dim,pseudoregion_t> ipfulltree;
-
+typedef vect<bool,dim>      bvect;
+typedef vect<int,dim>       ivect;
+typedef vect<CCTK_INT,dim>  jvect;
+typedef vect<CCTK_REAL,dim> rvect;
+typedef bbox<int,dim>       ibbox;
+typedef bbox<CCTK_INT,dim>  jbbox;
+typedef bbox<CCTK_REAL,dim> rbbox;
+typedef bboxset<int,dim>    ibset;
+  
 // (Try to replace these by b2vect and i2vect)
 typedef vect<vect<bool,2>,dim> bbvect;
 typedef vect<vect<int,2>,dim>  iivect;
+typedef vect<vect<CCTK_INT,2>,dim> jjvect;
 
 typedef vect<vect<bool,dim>,2> b2vect;
 typedef vect<vect<int,dim>,2>  i2vect;
 
 
 
+struct pseudoregion_t;
+struct region_t;
+
+typedef fulltree<int,dim,pseudoregion_t> ipfulltree;
+
+
+
 // A general type
 enum centering { error_centered, vertex_centered, cell_centered };
 
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index fa990a2b0..a3d50d08a 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -71,7 +71,6 @@ inline
 int
 dh::this_proc (int const rl, int const c) const
 {
-  // return c % dist::size();
   return h.processor (rl, c);
 }
 
@@ -1006,6 +1005,8 @@ void
 dh::
 recompose (int const rl, bool const do_prolongate)
 {
+  DECLARE_CCTK_PARAMETERS;
+  
   assert (rl>=0 and rl<h.reflevels());
   
   static Timer timer ("dh::recompose");
@@ -1015,13 +1016,31 @@ recompose (int const rl, bool const do_prolongate)
     (*f)->recompose_crop ();
   }
   
-  for (list<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
-    (*f)->recompose_allocate (rl);
+  if (combine_recompose) {
+    // Recompose all grid functions of this refinement levels at once.
+    // This may be faster, but requires more memory.
+    for (list<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+      (*f)->recompose_allocate (rl);
+    }
     for (comm_state state; not state.done(); state.step()) {
-      (*f)->recompose_fill (state, rl, do_prolongate);
+      for (list<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+        (*f)->recompose_fill (state, rl, do_prolongate);
+      }
     }
-    (*f)->recompose_free_old (rl);
-  } // for all grid functions of same vartype
+    for (list<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+      (*f)->recompose_free_old (rl);
+    }
+  } else {
+    // Recompose the grid functions sequentially.  This may be slower,
+    // but requires less memory.  This is the default.
+    for (list<ggf*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) {
+      (*f)->recompose_allocate (rl);
+      for (comm_state state; not state.done(); state.step()) {
+        (*f)->recompose_fill (state, rl, do_prolongate);
+      }
+      (*f)->recompose_free_old (rl);
+    }
+  }
   
   timer.stop (0);
 }
diff --git a/Carpet/CarpetLib/src/dist.cc b/Carpet/CarpetLib/src/dist.cc
index 7ac6f73f6..c870990fb 100644
--- a/Carpet/CarpetLib/src/dist.cc
+++ b/Carpet/CarpetLib/src/dist.cc
@@ -53,6 +53,27 @@ namespace dist {
     MPI_Finalize ();
   }
   
+  
+  
+  // Create an MPI datatype from a C datatype description
+  void create_mpi_datatype (size_t const count,
+                            mpi_struct_descr_t const descr[],
+                            MPI_Datatype & newtype)
+  {
+    int blocklengths[count];
+    MPI_Aint displacements[count];
+    MPI_Datatype types[count];
+    for (size_t n=0; n<count; ++n) {
+      blocklengths [n] = descr[n].blocklength;
+      displacements[n] = descr[n].displacement;
+      types        [n] = descr[n].type;
+    }
+    MPI_Type_struct (count, blocklengths, displacements, types, &newtype);
+    MPI_Type_commit (&newtype);
+  }
+  
+  
+  
   void checkpoint (const char* file, int line) {
     DECLARE_CCTK_PARAMETERS;
     if (verbose) {
diff --git a/Carpet/CarpetLib/src/dist.hh b/Carpet/CarpetLib/src/dist.hh
index ae7953f99..cf4cbd70c 100644
--- a/Carpet/CarpetLib/src/dist.hh
+++ b/Carpet/CarpetLib/src/dist.hh
@@ -30,6 +30,19 @@ namespace dist {
   void pseudoinit (MPI_Comm const c);
   void finalize ();
   
+  // Create MPI datatypes from C structures
+  struct mpi_struct_descr_t {
+    int          blocklength;
+    MPI_Aint     displacement;
+    MPI_Datatype type;
+  };
+  
+  void create_mpi_datatype (size_t const count,
+                            mpi_struct_descr_t const descr[],
+                            MPI_Datatype & newtype);
+  
+  
+  
   // Debugging output
 #define CHECKPOINT dist::checkpoint(__FILE__, __LINE__)
   void checkpoint (const char* file, int line);
diff --git a/Carpet/CarpetLib/src/gf.hh b/Carpet/CarpetLib/src/gf.hh
index 0ca125631..d5feb0a63 100644
--- a/Carpet/CarpetLib/src/gf.hh
+++ b/Carpet/CarpetLib/src/gf.hh
@@ -50,12 +50,14 @@ protected:
   
   virtual gdata* typed_data (int tl, int rl, int c, int ml)
   {
+    data<T>* const vl =
+      this->vectorleader
+      ? (data<T>*)(*this->vectorleader)(tl,rl,c,ml)
+      : NULL;
     return new data<T>(this->varindex,
                        h.refcent, this->transport_operator,
                        this->vectorlength, this->vectorindex,
-                       this->vectorleader
-                       ? (data<T>*)(*this->vectorleader)(tl,rl,c,ml)
-                       : NULL);
+                       vl);
   }
   
   
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 9b4a5f443..9a04af1df 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -272,6 +272,106 @@ local_components (int const rl)
 
 
 
+// Find the refinement level and component to which a grid point
+// belongs.  This uses a tree search over the superregions in the grid
+// struction, which should scale reasonably (O(n log n)) better with
+// the number of componets components.
+void
+gh::
+locate_position (rvect const & rpos,
+                 int const ml,
+                 int const minrl, int const maxrl,
+                 int & rl, int & c, ivect & aligned_ipos) const
+{
+  assert (ml>=0 and ml<mglevels());
+  assert (minrl>=0 and minrl<=maxrl and maxrl<=reflevels());
+  
+  // Try finer levels first
+  for (rl = maxrl-1; rl >= minrl; --rl) {
+    
+    // Align (round) the position to the nearest existing grid point
+    // on this refinement level
+    ivect const str = baseextent(ml,rl).stride();
+    aligned_ipos = ivect(floor(rpos / rvect(str) + rvect(0.5))) * str;
+    
+    gh::cregs const & regs = superregions.AT(rl);
+    
+    // Search all superregions linearly.  Each superregion corresponds
+    // to a "refined region", and the number of superregions is thus
+    // presumably independent of the number of processors.
+    for (size_t r = 0; r < regs.size(); ++r) {
+      region_t const & reg = regs.AT(r);
+      if (reg.extent.contains(aligned_ipos)) {
+        // We found the superregion to which this grid point belongs
+        
+        // Search the superregion hierarchically
+        pseudoregion_t const * const preg =
+          reg.processors->search(aligned_ipos);
+        assert (preg);
+        
+        // We now know the refinement level, component, and index to
+        // which this grid point belongs
+        c = preg->component;
+        return;
+      }
+    }
+  } // for rl
+  
+  // The point does not belong to any component on any refinement
+  // level
+  rl = -1;
+  c = -1;
+}
+
+void
+gh::
+locate_position (ivect const & ipos,
+                 int const ml,
+                 int const minrl, int const maxrl,
+                 int & rl, int & c, ivect & aligned_ipos) const
+{
+  assert (ml>=0 and ml<mglevels());
+  assert (minrl>=0 and minrl<=maxrl and maxrl<=reflevels());
+  
+  // Try finer levels first
+  for (rl = maxrl-1; rl >= minrl; --rl) {
+    
+    // Align (round) the position to the nearest existing grid point
+    // on this refinement level
+    ivect const str = baseextent(ml, rl).stride();
+    aligned_ipos = ivect(floor(rvect(ipos) / rvect(str) + rvect(0.5))) * str;
+    
+    gh::cregs const & regs = superregions.AT(rl);
+    
+    // Search all superregions linearly.  Each superregion corresponds
+    // to a "refined region", and the number of superregions is thus
+    // presumably independent of the number of processors.
+    for (size_t r = 0; r < regs.size(); ++r) {
+      region_t const & reg = regs.AT(r);
+      if (reg.extent.contains(aligned_ipos)) {
+        // We found the superregion to which this grid point belongs
+        
+        // Search the superregion hierarchically
+        pseudoregion_t const * const preg =
+          reg.processors->search(aligned_ipos);
+        assert (preg);
+        
+        // We now know the refinement level, component, and index to
+        // which this grid point belongs
+        c = preg->component;
+        return;
+      }
+    }
+  } // for rl
+  
+  // The point does not belong to any component on any refinement
+  // level
+  rl = -1;
+  c = -1;
+}
+
+
+
 // Time hierarchy management
 
 void
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index ffac38804..b80d71ca3 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -128,6 +128,16 @@ public:
   
   int local_components (const int rl) const;
   
+  void locate_position (rvect const & rpos,
+                        int const ml,
+                        int const minrl, int const maxrl,
+                        int & rl, int & c, ivect & aligned_ipos) const;
+  
+  void locate_position (ivect const & ipos,
+                        int const ml,
+                        int const minrl, int const maxrl,
+                        int & rl, int & c, ivect & aligned_ipos) const;
+  
   // Time hierarchy management
   void add (th * t);
   void remove (th * t);
diff --git a/Carpet/CarpetLib/src/prolongate_3d_o11_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_o11_rf2.cc
index 020158d4b..d731d524d 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_o11_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_o11_rf2.cc
@@ -35,20 +35,19 @@ namespace CarpetLib {
   RT
   coeff (int const i)
   {
-    RT const one = 1;
     static const RT coeffs[ncoeffs] = {
-           63/one*524288,
-      -   819/one*524288,
-         5005/one*524288,
-      - 19305/one*524288,
-        27027/one*262144,
-      - 63063/one*262144,
-       189189/one*262144,
-       135135/one*262144,
-      - 45045/one*524288,
-         9009/one*524288,
-      -  1287/one*524288,
-           91/one*524288
+           63/RT(524288.0),
+      -   819/RT(524288.0),
+         5005/RT(524288.0),
+      - 19305/RT(524288.0),
+        27027/RT(262144.0),
+      - 63063/RT(262144.0),
+       189189/RT(262144.0),
+       135135/RT(262144.0),
+      - 45045/RT(524288.0),
+         9009/RT(524288.0),
+      -  1287/RT(524288.0),
+           91/RT(524288.0)
     };
     return coeffs[i];
   }
@@ -74,7 +73,7 @@ namespace CarpetLib {
     typedef typename typeprops<T>::real RT;
     T res = typeprops<T>::fromreal (0);
     for (int i=0; i<ncoeffs; ++i) {
-      res += coeff<RT>(i) * interp0 (p + i*d1);
+      res += coeff<RT>(i) * interp0<T> (p + i*d1);
     }
     return res;
   }
@@ -90,7 +89,7 @@ namespace CarpetLib {
     typedef typename typeprops<T>::real RT;
     T res = typeprops<T>::fromreal (0);
     for (int i=0; i<ncoeffs; ++i) {
-      res += coeff<RT>(i) * interp1 (p + i*d2, d1);
+      res += coeff<RT>(i) * interp1<T> (p + i*d2, d1);
     }
     return res;
   }
@@ -107,7 +106,7 @@ namespace CarpetLib {
     typedef typename typeprops<T>::real RT;
     T res = typeprops<T>::fromreal (0);
     for (int i=0; i<ncoeffs; ++i) {
-      res += coeff<RT>(i) * interp2 (p + i*d3, d1, d2);
+      res += coeff<RT>(i) * interp2<T> (p + i*d3, d1, d2);
     }
     return res;
   }
@@ -245,7 +244,7 @@ namespace CarpetLib {
     
     // kernel
    l8000:
-    dst[DSTIND3(id,jd,kd)] = interp0 (& src[SRCIND3(is,js,ks)]);
+    dst[DSTIND3(id,jd,kd)] = interp0<T> (& src[SRCIND3(is,js,ks)]);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8001;
@@ -253,7 +252,7 @@ namespace CarpetLib {
     
     // kernel
    l8001:
-    dst[DSTIND3(id,jd,kd)] = interp1 (& src[SRCIND3(is-3,js,ks)], srcdi);
+    dst[DSTIND3(id,jd,kd)] = interp1<T> (& src[SRCIND3(is-3,js,ks)], srcdi);
     i = i+1;
     id = id+1;
     is = is+1;
@@ -277,7 +276,7 @@ namespace CarpetLib {
     
     // kernel
    l8010:
-    dst[DSTIND3(id,jd,kd)] = interp1 (& src[SRCIND3(is,js-3,ks)], srcdj);
+    dst[DSTIND3(id,jd,kd)] = interp1<T> (& src[SRCIND3(is,js-3,ks)], srcdj);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8011;
@@ -286,7 +285,7 @@ namespace CarpetLib {
     // kernel
    l8011:
     dst[DSTIND3(id,jd,kd)] =
-      interp2 (& src[SRCIND3(is-3,js-3,ks)], srcdi, srcdj);
+      interp2<T> (& src[SRCIND3(is-3,js-3,ks)], srcdi, srcdj);
     i = i+1;
     id = id+1;
     is = is+1;
@@ -326,7 +325,7 @@ namespace CarpetLib {
     
     // kernel
    l8100:
-    dst[DSTIND3(id,jd,kd)] = interp1 (& src[SRCIND3(is,js,ks-3)], srcdk);
+    dst[DSTIND3(id,jd,kd)] = interp1<T> (& src[SRCIND3(is,js,ks-3)], srcdk);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8101;
@@ -335,7 +334,7 @@ namespace CarpetLib {
     // kernel
    l8101:
     dst[DSTIND3(id,jd,kd)] =
-      interp2 (& src[SRCIND3(is-3,js,ks-3)], srcdi, srcdj);
+      interp2<T> (& src[SRCIND3(is-3,js,ks-3)], srcdi, srcdj);
     i = i+1;
     id = id+1;
     is = is+1;
@@ -360,7 +359,7 @@ namespace CarpetLib {
     // kernel
    l8110:
     dst[DSTIND3(id,jd,kd)] =
-      interp2 (& src[SRCIND3(is,js-3,ks-3)], srcdj, srcdk);
+      interp2<T> (& src[SRCIND3(is,js-3,ks-3)], srcdj, srcdk);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8111;
@@ -370,7 +369,7 @@ namespace CarpetLib {
    l8111:
     {
       dst[DSTIND3(id,jd,kd)] =
-        interp3 (& src[SRCIND3(is-3,js-3,ks-3)], srcdi, srcdj, srcdk);
+        interp3<T> (& src[SRCIND3(is-3,js-3,ks-3)], srcdi, srcdj, srcdk);
     }
     i = i+1;
     id = id+1;
diff --git a/Carpet/CarpetLib/src/prolongate_3d_o7_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_o7_rf2.cc
index fcf4d710f..a7341139f 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_o7_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_o7_rf2.cc
@@ -35,16 +35,15 @@ namespace CarpetLib {
   RT
   coeff (int const i)
   {
-    RT const one = 1;
     static const RT coeffs[ncoeffs] = {
-      -   5*one/2048,
-         49*one/2048,
-      - 245*one/2048,
-       1225*one/2048,
-       1225*one/2048,
-      - 245*one/2048,
-         49*one/2048,
-      -   5*one/2048
+      -   5/RT(2048.0),
+         49/RT(2048.0),
+      - 245/RT(2048.0),
+       1225/RT(2048.0),
+       1225/RT(2048.0),
+      - 245/RT(2048.0),
+         49/RT(2048.0),
+      -   5/RT(2048.0)
     };
     return coeffs[i];
   }
@@ -70,7 +69,7 @@ namespace CarpetLib {
     typedef typename typeprops<T>::real RT;
     T res = typeprops<T>::fromreal (0);
     for (int i=0; i<ncoeffs; ++i) {
-      res += coeff<RT>(i) * interp0 (p + i*d1);
+      res += coeff<RT>(i) * interp0<T> (p + i*d1);
     }
     return res;
   }
@@ -86,7 +85,7 @@ namespace CarpetLib {
     typedef typename typeprops<T>::real RT;
     T res = typeprops<T>::fromreal (0);
     for (int i=0; i<ncoeffs; ++i) {
-      res += coeff<RT>(i) * interp1 (p + i*d2, d1);
+      res += coeff<RT>(i) * interp1<T> (p + i*d2, d1);
     }
     return res;
   }
@@ -103,7 +102,7 @@ namespace CarpetLib {
     typedef typename typeprops<T>::real RT;
     T res = typeprops<T>::fromreal (0);
     for (int i=0; i<ncoeffs; ++i) {
-      res += coeff<RT>(i) * interp2 (p + i*d3, d1, d2);
+      res += coeff<RT>(i) * interp2<T> (p + i*d3, d1, d2);
     }
     return res;
   }
@@ -241,7 +240,7 @@ namespace CarpetLib {
     
     // kernel
    l8000:
-    dst[DSTIND3(id,jd,kd)] = interp0 (& src[SRCIND3(is,js,ks)]);
+    dst[DSTIND3(id,jd,kd)] = interp0<T> (& src[SRCIND3(is,js,ks)]);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8001;
@@ -249,7 +248,7 @@ namespace CarpetLib {
     
     // kernel
    l8001:
-    dst[DSTIND3(id,jd,kd)] = interp1 (& src[SRCIND3(is-3,js,ks)], srcdi);
+    dst[DSTIND3(id,jd,kd)] = interp1<T> (& src[SRCIND3(is-3,js,ks)], srcdi);
     i = i+1;
     id = id+1;
     is = is+1;
@@ -273,7 +272,7 @@ namespace CarpetLib {
     
     // kernel
    l8010:
-    dst[DSTIND3(id,jd,kd)] = interp1 (& src[SRCIND3(is,js-3,ks)], srcdj);
+    dst[DSTIND3(id,jd,kd)] = interp1<T> (& src[SRCIND3(is,js-3,ks)], srcdj);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8011;
@@ -282,7 +281,7 @@ namespace CarpetLib {
     // kernel
    l8011:
     dst[DSTIND3(id,jd,kd)] =
-      interp2 (& src[SRCIND3(is-3,js-3,ks)], srcdi, srcdj);
+      interp2<T> (& src[SRCIND3(is-3,js-3,ks)], srcdi, srcdj);
     i = i+1;
     id = id+1;
     is = is+1;
@@ -322,7 +321,7 @@ namespace CarpetLib {
     
     // kernel
    l8100:
-    dst[DSTIND3(id,jd,kd)] = interp1 (& src[SRCIND3(is,js,ks-3)], srcdk);
+    dst[DSTIND3(id,jd,kd)] = interp1<T> (& src[SRCIND3(is,js,ks-3)], srcdk);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8101;
@@ -331,7 +330,7 @@ namespace CarpetLib {
     // kernel
    l8101:
     dst[DSTIND3(id,jd,kd)] =
-      interp2 (& src[SRCIND3(is-3,js,ks-3)], srcdi, srcdj);
+      interp2<T> (& src[SRCIND3(is-3,js,ks-3)], srcdi, srcdj);
     i = i+1;
     id = id+1;
     is = is+1;
@@ -356,7 +355,7 @@ namespace CarpetLib {
     // kernel
    l8110:
     dst[DSTIND3(id,jd,kd)] =
-      interp2 (& src[SRCIND3(is,js-3,ks-3)], srcdj, srcdk);
+      interp2<T> (& src[SRCIND3(is,js-3,ks-3)], srcdj, srcdk);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8111;
@@ -366,7 +365,7 @@ namespace CarpetLib {
    l8111:
     {
       dst[DSTIND3(id,jd,kd)] =
-        interp3 (& src[SRCIND3(is-3,js-3,ks-3)], srcdi, srcdj, srcdk);
+        interp3<T> (& src[SRCIND3(is-3,js-3,ks-3)], srcdi, srcdj, srcdk);
     }
     i = i+1;
     id = id+1;
diff --git a/Carpet/CarpetLib/src/prolongate_3d_o9_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_o9_rf2.cc
index 8e2d6fc18..cd0d6038b 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_o9_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_o9_rf2.cc
@@ -35,18 +35,17 @@ namespace CarpetLib {
   RT
   coeff (int const i)
   {
-    RT const one = 1;
     static const RT coeffs[ncoeffs] = {
-      -   35*one/65536,
-         385*one/65536,
-      -  495*one/16384,
-        1617*one/16384,
-      - 8085*one/32768,
-       24255*one/32768,
-        8085*one/16384,
-      - 1155*one/16384,
-         693*one/65536,
-      -   55*one/65536
+      -   35/RT(65536.0),
+         385/RT(65536.0),
+      -  495/RT(16384.0),
+        1617/RT(16384.0),
+      - 8085/RT(32768.0),
+       24255/RT(32768.0),
+        8085/RT(16384.0),
+      - 1155/RT(16384.0),
+         693/RT(65536.0),
+      -   55/RT(65536.0)
     };
     return coeffs[i];
   }
@@ -72,7 +71,7 @@ namespace CarpetLib {
     typedef typename typeprops<T>::real RT;
     T res = typeprops<T>::fromreal (0);
     for (int i=0; i<ncoeffs; ++i) {
-      res += coeff<RT>(i) * interp0 (p + i*d1);
+      res += coeff<RT>(i) * interp0<T> (p + i*d1);
     }
     return res;
   }
@@ -88,7 +87,7 @@ namespace CarpetLib {
     typedef typename typeprops<T>::real RT;
     T res = typeprops<T>::fromreal (0);
     for (int i=0; i<ncoeffs; ++i) {
-      res += coeff<RT>(i) * interp1 (p + i*d2, d1);
+      res += coeff<RT>(i) * interp1<T> (p + i*d2, d1);
     }
     return res;
   }
@@ -105,7 +104,7 @@ namespace CarpetLib {
     typedef typename typeprops<T>::real RT;
     T res = typeprops<T>::fromreal (0);
     for (int i=0; i<ncoeffs; ++i) {
-      res += coeff<RT>(i) * interp2 (p + i*d3, d1, d2);
+      res += coeff<RT>(i) * interp2<T> (p + i*d3, d1, d2);
     }
     return res;
   }
@@ -243,7 +242,7 @@ namespace CarpetLib {
     
     // kernel
    l8000:
-    dst[DSTIND3(id,jd,kd)] = interp0 (& src[SRCIND3(is,js,ks)]);
+    dst[DSTIND3(id,jd,kd)] = interp0<T> (& src[SRCIND3(is,js,ks)]);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8001;
@@ -251,7 +250,7 @@ namespace CarpetLib {
     
     // kernel
    l8001:
-    dst[DSTIND3(id,jd,kd)] = interp1 (& src[SRCIND3(is-3,js,ks)], srcdi);
+    dst[DSTIND3(id,jd,kd)] = interp1<T> (& src[SRCIND3(is-3,js,ks)], srcdi);
     i = i+1;
     id = id+1;
     is = is+1;
@@ -275,7 +274,7 @@ namespace CarpetLib {
     
     // kernel
    l8010:
-    dst[DSTIND3(id,jd,kd)] = interp1 (& src[SRCIND3(is,js-3,ks)], srcdj);
+    dst[DSTIND3(id,jd,kd)] = interp1<T> (& src[SRCIND3(is,js-3,ks)], srcdj);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8011;
@@ -284,7 +283,7 @@ namespace CarpetLib {
     // kernel
    l8011:
     dst[DSTIND3(id,jd,kd)] =
-      interp2 (& src[SRCIND3(is-3,js-3,ks)], srcdi, srcdj);
+      interp2<T> (& src[SRCIND3(is-3,js-3,ks)], srcdi, srcdj);
     i = i+1;
     id = id+1;
     is = is+1;
@@ -324,7 +323,7 @@ namespace CarpetLib {
     
     // kernel
    l8100:
-    dst[DSTIND3(id,jd,kd)] = interp1 (& src[SRCIND3(is,js,ks-3)], srcdk);
+    dst[DSTIND3(id,jd,kd)] = interp1<T> (& src[SRCIND3(is,js,ks-3)], srcdk);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8101;
@@ -333,7 +332,7 @@ namespace CarpetLib {
     // kernel
    l8101:
     dst[DSTIND3(id,jd,kd)] =
-      interp2 (& src[SRCIND3(is-3,js,ks-3)], srcdi, srcdj);
+      interp2<T> (& src[SRCIND3(is-3,js,ks-3)], srcdi, srcdj);
     i = i+1;
     id = id+1;
     is = is+1;
@@ -358,7 +357,7 @@ namespace CarpetLib {
     // kernel
    l8110:
     dst[DSTIND3(id,jd,kd)] =
-      interp2 (& src[SRCIND3(is,js-3,ks-3)], srcdj, srcdk);
+      interp2<T> (& src[SRCIND3(is,js-3,ks-3)], srcdj, srcdk);
     i = i+1;
     id = id+1;
     if (i < regiext) goto l8111;
@@ -368,7 +367,7 @@ namespace CarpetLib {
    l8111:
     {
       dst[DSTIND3(id,jd,kd)] =
-        interp3 (& src[SRCIND3(is-3,js-3,ks-3)], srcdi, srcdj, srcdk);
+        interp3<T> (& src[SRCIND3(is-3,js-3,ks-3)], srcdi, srcdj, srcdk);
     }
     i = i+1;
     id = id+1;