this should be a working branch now

kharma/flux/flux.cpp

@@ -81,10 +81,9 @@ std::shared_ptr<KHARMAPackage> Flux::Initialize(ParameterInput *pin, std::shared
         throw std::runtime_error("Cannot enable lowered reconstruction on edges and poles!");
     if ((lower_edges || lower_poles) && recon != "weno5")
         throw std::runtime_error("Lowered reconstructions can only be enabled with weno5!");
-    bool use_ismr = packages.AllPackages().count("ISMR");
 
     int stencil = 0;
-    if (use_ismr) {
+    if (pin->GetOrAddBoolean("ismr", "on", false)) {
         // Override for ismr
         params.Add("recon", KReconstruction::Type::weno5_ismr);
         pin->SetString("flux", "reconstruction", "weno5_ismr");

kharma/flux/reconstruction.hpp

@@ -772,11 +772,7 @@ KOKKOS_INLINE_FUNCTION void ReconstructRowIsmr<Type::weno5_ismr, X1DIR>(partheno
                                         const int& k, const int& j, const int& is_l, const int& ie_l, 
                                         ScratchPad2D<Real> ql, ScratchPad2D<Real> qr, uint ismr_nlevels, int ng)
 {
-    if (j < ng + ismr_nlevels || j > P.GetDim(2) - 1 - ng - ismr_nlevels) {
-        KReconstruction::ReconstructX1<Type::linear_mc>(member, k, j, is_l, ie_l, P, ql, qr);
-    } else {
-        KReconstruction::ReconstructX1<Type::weno5>(member, k, j, is_l, ie_l, P, ql, qr);
-    }
+    KReconstruction::ReconstructX1<Type::weno5>(member, k, j, is_l, ie_l, P, ql, qr);
 }
 template <>
 KOKKOS_INLINE_FUNCTION void ReconstructRowIsmr<Type::weno5_ismr, X2DIR>(parthenon::team_mbr_t& member,

kharma/grmhd/grmhd.cpp

@@ -185,23 +185,6 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
     // No magnetic fields here. KHARMA should operate fine in GRHD without them,
     // so they are allocated only by B field packages.
 
-    // INTERNAL SMR
-    // TODO own package
-    // internal SMR :Added by Hyerin (03/07/24)
-	const bool ismr_poles = pmesh->packages.AllPackages().count("ISMR");
-    const uint ismr_nlevels = (ismr_poles) ? pmesh->packages.Get("ISMR")->Param<uint>("nlevels") : 0;
-    if (ismr_poles) {
-        uint ismr_nlevels = (uint) pin->GetOrAddInteger("GRMHD", "ismr_nlevels", 1);
-        params.Add("ismr_nlevels", ismr_nlevels);
-
-        // ISMR caches: not evolved, immediately copied to fluid state after averaging
-        m = Metadata({Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::OneCopy});
-        pkg->AddField("ismr_rho_avg", m);
-        pkg->AddField("ismr_u_avg", m);
-        m = Metadata({Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::OneCopy, Metadata::Vector}, s_vector);
-        pkg->AddField("ismr_uvec_avg", m);
-    }
-
     // A KHARMAPackage also contains quite a few "callbacks," or functions called at
     // specific points in a step if the package is loaded.
     // Generally, see the headers for function descriptions.
@@ -249,8 +232,8 @@ Real EstimateTimestep(MeshBlockData<Real> *rc)
     const auto& driver_pars = pmb->packages.Get("Driver")->AllParams();
 
     // Added by Hyerin (03/07/24)
-	const bool ismr_poles = pmesh->packages.AllPackages().count("ISMR");
-    const uint ismr_nlevels = (ismr_poles) ? pmesh->packages.Get("ISMR")->Param<uint>("nlevels") : 0;
+	const bool ismr_poles = pmb->packages.AllPackages().count("ISMR");
+    const uint ismr_nlevels = (ismr_poles) ? pmb->packages.Get("ISMR")->Param<uint>("nlevels") : 0;
     const bool polar_inner_x2 = pmb->boundary_flag[BoundaryFace::inner_x2] == BoundaryFlag::user;
     const bool polar_outer_x2 = pmb->boundary_flag[BoundaryFace::outer_x2] == BoundaryFlag::user;
 

kharma/ismr/ismr.hpp

@@ -0,0 +1,69 @@
+/* 
+ *  File: ismr.hpp
+ *  
+ *  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
+ *  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+ *  FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ *  DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ *  SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ *  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+ *  OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+ *  OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+#pragma once
+
+#include "decs.hpp"
+#include "types.hpp"
+
+#include <parthenon/parthenon.hpp>
+
+/**
+ * This package implements internal static mesh de-refinement by averaging variables next to the coordinate
+ * poles, creating effective "larger zones" and allowing increases to the timestep, without any modification
+ * to the existing data structures or block layout.
+ * Currently it is limited to spherical coordinates: use Parthenon's block-based refinement otherwise.
+ * Also note that the averaging operation corresponds to a first-order refinement/derefinement, possibly
+ * compromising overall second-order convergence in the very near-polar region.
+ * 
+ * The operator averages variables only in the phi-direction nearing the pole, mapping 1:2 zones
+ * for each of several levels.  In this way, zones near the pole can be much wider in phi without losing
+ * resolution in r, theta.
+ * 
+ * Idea is taken from Matthew Liska/H-AMR (Liska+ 2022)
+ * Implementation credit Hyerin Cho, please cite Cho+ 2024b (in prep) with description and first
+ * use of this implementation.
+ */
+namespace ISMR {
+
+/**
+ * Initialize ISMR parameters
+ */
+std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<Packages_t>& packages);
+
+/**
+ * Derefinement operation for fluid/cell-centered variables 
+ */
+TaskStatus DerefinePoles(MeshData<Real> *md);
+
+}

kharma/kharma.cpp

@@ -51,6 +51,7 @@
 #include "electrons.hpp"
 #include "implicit.hpp"
 #include "inverter.hpp"
+#include "ismr.hpp"
 #include "floors.hpp"
 #include "grmhd.hpp"
 #include "reductions.hpp"
@@ -419,6 +420,11 @@ Packages_t KHARMA::ProcessPackages(std::unique_ptr<ParameterInput> &pin)
     // Flux temporaries must be full size
     KHARMA::AddPackage(packages, Flux::Initialize, pin.get());
 
+	// ISMR temporaries must be full size
+    if (pin->GetOrAddBoolean("ismr", "on", false)) {
+        KHARMA::AddPackage(packages, ISMR::Initialize, pin.get());
+    }
+
     // And any dirichlet/constant boundaries
     // TODO avoid init if Parthenon will be handling all boundaries?
     KHARMA::AddPackage(packages, KBoundaries::Initialize, pin.get());