Change SetInterpolators and so on, before changes to slicedata

src/Evolution/DgSubcell/Actions/Initialize.hpp

@@ -100,7 +100,8 @@ struct SetSubcellGrid {
       subcell::Tags::ReconstructionOrder<Dim>,
       evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<Dim>,
       evolution::dg::subcell::Tags::InterpolatorsFromDgToNeighborFd<Dim>,
-      evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd<Dim>>;
+      evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd<Dim>,
+      evolution::dg::subcell::Tags::ProblemoChest<Dim>>;
   using compute_tags =
       tmpl::list<Tags::MeshCompute<Dim>, Tags::LogicalCoordinatesCompute<Dim>,
                  ::domain::Tags::MappedCoordinates<

src/Evolution/DgSubcell/Actions/ReconstructionCommunication.hpp

@@ -57,6 +57,7 @@
 #include "Parallel/GlobalCache.hpp"
 #include "Time/TimeStepId.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
+#include "Utilities/ErrorHandling/CaptureForError.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/Literals.hpp"
 #include "Utilities/MakeArray.hpp"
@@ -169,6 +170,23 @@ struct SendDataForReconstruction {
               static_cast<std::ptrdiff_t>(volume_data_to_slice.size() -
                                           cell_centered_flux.value().size())));
     }
+    // for (const auto [k]) {}
+    //   auto yolo =
+    //       db::get<evolution::dg::subcell::Tags::ProblemoChest<Dim>>(box);
+    // if (!yolo.empty() &&
+    //     element.id() == ElementId<Dim>{"[B0,(L1I0,L1I0,L1I1)]"}) {
+    //   // CAPTURE_FOR_ERROR(yolo);
+    //   CAPTURE_FOR_ERROR(
+    //       yolo.at(DirectionalId<Dim>{Direction<Dim>{0, Side::Lower},
+    //                             ElementId<Dim>{"[B5,(L1I1,L1I1,L1I1)]"}})
+    //           .value()
+    //           .problemo_dim);
+    //   CAPTURE_FOR_ERROR(yolo);
+    //   ERROR("wow");
+    // }
+
+    // modify slice data to not do anyhting for non empty interpolator
+    // not have the direction as key
     const DirectionMap<Dim, DataVector> all_sliced_data = slice_data(
         volume_data_to_slice, subcell_mesh.extents(), ghost_zone_size,
         element.internal_boundaries(), 0,
@@ -276,6 +294,22 @@ struct SendDataForReconstruction {
   }
 };
 
+// send data (modified sliced data where I don't send anyhting for problematic
+// direction) new action (receive data -> inbox -> checks to make sure it got
+// everyhitng all direction except problematic -> interpolation using the ghost
+// data plus yourself (volume subcell plus the neighbordat)) receive data
+// (untouched)
+
+/*!
+ * \brief Receive the subcell data from our neighbor for all non-problematic
+ * directions. Interpolate and send data for all problematic directions.
+ *
+ */
+template <size_t Dim, typename GhostDataMutator, bool LocalTimeStepping>
+struct ReceivedAndSendDataForReconstruction {
+  using inbox_tags = tmpl::list<
+      evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<Dim>>;
+};
 /*!
  * \brief Receive the subcell data from our neighbor, and accumulate the data
  * from the relaxed discrete maximum principle troubled-cell indicator.

src/Evolution/DgSubcell/SetInterpolators.hpp

@@ -26,6 +26,7 @@
 #include "Evolution/DgSubcell/Tags/Reconstructor.hpp"
 #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "Utilities/ErrorHandling/CaptureForError.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/TMPL.hpp"
@@ -54,7 +55,8 @@ struct SetInterpolators {
   using return_tags = tmpl::list<
       evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<Dim>,
       evolution::dg::subcell::Tags::InterpolatorsFromDgToNeighborFd<Dim>,
-      evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd<Dim>>;
+      evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd<Dim>,
+      evolution::dg::subcell::Tags::ProblemoChest<Dim>>;
   using argument_tags =
       tmpl::list<::domain::Tags::Element<Dim>, ::domain::Tags::Domain<Dim>,
                  domain::Tags::Mesh<Dim>, domain::Tags::Mesh<Dim>,
@@ -74,6 +76,9 @@ struct SetInterpolators {
       const gsl::not_null<
           DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>*>
           interpolators_neighbor_dg_to_fd_ptr,
+      const gsl::not_null<DirectionalIdMap<
+          Dim, std::optional<interpolators_detail::ProblemoBook>>*>
+          proeblemo_book_ptr,
       const Element<Dim>& element, const Domain<Dim>& domain,
       const Mesh<Dim>& my_dg_mesh,
       // Needs to be updated to support non-uniform h/p-refinement
@@ -167,11 +172,69 @@ struct SetInterpolators {
             neighbor_logical_ghost_zone_coords = get_logical_coords(
                 element_map, neighbor_grid_ghost_zone_coords);
 
+        // Here I want to set up the conditional that checks the
+        // logical coordinates (in the element) of the neighbors ghost
+        // if they exceed past the extremal grid points of the element
+        // aka leading to extrapolation.
+        bool problemo = false;
+        size_t max = std::numeric_limits<size_t>::max();
+        size_t problematic_dim = max;
+
+        for (size_t i = 0; i < neighbor_logical_ghost_zone_coords[0].size();
+             i++) {
+          for (size_t d = 0; d < Dim; d++) {
+            double ext = 1 - (1 / my_fd_mesh.extents(d));
+            if (std::abs(neighbor_logical_ghost_zone_coords.get(d)[i]) > ext) {
+              problemo = true;
+              // until we are confident, have this assert and comment the break
+              // below
+              ASSERT(problematic_dim == max || problematic_dim == d,
+                     "problem exists in more than one dimension");
+              problematic_dim = d;
+              // break;
+            }
+          }
+        }
+        if (problemo) {
+          auto new_neighbor_logical_ghost_zone_coords =
+              neighbor_logical_ghost_zone_coords;
+          auto new_basis = make_array<Dim>(my_fd_mesh.basis(0));
+          auto new_extents = make_array<Dim>(my_fd_mesh.extents(0));
+          auto new_quads = make_array<Dim>(my_fd_mesh.quadrature(0));
+          const double mf =
+              (my_fd_mesh.extents(problematic_dim)) /
+              (my_fd_mesh.extents(problematic_dim) + number_of_ghost_zones);
+          const double af =
+              (number_of_ghost_zones) /
+              (my_fd_mesh.extents(problematic_dim) + number_of_ghost_zones);
+          for (size_t d = 0; d < Dim; d++) {
+            gsl::at(new_basis, d) = my_fd_mesh.basis(d);
+            gsl::at(new_quads, d) = my_fd_mesh.quadrature(d);
+            if (d == problematic_dim) {
+              gsl::at(new_extents, d) =
+                  my_fd_mesh.extents(d) + number_of_ghost_zones;
+              for (size_t i = 0;
+                   i < new_neighbor_logical_ghost_zone_coords[0].size(); i++) {
+                new_neighbor_logical_ghost_zone_coords.get(d)[i] *= mf;
+                new_neighbor_logical_ghost_zone_coords.get(d)[i] -= af;
+              }
+            } else {
+              gsl::at(new_extents, d) = my_fd_mesh.extents(d);
+            }
+          }
+          const Mesh<Dim> new_mesh{new_extents, new_basis, new_quads};
+          (*interpolators_fd_to_neighbor_fd_ptr)[DirectionalId<Dim>{
+              direction, neighbor_id}] = intrp::Irregular<Dim>{
+              new_mesh, new_neighbor_logical_ghost_zone_coords};
+          (*proeblemo_book_ptr)[DirectionalId<Dim>{direction, neighbor_id}] =
+              interpolators_detail::ProblemoBook{problematic_dim};
+        } else {
+          (*interpolators_fd_to_neighbor_fd_ptr)[DirectionalId<Dim>{
+              direction, neighbor_id}] = intrp::Irregular<Dim>{
+              my_fd_mesh, neighbor_logical_ghost_zone_coords};
+        }
         // Set up interpolators for our local element to our neighbor's
         // ghost zones.
-        (*interpolators_fd_to_neighbor_fd_ptr)[DirectionalId<Dim>{
-            direction, neighbor_id}] = intrp::Irregular<Dim>{
-            my_fd_mesh, neighbor_logical_ghost_zone_coords};
         (*interpolators_dg_to_neighbor_fd_ptr)[DirectionalId<Dim>{
             direction, neighbor_id}] = intrp::Irregular<Dim>{
             my_dg_mesh, neighbor_logical_ghost_zone_coords};

src/Evolution/DgSubcell/Tags/Interpolators.hpp

@@ -16,12 +16,24 @@ template <size_t Dim>
 class Direction;
 template <size_t Dim>
 class ElementId;
+namespace PUP {
+class er;
+}  // namespace PUP
 namespace intrp {
 template <size_t Dim>
 class Irregular;
 }  // namespace intrp
 /// \endcond
 
+namespace interpolators_detail {
+struct ProblemoBook {
+  size_t problemo_dim;
+  /// Serialization for Charm++
+  // NOLINTNEXTLINE(google-runtime-references)
+  void pup(PUP::er& p) { p | problemo_dim; };
+};
+}  // namespace interpolators_detail
+
 namespace evolution::dg::subcell::Tags {
 /*!
  * \brief An `intrp::Irregular` from our FD grid to our neighbors' FD grid.
@@ -55,4 +67,15 @@ template <size_t Dim>
 struct InterpolatorsFromNeighborDgToFd : db::SimpleTag {
   using type = DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>;
 };
+
+/*!
+ * \brief This will be a book-keeping that records whether there the neighbor's
+ * ghost goes out of bound and the direction in which it does.
+ *
+ */
+template <size_t Dim>
+struct ProblemoChest : db::SimpleTag {
+  using type =
+      DirectionalIdMap<Dim, std::optional<interpolators_detail::ProblemoBook>>;
+};
 }  // namespace evolution::dg::subcell::Tags