Import Carpet

Ignore-this: 309b4dd613f4af2b84aa5d6743fdb6b3

5 files changed, 276 insertions, 24 deletions

diff --git a/Carpet/CarpetInterp2/README b/Carpet/CarpetInterp2/README
index 6480393c9..de19cb6e4 100644
--- a/Carpet/CarpetInterp2/README
+++ b/Carpet/CarpetInterp2/README
@@ -1,8 +1,9 @@
 Cactus Code Thorn CarpetInterp2
-Thorn Author(s)     : Erik Schnetter <schnetter@cct.lsu.edu>
-Thorn Maintainer(s) : Erik Schnetter <schnetter@cct.lsu.edu>
+Author(s)    : Erik Schnetter <schnetter@cct.lsu.edu>
+Maintainer(s): Erik Schnetter <schnetter@cct.lsu.edu>
+Licence      : GPLv2+
 --------------------------------------------------------------------------
 
-Purpose of the thorn:
+1. Purpose
 
 This thorn provides a parallel interpolator for Carpet.
diff --git a/Carpet/CarpetInterp2/interface.ccl b/Carpet/CarpetInterp2/interface.ccl
index 1a8e1c177..15b3919c7 100644
--- a/Carpet/CarpetInterp2/interface.ccl
+++ b/Carpet/CarpetInterp2/interface.ccl
@@ -46,8 +46,6 @@ CCTK_INT FUNCTION                                \
          CCTK_POINTER          IN  ddadxdx)
 USES FUNCTION MultiPatch_GlobalToLocal
 
-
-
 CCTK_INT FUNCTION                                      \
     InterpGridArrays                                   \
         (CCTK_POINTER_TO_CONST IN cctkGH,              \
@@ -62,3 +60,52 @@ CCTK_INT FUNCTION                                      \
 PROVIDES FUNCTION InterpGridArrays      \
     WITH CarpetInterp2_InterpGridArrays \
     LANGUAGE C
+
+CCTK_INT FUNCTION                                      \
+    MultiPatchInterpGridArrays                         \
+        (CCTK_POINTER_TO_CONST IN cctkGH,              \
+         CCTK_INT              IN N_dims,              \
+         CCTK_INT              IN order,               \
+         CCTK_INT              IN N_interp_points,     \
+         CCTK_INT        ARRAY IN interp_maps,         \
+         CCTK_POINTER_TO_CONST IN interp_coords,       \
+         CCTK_INT              IN N_input_arrays,      \
+         CCTK_INT        ARRAY IN input_array_indices, \
+         CCTK_INT              IN N_output_arrays,     \
+         CCTK_POINTER          IN output_arrays)
+PROVIDES FUNCTION MultiPatchInterpGridArrays      \
+    WITH CarpetInterp2_MultiPatchInterpGridArrays \
+    LANGUAGE C
+
+
+
+CCTK_POINTER FUNCTION                                  \
+    Interp2GridArraysSetup                             \
+        (CCTK_POINTER_TO_CONST IN cctkGH,              \
+         CCTK_INT              IN N_dims,              \
+         CCTK_INT              IN order,               \
+         CCTK_INT              IN N_interp_points,     \
+         CCTK_POINTER_TO_CONST IN interp_coords)
+PROVIDES FUNCTION Interp2GridArraysSetup      \
+    WITH CarpetInterp2_Interp2GridArraysSetup \
+    LANGUAGE C
+
+CCTK_INT FUNCTION                                      \
+    Interp2GridArrays                                  \
+        (CCTK_POINTER_TO_CONST IN cctkGH,              \
+         CCTK_POINTER_TO_CONST IN interp_handle,       \
+         CCTK_INT              IN N_input_arrays,      \
+         CCTK_INT        ARRAY IN input_array_indices, \
+         CCTK_INT              IN N_output_arrays,     \
+         CCTK_POINTER          IN output_arrays)
+PROVIDES FUNCTION Interp2GridArrays      \
+    WITH CarpetInterp2_Interp2GridArrays \
+    LANGUAGE C
+
+CCTK_INT FUNCTION                                      \
+    Interp2GridArraysFree                              \
+        (CCTK_POINTER_TO_CONST IN cctkGH,              \
+         CCTK_POINTER          IN interp_handle)
+PROVIDES FUNCTION Interp2GridArraysFree      \
+    WITH CarpetInterp2_Interp2GridArraysFree \
+    LANGUAGE C
diff --git a/Carpet/CarpetInterp2/src/fasterp.cc b/Carpet/CarpetInterp2/src/fasterp.cc
index 039c01a15..4b5126b78 100644
--- a/Carpet/CarpetInterp2/src/fasterp.cc
+++ b/Carpet/CarpetInterp2/src/fasterp.cc
@@ -25,7 +25,9 @@ namespace CarpetInterp2 {
   int const ipoison = -1234567890;
   CCTK_REAL const poison = -1.0e+12;
   
+  int get_poison (int const &) CCTK_ATTRIBUTE_CONST;
   int get_poison (int const &) { return ipoison; }
+  CCTK_REAL get_poison (CCTK_REAL const &) CCTK_ATTRIBUTE_CONST;
   CCTK_REAL get_poison (CCTK_REAL const &) { return poison; }
   
   template <typename T>
@@ -72,7 +74,9 @@ namespace CarpetInterp2 {
       {                                                                 \
         sizeof s.name / sizeof(type), /* count elements */              \
           (char*)&s.name - (char*)&s, /* offsetof doesn't work (why?) */ \
-          dist::datatype<type>()      /* find MPI datatype */           \
+          dist::mpi_datatype<type>(), /* find MPI datatype */           \
+          STRINGIFY(name), /* field name */                             \
+          STRINGIFY(type), /* type name */                              \
           }
       dist::mpi_struct_descr_t const descr[] = {
         ENTRY(int, mrc),
@@ -83,11 +87,12 @@ namespace CarpetInterp2 {
 #endif
         ENTRY(int, ind3d),
         ENTRY(CCTK_REAL, offset),
-        {1, sizeof(s), MPI_UB}
+        {1, sizeof(s), MPI_UB, "MPI_UB", "MPI_UB"}
       };
 #undef ENTRY
-      dist::create_mpi_datatype
-        (sizeof(descr) / sizeof(descr[0]), descr, newtype);
+      newtype =
+        dist::create_mpi_datatype (sizeof descr / sizeof descr[0], descr,
+                                   "fasterp_iloc_t::mpi_datatype", sizeof s);
       initialised = true;
     }
     return newtype;
@@ -240,6 +245,7 @@ namespace CarpetInterp2 {
       // round is not available with PGI compilers
       // CCTK_REAL const rx = round(x);
       CCTK_REAL const rx = floor(x+0.5);
+#warning "TODO: make eps a relative error my taking max(abs(xmin),abs(xmax)) into account"
       if (abs(x - rx) < eps) {
         // The interpolation point coincides with a grid point; no
         // interpolation is necessary (this is a special case)
@@ -478,7 +484,8 @@ namespace CarpetInterp2 {
   fasterp_setup_t (cGH const * restrict const cctkGH,
                    fasterp_glocs_t const & locations,
                    int const order_)
-    : order (order_)
+    : order (order_),
+      regridding_epoch (Carpet::regridding_epoch)
   {
     // Some global properties
     int const npoints = locations.size();
@@ -1026,6 +1033,12 @@ namespace CarpetInterp2 {
   {
     DECLARE_CCTK_PARAMETERS;
     
+    // Check regridding epoch
+    if (regridding_epoch != Carpet::regridding_epoch) {
+      CCTK_WARN (CCTK_WARN_ALERT,
+                 "The Carpet grid structure was changed since this fasterp_setup was created");
+    }
+    
     // Desired time level
     int const tl = 0;           // current time
     
@@ -1072,19 +1085,19 @@ namespace CarpetInterp2 {
       
       MPI_Irecv (& recv_points.AT(recv_proc.offset * nvars),
                  recv_proc.npoints * nvars,
-                 dist::datatype<CCTK_REAL>(), recv_proc.p, mpi_tag,
+                 dist::mpi_datatype<CCTK_REAL>(), recv_proc.p, mpi_tag,
                  comm_world, & recv_reqs.AT(pp));
 #ifdef CARPET_DEBUG
       MPI_Irecv (& recv_pn.AT(recv_proc.offset),
                  recv_proc.npoints * 2,
-                 dist::datatype<int>(), recv_proc.p, mpi_tag,
+                 dist::mpi_datatype<int>(), recv_proc.p, mpi_tag,
                  comm_world, & recv_reqs_pn.AT(pp));
 #endif
     }
     
     // Interpolate data and post Isends
     if (verbose) CCTK_INFO ("Interpolating and posting MPI_Isends");
-    // TODO: Use one arrays per processor?
+    // TODO: Use one array per processor?
     vector<CCTK_REAL> send_points (send_descr.npoints * nvars);
     fill_with_poison (send_points);
     vector<MPI_Request> send_reqs (send_descr.procs.size());
@@ -1107,6 +1120,8 @@ namespace CarpetInterp2 {
         int const rl = send_comp.mrc.rl;
         int const c  = send_comp.mrc.c;
         
+        int const lc = Carpet::vhh.AT(m)->get_local_component(rl,c);
+        
         // Get pointers to 3D data
         vector<CCTK_REAL const *> varptrs (nvars);
         for (size_t v=0; v<nvars; ++v) {
@@ -1117,11 +1132,13 @@ namespace CarpetInterp2 {
           varptrs.AT(v) =
             (CCTK_REAL const *)
             (* Carpet::arrdata.AT(gi).AT(m).data.AT(vi))
-            (tl, rl, c, Carpet::mglevel)->storage();
+            (tl, rl, lc, Carpet::mglevel)->storage();
           assert (varptrs.AT(v));
         }
         
-#pragma omp parallel for
+        // TODO: This loops seems unbalanced.  Maybe the different
+        // interpolations have different costs.
+#pragma omp parallel for schedule (dynamic, 1000)
         for (int n=0; n<int(send_comp.locs.size()); ++n) {
           size_t const ind = (send_comp.offset + n) * nvars;
           send_comp.locs.AT(n).interpolate
@@ -1161,12 +1178,12 @@ namespace CarpetInterp2 {
       
       MPI_Isend (& send_points.AT(send_proc.offset * nvars),
                  send_proc.npoints * nvars,
-                 dist::datatype<CCTK_REAL>(), send_proc.p, mpi_tag,
+                 dist::mpi_datatype<CCTK_REAL>(), send_proc.p, mpi_tag,
                  comm_world, & send_reqs.AT(pp));
 #ifdef CARPET_DEBUG
       MPI_Isend (& send_pn.AT(send_proc.offset),
                  send_proc.npoints * 2,
-                 dist::datatype<int>(), send_proc.p, mpi_tag,
+                 dist::mpi_datatype<int>(), send_proc.p, mpi_tag,
                  comm_world, & send_reqs_pn.AT(pp));
 #endif
     } // for pp
diff --git a/Carpet/CarpetInterp2/src/fasterp.hh b/Carpet/CarpetInterp2/src/fasterp.hh
index 393fd5144..1a58e8222 100644
--- a/Carpet/CarpetInterp2/src/fasterp.hh
+++ b/Carpet/CarpetInterp2/src/fasterp.hh
@@ -38,7 +38,7 @@ namespace CarpetInterp2 {
     static int components;
     
     static void determine_mrc_info ();
-    static int get_max_ind ()
+    static int get_max_ind () CCTK_ATTRIBUTE_PURE
     {
       return maps * reflevels * components;
     }
@@ -63,12 +63,12 @@ namespace CarpetInterp2 {
     mrc_t (int ind);
     
     // Convert a mrc into an integer index
-    int get_ind () const
+    int get_ind () const CCTK_ATTRIBUTE_PURE
     {
       return c + components * (rl + reflevels * m);
     }
     
-    bool operator== (mrc_t const & a) const
+    bool operator== (mrc_t const & a) const CCTK_ATTRIBUTE_PURE
     {
       return m == a.m and rl == a.rl and c == a.c;
     }
@@ -116,7 +116,7 @@ namespace CarpetInterp2 {
         coords[d].resize(n);
       }
     }
-    size_t size () const { return coords[0].size(); }
+    size_t size () const CCTK_ATTRIBUTE_PURE { return coords[0].size(); }
   };
   
   // A local location, given by map and local coordinates
@@ -130,7 +130,7 @@ namespace CarpetInterp2 {
         coords[d].resize(n);
       }
     }
-    size_t size () const { return maps.size(); }
+    size_t size () const CCTK_ATTRIBUTE_PURE { return maps.size(); }
   };
   
   // An integer location, given by map, refinement level, and
@@ -146,7 +146,7 @@ namespace CarpetInterp2 {
     int ind3d;                  // closest grid point
     rvect offset;               // in terms of grid points
     
-    static MPI_Datatype mpi_datatype ();
+    static MPI_Datatype mpi_datatype () CCTK_ATTRIBUTE_CONST;
     
     void output (ostream& os) const;
   };
@@ -277,6 +277,8 @@ namespace CarpetInterp2 {
     send_descr_t send_descr;
     int order;
     
+    int regridding_epoch;
+    
     void
     setup (cGH const * restrict cctkGH,
            fasterp_llocs_t const & locations);
@@ -300,7 +302,7 @@ namespace CarpetInterp2 {
     
     size_t
     get_npoints ()
-      const
+      const CCTK_ATTRIBUTE_PURE
     {
       return recv_descr.npoints;
     }
diff --git a/Carpet/CarpetInterp2/src/interp2.cc b/Carpet/CarpetInterp2/src/interp2.cc
index f8a17df99..ee026744d 100644
--- a/Carpet/CarpetInterp2/src/interp2.cc
+++ b/Carpet/CarpetInterp2/src/interp2.cc
@@ -75,4 +75,189 @@ namespace CarpetInterp2 {
     return 0;
   }
   
+  extern "C"
+  CCTK_INT
+  CarpetInterp2_MultiPatchInterpGridArrays (CCTK_POINTER_TO_CONST const cctkGH_,
+                                            CCTK_INT const N_dims,
+                                            CCTK_INT const order,
+                                            CCTK_INT const N_interp_points,
+                                            CCTK_INT const * const interp_maps,
+                                            CCTK_POINTER_TO_CONST const interp_coords_,
+                                            CCTK_INT const N_input_arrays,
+                                            CCTK_INT const * const input_array_indices,
+                                            CCTK_INT const N_output_arrays,
+                                            CCTK_POINTER const output_arrays_)
+  {
+    // Check input values and convert types
+    cGH const * const cctkGH = static_cast<cGH const *> (cctkGH_);
+    assert (cctkGH);
+    
+    assert (N_dims == dim);
+    
+    assert (N_interp_points >= 0);
+    
+    assert (interp_maps);
+    CCTK_REAL const * const * const interp_coords =
+      static_cast<CCTK_REAL const * const *> (interp_coords_);
+    assert (interp_coords);
+    for (int d=0; d<dim; ++d) {
+      assert (N_interp_points==0 or interp_coords[d]);
+    }
+    
+    assert (N_input_arrays >= 0);
+    
+    assert (input_array_indices);
+    for (int n=0; n<N_input_arrays; ++n) {
+      assert (input_array_indices[n]>=0 and
+              input_array_indices[n]<CCTK_NumVars());
+    }
+    
+    assert (N_output_arrays >= 0);
+    assert (N_output_arrays == N_input_arrays);
+    
+    CCTK_REAL * const * const output_arrays =
+      static_cast<CCTK_REAL * const *> (output_arrays_);
+    assert (output_arrays);
+    for (int n=0; n<N_output_arrays; ++n) {
+      assert (N_output_arrays==0 or output_arrays[n]);
+    }
+    
+    // Set up interpolation
+    fasterp_llocs_t locations (N_interp_points);
+    for (int i=0; i<N_interp_points; ++i) {
+      locations.maps.AT(i) = interp_maps[i];
+    }
+    for (int d=0; d<dim; ++d) {
+      for (int i=0; i<N_interp_points; ++i) {
+        locations.coords[d].AT(i) = interp_coords[d][i];
+      }
+    }
+    fasterp_setup_t const setup (cctkGH, locations, order);
+    
+    // Interpolate
+    vector<int> varinds (N_input_arrays);
+    for (int n=0; n<N_input_arrays; ++n) {
+      varinds.AT(n) = input_array_indices[n];
+    }
+    vector<CCTK_REAL *> values (N_input_arrays);
+    for (int n=0; n<N_input_arrays; ++n) {
+      values.AT(n) = output_arrays[n];
+    }
+    setup.interpolate (cctkGH, varinds, values);
+    
+    // Done
+    return 0;
+  }
+  
+  
+  
+  extern "C"
+  CCTK_POINTER
+  CarpetInterp2_Interp2GridArraysSetup
+  (CCTK_POINTER_TO_CONST const cctkGH_,
+   CCTK_INT const N_dims,
+   CCTK_INT const order,
+   CCTK_INT const N_interp_points,
+   CCTK_POINTER_TO_CONST const interp_coords_)
+  {
+    // Check input values and convert types
+    cGH const * const cctkGH = static_cast<cGH const *> (cctkGH_);
+    assert (cctkGH);
+    
+    assert (N_dims == dim);
+    
+    assert (N_interp_points >= 0);
+    
+    CCTK_REAL const * const * const interp_coords =
+      static_cast<CCTK_REAL const * const *> (interp_coords_);
+    assert (interp_coords);
+    for (int d=0; d<dim; ++d) {
+      assert (N_interp_points==0 or interp_coords[d]);
+    }
+    
+    // Set up interpolation
+    fasterp_glocs_t locations (N_interp_points);
+    for (int d=0; d<dim; ++d) {
+      for (int i=0; i<N_interp_points; ++i) {
+        locations.coords[d].AT(i) = interp_coords[d][i];
+      }
+    }
+    fasterp_setup_t * const setup =
+      new fasterp_setup_t (cctkGH, locations, order);
+    
+    // Done
+    return setup;
+  }
+  
+  extern "C"
+  CCTK_INT
+  CarpetInterp2_Interp2GridArrays
+  (CCTK_POINTER_TO_CONST const cctkGH_,
+   CCTK_POINTER_TO_CONST const setup_,
+   CCTK_INT const N_input_arrays,
+   CCTK_INT const * const input_array_indices,
+   CCTK_INT const N_output_arrays,
+   CCTK_POINTER const output_arrays_)
+  {
+    // Check input values and convert types
+    cGH const * const cctkGH = static_cast<cGH const *> (cctkGH_);
+    assert (cctkGH);
+    
+    fasterp_setup_t const * const setup =
+      static_cast<fasterp_setup_t const *> (setup_);
+    assert (setup);
+    
+    assert (N_input_arrays >= 0);
+    
+    assert (input_array_indices);
+    for (int n=0; n<N_input_arrays; ++n) {
+      assert (input_array_indices[n]>=0 and
+              input_array_indices[n]<CCTK_NumVars());
+    }
+    
+    assert (N_output_arrays >= 0);
+    assert (N_output_arrays == N_input_arrays);
+    
+    CCTK_REAL * const * const output_arrays =
+      static_cast<CCTK_REAL * const *> (output_arrays_);
+    assert (output_arrays);
+    for (int n=0; n<N_output_arrays; ++n) {
+      assert (N_output_arrays==0 or output_arrays[n]);
+    }
+    
+    // Interpolate
+    vector<int> varinds (N_input_arrays);
+    for (int n=0; n<N_input_arrays; ++n) {
+      varinds.AT(n) = input_array_indices[n];
+    }
+    vector<CCTK_REAL *> values (N_input_arrays);
+    for (int n=0; n<N_input_arrays; ++n) {
+      values.AT(n) = output_arrays[n];
+    }
+    setup->interpolate (cctkGH, varinds, values);
+    
+    // Done
+    return 0;
+  }
+  
+  extern "C"
+  CCTK_INT
+  CarpetInterp2_Interp2GridArraysFree
+  (CCTK_POINTER_TO_CONST const cctkGH_,
+   CCTK_POINTER          const setup_)
+  {
+    // Check input values and convert types
+    cGH const * const cctkGH = static_cast<cGH const *> (cctkGH_);
+    assert (cctkGH);
+    
+    fasterp_setup_t const * const setup =
+      static_cast<fasterp_setup_t const *> (setup_);
+    assert (setup);
+    
+    delete setup;
+    
+    // Done
+    return 0;
+  }
+  
 } // namespace CarpetInterp2