CarpetLib: Add new transport operator type op_copy

Add a new transport operator type op_copy.  This "prolongation"
operator does not interpolate in time, but rather copies the from the
newest time level instead.  Such grid functions need only one time
level.

This is intended for prolongating or restricing grid functions like
the ADM constraints; if done properly; they will have the same values
on the coarse and fine grids.  (However, this does not work for the
ADM constraints, because such grid functions still need to be set in
EVOL, not in ANALYSIS.)

darcs-hash:20050607160833-891bb-cfd1c7630f8996606328d7c7e9fe326561106aba.gz

3 files changed, 98 insertions, 14 deletions

diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index c8dca0d70..1d7de42e0 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -992,6 +992,7 @@ void data<T>
   
   switch (transport_operator) {
       
+  case op_copy:
   case op_Lagrange:
   case op_TVD:
   case op_ENO:
@@ -1036,6 +1037,60 @@ void data<T>
                     box, sext, dext);
   switch (transport_operator) {
     
+  case op_copy:
+    assert (times.size() == 1);
+    assert (srcs.size()>=1);
+    switch (order_space) {
+    case 0:
+    case 1:
+      if (all (sext.stride() == dext.stride() * 2)) {
+        CCTK_FNAME(prolongate_3d_real8_rf2)
+          ((const CCTK_REAL8*)srcs[0]->storage(),
+           srcshp[0], srcshp[1], srcshp[2],
+           (CCTK_REAL8*)storage(),
+           dstshp[0], dstshp[1], dstshp[2],
+           srcbbox, dstbbox, regbbox);
+      } else {
+        CCTK_FNAME(prolongate_3d_real8)
+          ((const CCTK_REAL8*)srcs[0]->storage(),
+           srcshp[0], srcshp[1], srcshp[2],
+           (CCTK_REAL8*)storage(),
+           dstshp[0], dstshp[1], dstshp[2],
+           srcbbox, dstbbox, regbbox);
+      }
+      break;
+    case 2:
+    case 3:
+      if (all (sext.stride() == dext.stride() * 2)) {
+        CCTK_FNAME(prolongate_3d_real8_o3_rf2)
+          ((const CCTK_REAL8*)srcs[0]->storage(),
+           srcshp[0], srcshp[1], srcshp[2],
+           (CCTK_REAL8*)storage(),
+           dstshp[0], dstshp[1], dstshp[2],
+           srcbbox, dstbbox, regbbox);
+      } else {
+        CCTK_FNAME(prolongate_3d_real8_o3)
+          ((const CCTK_REAL8*)srcs[0]->storage(),
+           srcshp[0], srcshp[1], srcshp[2],
+           (CCTK_REAL8*)storage(),
+           dstshp[0], dstshp[1], dstshp[2],
+           srcbbox, dstbbox, regbbox);
+      }
+      break;
+    case 4:
+    case 5:
+      CCTK_FNAME(prolongate_3d_real8_o5)
+        ((const CCTK_REAL8*)srcs[0]->storage(),
+         srcshp[0], srcshp[1], srcshp[2],
+         (CCTK_REAL8*)storage(),
+         dstshp[0], dstshp[1], dstshp[2],
+         srcbbox, dstbbox, regbbox);
+      break;
+    default:
+      assert (0);
+    }
+    break;
+    
   case op_Lagrange:
     switch (order_time) {
       
@@ -1396,12 +1451,22 @@ void data<CCTK_REAL8>
     min_time = min(min_time, times[tl]);
     max_time = max(max_time, times[tl]);
   }
-  if (time < min_time - eps || time > max_time + eps) {
-    ostringstream buf;
-    buf << "Internal error: extrapolation in time."
-        << "  time=" << time
-        << "  times=" << times;
-    CCTK_WARN (0, buf.str().c_str());
+  if (transport_operator != op_copy) {
+    if (time < min_time - eps || time > max_time + eps) {
+      ostringstream buf;
+      buf << "Internal error: extrapolation in time."
+          << "  time=" << time
+          << "  times=" << times;
+      CCTK_WARN (0, buf.str().c_str());
+    }
+  } else {
+    if (time > max_time + eps) {
+      ostringstream buf;
+      buf << "Internal error: extrapolation into the future."
+          << "  time=" << time
+          << "  times=" << times;
+      CCTK_WARN (0, buf.str().c_str());
+    }
   }
 }
 
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 1d7f7fccb..6989ba28d 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -477,9 +477,11 @@ void ggf::intercat (comm_state& state,
     {
       // (use the send boxes for communication)
       // interpolate the content
+      int const pos = d.prolongation_order_space;
+      int const pot
+        = transport_operator == op_copy ? 0 : prolongation_order_time;
       storage.at(ml1).at(rl1).at(c1).at(tl1)->interpolate_from
-      	(state, gsrcs, times, *r, time,
-	 d.prolongation_order_space, prolongation_order_time);
+      	(state, gsrcs, times, *r, time, pos, pot);
     }
   }
 }
@@ -513,13 +515,19 @@ void ggf::ref_bnd_prolongate (comm_state& state,
   assert (rl>=1);
   if (transport_operator == op_none) return;
   vector<int> tl2s;
-  // Interpolation in time
-  assert (timelevels() >= prolongation_order_time+1);
-  tl2s.resize(prolongation_order_time+1);
-  for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = i;
+  if (transport_operator != op_copy) {
+    // Interpolation in time
+    assert (timelevels() >= prolongation_order_time+1);
+    tl2s.resize(prolongation_order_time+1);
+    for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = i;
+  } else {
+    assert (timelevels() >= 1);
+    tl2s.resize(1);
+    tl2s.at(0) = 0;
+  }
   intercat (state,
             tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_bnd_coarse,
-	    tl2s,rl-1,  ml, time);
+            tl2s,rl-1,  ml, time);
 }
 
 // Restrict a multigrid level
diff --git a/Carpet/CarpetLib/src/operators.hh b/Carpet/CarpetLib/src/operators.hh
index e884415fb..d225feef5 100644
--- a/Carpet/CarpetLib/src/operators.hh
+++ b/Carpet/CarpetLib/src/operators.hh
@@ -1,6 +1,17 @@
 #ifndef OPERATORS_HH
 #define OPERATORS_HH
 
-enum operator_type { op_error, op_none, op_Lagrange, op_TVD, op_ENO };
+// Transport (i.e., prolongation and restriction) operator types
+
+enum operator_type
+{
+  op_error,                     // illegal operator type
+  op_none,                      // do not transport
+  op_copy,                      // use simple copying for prolongation
+                                // (needs only one time level)
+  op_Lagrange,                  // Lagrange interpolation (standard)
+  op_TVD,                       // use TVD stencils (for hydro)
+  op_ENO                        // use ENO stencils (for hydro)
+};
 
 #endif // OPERATORS_HH