Disable B3 reconnection near outer boundary, unstable.

Also add enough density to always recover momenta after Kastaun
inversion.  Should eliminate "bad velocity" failures.

Author: Ben Prather

kharma/b_ct/b_ct.cpp

@@ -73,6 +73,10 @@ std::shared_ptr<KHARMAPackage> B_CT::Initialize(ParameterInput *pin, std::shared
     if (lazy_prolongation && pin->GetString("parthenon/mesh", "refinement") == "adaptive")
         throw std::runtime_error("Cannot use non-divergence-preserving prolongation in AMR!");
 
+    // TODO don't set this unless we're reconnecting at boundaries (can't just check, we load Boundaries pkg later)
+    int reconnection_outer_buffer = pin->GetOrAddInteger("b_field", "reconnection_outer_buffer", 5);
+    params.Add("reconnection_outer_buffer", reconnection_outer_buffer);
+
     // FIELDS
 
     // Flags for B fields on faces.

kharma/b_ct/b_ct_boundaries.cpp

@@ -299,6 +299,8 @@ void B_CT::ReconnectBoundaryB3(MeshBlockData<Real> *rc, IndexDomain domain, cons
 
     const Real gam = pmb->packages.Get("GRMHD")->Param<Real>("gamma");
 
+    const int reconnection_outer_buffer = pmb->packages.Get("B_CT")->Param<int>("reconnection_outer_buffer");
+
     const Floors::Prescription floors = pmb->packages.Get("Floors")->Param<Floors::Prescription>("prescription");
     // Don't be fooled, this function does *not* support/preserve EMHD values
     const EMHD::EMHD_parameters& emhd_params = EMHD::GetEMHDParameters(pmb->packages);
@@ -308,7 +310,7 @@ void B_CT::ReconnectBoundaryB3(MeshBlockData<Real> *rc, IndexDomain domain, cons
     IndexRange3 bi = KDomain::GetRange(rc, IndexDomain::interior, F3, coarse);
     const int jf = (binner) ? bi.js : bi.je; // j index of last zone next to pole
     parthenon::par_for_outer(DEFAULT_OUTER_LOOP_PATTERN, "reduce_B3_" + bname, pmb->exec_space,
-        0, 1, 0, fpack.GetDim(4)-1, b.is, b.ie,
+        0, 1, 0, fpack.GetDim(4)-1, b.is, b.ie - reconnection_outer_buffer,
         KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int &v, const int& i) {
             // Sum the first rank of B3
             double B3_sum = 0.;

kharma/boundaries/boundaries.cpp

@@ -67,19 +67,13 @@ std::shared_ptr<KHARMAPackage> KBoundaries::Initialize(ParameterInput *pin, std:
     // Option to excise a bit at the poles when calculating fluxes
     bool excise_polar_flux = pin->GetOrAddBoolean("boundaries", "excise_polar_flux", false);
     params.Add("excise_polar_flux", excise_polar_flux);
-    if (excise_polar_flux) {
-        // These options are opposites
-        pin->SetBoolean("boundaries", "zero_polar_flux", false);
-        // TODO check whether the user explicitly set these false and yell instead of silent override
-        pin->SetBoolean("boundaries", "reconnect_B3_inner_x2", true);
-        pin->SetBoolean("boundaries", "reconnect_B3_outer_x2", true);
-        // These are themselves unstable, and don't help the wake much
-        //pin->SetBoolean("boundaries", "cancel_T3_inner_x2", true);
-        //pin->SetBoolean("boundaries", "cancel_T3_outer_x2", true);
-    }
     // Otherwise, those fluxes should be zero
-    bool zero_polar_flux = pin->GetOrAddBoolean("boundaries", "zero_polar_flux", spherical);
+    bool zero_polar_flux = pin->GetOrAddBoolean("boundaries", "zero_polar_flux", spherical && !excise_polar_flux);
     params.Add("zero_polar_flux", zero_polar_flux);
+    // Throw an error if both are set
+    if (excise_polar_flux && zero_polar_flux) {
+        throw std::runtime_error("Cannot set both boundaries/excise_polar_flux and boundaries/zero_polar_flux!");
+    }
 
     // Apply physical boundaries to conserved GRMHD variables rho u^r, T^mu_nu
     // Probably inadvisable?
@@ -199,7 +193,8 @@ std::shared_ptr<KHARMAPackage> KBoundaries::Initialize(ParameterInput *pin, std:
             bool clean_face_B = pin->GetOrAddBoolean("boundaries", "clean_face_B_" + bname, (btype == "outflow"));
             params.Add("clean_face_B_"+bname, clean_face_B);
             // Forcibly reconnect field loops that get trapped around the polar boundary
-            bool reconnect_B3 = pin->GetOrAddBoolean("boundaries", "reconnect_B3_" + bname, false);
+            // Needed to keep excised-flux transmitting boundaries stable
+            bool reconnect_B3 = pin->GetOrAddBoolean("boundaries", "reconnect_B3_" + bname, excise_flux);
             params.Add("reconnect_B3_"+bname, reconnect_B3);
 
             // Special EMF averaging.  Allows B3 to "slip" around the pole

kharma/floors/floors.hpp

@@ -102,7 +102,7 @@ class Prescription {
         // Constant sanity limits
         Real rho_min_const, u_min_const;
         // Purely geometric limits
-        Real rho_min_geom, u_min_geom, r_char;
+        Real rho_min_geom, u_min_geom, r_char, floors_switch_r;
         // Dynamic limits on magnetization/temperature
         Real bsq_over_rho_max, bsq_over_u_max, u_over_rho_max;
         // Limit entropy
@@ -113,7 +113,8 @@ class Prescription {
         bool use_r_char, temp_adjust_u, adjust_k;
         // Radius dependent floors?
         bool radius_dependent_floors;
-        Real floors_switch_r;
+        // Add density to respect the gamma ceiling?
+        bool use_rho_to_slow;
 };
 
 inline Prescription MakePrescription(parthenon::ParameterInput *pin, std::string block="floors")
@@ -163,6 +164,8 @@ inline Prescription MakePrescription(parthenon::ParameterInput *pin, std::string
     p.radius_dependent_floors = pin->GetOrAddBoolean("floors", "radius_dependent_floors", false); 
     p.floors_switch_r = pin->GetOrAddReal("floors", "floors_switch_r", 50.);
 
+    p.use_rho_to_slow = pin->GetOrAddBoolean("floors", "use_rho_to_slow", true);
+
     return p;
 }
 

kharma/inverter/inverter.cpp

@@ -203,6 +203,24 @@ inline void BlockPerformInversion(MeshBlockData<Real> *rc, IndexDomain domain, b
                 if (normal_frame_floors) {
                     fflagl |= Floors::determine_floors(G, P, m_p, gam, k, j, i, inverter_floors, inverter_floors_inner,
                         rhoflr_max, uflr_max);
+                    // Add a floor to density which controls wayward velocities
+                    if (inverter_floors.use_rho_to_slow) {
+                        // Calculate necessary rho
+                        const Real rho = std::max(P(m_p.RHO, k, j, i), rhoflr_max);
+                        const Real u = std::max(P(m_p.UU, k, j, i), uflr_max);
+                        const Real rhoh = rho + (gam + 1) * u;
+                        const Real S = m::sqrt(U(m_u.U1, k, j, i) * U(m_u.U1, k, j, i) +
+                                               U(m_u.U2, k, j, i) * U(m_u.U2, k, j, i) +
+                                               U(m_u.U3, k, j, i) * U(m_u.U3, k, j, i));
+                        const Real gamma = GRMHD::lorentz_calc(G, P, m_p, k, j, i, Loci::center);
+                        const Real rhoh_min = S / gamma * gamma * m::sqrt(1 - 1 / (gamma * gamma));
+                        if (rhoh < rhoh_min) {
+                            rhoflr_max = rho * rhoh_min/rhoh;
+                            uflr_max = u * rhoh_min/rhoh;
+                        }
+                    } // TODO try adding at different temperatures/only density or energy
+                    // TODO after adding here, we can use a much simpler solve:
+                    // Spar^2/(rhoh^2 W^4) + Sperp^2/(rhoh W^2 + B^2)^2 + 1/W^2 - 1 = 0
                 } else {
                     // Bare minimum floors for numerics, before applying the rest in user-selected frame
                     rhoflr_max = inverter_floors.rho_min_const;