Prolongation operator.

Arrange flags to split physics conserved vs. code conserved vars.
(Better flag organization forthcoming which will hide some
silliness here.)
In particular, we want the code to conserve the face B,
but the cell-centered version is still "conserved" in the
physical sense.

Add the prolongation operator from Olivares et al. 2019 in
generality, but currently hardcoded to alpha values matching
the operator in Toth & Roe.

Boundaries between refinement levels are still a problem,
but overall this is pretty close for a day of fiddling.

kharma/b_cd/b_cd.cpp

@@ -62,20 +62,21 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
     // B field as usual
     // TODO allow for implicit B here
     Metadata m = Metadata({Metadata::Real, Metadata::Cell, Metadata::Independent, Metadata::FillGhost,
-                 Metadata::Restart, Metadata::Conserved, Metadata::WithFluxes, Metadata::Vector}, s_vector);
+                 Metadata::Restart, Metadata::GetUserFlag("GRConserved"), Metadata::Conserved,
+                 Metadata::WithFluxes, Metadata::Vector}, s_vector);
     pkg->AddField("cons.B", m);
     m = Metadata({Metadata::Real, Metadata::Cell, Metadata::Derived,
-                  Metadata::Restart, Metadata::GetUserFlag("Primitive"), Metadata::Vector}, s_vector);
+                  Metadata::Restart, Metadata::GetUserFlag("GRPrimitive"), Metadata::Vector}, s_vector);
     pkg->AddField("prims.B", m);
 
     // Constraint damping scalar field psi.  Prim and cons forms correspond to B field forms,
     // i.e. differ by a factor of gdet.  This is apparently marginally more stable in some
     // circumstances.
     m = Metadata({Metadata::Real, Metadata::Cell, Metadata::Independent, Metadata::FillGhost,
-                  Metadata::Restart, Metadata::Conserved, Metadata::WithFluxes});
+                  Metadata::Restart, Metadata::GetUserFlag("GRConserved"), Metadata::Conserved, Metadata::WithFluxes});
     pkg->AddField("cons.psi_cd", m);
     m = Metadata({Metadata::Real, Metadata::Cell, Metadata::Derived,
-                  Metadata::Restart, Metadata::GetUserFlag("Primitive")});
+                  Metadata::Restart, Metadata::GetUserFlag("GRPrimitive")});
     pkg->AddField("prims.psi_cd", m);
 
     // We only update the divB field for output

kharma/b_cd/b_cd.hpp

@@ -47,7 +47,7 @@ using namespace parthenon;
  *
  * This requires only the values at cell centers, and preserves a cell-centered divergence representation
  * 
- * This implementation includes conversion from "primitive" to "conserved" B and back,
+ * This implementation includes conversion from "GRPrimitive" to "conserved" B and back,
  * i.e. between field strength and flux via multiplying by gdet.
  */
 namespace B_CD {

kharma/b_cleanup/b_cleanup.cpp

@@ -146,10 +146,10 @@ std::shared_ptr<KHARMAPackage> B_Cleanup::Initialize(ParameterInput *pin, std::s
         MetadataFlag areWeImplicit = (implicit_b) ? Metadata::GetUserFlag("Implicit")
                                                     : Metadata::GetUserFlag("Explicit");
 
-        // Flags for B fields.  "Primitive" form is field, "conserved" is flux
-        std::vector<MetadataFlag> flags_prim = {Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::GetUserFlag("Primitive"),
+        // Flags for B fields.  "GRPrimitive" form is field, "conserved" is flux
+        std::vector<MetadataFlag> flags_prim = {Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::GetUserFlag("GRPrimitive"),
                                                 Metadata::Restart, Metadata::GetUserFlag("MHD"), areWeImplicit, Metadata::Vector};
-        std::vector<MetadataFlag> flags_cons = {Metadata::Real, Metadata::Cell, Metadata::Independent, Metadata::Conserved,
+        std::vector<MetadataFlag> flags_cons = {Metadata::Real, Metadata::Cell, Metadata::Independent, Metadata::GetUserFlag("GRConserved"), Metadata::Conserved,
                                                 Metadata::WithFluxes, Metadata::FillGhost, Metadata::GetUserFlag("MHD"), areWeImplicit, Metadata::Vector};
 
         auto m = Metadata(flags_prim, s_vector);

kharma/b_ct/b_ct.cpp

@@ -41,6 +41,7 @@
 #include "kharma_driver.hpp"
 
 #include <parthenon/parthenon.hpp>
+#include <prolong_restrict/pr_ops.hpp>
 
 using namespace parthenon;
 
@@ -77,22 +78,23 @@ std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared
     // TODO maybe one day implicit?
 
     // Flags for B fields on faces.
-    // We don't mark these as "Primitive" and "Conserved" else they'd be bundled
+    // We don't mark these as "GRPrimitive" and "GRConserved" else they'd be bundled
     // with all the cell vars in a bunch of places we don't want
     std::vector<MetadataFlag> flags_prim_f = {Metadata::Real, Metadata::Face, Metadata::Derived,
                                             Metadata::GetUserFlag("Explicit")};
-    std::vector<MetadataFlag> flags_cons_f = {Metadata::Real, Metadata::Face, Metadata::Independent,
+    std::vector<MetadataFlag> flags_cons_f = {Metadata::Real, Metadata::Face, Metadata::Independent, Metadata::Conserved,
                                               Metadata::GetUserFlag("Explicit"), Metadata::FillGhost}; // TODO TODO Restart
     auto m = Metadata(flags_prim_f);
     pkg->AddField("prims.fB", m);
     m = Metadata(flags_cons_f);
+    m.RegisterRefinementOps<parthenon::refinement_ops::ProlongateSharedMinMod, parthenon::refinement_ops::RestrictAverage, ProlongateInternalOlivares>();
     pkg->AddField("cons.fB", m);
 
     // Cell-centered versions.  Needed for BS, not for other schemes.
     // Probably will want to keep primitives for e.g. correct PtoU of MHD vars, but cons maybe can go
-    std::vector<MetadataFlag> flags_prim = {Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::GetUserFlag("Primitive"),
+    std::vector<MetadataFlag> flags_prim = {Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::GetUserFlag("GRPrimitive"),
                                             Metadata::GetUserFlag("MHD"), Metadata::GetUserFlag("Explicit"), Metadata::Vector};
-    std::vector<MetadataFlag> flags_cons = {Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::Conserved, Metadata::WithFluxes,
+    std::vector<MetadataFlag> flags_cons = {Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::GetUserFlag("GRConserved"), Metadata::WithFluxes,
                                             Metadata::GetUserFlag("MHD"), Metadata::GetUserFlag("Explicit"), Metadata::Vector};
     std::vector<int> s_vector({NVEC});
     m = Metadata(flags_prim, s_vector);
@@ -102,7 +104,7 @@ std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared
 
     // EMF on edges.
     // TODO only sync when needed
-    std::vector<MetadataFlag> flags_emf = {Metadata::Real, Metadata::Edge, Metadata::Derived, Metadata::OneCopy, Metadata::FillGhost};
+    std::vector<MetadataFlag> flags_emf = {Metadata::Real, Metadata::Edge, Metadata::Derived, Metadata::OneCopy};
     m = Metadata(flags_emf);
     pkg->AddField("B_CT.emf", m);
 
@@ -125,7 +127,6 @@ std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared
 
     // Register the other callbacks
     pkg->PostStepDiagnosticsMesh = B_CT::PostStepDiagnostics;
-    // TODO TODO prolongation/restriction will be registered here too
 
     // The definition of MaxDivB we care about actually changes per-transport,
     // so calculating it is handled by the transport package
@@ -138,15 +139,15 @@ std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared
 
     // List (vector) of HistoryOutputVars that will all be enrolled as output variables
     // LATER
-    // parthenon::HstVar_list hst_vars = {};
-    // hst_vars.emplace_back(parthenon::HistoryOutputVar(UserHistoryOperation::max, B_CT::MaxDivB, "MaxDivB"));
-    // // Event horizon magnetization.  Might be the same or different for different representations?
-    // if (pin->GetBoolean("coordinates", "spherical")) {
-    //     hst_vars.emplace_back(parthenon::HistoryOutputVar(UserHistoryOperation::sum, ReducePhi0, "Phi_0"));
-    //     hst_vars.emplace_back(parthenon::HistoryOutputVar(UserHistoryOperation::sum, ReducePhi5, "Phi_EH"));
-    // }
-    // // add callbacks for HST output to the Params struct, identified by the `hist_param_key`
-    // pkg->AddParam<>(parthenon::hist_param_key, hst_vars);
+    parthenon::HstVar_list hst_vars = {};
+    hst_vars.emplace_back(parthenon::HistoryOutputVar(UserHistoryOperation::max, B_CT::MaxDivB, "MaxDivB"));
+    // Event horizon magnetization.  Might be the same or different for different representations?
+    if (pin->GetBoolean("coordinates", "spherical")) {
+        // hst_vars.emplace_back(parthenon::HistoryOutputVar(UserHistoryOperation::sum, ReducePhi0, "Phi_0"));
+        // hst_vars.emplace_back(parthenon::HistoryOutputVar(UserHistoryOperation::sum, ReducePhi5, "Phi_EH"));
+    }
+    // add callbacks for HST output to the Params struct, identified by the `hist_param_key`
+    pkg->AddParam<>(parthenon::hist_param_key, hst_vars);
 
     return pkg;
 }

kharma/b_ct/b_ct.hpp

@@ -210,4 +210,122 @@ KOKKOS_INLINE_FUNCTION Real upwind_diff(const VariableFluxPack<Real>& B_U, const
     }
 }
 
+// Only by formatting may the following be made even a little comprehensible.
+
+template<int diff_face, int diff_side, int offset, int DIM>
+KOKKOS_FORCEINLINE_FUNCTION Real F(const ParArrayND<Real, VariableState> &fine, int l, int m, int n, int fk, int fj, int fi)
+{
+    // TODO compile-time error on misuse?
+    constexpr int df_is_k = 2*(diff_face == V3 && DIM > 2);
+    constexpr int df_is_j = 2*(diff_face == V2 && DIM > 1);
+    constexpr int df_is_i = 2*(diff_face == V1 && DIM > 0);
+    constexpr int ds_is_k = (diff_side == V3 && DIM > 2);
+    constexpr int ds_is_j = (diff_side == V2 && DIM > 1);
+    constexpr int ds_is_i = (diff_side == V1 && DIM > 0);
+    constexpr int of_is_k = (offset == V3 && DIM > 2);
+    constexpr int of_is_j = (offset == V2 && DIM > 1);
+    constexpr int of_is_i = (offset == V1 && DIM > 0);
+    // if(fi == 10 && fj == 10 && fk == 0) {
+    //     fprintf(stderr, "F facediff dir %d sidediff dirr %d off dir %d\nadding terms %d %d %d, -%d %d %d, -%d %d %d, %d %d %d\n",
+    //             diff_face, diff_side, offset,
+    //             df_is_i+ds_is_i+of_is_i, df_is_j+ds_is_j+of_is_j, df_is_k+ds_is_k+of_is_k,
+    //             ds_is_i+of_is_i,         ds_is_j+of_is_j,         ds_is_k+of_is_k,
+    //             df_is_i+of_is_i,         df_is_j+of_is_j,         df_is_k+of_is_k,
+    //             of_is_i                , of_is_j                , of_is_k);
+    // }
+    return fine(diff_face, l, m, n, fk+df_is_k+ds_is_k+of_is_k, fj+df_is_j+ds_is_j+of_is_j, fi+df_is_i+ds_is_i+of_is_i)
+         - fine(diff_face, l, m, n, fk+ds_is_k+of_is_k        , fj+ds_is_j+of_is_j        , fi+ds_is_i+of_is_i)
+         - fine(diff_face, l, m, n, fk+df_is_k+of_is_k        , fj+df_is_j+of_is_j        , fi+df_is_i+of_is_i)
+         + fine(diff_face, l, m, n, fk+of_is_k                , fj+of_is_j                , fi+of_is_i);
+}
+
+struct ProlongateInternalOlivares {
+  static constexpr bool OperationRequired(TopologicalElement fel,
+                                          TopologicalElement cel) {
+    return fel == cel;
+  }
+
+  template <int DIM, TopologicalElement el = TopologicalElement::CC,
+            TopologicalElement cel = TopologicalElement::CC>
+  KOKKOS_FORCEINLINE_FUNCTION static void
+  Do(const int l, const int m, const int n, const int k, const int j, const int i,
+     const IndexRange &ckb, const IndexRange &cjb, const IndexRange &cib,
+     const IndexRange &kb, const IndexRange &jb, const IndexRange &ib,
+     const Coordinates_t &coords, const Coordinates_t &coarse_coords,
+     const ParArrayND<Real, VariableState> *,
+     const ParArrayND<Real, VariableState> *pfine) {
+
+    // Definitely exit on what we can't handle
+    if constexpr (el != TE::F1 && el != TE::F2 && el != TE::F3)
+        return;
+
+    // Handle permutations "naturally."
+    // Olivares et al. is fond of listing x1 versions which permute,
+    // this makes translating/checking those easier
+    constexpr int me = static_cast<int>(el) % 3;
+    constexpr int next = (me+1) % 3;
+    constexpr int third = (me+2) % 3;
+
+    // Exit if we're computing a trivial direction
+    if constexpr ((me == V3 && !(DIM > 2)) || (me == V2 && !(DIM > 1)) || (me == V1 && !(DIM > 0)))
+        return;
+
+    // Fine array, indices
+    auto &fine = *pfine;
+    const int fi = (DIM > 0) ? (i - cib.s) * 2 + ib.s : ib.s;
+    const int fj = (DIM > 1) ? (j - cjb.s) * 2 + jb.s : jb.s;
+    const int fk = (DIM > 2) ? (k - ckb.s) * 2 + kb.s : kb.s;
+
+    // Coefficients selecting a particular formula (see Olivares et al. 2019)
+    // TODO options here. This corresponds to Cunningham, but we could have:
+    // 1. differences of squares of zone dimesnions (Toth)
+    // 2. heuristic based on flux difference of top vs bottom halves (Olivares)
+    //constexpr Real a[3] = {0., 0., 0.};
+    const Real a[3] = {(SQR(coords.Dxc<2>(j)) - SQR(coords.Dxc<3>(k)))/(SQR(coords.Dxc<2>(j)) + SQR(coords.Dxc<3>(k))),
+                       (SQR(coords.Dxc<3>(k)) - SQR(coords.Dxc<1>(i)))/(SQR(coords.Dxc<3>(k)) + SQR(coords.Dxc<1>(i))),
+                       (SQR(coords.Dxc<1>(i)) - SQR(coords.Dxc<2>(j)))/(SQR(coords.Dxc<1>(i)) + SQR(coords.Dxc<2>(j)))};
+
+    // Coefficients for each term evaluating the four sub-faces
+    const Real coeff[4][4] = {{3 + a[next], 1 - a[next], 3 - a[third], 1 + a[third]},
+                              {3 + a[next], 1 - a[next], 1 + a[third], 3 - a[third]},
+                              {1 - a[next], 3 + a[next], 3 - a[third], 1 + a[third]},
+                              {1 - a[next], 3 + a[next], 1 + a[third], 3 - a[third]}};
+
+    constexpr int diff_k = (me == V3), diff_j = (me == V2), diff_i = (me == V1);
+    // if(fi == 10 && fj == 10 && fk == 0) {
+    //     fprintf(stderr, "Prolongating %d %d %d EL %d, DIM %d\n", fi, fj, fk, static_cast<int>(el), DIM);
+    //     fprintf(stderr, "Differencing %d %d %d\n", diff_i, diff_j, diff_k);
+    // }
+
+    // Iterate through the 4 sub-faces
+    for (int elem=0; elem < 4; elem++) {
+        // Make sure we can offset in other directions before doing so, though
+        // TODO eliminate redundant work or template these so the compiler can?
+        const int off_i = (DIM > 0) ? elem%2*(me == V2) + elem/2*(me == V3) + (me == V1) : 0;
+        const int off_j = (DIM > 1) ? elem%2*(me == V3) + elem/2*(me == V1) + (me == V2) : 0;
+        const int off_k = (DIM > 2) ? elem%2*(me == V1) + elem/2*(me == V2) + (me == V3) : 0;
+
+        fine(me, l, m, n, fk+off_k, fj+off_j, fi+off_i) =
+            // Average faces on either side of us in selected direction (diff), on each of the 4 sub-faces (off)
+            0.5*(fine(me, l, m, n, fk+off_k-diff_k, fj+off_j-diff_j, fi+off_i-diff_i) +
+                 fine(me, l, m, n, fk+off_k+diff_k, fj+off_j+diff_j, fi+off_i+diff_i)) +
+            1./16*(coeff[elem][0]*F<next ,me,-1,DIM>(fine, l, m, n, fk, fj, fi) + coeff[elem][1]*F<next,me,third,DIM>(fine, l, m, n, fk, fj, fi)
+                 + coeff[elem][2]*F<third,me,-1,DIM>(fine, l, m, n, fk, fj, fi) + coeff[elem][3]*F<third,me,next,DIM>(fine, l, m, n, fk, fj, fi));
+
+        // if(fi == 10 && fj == 10 && fk == 0 && me == V1) {
+        //     fprintf(stderr, "Elem %d Offset %d %d %d set %g\n", elem, off_i, off_j, off_k, fine(me, l, m, n, fk+off_k, fj+off_j, fi+off_i));
+        //     fprintf(stderr, "Averaging faces %d %d %d and %d %d %d (%g & %g)\n", fi+off_i-diff_i, fj+off_j-diff_j, fk+off_k-diff_k, 
+        //         fi+off_i+diff_i, fj+off_j+diff_j, fk+off_k+diff_k, 
+        //         fine(me, l, m, n, fk+off_k-diff_k, fj+off_j-diff_j, fi+off_i-diff_i),
+        //         fine(me, l, m, n, fk+off_k+diff_k, fj+off_j+diff_j, fi+off_i+diff_i));
+        //     fprintf(stderr, "Coeffs %g %g %g %g\n", coeff[elem][0]*F<next,me,-1,DIM>(fine, l, m, n, fk, fj, fi),
+        //                                             coeff[elem][1]*F<next,me,third,DIM>(fine, l, m, n, fk, fj, fi),
+        //                                             coeff[elem][2]*F<third,me,-1,DIM>(fine, l, m, n, fk, fj, fi),
+        //                                             coeff[elem][3]*F<third,me,next,DIM>(fine, l, m, n, fk, fj, fi));
+        // }
+
+    }
+  }
+};
+
 }

kharma/b_flux_ct/b_flux_ct.cpp

@@ -102,10 +102,10 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
     MetadataFlag areWeImplicit = (implicit_b) ? Metadata::GetUserFlag("Implicit")
                                               : Metadata::GetUserFlag("Explicit");
 
-    // Flags for B fields.  "Primitive" form is field, "conserved" is flux
-    std::vector<MetadataFlag> flags_prim = {Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::GetUserFlag("Primitive"),
+    // Flags for B fields.  "GRPrimitive" form is field, "conserved" is flux
+    std::vector<MetadataFlag> flags_prim = {Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::GetUserFlag("GRPrimitive"),
                                             Metadata::Restart, Metadata::GetUserFlag("MHD"), areWeImplicit, Metadata::Vector};
-    std::vector<MetadataFlag> flags_cons = {Metadata::Real, Metadata::Cell, Metadata::Independent, Metadata::Conserved,
+    std::vector<MetadataFlag> flags_cons = {Metadata::Real, Metadata::Cell, Metadata::Independent, Metadata::GetUserFlag("GRConserved"), Metadata::Conserved,
                                             Metadata::WithFluxes, Metadata::FillGhost, Metadata::GetUserFlag("MHD"), areWeImplicit, Metadata::Vector};
 
     auto m = Metadata(flags_prim, s_vector);

kharma/b_flux_ct/b_flux_ct.hpp

@@ -46,7 +46,7 @@
  *
  * This requires only the values at cell centers
  * 
- * This implementation includes conversion from "primitive" to "conserved" B and back
+ * This implementation includes conversion from "GRPrimitive" to "conserved" B and back
  */
 namespace B_FluxCT {
 /**

kharma/domain.hpp

@@ -42,7 +42,6 @@ namespace KDomain {
 
 /**
  * Functions for checking boundaries in 3D.
- * Uses IndexRange objects, or this would be in kharma_utils.hpp
  */
 KOKKOS_INLINE_FUNCTION bool outside(const int& k, const int& j, const int& i,
                                     const IndexRange& kb, const IndexRange& jb, const IndexRange& ib)
@@ -66,6 +65,7 @@ KOKKOS_INLINE_FUNCTION bool inside(const int& k, const int& j, const int& i, con
 }
 
 // TODO(BSP) these really should be in Parthenon
+// There's a templated way to do it I forget, but this would be easier
 template<typename T>
 inline const int& GetNDim(MeshBlockData<T>* rc)
 { return rc->GetBlockPointer()->pmy_mesh->ndim; }

kharma/driver/imex_step.cpp

@@ -171,7 +171,7 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta
         // If evolving GRMHD explicitly, UtoP needs a guess in order to converge, so we copy in md_sub_step_init
         auto t_copy_prims = t_none;
         if (!pkgs.at("GRMHD")->Param<bool>("implicit")) {
-            t_copy_prims = tl.AddTask(t_none, Copy<MeshData<Real>>, std::vector<MetadataFlag>({Metadata::GetUserFlag("HD"), Metadata::GetUserFlag("Primitive")}),
+            t_copy_prims = tl.AddTask(t_none, Copy<MeshData<Real>>, std::vector<MetadataFlag>({Metadata::GetUserFlag("HD"), Metadata::GetUserFlag("GRPrimitive")}),
                                       md_sub_step_init.get(), md_solver.get());
         }
 
@@ -199,7 +199,7 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta
             // Copy the primitives to the `linesearch` MeshData object if linesearch was enabled.
             auto t_copy_linesearch = t_guess_ready;
             if (use_linesearch) {
-                t_copy_linesearch = tl.AddTask(t_guess_ready, Copy<MeshData<Real>>, std::vector<MetadataFlag>({Metadata::GetUserFlag("Primitive")}),
+                t_copy_linesearch = tl.AddTask(t_guess_ready, Copy<MeshData<Real>>, std::vector<MetadataFlag>({Metadata::GetUserFlag("GRPrimitive")}),
                                                 md_solver.get(), md_linesearch.get());
             }
 

kharma/driver/kharma_driver.hpp

@@ -68,7 +68,7 @@ class KHARMADriver : public MultiStageDriver {
          * so that the driver can repeat calls to create a predictor-corrector, RK2/4, etc.
          * 
          * Unlike MHD, GRMHD must keep two forms of the variables: the conserved variables, and a set of
-         * "primitive" variables more amenable to reconstruction.  To evolve the fluid, the code must:
+         * "GRPrimitive" variables more amenable to reconstruction.  To evolve the fluid, the code must:
          * 1. Reconstruct the right- and left-going components at zone faces, given the primitive variables
          * 2. Calculate the fluxes of conserved quantities through the faces
          * 2a. Apply any fixes to fluxes (e.g., for the magnetic field)