Merge pull request #144 from AFD-Illinois/fix/excised-timestep

c-prather · web-flow · commit b15adfc5fb49 · 2025-09-17T15:03:12.000-06:00
Fix the timestep when `excise_polar_flux=true`
diff --git a/kharma/b_cleanup/b_cleanup.cpp b/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());
     }
diff --git a/kharma/b_ct/b_ct.cpp b/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.
@@ -519,6 +519,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) {
@@ -545,7 +550,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
diff --git a/kharma/b_ct/b_ct_boundaries.cpp b/kharma/b_ct/b_ct_boundaries.cpp
@@ -341,6 +341,7 @@ void B_CT::ReconnectBoundaryB3(MeshBlockData<Real> *rc, IndexDomain domain, cons
                     Inverter::u_to_p<Inverter::Type::kastaun>(G, U, m_u, gam, k, jf, i, P, m_p, Loci::center,
                                                               floors, 25, 1e-12);
                     // 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)
diff --git a/kharma/boundaries/boundaries.cpp b/kharma/boundaries/boundaries.cpp
@@ -173,8 +173,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);
 
@@ -196,8 +195,8 @@ std::shared_ptr<KHARMAPackage> KBoundaries::Initialize(ParameterInput *pin, std:
             bool reconnect_B3 = pin->GetOrAddBoolean("boundaries", "reconnect_B3_" + bname, false);
             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
@@ -704,6 +703,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};
@@ -770,7 +773,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;
                             }
@@ -817,21 +821,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
@@ -903,7 +905,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;
 
diff --git a/kharma/grmhd/grmhd.cpp b/kharma/grmhd/grmhd.cpp
@@ -39,6 +39,7 @@
 // TODO eliminate when Parthenon gets reduction types
 #include "Kokkos_Core.hpp"
 
+#include "b_ct.hpp"
 #include "boundaries.hpp"
 #include "current.hpp"
 #include "floors.hpp"
@@ -257,6 +258,7 @@ Real EstimateTimestep(MeshData<Real> *md)
     const auto& grmhd_pars = pmesh->packages.Get("GRMHD")->AllParams();
 
     // If we have to recompute ctop anywhere, we do it now
+    // TODO only call if one of the reconnect_ or excised polar is enabled
     UpdateAveragedCtop(md);
 
     // Other things we might have to return (light-crossing, pre-set timestep, etc.)
@@ -312,12 +314,17 @@ Real EstimateTimestep(MeshData<Real> *md)
         const auto& cmax  = rc->PackVariables(std::vector<std::string>{"Flux.cmax"});
         const auto& cmin  = rc->PackVariables(std::vector<std::string>{"Flux.cmin"});
 
+        auto& boundaries = pmesh->packages.Get<KHARMAPackage>("Boundaries")->AllParams();
+        const bool excise_inner_x2 = boundaries.Get<bool>("excise_flux_inner_x2");
+        const bool excise_outer_x2 = boundaries.Get<bool>("excise_flux_outer_x2");
+
         double block_min_ndt = 0.;
         pmb->par_reduce("ndt_min", b.ks, b.ke, b.js, b.je, b.is, b.ie,
             KOKKOS_LAMBDA (const int k, const int j, const int i,
                         double &local_result) {
                 const auto& G = cmax.GetCoords();
                 int ismr_factor = 1;
+                double excise_factor = 1.0;
                 double courant_limit = 1.0;
                 if (ismr_poles && polar_inner_x2 && j < (b.js + ismr_nlevels)) {
                     ismr_factor = m::pow(2, ismr_nlevels - (j - b.js));
@@ -328,8 +335,14 @@ Real EstimateTimestep(MeshData<Real> *md)
                     courant_limit = 0.5;
                 }
 
+                if (excise_inner_x2 && polar_inner_x2 && j == b.js) {
+                    excise_factor = 0.5;
+                } else if (excise_outer_x2 && polar_outer_x2 && j == b.je) {
+                    excise_factor = 0.5;
+                }
+
                 double ndt_zone = courant_limit / (1 / (G.Dxc<1>(i) /  m::max(cmax(V1, k, j, i), cmin(V1, k, j, i))) +
-                                    1 / (G.Dxc<2>(j) /  m::max(cmax(V2, k, j, i), cmin(V2, k, j, i))) +
+                                    1 / (G.Dxc<2>(j) * excise_factor /  m::max(cmax(V2, k, j, i), cmin(V2, k, j, i))) +
                                     1 / (G.Dxc<3>(k) * ismr_factor /  m::max(cmax(V3, k, j, i), cmin(V3, k, j, i))));
 
                 if (!m::isnan(ndt_zone) && (ndt_zone < local_result)) {
@@ -642,6 +655,8 @@ void CancelBoundaryT3(MeshBlockData<Real> *rc, IndexDomain domain, bool coarse)
 void UpdateAveragedCtop(MeshData<Real> *md)
 {
     auto pmesh = md->GetMeshPointer();
+    if (pmesh->packages.AllPackages().count("B_CT"))
+        B_CT::MeshUtoP(md, IndexDomain::interior);
     auto& params = pmesh->packages.Get<KHARMAPackage>("Boundaries")->AllParams();
     for (auto &pmb : pmesh->block_list) {
         auto &rc = pmb->meshblock_data.Get(md->StageName());
@@ -661,7 +676,8 @@ void UpdateAveragedCtop(MeshData<Real> *md)
             // If we've modified values on the pole...
             if (params.Get<bool>("cancel_T3_" + bname) ||
                 params.Get<bool>("cancel_U3_" + bname) ||
-                b3_is_reconnected) {
+                b3_is_reconnected ||
+                params.Get<bool>("excise_flux_" + bname)) {
                 // ...and if this face of the block corresponds to a global boundary...
                 if (pmb->boundary_flag[bface] == BoundaryFlag::user) {
                     PackIndexMap prims_map, cons_map;
@@ -699,10 +715,6 @@ void UpdateAveragedCtop(MeshData<Real> *md)
                             Flux::vchar_global(G, P, m_p, Dtmp, gam, emhd_params, k, jf, i, Loci::center, X3DIR,
                                         cmax(V3, k, jf, i), cmin_minus);
                             cmin(V3, k, jf, i) = -cmin_minus;
-                            if (half_cells) {
-                                cmin(bdir-1, k, jf, i) *= 0.5;
-                                cmax(bdir-1, k, jf, i) *= 0.5;
-                            }
                         }
                     );
                 }
diff --git a/kharma/prob/problem.cpp b/kharma/prob/problem.cpp
@@ -154,12 +154,11 @@ void KHARMA::ProblemGenerator(MeshBlock *pmb, ParameterInput *pin)
     // If needed, they are applied within the problem-specific call.
     // See InitializeFMTorus in fm_torus.cpp for the details for torus problems.
 
-    // Note we no longer call PtoU here either, as GRMHD variables' PtoU requires
-    // the magnetic field, which is added in PostInitialize, after all blocks
-    // are filled with other variables (it can be related to density averages which
-    // require correct ghost zones)
-    // If the B field will depend on the conserved variables (for some reason?)
-    // they must be computed by the particular problem.
+    // This is a temporary PtoU call.  It will underestimate T^0_0 for magnetized
+    // problems, since the magnetic field is not yet initialized.
+    // However, the polar mitigations expect P,U in a consistent state,
+    // so we have to give them something.
+    Flux::BlockPtoU(rc.get(), IndexDomain::interior);
 
     EndFlag();
 }
