CarpetLib: Support prolongating faces

Add infrastructure to prolongate faces of grid functions. This may be used during refluxing to apply part of the correction to the fine grid.

5 files changed, 119 insertions, 1 deletions

diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 1892336c1..799bae113 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -66,6 +66,31 @@ init_fast_ref_refl_sendrecv ()
   fast_ref_refl_sendrecv[2][1] = & fast_dboxes::fast_ref_refl_sendrecv_2_1;
 }
 
+vect<vect<dh::srpvect dh::fast_dboxes::*,2>,dim>
+dh::fast_dboxes::
+fast_ref_refl_prol_sendrecv;
+
+void
+dh::fast_dboxes::
+init_fast_ref_refl_prol_sendrecv ()
+{
+  static bool initialised = false;
+  if (initialised) return;
+  initialised = true;
+  fast_ref_refl_prol_sendrecv[0][0] =
+    & fast_dboxes::fast_ref_refl_prol_sendrecv_0_0;
+  fast_ref_refl_prol_sendrecv[0][1] =
+    & fast_dboxes::fast_ref_refl_prol_sendrecv_0_1;
+  fast_ref_refl_prol_sendrecv[1][0] =
+    & fast_dboxes::fast_ref_refl_prol_sendrecv_1_0;
+  fast_ref_refl_prol_sendrecv[1][1] =
+    & fast_dboxes::fast_ref_refl_prol_sendrecv_1_1;
+  fast_ref_refl_prol_sendrecv[2][0] =
+    & fast_dboxes::fast_ref_refl_prol_sendrecv_2_0;
+  fast_ref_refl_prol_sendrecv[2][1] =
+    & fast_dboxes::fast_ref_refl_prol_sendrecv_2_1;
+}
+
 
 
 // Constructors
@@ -82,6 +107,7 @@ dh (gh & h_,
     prolongation_orders_space(prolongation_orders_space_)
 {
   fast_dboxes::init_fast_ref_refl_sendrecv();
+  fast_dboxes::init_fast_ref_refl_prol_sendrecv();
   size_t const maxreflevels = h.reffacts.size();
   assert (ghost_widths.size() >= maxreflevels);
   assert (buffer_widths.size() >= maxreflevels);
@@ -1263,6 +1289,8 @@ regrid (bool const do_init)
                 
                 srpvect fast_dboxes::* const fast_ref_refl_sendrecv =
                   fast_dboxes::fast_ref_refl_sendrecv[dir][face];
+                srpvect fast_dboxes::* const fast_ref_refl_prol_sendrecv =
+                  fast_dboxes::fast_ref_refl_prol_sendrecv[dir][face];
                 
                 // Refluxing must fill all coarse refluxing boundary
                 // points, and may use all fine points
@@ -1309,6 +1337,22 @@ regrid (bool const do_init)
                       (fast_level_otherproc.*fast_ref_refl_sendrecv).
                         push_back (preg);
                     }
+                    
+                    // reflux-prolongation
+                    sendrecv_pseudoregion_t const preg_prol
+                      (shifted_recv, oc, shifted_send, c);
+                    if (verbose) {
+                      cout << "REF ref_refL_prol ml=" << ml << " rl=" << rl << " olc=" << olc << " c=" << c << " oc=" << oc << " dir=" << dir << " face=" << face << "\n"
+                           << "   preg=" << preg_prol << "\n";
+                    }
+                    (fast_level.*fast_ref_refl_prol_sendrecv).push_back
+                      (preg_prol);
+                    if (not on_this_proc (orl, oc)) {
+                      fast_dboxes & fast_level_otherproc =
+                        fast_level_otherprocs.AT(this_proc(orl, oc));
+                      (fast_level_otherproc.*fast_ref_refl_prol_sendrecv).
+                        push_back (preg_prol);
+                    }
                   }
                   
                   needrecv -= ovlp;
@@ -1881,7 +1925,6 @@ regrid (bool const do_init)
             for (int face = 0; face < 2; ++ face) {
               srpvect fast_dboxes::* const fast_ref_refl_sendrecv =
                 fast_dboxes::fast_ref_refl_sendrecv[dir][face];
-              
               broadcast_schedule (fast_level_otherprocs, fast_olevel,
                                   fast_ref_refl_sendrecv);
             }
@@ -1889,6 +1932,19 @@ regrid (bool const do_init)
           timer_bcast_comm_ref_refl.stop();
         }
         
+        static Carpet::Timer timer_bcast_comm_ref_refl_prol
+          ("ref_refl_prol");
+        timer_bcast_comm_ref_refl_prol.start();
+        for (int dir = 0; dir < dim; ++ dir) {
+          for (int face = 0; face < 2; ++ face) {
+            srpvect fast_dboxes::* const fast_ref_refl_prol_sendrecv =
+              fast_dboxes::fast_ref_refl_prol_sendrecv[dir][face];
+            broadcast_schedule (fast_level_otherprocs, fast_level,
+                                fast_ref_refl_prol_sendrecv);
+          }
+        }
+        timer_bcast_comm_ref_refl_prol.stop();
+        
         // TODO: Maybe broadcast old2new schedule only if do_init is
         // set
         static Carpet::Timer timer_bcast_comm_old2new_sync
@@ -2283,6 +2339,12 @@ mpi_datatype (dh::fast_dboxes const &)
       ENTRY (dh::srpvect, fast_ref_refl_sendrecv_1_1),
       ENTRY (dh::srpvect, fast_ref_refl_sendrecv_2_0),
       ENTRY (dh::srpvect, fast_ref_refl_sendrecv_2_1),
+      ENTRY (dh::srpvect, fast_ref_refl_prol_sendrecv_0_0),
+      ENTRY (dh::srpvect, fast_ref_refl_prol_sendrecv_0_1),
+      ENTRY (dh::srpvect, fast_ref_refl_prol_sendrecv_1_0),
+      ENTRY (dh::srpvect, fast_ref_refl_prol_sendrecv_1_1),
+      ENTRY (dh::srpvect, fast_ref_refl_prol_sendrecv_2_0),
+      ENTRY (dh::srpvect, fast_ref_refl_prol_sendrecv_2_1),
       {1, sizeof s, MPI_UB, "MPI_UB", "MPI_UB"}
     };
 #undef ENTRY
@@ -2405,6 +2467,12 @@ memory ()
     memoryof (fast_ref_refl_sendrecv_1_1) +
     memoryof (fast_ref_refl_sendrecv_2_0) +
     memoryof (fast_ref_refl_sendrecv_2_1) +
+    memoryof (fast_ref_refl_prol_sendrecv_0_0) +
+    memoryof (fast_ref_refl_prol_sendrecv_0_1) +
+    memoryof (fast_ref_refl_prol_sendrecv_1_0) +
+    memoryof (fast_ref_refl_prol_sendrecv_1_1) +
+    memoryof (fast_ref_refl_prol_sendrecv_2_0) +
+    memoryof (fast_ref_refl_prol_sendrecv_2_1) +
     memoryof (do_init) +
     memoryof (fast_old2new_sync_sendrecv) +
     memoryof (fast_old2new_ref_prol_sendrecv);
@@ -2491,6 +2559,7 @@ input (istream & is)
     consume (is, "fine_boundary:");
     is >> fine_boundary;
     skipws (is);
+    // TODO: read boundary sizes and boundary offsets
     consume (is, "}");
   } catch (input_error & err) {
     cout << "Input error while reading a dh::local_dboxes" << endl;
@@ -2710,6 +2779,12 @@ output (ostream & os)
      << "   fast_ref_refl_sendrecv_1_1: " << fast_ref_refl_sendrecv_1_1 << eol
      << "   fast_ref_refl_sendrecv_2_0: " << fast_ref_refl_sendrecv_2_0 << eol
      << "   fast_ref_refl_sendrecv_2_1: " << fast_ref_refl_sendrecv_2_1 << eol
+     << "   fast_ref_refl_prol_sendrecv_0_0: " << fast_ref_refl_prol_sendrecv_0_0 << eol
+     << "   fast_ref_refl_prol_sendrecv_0_1: " << fast_ref_refl_prol_sendrecv_0_1 << eol
+     << "   fast_ref_refl_prol_sendrecv_1_0: " << fast_ref_refl_prol_sendrecv_1_0 << eol
+     << "   fast_ref_refl_prol_sendrecv_1_1: " << fast_ref_refl_prol_sendrecv_1_1 << eol
+     << "   fast_ref_refl_prol_sendrecv_2_0: " << fast_ref_refl_prol_sendrecv_2_0 << eol
+     << "   fast_ref_refl_prol_sendrecv_2_1: " << fast_ref_refl_prol_sendrecv_2_1 << eol
      << "   do_init: " << do_init << eol
      << "   fast_old2new_sync_sendrecv: " << fast_old2new_sync_sendrecv << eol
      << "   fast_old2new_ref_prol_sendrecv: " << fast_old2new_ref_prol_sendrecv << eol
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index 2ec339541..e0c78fd2a 100644
--- a/Carpet/CarpetLib/src/dh.hh
+++ b/Carpet/CarpetLib/src/dh.hh
@@ -166,6 +166,18 @@ public:
     static
     void init_fast_ref_refl_sendrecv ();
     
+    // vect<vect<srpvect,2>,dim> fast_ref_refl_prol_sendrecv;
+    srpvect fast_ref_refl_prol_sendrecv_0_0;
+    srpvect fast_ref_refl_prol_sendrecv_0_1;
+    srpvect fast_ref_refl_prol_sendrecv_1_0;
+    srpvect fast_ref_refl_prol_sendrecv_1_1;
+    srpvect fast_ref_refl_prol_sendrecv_2_0;
+    srpvect fast_ref_refl_prol_sendrecv_2_1;
+    static
+    vect<vect<srpvect fast_dboxes::*,2>,dim> fast_ref_refl_prol_sendrecv;
+    static
+    void init_fast_ref_refl_prol_sendrecv ();
+    
     // Regridding schedule:
     
     bool do_init;               // the srpvects below are only defined
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index dbabfe1c7..17f9d9406 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -516,6 +516,32 @@ ref_reflux_all (comm_state & state,
 
 
 
+// Reflux-prolongate a refinement level
+void
+ggf::
+ref_reflux_prolongate_all (comm_state & state,
+                           int const tl, int const rl, int const ml,
+                           int const dir, int const face)
+{
+  static Timer timer ("ref_reflux_prolongate_all");
+  timer.start ();
+  // Require same times
+  static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
+  assert (abs(t.get_time(ml,rl,tl) - t.get_time(ml,rl-1,tl)) <=
+          1.0e-8 * (1.0 + abs(t.get_time(ml,rl,tl))));
+  islab slabinfo;
+  slabinfo.is_centered = 1 - ivect::dir(dir);
+  transfer_from_all (state,
+                     tl,rl,ml,
+                     dh::fast_dboxes::fast_ref_refl_prol_sendrecv[dir][face],
+                     tl,rl-1,ml,
+                     false, false,
+                     &slabinfo);
+  timer.stop (0);
+}
+
+
+
 // Transfer regions of all components
 void
 ggf::
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 93a5a78bd..2e8b9f4ad 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -144,6 +144,10 @@ public:
   // Reflux a refinement level
   void ref_reflux_all (comm_state& state,
                        int tl, int rl, int ml, int dir, int face);
+
+  // Reflux a refinement level
+  void ref_reflux_prolongate_all (comm_state& state,
+                                  int tl, int rl, int ml, int dir, int face);
   
   
   
diff --git a/Carpet/CarpetLib/src/vect.cc b/Carpet/CarpetLib/src/vect.cc
index 661fd0190..9fcaf0e4f 100644
--- a/Carpet/CarpetLib/src/vect.cc
+++ b/Carpet/CarpetLib/src/vect.cc
@@ -121,6 +121,7 @@ template void vect<vect<int,2>,dim>::output (ostream& os) const;
 template void vect<vect<bool,dim>,2>::output (ostream& os) const;
 template void vect<vect<int,dim>,2>::output (ostream& os) const;
 template void vect<vect<CCTK_REAL,dim>,2>::output (ostream& os) const;
+template void vect<vect<vector<int>,2>,dim>::output (ostream& os) const;