Merge branch 'fix/papercuts' into parthenon-bump

Author: Ben Prather

external/patches/parthenon-remove-hangdetect.patch

@@ -0,0 +1,15 @@
+diff --git a/src/utils/communication_buffer.hpp b/src/utils/communication_buffer.hpp
+index edfc04f1..b2c44091 100644
+--- a/src/utils/communication_buffer.hpp
++++ b/src/utils/communication_buffer.hpp
+@@ -327,8 +327,8 @@ bool CommBuffer<T>::TryReceive() noexcept {
+       *comm_type_ == BuffCommType::sparse_receiver) {
+ #ifdef MPI_PARALLEL
+     (*nrecv_tries_)++;
+-    PARTHENON_REQUIRE(*nrecv_tries_ < 1e8,
+-                      "MPI probably hanging after 1e8 receive tries.");
++    //PARTHENON_REQUIRE(*nrecv_tries_ < 1e8,
++    //                  "MPI probably hanging after 1e8 receive tries.");
+ 
+     TryStartReceive();
+ 

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
@@ -210,7 +210,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 {
@@ -282,7 +282,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 {
@@ -291,7 +291,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

@@ -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", 10);
+    params.Add("reconnection_outer_buffer", reconnection_outer_buffer);
+
     // FIELDS
 
     // Flags for B fields on faces.
@@ -385,7 +389,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.
@@ -528,6 +532,11 @@ TaskStatus B_CT::DerefinePoles(MeshData<Real> *md)
                 const int offset = (binner) ? 1 : -1; // offset to read the physical face values
                 const int point_out = offset; // if F2 B field at j_f + offset face is positive when pointing out of the cell, +1.
 
+                // Should we allow flux through the pole?
+                auto &bpars = pmesh->packages.Get("Boundaries")->AllParams();
+                const bool allow_flux = binner ? bpars.Get<bool>("excise_flux_inner_x2"):
+                                                 bpars.Get<bool>("excise_flux_outer_x2");
+
                 // F1 average
                 pmb->par_for("B_CT_derefine_poles_avg_F1", bCC.ks, bCC.ke, j_p.s, j_p.e, bF1.is, bF1.ie,
                     KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
@@ -554,7 +563,7 @@ TaskStatus B_CT::DerefinePoles(MeshData<Real> *md)
                         // starting k-index of the coarse cell
                         const int k_start = k - k_fine;
 
-                        if (j == j_f) {
+                        if (!allow_flux && j == j_f) {
                             // The fine cells have 0 fluxes through the physical-ghost boundaries.
                             B_avg(F2, 0, k, j, i) = 0.;
                         } else { // average the fine cells

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 + 6, 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.;
@@ -339,8 +341,9 @@ 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,
-                                                              25, 1e-12);
+                                                              25, 1e-12, false);
                     // Floor them
+                    // TODO THIS IS IN FLUID FRAME
                     int fflag = Floors::apply_geo_floors(G, P, m_p, gam, k, jf, i, floors, floors, Loci::center);
                     // Recalculate U on anything we floored
                     if (fflag)

kharma/boundaries/boundaries.cpp

@@ -67,12 +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 *completely* incompatible
-        pin->SetBoolean("boundaries", "zero_polar_flux", false);
-    }
     // 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?
@@ -174,8 +175,7 @@ std::shared_ptr<KHARMAPackage> KBoundaries::Initialize(ParameterInput *pin, std:
         // Ensure fluxes through the zero-size face at the pole are zero
         bool zero_flux = pin->GetOrAddBoolean("boundaries", "zero_flux_" + bname, zero_polar_flux && bdir == X2DIR);
         params.Add("zero_flux_" + bname, zero_flux);
-
-        // Ensure fluxes through the zero-size face at the pole are zero
+        // OR allow them via faux-excision
         bool excise_flux = pin->GetOrAddBoolean("boundaries", "excise_flux_" + bname, excise_polar_flux && bdir == X2DIR);
         params.Add("excise_flux_" + bname, excise_flux);
 
@@ -194,11 +194,12 @@ 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 B slippage, e.g. around pole for transmitting conditions
-            // Useful for certain dirichlet conditions e.g. multizone
+            // Special EMF averaging.  Allows B3 to "slip" around the pole
+            // Also useful to allow coherent motion even with Dirichlet boundaries, for e.g. multizone
             bool average_EMF = pin->GetOrAddBoolean("boundaries", "average_EMF_" + bname, (btype == "transmitting"));
             params.Add("average_EMF_"+bname, average_EMF);
             // Otherwise, always zero EMFs to prevent B field escaping the domain in polar/dirichlet bounds
@@ -705,6 +706,10 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
                 if (pmb->boundary_flag[bface] == BoundaryFlag::user) {
                     if (bdir != 2) throw std::runtime_error("Excised polar fluxes only fully implemented in X2!");
 
+                    // Pack w/B to match indices with the `Flux.X` below
+                    // We won't *update* B field though
+                    auto &F = rc->PackVariablesAndFluxes({Metadata::WithFluxes}, cons_map);
+
                     // Going to need the primitive vars
                     PackIndexMap prims_map;
                     std::vector<MetadataFlag> prims_flags = {Metadata::GetUserFlag("Primitive"), Metadata::Cell};
@@ -771,7 +776,8 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
 
                             // Use LLF flux
                             PLOOP {
-                                F.flux(dir, ip, k, j, i) = Flux::llf(Fl_all(ip, k, j, i), Fr_all(ip, k, j, i),
+                                if (ip != m_u.B1 && ip != m_u.B2 && ip != m_u.B3)
+                                    F.flux(dir, ip, k, j, i) = Flux::llf(Fl_all(ip, k, j, i), Fr_all(ip, k, j, i),
                                                                     cmax(dir-1, k, j, i), cmin(dir-1, k, j, i),
                                                                     Ul_all(ip, k, j, i), Ur_all(ip, k, j, i)) * 0.5;
                             }
@@ -818,21 +824,19 @@ TaskStatus KBoundaries::FixFlux(MeshData<Real> *md)
 
                             // Use LLF flux
                             PLOOP {
-                                F.flux(bdir, ip, k, j, i) = Flux::llf(Fl_all(ip, k, j, i), Fr_all(ip, k, j, i),
-                                                                    cmax(bdir-1, k, j, i), cmin(bdir-1, k, j, i),
-                                                                    Ul_all(ip, k, j, i), Ur_all(ip, k, j, i));
-                                // Reduce the X1 flux in a semi-consistent way
-                                const int jc = (binner) ? j_cell + 1 : j_cell;
-                                F.flux(X1DIR, ip, k, j_cell, i) *= 0.5
-                                    * (G.gdet(Loci::face1, j_cell, i) + G.gdet(Loci::corner, jc, i)) / 2 / G.gdet(Loci::face1, j_cell, i);
-                                // This is also a decent guess, but less accurate than recalculating as above
-                                // F.flux(X3DIR, ip, k, j_cell, i) *= 0.5
-                                //     * G.gdet(loc, j_cell, i) / G.gdet(Loci::center, j_cell, i);
+                                if (ip != m_u.B1 && ip != m_u.B2 && ip != m_u.B3) {
+                                    F.flux(bdir, ip, k, j, i) = Flux::llf(Fl_all(ip, k, j, i), Fr_all(ip, k, j, i),
+                                                                        cmax(bdir-1, k, j, i), cmin(bdir-1, k, j, i),
+                                                                        Ul_all(ip, k, j, i), Ur_all(ip, k, j, i));
+                                    // Reduce the X1 flux in a semi-consistent way
+                                    const int jc = (binner) ? j_cell + 1 : j_cell;
+                                    F.flux(X1DIR, ip, k, j_cell, i) *= 0.5
+                                        * (G.gdet(Loci::face1, j_cell, i) + G.gdet(Loci::corner, jc, i)) / 2 / G.gdet(Loci::face1, j_cell, i);
+                                    // This is also a decent guess, but less accurate than recalculating as above
+                                    // F.flux(X3DIR, ip, k, j_cell, i) *= 0.5
+                                    //     * G.gdet(loc, j_cell, i) / G.gdet(Loci::center, j_cell, i);
+                                }
                             }
-
-                            // Account for the half-size in the timestep later
-                            cmax(bdir-1, k, j, i) *= 2;
-                            cmin(bdir-1, k, j, i) *= 2;
                         }
                     );
                     // Then average to make absolutely sure fluxes match
@@ -904,7 +908,9 @@ void KBoundaries::AddSource(MeshData<Real> *md, MeshData<Real> *mdudt, IndexDoma
                         b.ks = b.ke = (binner) ? bi.ks : bi.ke;
                     }
 
-                    auto &dUdt = rc->PackVariables({Metadata::WithFluxes});
+                    // The magnetic field is probably defined at faces; even if it's defined in cells,
+                    // we shouldn't be monkeying with it.  We just do not adjust it here.
+                    auto &dUdt = rc->PackVariables({Metadata::GetUserFlag("HD"), Metadata::WithFluxes});
                     const auto& G = pmb->coords;
                     const Loci loc = (binner) ? Loci::outer_half : Loci::inner_half;
 

kharma/driver/imex_step.cpp

@@ -83,7 +83,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", base);
         for (int i = 1; i < integrator->nstages; i++)

kharma/driver/kharma_driver.cpp

@@ -378,7 +378,9 @@ TaskID KHARMADriver::AddStateUpdate(TaskID& t_start, TaskList& tl, MeshData<Real
     auto t_copy_prims = t_update;
     auto pmb0  = md_full_step_init->GetBlockData(0)->GetBlockPointer();
     auto& pkgs = pmb0->packages.AllPackages();
-    if (!pkgs.at("GRMHD")->Param<bool>("implicit")) {
+    // If we're explicitly evolving, UtoP needs a guess (except Kastaun inverter)
+    if (!pkgs.at("GRMHD")->Param<bool>("implicit") &&
+        pkgs.at("Inverter")->Param<Inverter::Type>("inverter_type") != Inverter::Type::kastaun) {
         t_copy_prims = tl.AddTask(t_start, Copy<MeshData<Real>>,
                                     std::vector<MetadataFlag>({Metadata::GetUserFlag("HD"), Metadata::GetUserFlag("Primitive")}),
                                     md_sub_step_init, md_update);

kharma/driver/simple_step.cpp

@@ -44,6 +44,7 @@ TaskCollection KHARMADriver::MakeSimpleTaskCollection(BlockList_t &blocks, int s
     auto& flux_pkg = pkgs.at("Fluxes")->AllParams();
     auto& inverter_pkg = pkgs.at("Inverter")->AllParams();
     PARTHENON_REQUIRE(!pkgs.count("B_CT"), "Face-centered B not compatible with simple driver!");
+    // TODO require we're not using B_Cleanup as transport too
     PARTHENON_REQUIRE(!pkgs.count("Electrons"), "Electrons not compatible with simple driver!");
     PARTHENON_REQUIRE(!flux_pkg.Get<bool>("use_fofc"), "Flux corrections not compatible with simple driver!");
     PARTHENON_REQUIRE(inverter_pkg.Get<Inverter::Type>("inverter_type") == Inverter::Type::kastaun,

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;
 }