B Cleanup for Faces

Mainly this PR adds a laplacian for face-centered B to `CleanupDivergence`,
which is automatically used if the B_CT package is loaded.  It doesn't actually
seem faster than the cell-centered version despite having a smaller stencil.
It also adds a "PtoU" for face-centered fields, which takes a cell-centered B
and interpolates it to faces.  Yes, this means interpolating interpolated values,
which is not ideal but doesn't seem egregious.
Finally, it adds a version of the `resize` test with face-centered fields.

I also rephrased the face divB in terms of derivatives for clarity, spit out Flux-CT
device functions into a header like Face-CT has, and modified some uncommon
device-side B field functions a bit.

Finally, this adds an unrelated feature allowing re-initialization of electron temps
at user option, designed for restarting ideal simulations w/electrons enabled.

kharma/b_cleanup/b_cleanup.hpp

@@ -71,10 +71,18 @@ bool CleanupThisStep(Mesh* pmesh, int nstep);
  * Extra MeshData arg is just to satisfy Parthenon solver calling convention
  */
 TaskStatus CornerLaplacian(MeshData<Real>* md, const std::string& p_var, MeshData<Real>* md_again, const std::string& lap_var);
+/**
+ * Calculate the laplacian using divergence at centers.
+ */
+TaskStatus CenterLaplacian(MeshData<Real>* md, const std::string& p_var, MeshData<Real>* md_again, const std::string& lap_var);
 
 /**
- * Apply B -= grad(P) to subtract divergence from the magnetic field
+ * Apply B -= grad(P) on cell centers to subtract divergence from the magnetic field
+ */
+TaskStatus ApplyPCenter(MeshData<Real> *msolve, MeshData<Real> *md);
+/**
+ * Apply B -= grad(P) on faces to subtract divergence from the magnetic field
  */
-TaskStatus ApplyP(MeshData<Real> *msolve, MeshData<Real> *md);
+TaskStatus ApplyPFace(MeshData<Real> *msolve, MeshData<Real> *md);
 
 } // namespace B_Cleanup

kharma/b_ct/b_ct.cpp

@@ -202,6 +202,74 @@ TaskStatus B_CT::BlockUtoP(MeshBlockData<Real> *rc, IndexDomain domain, bool coa
     return TaskStatus::complete;
 }
 
+TaskStatus B_CT::DangerousPtoU(MeshData<Real> *md, IndexDomain domain, bool coarse)
+{
+    auto B_Uf = md->PackVariables(std::vector<std::string>{"cons.fB"});
+    auto B_U = md->PackVariables(std::vector<std::string>{"cons.B"});
+    auto B_P = md->PackVariables(std::vector<std::string>{"prims.B"});
+
+    // Figure out indices
+    const IndexRange block = IndexRange{0, B_Uf.GetDim(5)-1};
+
+    auto pmb0 = md->GetBlockData(0)->GetBlockPointer();
+    // Average the primitive vals to faces and multiply by gdet
+    const IndexRange3 bf1 = KDomain::GetRange(md, domain, F1, coarse);
+    pmb0->par_for("PtoU_B_F1", block.s, block.e, bf1.ks, bf1.ke, bf1.js, bf1.je, bf1.is, bf1.ie,
+        KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i) {
+            const auto& G = B_Uf.GetCoords(b);
+            B_Uf(b, F1, 0, k, j, i) = G.gdet(Loci::face1, j, i) * (B_P(b, V1, k, j, i-1) + B_P(b, V1, k, j, i)) / 2;
+        }
+    );
+    const IndexRange3 bf2 = KDomain::GetRange(md, domain, F2, coarse);
+    pmb0->par_for("PtoU_B_F2", block.s, block.e, bf2.ks, bf2.ke, bf2.js, bf2.je, bf2.is, bf2.ie,
+        KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i) {
+            const auto& G = B_Uf.GetCoords(b);
+            B_Uf(b, F2, 0, k, j, i) = G.gdet(Loci::face2, j, i) * (B_P(b, V2, k, j-1, i) + B_P(b, V2, k, j, i)) / 2;
+        }
+    );
+    const IndexRange3 bf3 = KDomain::GetRange(md, domain, F3, coarse);
+    pmb0->par_for("PtoU_B_F3", block.s, block.e, bf3.ks, bf3.ke, bf3.js, bf3.je, bf3.is, bf3.ie,
+        KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i) {
+            const auto& G = B_Uf.GetCoords(b);
+            B_Uf(b, F3, 0, k, j, i) = G.gdet(Loci::face3, j, i) * (B_P(b, V3, k-1, j, i) + B_P(b, V3, k, j, i)) / 2;
+        }
+    );
+
+    // Make sure B on poles is still zero, even though we've interpolated
+    if (pmb0->coords.coords.is_spherical()) {
+        for (int i=0; i < md->GetMeshPointer()->GetNumMeshBlocksThisRank(); i++) {
+            auto rc = md->GetBlockData(i);
+            auto pmb = rc->GetBlockPointer();
+            auto dB_block = rc->PackVariables(std::vector<std::string>{"dB"});
+            if (pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user) {
+                pmb->par_for("dB_boundary", bf2.ks, bf2.ke, bf2.is, bf2.ie,
+                    KOKKOS_LAMBDA (const int &k, const int &i) {
+                        dB_block(F2, 0, k, bf2.js, i) = 0.;
+                    }
+                );
+            }
+            if (pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user) {
+                pmb->par_for("dB_boundary", bf2.ks, bf2.ke, bf2.is, bf2.ie,
+                    KOKKOS_LAMBDA (const int &k, const int &i) {
+                        dB_block(F2, 0, k, bf2.je, i) = 0.;
+                    }
+                );
+            }
+        }
+    }
+
+    // Also recover conserved B at centers, just in case
+    const IndexRange3 bc = KDomain::GetRange(md, domain, CC, coarse);
+    pmb0->par_for("UtoP_B_centerPtoU", block.s, block.e, 0, NVEC-1, bc.ks, bc.ke, bc.js, bc.je, bc.is, bc.ie,
+        KOKKOS_LAMBDA (const int &b, const int &v, const int &k, const int &j, const int &i) {
+            const auto& G = B_U.GetCoords(b);
+            B_U(b, v, k, j, i) = B_P(b, v, k, j, i) * G.gdet(Loci::center, j, i);
+        }
+    );
+
+    return TaskStatus::complete;
+}
+
 TaskStatus B_CT::CalculateEMF(MeshData<Real> *md)
 {
     auto pmesh = md->GetMeshPointer();

kharma/b_ct/b_ct.hpp

@@ -52,17 +52,23 @@ namespace B_CT {
 std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<Packages_t>& packages);
 
 /**
- * Get the primitive variables, which in Parthenon's nomenclature are "derived".
- * Also applies floors to the calculated primitives, and fixes up any inversion errors
+ * Get the primitive variables, which in Parthenon's nomenclature are "derived"
  * 
  * Defaults to entire domain, as the KHARMA algorithm relies on applying UtoP over ghost zones.
  * 
- * input: Conserved B = sqrt(-gdet) * B^i
+ * input: Conserved B = sqrt(-g) * B^i
  * output: Primitive B = B^i
  */
 TaskStatus BlockUtoP(MeshBlockData<Real> *mbd, IndexDomain domain, bool coarse=false);
 TaskStatus MeshUtoP(MeshData<Real> *md, IndexDomain domain, bool coarse=false);
 
+/**
+ * Interpolate primitive B to faces, then multiply by gdet
+ * DANGEROUS as it replaces conserved field state. Only used for resizing runs,
+ * and only with B_Cleanup::CleanupDivergence directly afterward!
+ */
+TaskStatus DangerousPtoU(MeshData<Real> *md, IndexDomain domain, bool coarse=false);
+
 /**
  * Calculate the EMF around edges of faces caused by the flux of B field
  * through each face.

kharma/b_ct/b_ct_functions.hpp

@@ -47,12 +47,31 @@ namespace B_CT {
 template<typename Global>
 KOKKOS_INLINE_FUNCTION Real face_div(const GRCoordinates &G, Global &v, const int &ndim, const int &k, const int &j, const int &i)
 {
-    Real du = (v(F1, 0, k, j, i + 1) * G.Volume<F1>(k, j, i + 1) - v(F1, 0, k, j, i) * G.Volume<F1>(k, j, i));
+    Real du = (v(F1, 0, k, j, i + 1) - v(F1, 0, k, j, i)) / G.Dxc<1>(k, j, i);
     if (ndim > 1)
-        du += (v(F2, 0, k, j + 1, i) * G.Volume<F2>(k, j + 1, i) - v(F2, 0, k, j, i) * G.Volume<F2>(k, j, i));
+        du += (v(F2, 0, k, j + 1, i) - v(F2, 0, k, j, i)) / G.Dxc<2>(k, j, i);
     if (ndim > 2)
-        du += (v(F3, 0, k + 1, j, i) * G.Volume<F3>(k + 1, j, i) - v(F3, 0, k, j, i) * G.Volume<F3>(k, j, i));
-    return du / G.Volume<CC>(k, j, i);
+        du += (v(F3, 0, k + 1, j, i) - v(F3, 0, k, j, i)) / G.Dxc<3>(k, j, i);
+    return du;
+}
+
+/**
+ * Gradient on one face, denoted by DIR. Split vs Flux-CT version because generally
+ * this will be needed in separate loops setting F1, F2, F3, with separate bounds
+ * Note this is backward-difference: face N borders cells N-1 & N
+ */
+template<CoordinateDirection DIR>
+KOKKOS_FORCEINLINE_FUNCTION double face_grad(const GRCoordinates& G, const VariablePack<Real>& P,
+                                          const int& k, const int& j, const int& i)
+{
+    // Backward gradient to cell faces
+    if constexpr (DIR == X1DIR) {
+        return (P(0, k, j, i) - P(0, k, j, i-1)) / G.Dxc<1>(k, j, i);
+    } else if constexpr (DIR == X2DIR) {
+        return (P(0, k, j, i) - P(0, k, j-1, i)) / G.Dxc<2>(k, j, i);
+    } else if constexpr (DIR == X3DIR) {
+        return (P(0, k, j, i) - P(0, k-1, j, i)) / G.Dxc<3>(k, j, i);
+    }
 }
 
 template<TE el, int NDIM>

kharma/b_flux_ct/b_flux_ct.cpp

@@ -594,7 +594,7 @@ double MaxDivB(MeshData<Real> *md)
         pmb->par_reduce("divB_max", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
             KOKKOS_LAMBDA (const int &k, const int &j, const int &i, double &local_result) {
                 const auto& G = B_U.GetCoords(b);
-                const double local_divb = m::abs(corner_div(G, B_U, b, k, j, i, ndim > 2));
+                const double local_divb = m::abs(corner_div(G, B_U(b), k, j, i, ndim > 2));
                 if (local_divb > local_result) local_result = local_divb;
             }
         , max_reducer);
@@ -669,7 +669,7 @@ void CalcDivB(MeshData<Real> *md, std::string divb_field_name)
         pmb->par_for("calc_divB", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
             KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
                 const auto& G = B_U.GetCoords(b);
-                divB(b, 0, k, j, i) = corner_div(G, B_U, b, k, j, i, ndim > 2);
+                divB(b, 0, k, j, i) = corner_div(G, B_U(b), k, j, i, ndim > 2);
             }
         );
     }
@@ -695,7 +695,7 @@ void FillOutput(MeshBlock *pmb, ParameterInput *pin)
     pmb->par_for("divB_output", kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
         KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
             const auto& G = B_U.GetCoords();
-            divB(0, k, j, i) = corner_div(G, B_U, 0, k, j, i, ndim > 2);
+            divB(0, k, j, i) = corner_div(G, B_U, k, j, i, ndim > 2);
         }
     );
 }

kharma/b_flux_ct/b_flux_ct.hpp

@@ -33,6 +33,7 @@
  */
 #pragma once
 
+#include "b_flux_ct_functions.hpp"
 #include "decs.hpp"
 #include "grmhd_functions.hpp"
 #include "reductions.hpp"
@@ -149,149 +150,4 @@ inline Real ReducePhi5(MeshData<Real> *md)
     return Reductions::EHReduction<Reductions::Var::phi, Real>(md, UserHistoryOperation::sum, 5);
 }
 
-/**
- * ND divergence, averaging to cell corners
- * TODO likely better templated, as with all ND stuff
- */
-template<typename Global>
-KOKKOS_FORCEINLINE_FUNCTION double corner_div(const GRCoordinates& G, const Global& B_U, const int& b,
-                                         const int& k, const int& j, const int& i, const bool& do_3D)
-{
-    const double norm = (do_3D) ? 0.25 : 0.5;
-    // 2D divergence, averaging to corners
-    double term1 = B_U(b, V1, k, j, i)   - B_U(b, V1, k, j, i-1) +
-                   B_U(b, V1, k, j-1, i) - B_U(b, V1, k, j-1, i-1);
-    double term2 = B_U(b, V2, k, j, i)   - B_U(b, V2, k, j-1, i) +
-                   B_U(b, V2, k, j, i-1) - B_U(b, V2, k, j-1, i-1);
-    double term3 = 0.;
-    if (do_3D) {
-        // Average to corners in 3D, add 3rd flux
-        term1 +=  B_U(b, V1, k-1, j, i)   + B_U(b, V1, k-1, j-1, i)
-                - B_U(b, V1, k-1, j, i-1) - B_U(b, V1, k-1, j-1, i-1);
-        term2 +=  B_U(b, V2, k-1, j, i)   + B_U(b, V2, k-1, j, i-1)
-                - B_U(b, V2, k-1, j-1, i) - B_U(b, V2, k-1, j-1, i-1);
-        term3 =   B_U(b, V3, k, j, i)     + B_U(b, V3, k, j-1, i)
-                + B_U(b, V3, k, j, i-1)   + B_U(b, V3, k, j-1, i-1)
-                - B_U(b, V3, k-1, j, i)   - B_U(b, V3, k-1, j-1, i)
-                - B_U(b, V3, k-1, j, i-1) - B_U(b, V3, k-1, j-1, i-1);
-    }
-    return norm*term1/G.Dxc<1>(i) + norm*term2/G.Dxc<2>(j) + norm*term3/G.Dxc<3>(k);
-}
-template<typename Global>
-KOKKOS_FORCEINLINE_FUNCTION double corner_div(const GRCoordinates& G, const Global& P, const VarMap& m_p, 
-                                         const int& b, const int& k, const int& j, const int& i,
-                                         const bool& do_3D)
-{
-    const double norm = (do_3D) ? 0.25 : 0.5;
-    // 2D divergence, averaging to corners
-    double term1 = P(b, m_p.B1, k, j, i)   - P(b, m_p.B1, k, j, i-1) +
-                   P(b, m_p.B1, k, j-1, i) - P(b, m_p.B1, k, j-1, i-1);
-    double term2 = P(b, m_p.B2, k, j, i)   - P(b, m_p.B2, k, j-1, i) +
-                   P(b, m_p.B2, k, j, i-1) - P(b, m_p.B2, k, j-1, i-1);
-    double term3 = 0.;
-    if (do_3D) {
-        // Average to corners in 3D, add 3rd flux
-        term1 +=  P(b, m_p.B1, k-1, j, i)   + P(b, m_p.B1, k-1, j-1, i)
-                - P(b, m_p.B1, k-1, j, i-1) - P(b, m_p.B1, k-1, j-1, i-1);
-        term2 +=  P(b, m_p.B2, k-1, j, i)   + P(b, m_p.B2, k-1, j, i-1)
-                - P(b, m_p.B2, k-1, j-1, i) - P(b, m_p.B2, k-1, j-1, i-1);
-        term3 =   P(b, m_p.B3, k, j, i)     + P(b, m_p.B3, k, j-1, i)
-                + P(b, m_p.B3, k, j, i-1)   + P(b, m_p.B3, k, j-1, i-1)
-                - P(b, m_p.B3, k-1, j, i)   - P(b, m_p.B3, k-1, j-1, i)
-                - P(b, m_p.B3, k-1, j, i-1) - P(b, m_p.B3, k-1, j-1, i-1);
-    }
-    return norm*term1/G.Dxc<1>(i) + norm*term2/G.Dxc<2>(j) + norm*term3/G.Dxc<3>(k);
-}
-
-/**
- * 2D or 3D gradient, averaging to cell centers from corners.
- * Note this is forward-difference, while previous def is backward
- */
-template<typename Global>
-KOKKOS_FORCEINLINE_FUNCTION void center_grad(const GRCoordinates& G, const Global& P, const int& b,
-                                          const int& k, const int& j, const int& i, const bool& do_3D,
-                                          double& B1, double& B2, double& B3)
-{
-    const double norm = (do_3D) ? 0.25 : 0.5;
-    // 2D gradient, averaging to centers
-    double term1 =  P(b, 0, k, j+1, i+1) + P(b, 0, k, j, i+1)
-                  - P(b, 0, k, j+1, i)   - P(b, 0, k, j, i);
-    double term2 =  P(b, 0, k, j+1, i+1) + P(b, 0, k, j+1, i)
-                  - P(b, 0, k, j, i+1)   - P(b, 0, k, j, i);
-    double term3 = 0.;
-    if (do_3D) {
-        // Average to centers in 3D, add 3rd flux
-        term1 += P(b, 0, k+1, j+1, i+1) + P(b, 0, k+1, j, i+1)
-               - P(b, 0, k+1, j+1, i)   - P(b, 0, k+1, j, i);
-        term2 += P(b, 0, k+1, j+1, i+1) + P(b, 0, k+1, j+1, i)
-               - P(b, 0, k+1, j, i+1)   - P(b, 0, k+1, j, i);
-        term3 =  P(b, 0, k+1, j+1, i+1) + P(b, 0, k+1, j, i+1)
-               + P(b, 0, k+1, j+1, i)   + P(b, 0, k+1, j, i)
-               - P(b, 0, k, j+1, i+1)   - P(b, 0, k, j, i+1)
-               - P(b, 0, k, j+1, i)     - P(b, 0, k, j, i);
-    }
-    B1 = norm*term1/G.Dxc<1>(i);
-    B2 = norm*term2/G.Dxc<2>(j);
-    B3 = norm*term3/G.Dxc<3>(k);
-}
-
-KOKKOS_FORCEINLINE_FUNCTION void averaged_curl_3D(const GRCoordinates& G, const GridVector& A, const GridVector& B_U,
-                                             const int& k, const int& j, const int& i)
-{
-    // Take a flux-ct step from the corner potentials.
-    // This needs to be 3D because post-tilt A may not point in the phi direction only
-
-    // A3,2 derivative
-    const Real A3c2f = (A(V3, k, j + 1, i)     + A(V3, k, j + 1, i + 1) + 
-                        A(V3, k + 1, j + 1, i) + A(V3, k + 1, j + 1, i + 1)) / 4;
-    const Real A3c2b = (A(V3, k, j, i)     + A(V3, k, j, i + 1) +
-                        A(V3, k + 1, j, i) + A(V3, k + 1, j, i + 1)) / 4;
-    // A2,3 derivative
-    const Real A2c3f = (A(V2, k + 1, j, i)     + A(V2, k + 1, j, i + 1) +
-                        A(V2, k + 1, j + 1, i) + A(V2, k + 1, j + 1, i + 1)) / 4;
-    const Real A2c3b = (A(V2, k, j, i)     + A(V2, k, j, i + 1) +
-                        A(V2, k, j + 1, i) + A(V2, k, j + 1, i + 1)) / 4;
-    B_U(V1, k, j, i) = (A3c2f - A3c2b) / G.Dxc<2>(j) - (A2c3f - A2c3b) / G.Dxc<3>(k);
-
-    // A1,3 derivative
-    const Real A1c3f = (A(V1, k + 1, j, i)     + A(V1, k + 1, j, i + 1) + 
-                        A(V1, k + 1, j + 1, i) + A(V1, k + 1, j + 1, i + 1)) / 4;
-    const Real A1c3b = (A(V1, k, j, i)     + A(V1, k, j, i + 1) +
-                        A(V1, k, j + 1, i) + A(V1, k, j + 1, i + 1)) / 4;
-    // A3,1 derivative
-    const Real A3c1f = (A(V3, k, j, i + 1)     + A(V3, k + 1, j, i + 1) +
-                        A(V3, k, j + 1, i + 1) + A(V3, k + 1, j + 1, i + 1)) / 4;
-    const Real A3c1b = (A(V3, k, j, i)     + A(V3, k + 1, j, i) +
-                        A(V3, k, j + 1, i) + A(V3, k + 1, j + 1, i)) / 4;
-    B_U(V2, k, j, i) = (A1c3f - A1c3b) / G.Dxc<3>(k) - (A3c1f - A3c1b) / G.Dxc<1>(i);
-
-    // A2,1 derivative
-    const Real A2c1f = (A(V2, k, j, i + 1)     + A(V2, k, j + 1, i + 1) + 
-                        A(V2, k + 1, j, i + 1) + A(V2, k + 1, j + 1, i + 1)) / 4;
-    const Real A2c1b = (A(V2, k, j, i)     + A(V2, k, j + 1, i) +
-                        A(V2, k + 1, j, i) + A(V2, k + 1, j + 1, i)) / 4;
-    // A1,2 derivative
-    const Real A1c2f = (A(V1, k, j + 1, i)     + A(V1, k, j + 1, i + 1) +
-                        A(V1, k + 1, j + 1, i) + A(V1, k + 1, j + 1, i + 1)) / 4;
-    const Real A1c2b = (A(V1, k, j, i)     + A(V1, k, j, i + 1) +
-                        A(V1, k + 1, j, i) + A(V1, k + 1, j, i + 1)) / 4;
-    B_U(V3, k, j, i) = (A2c1f - A2c1b) / G.Dxc<1>(i) - (A1c2f - A1c2b) / G.Dxc<2>(j);
-}
-
-KOKKOS_FORCEINLINE_FUNCTION void averaged_curl_2D(const GRCoordinates& G, const GridVector& A, const GridVector& B_U,
-                                             const int& k, const int& j, const int& i)
-{
-    // A3,2 derivative
-    const Real A3c2f = (A(V3, k, j + 1, i) + A(V3, k, j + 1, i + 1)) / 2;
-    const Real A3c2b = (A(V3, k, j, i)     + A(V3, k, j, i + 1)) / 2;
-    B_U(V1, k, j, i) = (A3c2f - A3c2b) / G.Dxc<2>(j);
-
-    // A3,1 derivative
-    const Real A3c1f = (A(V3, k, j, i + 1) + A(V3, k, j + 1, i + 1)) / 2;
-    const Real A3c1b = (A(V3, k, j, i)     + A(V3, k, j + 1, i)) / 2;
-    B_U(V2, k, j, i) = - (A3c1f - A3c1b) / G.Dxc<1>(i);
-
-    B_U(V3, k, j, i) = 0;
-}
-
 }