Reconstruct more fluxes when using excised transmitting poles

kharma/b_cleanup/b_cleanup.cpp

@@ -1,25 +1,25 @@
-/* 
+/*
  *  File: b_cleanup.cpp
- *  
+ *
  *  BSD 3-Clause License
- *  
+ *
  *  Copyright (c) 2020, AFD Group at UIUC
  *  All rights reserved.
- *  
+ *
  *  Redistribution and use in source and binary forms, with or without
  *  modification, are permitted provided that the following conditions are met:
- *  
+ *
  *  1. Redistributions of source code must retain the above copyright notice, this
  *     list of conditions and the following disclaimer.
- *  
+ *
  *  2. Redistributions in binary form must reproduce the above copyright notice,
  *     this list of conditions and the following disclaimer in the documentation
  *     and/or other materials provided with the distribution.
- *  
+ *
  *  3. Neither the name of the copyright holder nor the names of its
  *     contributors may be used to endorse or promote products derived from
  *     this software without specific prior written permission.
- *  
+ *
  *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
@@ -202,7 +202,7 @@ TaskStatus B_Cleanup::CleanupDivergence(std::shared_ptr<MeshData<Real>>& md)
     }
 
     // Calculate/print inital max divB exactly as we would during run
-    double divb_start; 
+    double divb_start;
     if (use_b_ct) {
         divb_start = B_CT::GlobalMaxDivB(md.get());
     } else {
@@ -274,7 +274,7 @@ TaskStatus B_Cleanup::CleanupDivergence(std::shared_ptr<MeshData<Real>>& md)
         // Synchronize to update cons.B's ghost zones
         KHARMADriver::SyncAllBounds(md);
         // Make sure prims.B reflects solution
-        B_CT::MeshUtoP(md.get(), IndexDomain::entire, false);
+        B_CT::MeshUtoP(md.get(), IndexDomain::entire);
         // Recalculate divB max for one last check
         divb_end = B_CT::GlobalMaxDivB(md.get());
     } else {
@@ -283,7 +283,7 @@ TaskStatus B_Cleanup::CleanupDivergence(std::shared_ptr<MeshData<Real>>& md)
         // Synchronize to update cons.B's ghost zones
         KHARMADriver::SyncAllBounds(md);
         // Make sure prims.B reflects solution
-        B_FluxCT::MeshUtoP(md.get(), IndexDomain::entire, false);
+        B_FluxCT::MeshUtoP(md.get(), IndexDomain::entire);
         // Recalculate divB max for one last check
         divb_end = B_FluxCT::GlobalMaxDivB(md.get());
     }

kharma/b_ct/b_ct.cpp

@@ -376,7 +376,7 @@ TaskStatus B_CT::CalculateEMF(MeshData<Real> *md)
             );
         } else if (scheme == "gs05_c" || scheme == "sg07") {
             auto& rho = md->PackVariablesAndFluxes(std::vector<std::string>{"cons.rho"});
-            pmb0->par_for("B_CT_emf_GS05_c", block.s, block.e, b1.ks, b1.ke, b1.js, b1.je, b1.is, b1.ie,
+            pmb0->par_for("B_CT_emf_SG07", block.s, block.e, b1.ks, b1.ke, b1.js, b1.je, b1.is, b1.ie,
                 KOKKOS_LAMBDA (const int &bl, const int &k, const int &j, const int &i) {
                     // Following adapted closely from AthenaK, including clever use of the mass flux for the
                     // sign of the contact mode.

kharma/b_ct/b_ct_boundaries.cpp

@@ -339,7 +339,7 @@ void B_CT::ReconnectBoundaryB3(MeshBlockData<Real> *rc, IndexDomain domain, cons
 
                     // Recover primitive GRMHD variables from our modified U
                     Inverter::u_to_p<Inverter::Type::kastaun>(G, U, m_u, gam, k, jf, i, P, m_p, Loci::center,
-                                                              floors, 25, 1e-12);
+                                                              25, 1e-12, false);
                     // Floor them
                     int fflag = Floors::apply_geo_floors(G, P, m_p, gam, k, jf, i, floors, floors, Loci::center);
                     // Recalculate U on anything we floored

kharma/driver/imex_step.cpp

@@ -81,7 +81,7 @@ TaskCollection KHARMADriver::MakeImExTaskCollection(BlockList_t &blocks, int sta
     // TODO these can now be reduced by including the var lists/flags which actually need to be allocated
     // TODO except the Copy they can be run on step 1 only
     if (stage == 1) {
-        auto &base = pmesh->mesh_data.Get();
+        auto &base = pmesh->mesh_data.Get("base");
         // Fluxes
         pmesh->mesh_data.Add("dUdt");
         for (int i = 1; i < integrator->nstages; i++)

kharma/driver/kharma_step.cpp

@@ -99,7 +99,7 @@ TaskCollection KHARMADriver::MakeDefaultTaskCollection(BlockList_t &blocks, int
     // TODO these can now be reduced by including the var lists/flags which actually need to be allocated
     // TODO except the Copy they can be run on step 1 only
     if (stage == 1) {
-        auto &base = pmesh->mesh_data.Get();
+        auto &base = pmesh->mesh_data.Get("base");
         // Fluxes
         pmesh->mesh_data.Add("dUdt");
         for (int i = 1; i < integrator->nstages; i++)

kharma/floors/floors.cpp

@@ -106,6 +106,10 @@ std::shared_ptr<KHARMAPackage> Floors::Initialize(ParameterInput *pin, std::shar
     else
         params.Add("prescription_inner", MakePrescriptionInner(pin, MakePrescription(pin)), "floors");
 
+    // Sometimes we want the floors package, but don't want to apply anything, for tests
+    bool disable_call = pin->GetOrAddBoolean("floors", "disable_call", false);
+    params.Add("disable_call", disable_call);
+
     // These preserve floor values between the "mark" pass and the actual floor application
     // We need them even if floors are disabled, to apply initial values based on some prescription
     // as a part of problem setup
@@ -120,9 +124,14 @@ std::shared_ptr<KHARMAPackage> Floors::Initialize(ParameterInput *pin, std::shar
     m = Metadata({Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::OneCopy, Metadata::Overridable});
     pkg->AddField("pflag", m);
 
-    // TODO(BSP) THIS IS THE ONLY MeshApplyFloors.  Any others will NOT BE CALLED.
-    // Use BlockApplyFloors in your packages or fix Packages::MeshApplyFloors
-    pkg->MeshApplyFloors = Floors::ApplyGRMHDFloors;
+    // Don't actually call the usual floor function if we're using normal frame w/Kastaun,
+    // floors will be applied during the inversion call.
+    // Also allow manually disabling the call, for testing
+    if (!disable_call && frame != InjectionFrame::normal_kastaun) {
+        // TODO(BSP) THIS IS THE ONLY MeshApplyFloors.  Any others will NOT BE CALLED.
+        // Use BlockApplyFloors in your packages or fix Packages::MeshApplyFloors
+        pkg->MeshApplyFloors = Floors::ApplyGRMHDFloors;
+    }
     pkg->PostStepDiagnosticsMesh = Floors::PostStepDiagnostics;
 
     // List (vector) of HistoryOutputVars that will all be enrolled as output variables
@@ -187,7 +196,7 @@ TaskStatus Floors::ApplyInitialFloors(ParameterInput *pin, MeshBlockData<Real> *
             // Initial floors, so the radius-dependence of floors don't matter that much. 
             int fflag = determine_floors(G, P, m_p, gam, k, j, i, floors, floors, rhoflr_max, uflr_max);
             if (fflag) {
-                apply_floors<InjectionFrame::fluid>(G, P, m_p, gam, k, j, i, rhoflr_max, uflr_max, floors, U, m_u);
+                apply_floors<InjectionFrame::fluid>(G, P, m_p, gam, k, j, i, rhoflr_max, uflr_max, U, m_u);
                 apply_ceilings(G, P, m_p, gam, k, j, i, floors, floors, U, m_u);
                 // P->U for any modified zones
                 Flux::p_to_u_mhd(G, P, m_p, emhd_params, gam, k, j, i, U, m_u, Loci::center);
@@ -223,6 +232,7 @@ TaskStatus Floors::DetermineGRMHDFloors(MeshData<Real> *md, IndexDomain domain,
     pmb0->par_for("determine_floors", block.s, block.e, b.ks, b.ke, b.js, b.je, b.is, b.ie,
         KOKKOS_LAMBDA (const int &b, const int &k, const int &j, const int &i) {
             const auto& G = P.GetCoords(b);
+            // The inverter might have set some floor flags, so we add to that non-destructively
             fflag(b, 0, k, j, i) = static_cast<int>(fflag(b, 0, k, j, i)) |
                                     determine_floors(G, P(b), m_p, gam, k, j, i, floors, floors_inner,
                                                      floor_vals(b, rhofi, k, j, i), floor_vals(b, ufi, k, j, i));
@@ -239,6 +249,12 @@ TaskStatus Floors::ApplyGRMHDFloors(MeshData<Real> *md, IndexDomain domain)
 {
     auto pmesh = md->GetMeshPointer();
     const auto& pars = pmesh->packages.Get("Floors")->AllParams();
+
+    // Determine floors
+    const Floors::Prescription floors = pars.Get<Floors::Prescription>("prescription");
+    const Floors::Prescription floors_inner = pars.Get<Floors::Prescription>("prescription_inner");
+    DetermineGRMHDFloors(md, domain, floors, floors_inner);
+
     if (pars.Get<InjectionFrame>("frame") == InjectionFrame::normal_kastaun) {
         return ApplyFloorsInFrame<InjectionFrame::normal_kastaun>(md, domain);
     } else if (pars.Get<InjectionFrame>("frame") == InjectionFrame::normal_onedw) {

kharma/floors/floors_functions.hpp

@@ -163,7 +163,6 @@ KOKKOS_INLINE_FUNCTION int determine_floors(const GRCoordinates& G, const Variab
 
 #define FLOOR_ONE_ARGS const GRCoordinates& G, const VariablePack<Real>& P, const VarMap& m_p, const Real& gam, \
                         const int& k, const int& j, const int& i, const Real& rhoflr_max, const Real& uflr_max, \
-                        const Floors::Prescription& floors, \
                         const VariablePack<Real>& U, const VarMap& m_u
 
 /**
@@ -266,59 +265,57 @@ KOKKOS_INLINE_FUNCTION int apply_floors<InjectionFrame::normal_onedw>(FLOOR_ONE_
 {
     // Add the material in the normal observer frame.
     // 1. Calculate how much material we're adding.
-    // This is an estimate, as it's what we'd have to do in fluid frame
+    // By using the existing velocities for the rho*u^0 and T^0_0 contributions,
+    // We produce a guaranteed overestimate (since NOF floors will slow the material,
+    // reducing the Lorentz factor)
     const Real rho_add    = m::max(0., rhoflr_max - P(m_p.RHO, k, j, i));
     const Real u_add      = m::max(0., uflr_max - P(m_p.UU, k, j, i));
-    const Real uvec[NVEC] = {0}, B[NVEC] = {0};
+    const Real uvec[NVEC] = {P(m_p.U1, k, j, i), P(m_p.U2, k, j, i), P(m_p.U3, k, j, i)};
+    const Real B[NVEC] = {0.};
 
     // 2. Calculate the increase in conserved mass/energy corresponding to the new material.
     Real rho_ut, T[GR_DIM];
     GRMHD::p_to_u_mhd(G, rho_add, u_add, uvec, B, gam, k, j, i, rho_ut, T, Loci::center);
 
     // 3. Add new conserved mass/energy to the current "conserved" state.
-    // Also add to the local primitives as a guess
+    U(m_u.RHO, k, j, i) += rho_ut;
+    U(m_u.UU, k, j, i)  += T[0];
+    // Also add to the local primitives to produce a better guess
     P(m_p.RHO, k, j, i) += rho_add;
     P(m_p.UU, k, j, i)  += u_add;
-    // Add any velocity here
-    U(m_u.RHO, k, j, i) += rho_ut;
-    U(m_u.UU, k, j, i)  += T[0]; // Note that m_u.U1 != m_u.UU + 1 necessarily
-    U(m_u.U1, k, j, i)  += T[1];
-    U(m_u.U2, k, j, i)  += T[2];
-    U(m_u.U3, k, j, i)  += T[3];
 
-    // Recover primitive variables from conserved versions.  Use Kastaun with safe parameters so we don't fail often
+    // Recover primitive variables from conserved versions
     return Inverter::u_to_p<Inverter::Type::onedw>(G, U, m_u, gam, k, j, i, P, m_p, Loci::center,
-                                                     floors, 8, 1e-8);
+                                                    8, 1e-8, false);
 }
 
 template<>
 KOKKOS_INLINE_FUNCTION int apply_floors<InjectionFrame::normal_kastaun>(FLOOR_ONE_ARGS)
 {
     // Add the material in the normal observer frame.
     // 1. Calculate how much material we're adding.
-    // This is an estimate, as it's what we'd have to do in fluid frame
+    // By using the existing velocities for the rho*u^0 and T^0_0 contributions,
+    // We produce a guaranteed overestimate (since NOF floors will slow the material,
+    // reducing the Lorentz factor)
     const Real rho_add    = m::max(0., rhoflr_max - P(m_p.RHO, k, j, i));
     const Real u_add      = m::max(0., uflr_max - P(m_p.UU, k, j, i));
-    const Real uvec[NVEC] = {0}, B[NVEC] = {0};
+    // TODO turn this into an option maybe? Or maybe set rhoflr/uflr using it
+    //const Real uvec[NVEC] = {P(m_p.U1, k, j, i), P(m_p.U2, k, j, i), P(m_p.U3, k, j, i)};
+    const Real uvec[NVEC] = {0.};
+    const Real B[NVEC] = {0.};
 
     // 2. Calculate the increase in conserved mass/energy corresponding to the new material.
     Real rho_ut, T[GR_DIM];
     GRMHD::p_to_u_mhd(G, rho_add, u_add, uvec, B, gam, k, j, i, rho_ut, T, Loci::center);
 
     // 3. Add new conserved mass/energy to the current "conserved" state.
-    // Also add to the local primitives as a guess
-    P(m_p.RHO, k, j, i) += rho_add;
-    P(m_p.UU, k, j, i)  += u_add;
-    // Add any velocity here
+    // (no need to modify the guess for Kastaun, esp once we sync mu)
     U(m_u.RHO, k, j, i) += rho_ut;
-    U(m_u.UU, k, j, i)  += T[0]; // Note that m_u.U1 != m_u.UU + 1 necessarily
-    U(m_u.U1, k, j, i)  += T[1];
-    U(m_u.U2, k, j, i)  += T[2];
-    U(m_u.U3, k, j, i)  += T[3];
+    U(m_u.UU, k, j, i)  += T[0];
 
-    // Recover primitive variables from conserved versions.  Use Kastaun with safe parameters so we don't fail often
+    // Recover new primitive variables.  Use Kastaun with safe parameters so we don't fail often
     return Inverter::u_to_p<Inverter::Type::kastaun>(G, U, m_u, gam, k, j, i, P, m_p, Loci::center,
-                                                     floors, 25, 1e-12);
+                                                     25, 1e-12, true);
 }
 
 /**

kharma/floors/floors_impl.hpp

@@ -69,9 +69,6 @@ TaskStatus ApplyFloorsInFrame(MeshData<Real> *md, IndexDomain domain)
     const Floors::Prescription floors = pmb0->packages.Get("Floors")->Param<Floors::Prescription>("prescription");
     const Floors::Prescription floors_inner = pmb0->packages.Get("Floors")->Param<Floors::Prescription>("prescription_inner");
 
-    // Determine floors
-    DetermineGRMHDFloors(md, domain, floors, floors_inner);
-
     const IndexRange3 b = KDomain::GetRange(md, domain);
     const IndexRange block = IndexRange{0, P.GetDim(5) - 1};
     pmb0->par_for("apply_floors", block.s, block.e, b.ks, b.ke, b.js, b.je, b.is, b.ie,
@@ -85,13 +82,13 @@ TaskStatus ApplyFloorsInFrame(MeshData<Real> *md, IndexDomain domain)
                     if (G.r(k, j, i) > switch_r) {
                         pflag_l = apply_floors<InjectionFrame::fluid>(G, P(b), m_p, gam, k, j, i,
                                             floor_vals(b, rhofi, k, j, i), floor_vals(b, ufi, k, j, i),
-                                            floors, U(b), m_u);
+                                            U(b), m_u);
                     } else {
                         // TODO should mixed frames respect Kastaun vs 1Dw?
                         // Since no prior simulations use mixed frames thus requiring back-compat, I said no
                         pflag_l = apply_floors<InjectionFrame::normal_kastaun>(G, P(b), m_p, gam, k, j, i,
                                             floor_vals(b, rhofi, k, j, i), floor_vals(b, ufi, k, j, i),
-                                            floors, U(b), m_u);
+                                            U(b), m_u);
                     }
                 } else if (frame == InjectionFrame::mixed_normal_drift) {
                     FourVectors Dtmp;
@@ -101,16 +98,16 @@ TaskStatus ApplyFloorsInFrame(MeshData<Real> *md, IndexDomain domain)
                     if (mag_switch < switch_beta) {
                         pflag_l = apply_floors<InjectionFrame::drift>(G, P(b), m_p, gam, k, j, i,
                                             floor_vals(b, rhofi, k, j, i), floor_vals(b, ufi, k, j, i),
-                                            floors, U(b), m_u);
+                                            U(b), m_u);
                     } else {
                         pflag_l = apply_floors<InjectionFrame::normal_kastaun>(G, P(b), m_p, gam, k, j, i,
                                             floor_vals(b, rhofi, k, j, i), floor_vals(b, ufi, k, j, i),
-                                            floors, U(b), m_u);
+                                            U(b), m_u);
                     }
                 } else {
                     pflag_l = apply_floors<frame>(G, P(b), m_p, gam, k, j, i,
                                         floor_vals(b, rhofi, k, j, i), floor_vals(b, ufi, k, j, i),
-                                        floors, U(b), m_u);
+                                        U(b), m_u);
                 }
 
                 // Record the pflag if nonzero, that is, if *either* the initial inversion or

kharma/flux/flux.cpp

@@ -180,6 +180,8 @@ std::shared_ptr<KHARMAPackage> Flux::Initialize(ParameterInput *pin, std::shared
     params.Add("use_fofc", use_fofc);
 
     if (use_fofc) {
+        // TODO check floors are enabled!  We can't do fofc without them
+
         // FOFC-specific options
         bool use_glf = pin->GetOrAddBoolean("fofc", "use_glf", false);
         params.Add("fofc_use_glf", use_glf);