Fix the boundaries of EMHD problems

All MHD variables in bondi_viscous now converge as expected, and
boundaries are applied to dP as expected. Source term seems to be
much, much too large for some reason.

Also Vbump Kokkos to fix a CUDA segfault (again?)

kharma/b_ct/b_ct.cpp

@@ -73,7 +73,7 @@ std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared
     params.Add("ct_scheme", ct_scheme);
     // Use the default Parthenon prolongation operator, rather than the divergence-preserving one
     // This relies entirely on the EMF communication for preserving the divergence
-    bool lazy_prolongation = pin->GetOrAddBoolean("b_field", "lazy_prolongation", true);
+    bool lazy_prolongation = pin->GetOrAddBoolean("b_field", "lazy_prolongation", false);
     // Need to preserve divergence if you refine/derefine during sim i.e. AMR
     if (lazy_prolongation && pin->GetString("parthenon/mesh", "refinement") == "adaptive")
         throw std::runtime_error("Cannot use non-preserving prolongation in AMR!");
@@ -83,6 +83,7 @@ std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared
     // Flags for B fields on faces.
     // We don't mark these as "Primitive" and "Conserved" else they'd be bundled
     // with all the cell vars in a bunch of places we don't want
+    // Also note we *always* sync B field conserved var
     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,
@@ -172,6 +173,8 @@ void B_CT::BlockUtoP(MeshBlockData<Real> *rc, IndexDomain domain, bool coarse)
     auto B_U = rc->PackVariables(std::vector<std::string>{"cons.B"});
     auto B_P = rc->PackVariables(std::vector<std::string>{"prims.B"});
     const auto& G = pmb->coords;
+    // Return if we're not syncing U & P at all (e.g. edges)
+    if (B_Uf.GetDim(4) == 0) return;
 
     // TODO get rid of prims on faces probably
 
@@ -213,7 +216,7 @@ TaskStatus B_CT::CalculateEMF(MeshData<Real> *md)
 
     // Figure out indices
     const IndexRange3 b = KDomain::GetRange(md, IndexDomain::interior, 0, 0);
-    const IndexRange3 b1 = KDomain::GetRange(md, IndexDomain::interior, -1, 1);
+    const IndexRange3 b1 = KDomain::GetRange(md, IndexDomain::interior, -1, 2);
     const IndexRange block = IndexRange{0, emf_pack.GetDim(5)-1};
     const int kd = ndim > 2 ? 1 : 0;
     const int jd = ndim > 1 ? 1 : 0;

kharma/b_flux_ct/b_flux_ct.cpp

@@ -109,6 +109,14 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
     auto flags_cons = packages->Get("Driver")->Param<std::vector<MetadataFlag>>("cons_flags");
     flags_cons.insert(flags_cons.end(), flags_b.begin(), flags_b.end());
 
+    // Always sync B field conserved var, for standardization with B_CT
+    if (!flags_cons.count(Metadata::FillGhost)) {
+        flags_cons.push_back(Metadata::FillGhost);
+    }
+    if (flags_prims.count(Metadata::FillGhost)) {
+        flags_prims.erase(std::remove(flags_prims.begin(), flags_prims.end(), Metadata::FillGhost), flags_prims.end()); 
+    }
+
     auto m = Metadata(flags_prim, s_vector);
     pkg->AddField("prims.B", m);
     m = Metadata(flags_cons, s_vector);

kharma/boundaries/boundaries.cpp

@@ -128,6 +128,11 @@ std::shared_ptr<KHARMAPackage> KBoundaries::Initialize(ParameterInput *pin, std:
         bool zero_flux = pin->GetOrAddBoolean("boundaries", "zero_flux_" + bname, zero_polar_flux && bdir == X2DIR);
         params.Add("zero_flux_" + bname, zero_flux);
 
+        // Allow specifically dP to outflow in otherwise Dirichlet conditions
+        // Only used for viscous_bondi problem
+        bool outflow_EMHD = pin->GetOrAddBoolean("boundaries", "outflow_EMHD_" + bname, false);
+        params.Add("outflow_EMHD_" + bname, outflow_EMHD);
+
         // BOUNDARY TYPES
         // Get the boundary type we specified in kharma
         auto btype = pin->GetString("boundaries", bname);
@@ -248,6 +253,13 @@ void KBoundaries::ApplyBoundary(std::shared_ptr<MeshBlockData<Real>> &rc, IndexD
     const auto btype_name = params.Get<std::string>(bname);
     const auto bdir = BoundaryDirection(bface);
 
+    // If we're pretending to sync primitives, but applying physical bounds
+    // to conserved variables, make sure we're up to date
+    if (pmb->packages.Get<KHARMAPackage>("Driver")->Param<bool>("prims_are_fundamental") &&
+        params.Get<bool>("domain_bounds_on_conserved")) {
+        Flux::BlockPtoU_Send(rc.get(), domain, coarse);
+    }
+
     Flag("Apply "+bname+" boundary: "+btype_name);
     pkg->KBoundaries[bface](rc, coarse);
     EndFlag();
@@ -271,6 +283,26 @@ void KBoundaries::ApplyBoundary(std::shared_ptr<MeshBlockData<Real>> &rc, IndexD
         EndFlag();
     }
 
+    // Allow specifically dP to outflow in otherwise Dirichlet conditions
+    // Only used for viscous_bondi problem
+    // TODO make this more general?
+    if (params.Get<bool>("outflow_EMHD_" + bname)) {
+        Flag("OutflowEMHD_"+bname);
+        auto EMHDg = rc->PackVariables({Metadata::GetUserFlag("EMHDVar"), Metadata::FillGhost});
+        const auto &bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds;
+        const auto &range = (bdir == 1) ? bounds.GetBoundsI(IndexDomain::interior)
+                                : (bdir == 2 ? bounds.GetBoundsJ(IndexDomain::interior)
+                                    : bounds.GetBoundsK(IndexDomain::interior));
+        const int ref = BoundaryIsInner(domain) ? range.s : range.e;
+        pmb->par_for_bndry(
+            "outflow_EMHD", IndexRange{0,EMHDg.GetDim(4)-1}, domain, CC, coarse,
+            KOKKOS_LAMBDA (const int &v, const int &k, const int &j, const int &i) {
+                EMHDg(v, k, j, i) = EMHDg(v, (bdir == 3) ? ref : k, (bdir == 2) ? ref : j, (bdir == 1) ? ref : i);
+            }
+        );
+        EndFlag();
+    }
+
     /*
     * KHARMA is very particular about corner boundaries.
     * In particular, we apply the outflow boundary over ALL X2 & X3.
@@ -320,7 +352,8 @@ void KBoundaries::ApplyBoundary(std::shared_ptr<MeshBlockData<Real>> &rc, IndexD
             Packages::BoundaryPtoUElseUtoP(rc.get(), domain, coarse);
         }
     } else {
-        Packages::BlockUtoP(rc.get(), domain, coarse);
+        // These get applied the same way regardless of driver
+        Packages::BoundaryUtoP(rc.get(), domain, coarse);
     }
 
     EndFlag();
@@ -350,8 +383,6 @@ void KBoundaries::CorrectBPrimitive(std::shared_ptr<MeshBlockData<Real>>& rc, In
 {
     Flag("CorrectBPrimitive");
     std::shared_ptr<MeshBlock> pmb = rc->GetBlockPointer();
-    const Real gam = pmb->packages.Get("GRMHD")->Param<Real>("gamma");
-
     auto B_P = rc->PackVariables(std::vector<std::string>{"prims.B"});
     // Return if no field to correct
     if (B_P.GetDim(4) == 0) return;

kharma/driver/imex_step.cpp

@@ -84,7 +84,7 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta
         if (use_jcon) {
             pmesh->mesh_data.Add("preserve");
             // Above only copies on allocate -- ensure we copy every step
-            Copy<MeshData<Real>>({}, base.get(), pmesh->mesh_data.Get("preserve").get());
+            Copy<MeshData<Real>>({Metadata::Cell}, base.get(), pmesh->mesh_data.Get("preserve").get());
         }
         if (use_implicit) {
             // When solving, we need a temporary copy with any explicit updates,

kharma/emhd/emhd.cpp

@@ -167,6 +167,8 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
 
     // UtoP function specifically for boundary sync (KHARMA must sync cons for AMR) and output
     pkg->BoundaryUtoP = EMHD::BlockUtoP;
+    // If we wanted to apply the domian boundaries to primitive EMHD variables
+    //pkg->DomainBoundaryPtoU = EMHD::BlockPtoU;
 
     // Add all explicit source terms -- implicit terms are called from Implicit::Step
     pkg->AddSource = EMHD::AddSource;

kharma/flux/flux.cpp

@@ -168,12 +168,19 @@ TaskStatus Flux::BlockPtoU_Send(MeshBlockData<Real> *rc, IndexDomain domain, boo
 
     const EMHD::EMHD_parameters& emhd_params = EMHD::GetEMHDParameters(pmb->packages);
 
-    // Pack variables
+    // Pack variables. We never want to run this on the B field
+    using FC = Metadata::FlagCollection;
+    auto cons_flags = FC(Metadata::Conserved, Metadata::Cell, Metadata::GetUserFlag("HD"));
+    if (pmb->packages.AllPackages().count("EMHD"))
+        cons_flags = cons_flags + FC(Metadata::Conserved, Metadata::Cell, Metadata::GetUserFlag("EMHDVar"));
     PackIndexMap prims_map, cons_map;
-    const auto& P = rc->PackVariables({Metadata::GetUserFlag("Primitive")}, prims_map);
-    const auto& U = rc->PackVariables({Metadata::Conserved}, cons_map);
+    const auto& P = rc->PackVariables({Metadata::GetUserFlag("Primitive"), Metadata::Cell}, prims_map);
+    const auto& U = rc->PackVariables(cons_flags, cons_map);
     const VarMap m_u(cons_map, true), m_p(prims_map, false);
 
+    // Return if we're not syncing U & P at all (e.g. edges)
+    if (P.GetDim(4) == 0) return TaskStatus::complete;
+
     auto bounds = coarse ? pmb->c_cellbounds : pmb->cellbounds;
     IndexRange ib = bounds.GetBoundsI(domain);
     IndexRange jb = bounds.GetBoundsJ(domain);

tests/bondi_viscous/check.py

@@ -23,8 +23,8 @@
     outputdir = './'
     kharmadir = '../../'
 
-    NVAR  = 3
-    VARS  = ['rho', 'u', 'dP']
+    NVAR  = 4
+    VARS  = ['rho', 'u', 'dP', 'B']
     RES   = [int(r) for r in sys.argv[1].split(",")]
     LONG  = sys.argv[2]
     SHORT = sys.argv[3]
@@ -44,12 +44,15 @@
         state.params['eta'] = eta
         state.params['tau'] = tau
         dP_check = bondi.compute_dP(mdot, rc, gam, dump.grid, eta, tau)
+        state.cache['dP'] = dP_check
 
         # load code data
         dump = pyharm.load_dump("emhd_2d_{}_end_emhd2d_weno.phdf".format(res))
 
-        rho, uu, dP_tilde = dump['RHO'], dump['UU'], dump['dP']
+        # TODO iterate on names here
+        #rho, uu, dP_tilde = dump['RHO'], dump['UU'], dump['dP']
         #rho, uu = dump['RHO'], dump['UU']
+        rho, uu, dP_tilde, B1 = dump['RHO'], dump['UU'], dump['dP'], dump['B1']
 
         # compute dP
         if dump['emhd/higher_order_terms'] == "true":
@@ -61,17 +64,30 @@
             dP = dP_tilde
 
         # Plot
-        fig = plt.figure(figsize=(6,6))
-        ax = fig.add_subplot(1,1,1)
-        pplt.plot_diff_xz(ax, dump, state, 'rho')
+        for var in ['rho', 'u', 'B1', 'dP']:
+            fig = plt.figure(figsize=(6,6))
+            ax = fig.add_subplot(1,1,1)
+            pplt.plot_diff_xz(ax, dump, state, var)
+            plt.legend()
+            fig.savefig("compare_{}_{}.png".format(var, res))
+            plt.close(fig)
+
+        r_start_ind = 1
+        radius = np.mean(dump.grid['r'][r_start_ind:], axis=(1,2))
+        plt.plot(radius, dP_check[r_start_ind:], label='dP ODE check')
+        plt.plot(radius, np.mean(dump['dP'][r_start_ind:], axis=(1,2)), label='dP0 ODE check')
+        plt.plot(radius, np.mean(state['ucon'][1][r_start_ind:], axis=(1,2)), label='ur')
+        #plt.plot(radius, np.mean(coeff[r_start_ind:], axis=(1,2)), label='coeff')
         plt.legend()
-        fig.savefig("compare_rho_{}.png".format(res))
-        plt.close(fig)
+        plt.savefig('dP_soln_new.png')
+        plt.close()
+
 
         # compute L1 norm
-        L1[r,0] = np.mean(np.fabs(rho[:,0,0] - state['rho'][:,0,0]))
-        L1[r,1] = np.mean(np.fabs(uu[:,0,0]  - state['u'][:,0,0]))
-        L1[r,2] = np.mean(np.fabs(dP[:,0,0]  - dP_check)[1:-1])
+        L1[r,0] = np.mean(np.fabs(rho - state['rho'])[1:-1])
+        L1[r,1] = np.mean(np.fabs(uu  - state['u']))
+        L1[r,2] = np.mean(np.fabs(dP  - dP_check)[1:-1])
+        L1[r,3] = np.mean(np.fabs(B1  - state['B1']))
 
     # MEASURE CONVERGENCE
     L1 = np.array(L1)