diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index 1475ae8e1db29..5ac54383ff4a9 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -22,5 +22,6 @@ add_subdirectory(particle_leapfrog)
 add_subdirectory(particle_tracers)
 add_subdirectory(poisson)
 add_subdirectory(poisson_gmg)
+add_subdirectory(linear_solvers)
 add_subdirectory(diffusion)
 add_subdirectory(sparse_advection)
diff --git a/example/linear_solvers/CMakeLists.txt b/example/linear_solvers/CMakeLists.txt
new file mode 100644
index 0000000000000..88913cc5db881
--- /dev/null
+++ b/example/linear_solvers/CMakeLists.txt
@@ -0,0 +1,36 @@
+#=========================================================================================
+# (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved.
+#
+# This program was produced under U.S. Government contract 89233218CNA000001 for Los
+# Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+# for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+# in the program are reserved by Triad National Security, LLC, and the U.S. Department
+# of Energy/National Nuclear Security Administration. The Government is granted for
+# itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+# license in this material to reproduce, prepare derivative works, distribute copies to
+# the public, perform publicly and display publicly, and to permit others to do so.
+#=========================================================================================
+
+get_property(DRIVER_LIST GLOBAL PROPERTY DRIVERS_USED_IN_TESTS)
+if( "linear-solvers-example" IN_LIST DRIVER_LIST OR NOT PARTHENON_DISABLE_EXAMPLES)
+  add_executable(
+      linear-solvers-example
+          linear_solver_driver.cpp
+          linear_solver_driver.hpp
+          poisson_cell_equation.hpp
+          poisson_cell_package.cpp
+          poisson_cell_package.hpp
+          poisson_nodal_equation.hpp
+          poisson_nodal_package.cpp
+          poisson_nodal_package.hpp
+          helmholtz_equation.cpp
+          helmholtz_equation.hpp
+          helmholtz_package.cpp
+          helmholtz_package.hpp
+          main.cpp
+          parthenon_app_inputs.cpp
+          variable_type.hpp
+  )
+  target_link_libraries(linear-solvers-example PRIVATE Parthenon::parthenon)
+  lint_target(linear-solvers-example)
+endif()
diff --git a/example/linear_solvers/helmholtz_equation.cpp b/example/linear_solvers/helmholtz_equation.cpp
new file mode 100644
index 0000000000000..0d5326087000f
--- /dev/null
+++ b/example/linear_solvers/helmholtz_equation.cpp
@@ -0,0 +1,169 @@
+//========================================================================================
+// (C) (or copyright) 2023-2024. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#include <memory>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <kokkos_abstraction.hpp>
+#include <parthenon/package.hpp>
+
+#include "helmholtz_equation.hpp"
+
+namespace helmholtz_package {
+using namespace parthenon::package::prelude;
+
+parthenon::TaskStatus
+HelmholtzEquation::AxImpl(std::shared_ptr<parthenon::MeshData<Real>> &md_in,
+                          std::shared_ptr<parthenon::MeshData<Real>> &md_out) {
+  using namespace parthenon;
+  using TE = TopologicalElement;
+  auto pkg = md_in->GetMeshPointer()->packages.Get("helmholtz_package");
+  const auto alpha = pkg->Param<Real>("diagonal_alpha");
+
+  const int ndim = md_in->GetMeshPointer()->ndim;
+  IndexRange ib = md_in->GetBoundsI(IndexDomain::interior);
+  IndexRange jb = md_in->GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = md_in->GetBoundsK(IndexDomain::interior);
+
+  auto desc = parthenon::MakePackDescriptorFromTypeList<IndependentVars>(md_in.get());
+  auto pack_in = desc.GetPack(md_in.get());
+  auto pack_out = desc.GetPack(md_out.get());
+
+  const int ioff = ndim > 0;
+  const int joff = ndim > 1;
+  const int koff = ndim > 2;
+  parthenon::par_for(
+      "HelmholtzEquation::Ax", 0, pack_in.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s,
+      ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        const auto &coords = pack_in.GetCoordinates(b);
+        const Real dx1 = coords.template Dxc<X1DIR>(k, j, i);
+        const Real dx2 = coords.template Dxc<X2DIR>(k, j, i);
+        const Real dx3 = coords.template Dxc<X3DIR>(k, j, i);
+
+        Real Ax = -alpha * pack_in(b, TE::CC, vcc_t(), k, j, i);
+        Ax -= (pack_in(b, TE::F1, vfc_t(), k, j, i + ioff) -
+               pack_in(b, TE::F1, vfc_t(), k, j, i)) /
+              dx1;
+        Ax -= (pack_in(b, TE::F2, vfc_t(), k, j + joff, i) -
+               pack_in(b, TE::F2, vfc_t(), k, j, i)) /
+              dx2;
+        Ax -= (pack_in(b, TE::F3, vfc_t(), k + koff, j, i) -
+               pack_in(b, TE::F3, vfc_t(), k, j, i)) /
+              dx3;
+
+        pack_out(b, TE::CC, vcc_t(), k, j, i) = Ax;
+      });
+  std::vector<TE> tes{TE::F1};
+  if (ndim > 1) tes.push_back(TE::F2);
+  if (ndim > 2) tes.push_back(TE::F3);
+  for (auto &&te : tes) {
+    IndexRange ib = md_in->GetBoundsI(IndexDomain::interior, te);
+    IndexRange jb = md_in->GetBoundsJ(IndexDomain::interior, te);
+    IndexRange kb = md_in->GetBoundsK(IndexDomain::interior, te);
+    const int ioff = TopologicalOffsetI(te) * (ndim > 0);
+    const int joff = TopologicalOffsetJ(te) * (ndim > 1);
+    const int koff = TopologicalOffsetK(te) * (ndim > 2);
+    parthenon::par_for(
+        "HelmholtzEquation::Ax", 0, pack_in.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e,
+        ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          const auto &coords = pack_in.GetCoordinates(b);
+          const Real dx1 = coords.template Dxc<X1DIR>(k, j, i);
+          const Real dx2 = coords.template Dxc<X2DIR>(k, j, i);
+          const Real dx3 = coords.template Dxc<X3DIR>(k, j, i);
+
+          Real Ax = alpha * pack_in(b, te, vfc_t(), k, j, i);
+          Ax += (pack_in(b, TE::CC, vcc_t(), k, j, i) -
+                 pack_in(b, TE::CC, vcc_t(), k - koff, j - joff, i - ioff)) /
+                dx1;
+          pack_out(b, te, vfc_t(), k, j, i) = Ax;
+        });
+  }
+  return TaskStatus::complete;
+}
+
+parthenon::TaskStatus
+HelmholtzEquation::SetBoundary(std::shared_ptr<parthenon::MeshData<Real>> &md,
+                               bool coarse) {
+  using namespace parthenon;
+
+  using TE = TopologicalElement;
+  const int ndim = md->GetMeshPointer()->ndim;
+  
+  CellLevel cl = coarse ? CellLevel::coarse : CellLevel::same;
+
+  std::set<PDOpt> opts{};
+  if (coarse) opts.emplace(PDOpt::Coarse);
+  auto desc = parthenon::MakePackDescriptor<vfc_t>(md.get(), {}, opts);
+  auto pack = desc.GetPack(md.get(), GetBlockSelector::OnPhysicalBoundary());
+
+  std::vector<TE> tes{TE::F1};
+  if (ndim > 1) tes.push_back(TE::F2);
+  if (ndim > 2) tes.push_back(TE::F3);
+  for (auto &&te : tes) {
+    IndexRange ib = md->GetBoundsI(cl, IndexDomain::interior, te);
+    IndexRange jb = md->GetBoundsJ(cl, IndexDomain::interior, te);
+    IndexRange kb = md->GetBoundsK(cl, IndexDomain::interior, te);
+
+    parthenon::par_for(
+        "PoissonNodal::SetBoundary", 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e,
+        ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          const int oi = TopologicalOffsetI(te) * ((ib.e == i) - (ib.s == i));
+          const int oj = TopologicalOffsetJ(te) * ((jb.e == j) - (jb.s == j));
+          const int ok = TopologicalOffsetK(te) * ((kb.e == k) - (kb.s == k));
+          if (pack.IsPhysicalBoundary(b, ok, oj, oi)) pack(b, te, vfc_t(), k, j, i) = 0.0;
+        });
+  }
+  return TaskStatus::complete;
+}
+
+parthenon::TaskStatus
+HelmholtzEquation::SetDiagonal(std::shared_ptr<parthenon::MeshData<Real>> & /*md_mat*/,
+                               std::shared_ptr<parthenon::MeshData<Real>> &md_diag) {
+  using namespace parthenon;
+  using TE = TopologicalElement;
+  auto pkg = md_diag->GetMeshPointer()->packages.Get("helmholtz_package");
+  const auto alpha = pkg->Param<Real>("diagonal_alpha");
+
+  const int ndim = md_diag->GetMeshPointer()->ndim;
+  IndexRange ib = md_diag->GetBoundsI(IndexDomain::interior);
+  IndexRange jb = md_diag->GetBoundsJ(IndexDomain::interior);
+  IndexRange kb = md_diag->GetBoundsK(IndexDomain::interior);
+
+  auto desc = parthenon::MakePackDescriptorFromTypeList<IndependentVars>(md_diag.get());
+  auto pack_diag = desc.GetPack(md_diag.get());
+
+  parthenon::par_for(
+      "HelmholtzEquation::Ax", 0, pack_diag.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e,
+      ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        pack_diag(b, TE::CC, vcc_t(), k, j, i) = -alpha;
+      });
+
+  std::vector<TE> tes{TE::F1};
+  if (ndim > 1) tes.push_back(TE::F2);
+  if (ndim > 2) tes.push_back(TE::F3);
+  for (auto &&te : tes) {
+    IndexRange ib = md_diag->GetBoundsI(IndexDomain::interior, te);
+    IndexRange jb = md_diag->GetBoundsJ(IndexDomain::interior, te);
+    IndexRange kb = md_diag->GetBoundsK(IndexDomain::interior, te);
+    parthenon::par_for(
+        "HelmholtzEquation::Ax", 0, pack_diag.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e,
+        ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          pack_diag(b, te, vfc_t(), k, j, i) = alpha;
+        });
+  }
+  return TaskStatus::complete;
+}
+
+} // namespace helmholtz_package
diff --git a/example/linear_solvers/helmholtz_equation.hpp b/example/linear_solvers/helmholtz_equation.hpp
new file mode 100644
index 0000000000000..58e858f8e745d
--- /dev/null
+++ b/example/linear_solvers/helmholtz_equation.hpp
@@ -0,0 +1,65 @@
+//========================================================================================
+// (C) (or copyright) 2023-2024. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef EXAMPLE_LINEAR_SOLVERS_HELMHOLTZ_EQUATION_HPP_
+#define EXAMPLE_LINEAR_SOLVERS_HELMHOLTZ_EQUATION_HPP_
+
+#include <memory>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <kokkos_abstraction.hpp>
+#include <parthenon/package.hpp>
+
+#include "helmholtz_package.hpp"
+
+namespace helmholtz_package {
+using namespace parthenon::package::prelude;
+
+// This class implement methods for calculating A.x = y and returning the diagonal of A,
+// where A is the the matrix representing the discretized Poisson equation on the grid.
+// Here we implement the Laplace operator in terms of a flux divergence to (potentially)
+// consistently deal with coarse fine boundaries on the grid. Only the routines Ax and
+// SetDiagonal need to be defined for interfacing this with solvers. The other methods
+// are internal, but can't be marked private or protected because they launch kernels
+// on device.
+class HelmholtzEquation {
+ public:
+  using vcc_t = u;
+  using vfc_t = F;
+  using IndependentVars = parthenon::TypeList<vcc_t, vfc_t>;
+
+  HelmholtzEquation(parthenon::ParameterInput *pin, const std::string &label) {}
+
+  parthenon::TaskID Ax(parthenon::TaskList &tl, parthenon::TaskID depends_on,
+                       std::shared_ptr<parthenon::MeshData<Real>> & /*md_mat*/,
+                       std::shared_ptr<parthenon::MeshData<Real>> &md_in,
+                       std::shared_ptr<parthenon::MeshData<Real>> &md_out) {
+    return tl.AddTask(depends_on, AxImpl, md_in, md_out);
+  }
+
+  static parthenon::TaskStatus AxImpl(std::shared_ptr<parthenon::MeshData<Real>> &md_in,
+                                      std::shared_ptr<parthenon::MeshData<Real>> &md_out);
+
+  static parthenon::TaskStatus SetBoundary(std::shared_ptr<parthenon::MeshData<Real>> &md,
+                                           bool coarse);
+
+  parthenon::TaskStatus
+  SetDiagonal(std::shared_ptr<parthenon::MeshData<Real>> & /*md_mat*/,
+              std::shared_ptr<parthenon::MeshData<Real>> &md_diag);
+};
+
+} // namespace helmholtz_package
+
+#endif // EXAMPLE_LINEAR_SOLVERS_HELMHOLTZ_EQUATION_HPP_
diff --git a/example/linear_solvers/helmholtz_package.cpp b/example/linear_solvers/helmholtz_package.cpp
new file mode 100644
index 0000000000000..81cbdf2de7eec
--- /dev/null
+++ b/example/linear_solvers/helmholtz_package.cpp
@@ -0,0 +1,196 @@
+//========================================================================================
+// (C) (or copyright) 2021-2024. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <bvals/boundary_conditions_generic.hpp>
+#include <coordinates/coordinates.hpp>
+#include <parthenon/driver.hpp>
+#include <parthenon/package.hpp>
+#include <solvers/bicgstab_solver.hpp>
+#include <solvers/cg_solver.hpp>
+#include <solvers/solver_utils.hpp>
+
+#include "defs.hpp"
+#include "helmholtz_equation.hpp"
+#include "helmholtz_package.hpp"
+#include "kokkos_abstraction.hpp"
+
+using namespace parthenon::package::prelude;
+using parthenon::HostArray1D;
+namespace helmholtz_package {
+
+using namespace parthenon;
+using namespace parthenon::BoundaryFunction;
+// We need to register FixedFace boundary conditions by hand since they can't
+// be chosen in the parameter input file. FixedFace boundary conditions assume
+// Dirichlet booundary conditions on the face of the domain and linearly extrapolate
+// into the ghosts to ensure the linear reconstruction on the block face obeys the
+// chosen boundary condition. Just setting the ghost zones of CC variables to a fixed
+// value results in poor MG convergence because the effective BC at the face
+// changes with MG level.
+
+// Build type that selects only variables within the helmholtz namespace. Internal solver
+// variables have the namespace of input variables prepended, so they will also be
+// selected by this type.
+struct any_helmholtz : public parthenon::variable_names::base_t<true> {
+  template <class... Ts>
+  KOKKOS_INLINE_FUNCTION any_helmholtz(Ts &&...args)
+      : base_t<true>(std::forward<Ts>(args)...) {}
+  static std::string name() { return "helmholtz[.].*"; }
+};
+
+template <CoordinateDirection DIR, BCSide SIDE>
+auto GetBC() {
+  return [](std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) -> void {
+    using namespace parthenon;
+    using namespace parthenon::BoundaryFunction;
+    GenericBC<DIR, SIDE, BCType::FixedFace, any_helmholtz>(rc, coarse, 0.0);
+  };
+}
+
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
+  auto pkg = std::make_shared<StateDescriptor>("helmholtz_package");
+
+  // Set boundary conditions for helmholtz variables
+  using BF = parthenon::BoundaryFace;
+  pkg->UserBoundaryFunctions[BF::inner_x1].push_back(GetBC<X1DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::inner_x2].push_back(GetBC<X2DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::inner_x3].push_back(GetBC<X3DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::outer_x1].push_back(GetBC<X1DIR, BCSide::Outer>());
+  pkg->UserBoundaryFunctions[BF::outer_x2].push_back(GetBC<X2DIR, BCSide::Outer>());
+  pkg->UserBoundaryFunctions[BF::outer_x3].push_back(GetBC<X3DIR, BCSide::Outer>());
+
+  Real diagonal_alpha = pin->GetOrAddReal("helmholtz", "diagonal_alpha", 0.0);
+  pkg->AddParam<>("diagonal_alpha", diagonal_alpha);
+
+  std::string solver = pin->GetOrAddString("helmholtz", "solver", "MG");
+  pkg->AddParam<>("solver", solver);
+
+  bool use_exact_rhs = pin->GetOrAddBoolean("helmholtz", "use_exact_rhs", false);
+  pkg->AddParam<>("use_exact_rhs", use_exact_rhs);
+
+  std::string prolong =
+      pin->GetOrAddString("helmholtz", "boundary_prolongation", "Linear");
+
+  HelmholtzEquation eq(pin, "helmholtz");
+  pkg->AddParam<>("helmholtz_equation", eq, parthenon::Params::Mutability::Mutable);
+
+  std::shared_ptr<parthenon::solvers::SolverBase> psolver;
+  using prolongator_t = parthenon::solvers::ProlongationBlockInteriorDefault;
+  using preconditioner_t = parthenon::solvers::MGSolver<HelmholtzEquation, prolongator_t>;
+  const std::string base_label = "base";
+  const std::string u_label = "helmholtz_u";
+  const std::string rhs_label = "helmholtz_rhs";
+  if (solver == "MG") {
+    psolver =
+        std::make_shared<parthenon::solvers::MGSolver<HelmholtzEquation, prolongator_t>>(
+            base_label, u_label, rhs_label, pin, "helmholtz/solver_params",
+            HelmholtzEquation(pin, "helmholtz"));
+  } else if (solver == "CG") {
+    psolver = std::make_shared<
+        parthenon::solvers::CGSolver<HelmholtzEquation, preconditioner_t>>(
+        base_label, u_label, rhs_label, pin, "helmholtz/solver_params",
+        HelmholtzEquation(pin, "helmholtz"));
+  } else if (solver == "BiCGSTAB") {
+    psolver = std::make_shared<
+        parthenon::solvers::BiCGSTABSolver<HelmholtzEquation, preconditioner_t>>(
+        base_label, u_label, rhs_label, pin, "helmholtz/solver_params",
+        HelmholtzEquation(pin, "helmholtz"));
+  } else {
+    PARTHENON_FAIL("Unknown solver type.");
+  }
+  pkg->AddParam<>("solver_pointer", psolver);
+
+  using namespace parthenon::refinement_ops;
+
+  std::vector<MetadataFlag> flags_cc{Metadata::Cell, Metadata::Independent,
+                                     Metadata::FillGhost, Metadata::GMGRestrict,
+                                     Metadata::GMGProlongate};
+  std::vector<MetadataFlag> flags_fc{Metadata::Face, Metadata::Independent,
+                                     Metadata::FillGhost, Metadata::GMGRestrict,
+                                     Metadata::GMGProlongate};
+  auto mflux_comm_cc = Metadata(flags_cc);
+  auto mflux_comm_fc = Metadata(flags_fc);
+  if (prolong == "Linear") {
+    mflux_comm_cc.RegisterRefinementOps<ProlongateSharedLinear, RestrictAverage>();
+    mflux_comm_fc.RegisterRefinementOps<ProlongateSharedLinear, RestrictAverage>();
+  } else if (prolong == "Constant") {
+    mflux_comm_cc.RegisterRefinementOps<ProlongatePiecewiseConstant, RestrictAverage>();
+    mflux_comm_fc.RegisterRefinementOps<ProlongatePiecewiseConstant, RestrictAverage>();
+  } else {
+    PARTHENON_FAIL("Unknown prolongation method for Helmholtz boundaries.");
+  }
+  // u is the solution vector that starts with an initial guess and then gets updated
+  // by the solver
+  pkg->AddField<u>(mflux_comm_cc);
+  pkg->AddField<F>(mflux_comm_fc);
+
+  return pkg;
+}
+
+parthenon::TaskStatus
+SetVector(parthenon::ParameterInput *pin, bool use_exponential,
+          std::shared_ptr<parthenon::MeshData<parthenon::Real>> md) {
+  using namespace parthenon;
+  Real x0 = pin->GetOrAddReal("helmholtz", "x0", 0.0);
+  Real y0 = pin->GetOrAddReal("helmholtz", "y0", 0.0);
+  Real z0 = pin->GetOrAddReal("helmholtz", "z0", 0.0);
+  Real radius0 = pin->GetOrAddReal("helmholtz", "radius", 0.1);
+  const int ndim = md->GetMeshPointer()->ndim;
+
+  auto desc = MakePackDescriptor<u, F>(md.get());
+  auto pack = desc.GetPack(md.get());
+
+  using TE = parthenon::TopologicalElement;
+  auto ib = md->GetBoundsI(IndexDomain::entire, TE::CC);
+  auto jb = md->GetBoundsJ(IndexDomain::entire, TE::CC);
+  auto kb = md->GetBoundsK(IndexDomain::entire, TE::CC);
+
+  parthenon::par_for(
+      "Helmholtz::rhs_u", 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        const auto &coords = pack.GetCoordinates(b);
+        Real x1 = coords.X<1, TE::CC>(i);
+        Real x2 = coords.X<2, TE::CC>(j);
+        Real x3 = coords.X<2, TE::CC>(k);
+        Real rad = (x1 - x0) * (x1 - x0);
+        if (ndim > 1) rad += (x2 - y0) * (x2 - y0);
+        if (ndim > 2) rad += (x3 - z0) * (x3 - z0);
+        rad = std::sqrt(rad);
+
+        pack(b, TE::CC, u(), k, j, i) = rad < radius0 ? 1.0 : 0.0;
+        if (use_exponential) pack(b, TE::CC, u(), k, j, i) = -exp(-10.0 * rad * rad);
+      });
+
+  for (auto te : {TE::F1, TE::F2, TE::F3}) {
+    auto ib = md->GetBoundsI(IndexDomain::entire, te);
+    auto jb = md->GetBoundsJ(IndexDomain::entire, te);
+    auto kb = md->GetBoundsK(IndexDomain::entire, te);
+    parthenon::par_for(
+        "Helmholtz::rhs_F", 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+        KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          pack(b, te, F(), k, j, i) = 0.0;
+        });
+  }
+
+  return TaskStatus::complete;
+}
+
+} // namespace helmholtz_package
diff --git a/example/linear_solvers/helmholtz_package.hpp b/example/linear_solvers/helmholtz_package.hpp
new file mode 100644
index 0000000000000..6ca07172fc191
--- /dev/null
+++ b/example/linear_solvers/helmholtz_package.hpp
@@ -0,0 +1,37 @@
+//========================================================================================
+// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef EXAMPLE_LINEAR_SOLVERS_HELMHOLTZ_PACKAGE_HPP_
+#define EXAMPLE_LINEAR_SOLVERS_HELMHOLTZ_PACKAGE_HPP_
+
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <kokkos_abstraction.hpp>
+#include <parthenon/package.hpp>
+
+#include "variable_type.hpp"
+
+namespace helmholtz_package {
+using namespace parthenon::package::prelude;
+
+VARIABLE(helmholtz, u);
+VARIABLE(helmholtz, F);
+
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
+TaskStatus SetVector(ParameterInput *pin, bool use_exponential,
+                     std::shared_ptr<MeshData<Real>> md);
+} // namespace helmholtz_package
+
+#endif // EXAMPLE_LINEAR_SOLVERS_HELMHOLTZ_PACKAGE_HPP_
diff --git a/example/linear_solvers/linear_solver_driver.cpp b/example/linear_solvers/linear_solver_driver.cpp
new file mode 100644
index 0000000000000..67cdb64ab2d6d
--- /dev/null
+++ b/example/linear_solvers/linear_solver_driver.cpp
@@ -0,0 +1,175 @@
+//========================================================================================
+// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+#include <algorithm>
+#include <cstdio>
+#include <memory>
+#include <string>
+#include <vector>
+
+// Local Includes
+#include <amr_criteria/refinement_package.hpp>
+#include <bvals/comms/bvals_in_one.hpp>
+#include <interface/metadata.hpp>
+#include <interface/update.hpp>
+#include <mesh/meshblock_pack.hpp>
+#include <parthenon/driver.hpp>
+#include <prolong_restrict/prolong_restrict.hpp>
+#include <solvers/solver_base.hpp>
+#include <solvers/solver_utils.hpp>
+
+#include "helmholtz_equation.hpp"
+#include "helmholtz_package.hpp"
+#include "linear_solver_driver.hpp"
+#include "poisson_cell_equation.hpp"
+#include "poisson_cell_package.hpp"
+#include "poisson_nodal_equation.hpp"
+#include "poisson_nodal_package.hpp"
+
+using namespace parthenon::driver::prelude;
+
+namespace linear_solver_example {
+
+parthenon::DriverStatus LinearSolverDriver::Execute() {
+  using namespace parthenon;
+
+  pouts->MakeOutputs(pmesh, pinput);
+  ConstructAndExecuteTaskLists<>(this);
+  pouts->MakeOutputs(pmesh, pinput);
+
+  return DriverStatus::complete;
+}
+
+TaskCollection LinearSolverDriver::MakeTaskCollection(BlockList_t &blocks) {
+  using namespace parthenon;
+  TaskCollection tc;
+
+  {
+    using namespace poisson_nodal_package;
+    initialize_vector_func_t Initialize = [](ParameterInput *,
+                                             std::shared_ptr<MeshData<Real>>) {
+      return TaskStatus::complete;
+    };
+    auto SetRHS = [](auto *pinput, auto pmd) { return SetVector(pinput, false, pmd); };
+    auto SetExact = [](auto *pinput, auto pmd) { return SetVector(pinput, true, pmd); };
+    AddSolverTaskRegion<PoissonEquation<u>::IndependentVars>(
+        tc, "poisson_nodal_package", Initialize, SetRHS, SetExact);
+  }
+
+  {
+    using namespace poisson_cell_package;
+    auto SetRHS = [](auto *pinput, auto pmd) { return SetVector(pinput, false, pmd); };
+    auto SetExact = [](auto *pinput, auto pmd) { return SetVector(pinput, true, pmd); };
+    AddSolverTaskRegion<PoissonEquation<u, D>::IndependentVars>(
+        tc, "poisson_cell_package", SetD, SetRHS, SetExact);
+  }
+
+  {
+    using namespace helmholtz_package;
+    initialize_vector_func_t Initialize = [](ParameterInput *,
+                                             std::shared_ptr<MeshData<Real>>) {
+      return TaskStatus::complete;
+    };
+    auto SetRHS = [](auto *pinput, auto pmd) { return SetVector(pinput, false, pmd); };
+    auto SetExact = [](auto *pinput, auto pmd) { return SetVector(pinput, true, pmd); };
+    AddSolverTaskRegion<HelmholtzEquation::IndependentVars>(tc, "helmholtz_package",
+                                                            Initialize, SetRHS, SetExact);
+  }
+
+  return tc;
+}
+
+template <class solver_TL>
+void LinearSolverDriver::AddSolverTaskRegion(parthenon::TaskCollection &tc,
+                                             std::string package_label,
+                                             initialize_vector_func_t Initialize,
+                                             initialize_vector_func_t SetRHS,
+                                             initialize_vector_func_t SetExact) {
+  using namespace parthenon;
+  TaskID none(0);
+
+  auto pkg = pmesh->packages.Get(package_label);
+  auto use_exact_rhs = pkg->Param<bool>("use_exact_rhs");
+  auto psolver =
+      pkg->Param<std::shared_ptr<parthenon::solvers::SolverBase>>("solver_pointer");
+
+  auto partitions = pmesh->GetDefaultBlockPartitions();
+  const int num_partitions = partitions.size();
+  TaskRegion &region = tc.AddRegion(num_partitions);
+  for (int i = 0; i < num_partitions; ++i) {
+    TaskList &tl = region[i];
+    auto &md = pmesh->mesh_data.Add("base", partitions[i]);
+    auto &field_labels = psolver->GetFieldLabels();
+    auto &md_u =
+        pmesh->mesh_data.Add(psolver->GetSolutionContainerLabel(), md, field_labels);
+    auto &md_rhs =
+        pmesh->mesh_data.Add(psolver->GetRHSContainerLabel(), md, field_labels);
+    auto &md_exact = pmesh->mesh_data.Add("exact_" + package_label, md, field_labels);
+
+    auto initialize = tl.AddTask(none, Initialize, pinput, md);
+    // set the rhs
+    auto set_rhs = tl.AddTask(initialize, SetRHS, pinput, md_rhs);
+
+    // Possibly set rhs <- A.u_exact for a given u_exact so that the exact solution is
+    // known when we solve A.u = rhs
+    if (use_exact_rhs) {
+      auto set_exact = tl.AddTask(set_rhs, SetExact, pinput, md_exact);
+      auto comm =
+          AddBoundaryExchangeTasks<BoundaryType::any>(set_exact, tl, md_exact, true);
+      set_rhs = psolver->Ax(tl, comm, md, md_exact, md_rhs);
+    }
+
+    // Set initial solution guess to zero
+    auto zero_u = tl.AddTask(set_rhs, TF(solvers::utils::SetToZero<solver_TL>), md_u);
+    auto setup = psolver->AddSetupTasks(tl, zero_u, i, pmesh);
+    auto solve = psolver->AddTasks(tl, setup, i, pmesh);
+
+    // If we are using a rhs to which we know the exact solution, compare our computed
+    // solution to the exact solution
+    if (use_exact_rhs) {
+      auto diff = tl.AddTask(solve, solvers::utils::AddFieldsAndStore<solver_TL>,
+                             md_exact, md_u, md_exact, 1.0, -1.0);
+      auto get_err =
+          solvers::utils::DotProduct<solver_TL>(diff, tl, &err, md_exact, md_exact);
+      tl.AddTask(
+          get_err,
+          [package_label](LinearSolverDriver *driver, int partition,
+                          std::shared_ptr<parthenon::solvers::SolverBase> psolver) {
+            if (partition != 0) return TaskStatus::complete;
+            driver->final_rms_error[package_label] =
+                std::sqrt(driver->err.val / driver->pmesh->GetTotalCells());
+            driver->final_rms_residual[package_label] = psolver->GetFinalResidual();
+            if (Globals::my_rank == 0)
+              printf("Final residual: %e\n", driver->final_rms_residual[package_label]);
+            printf("Final rms error: %e\n", driver->final_rms_error[package_label]);
+            return TaskStatus::complete;
+          },
+          this, i, psolver);
+    } else {
+      tl.AddTask(
+          solve,
+          [package_label](LinearSolverDriver *driver, int partition,
+                          std::shared_ptr<parthenon::solvers::SolverBase> psolver) {
+            if (partition != 0) return TaskStatus::complete;
+            driver->final_rms_error[package_label] = 0.0;
+            driver->final_rms_residual[package_label] = psolver->GetFinalResidual();
+            if (Globals::my_rank == 0)
+              printf("Final residual: %e\n", driver->final_rms_residual[package_label]);
+            return TaskStatus::complete;
+          },
+          this, i, psolver);
+    }
+  }
+}
+
+} // namespace linear_solver_example
diff --git a/example/linear_solvers/linear_solver_driver.hpp b/example/linear_solvers/linear_solver_driver.hpp
new file mode 100644
index 0000000000000..802b5722f0a3c
--- /dev/null
+++ b/example/linear_solvers/linear_solver_driver.hpp
@@ -0,0 +1,61 @@
+//========================================================================================
+// (C) (or copyright) 2021-2024. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+#ifndef EXAMPLE_LINEAR_SOLVERS_LINEAR_SOLVER_DRIVER_HPP_
+#define EXAMPLE_LINEAR_SOLVERS_LINEAR_SOLVER_DRIVER_HPP_
+
+#include <map>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include <kokkos_abstraction.hpp>
+#include <parthenon/driver.hpp>
+#include <parthenon/package.hpp>
+
+namespace linear_solver_example {
+using namespace parthenon::driver::prelude;
+
+class LinearSolverDriver : public Driver {
+ public:
+  LinearSolverDriver(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm)
+      : Driver(pin, app_in, pm) {
+    InitializeOutputs();
+  }
+  // This next function essentially defines the driver.
+  TaskCollection MakeTaskCollection(BlockList_t &blocks);
+
+  DriverStatus Execute() override;
+
+  std::map<std::string, Real> final_rms_error, final_rms_residual;
+
+  // Necessary reductions for checking error from exact solution
+  AllReduce<Real> err;
+
+ private:
+  using initialize_vector_func_t = std::function<parthenon::TaskStatus(
+      parthenon::ParameterInput *,
+      std::shared_ptr<parthenon::MeshData<parthenon::Real>>)>;
+  template <class solver_TypeList>
+  void AddSolverTaskRegion(parthenon::TaskCollection &tc, std::string pacakge_label,
+                           initialize_vector_func_t Initialize,
+                           initialize_vector_func_t SetRHS,
+                           initialize_vector_func_t SetExact);
+};
+
+void ProblemGenerator(Mesh *pm, parthenon::ParameterInput *pin, MeshData<Real> *md);
+parthenon::Packages_t ProcessPackages(std::unique_ptr<parthenon::ParameterInput> &pin);
+
+} // namespace linear_solver_example
+
+#endif // EXAMPLE_LINEAR_SOLVERS_LINEAR_SOLVER_DRIVER_HPP_
diff --git a/example/linear_solvers/main.cpp b/example/linear_solvers/main.cpp
new file mode 100644
index 0000000000000..0493c6ac0d900
--- /dev/null
+++ b/example/linear_solvers/main.cpp
@@ -0,0 +1,64 @@
+//========================================================================================
+// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+#include "parthenon_manager.hpp"
+
+#include "linear_solver_driver.hpp"
+
+int main(int argc, char *argv[]) {
+  using parthenon::ParthenonManager;
+  using parthenon::ParthenonStatus;
+  ParthenonManager pman;
+
+  // Redefine parthenon defaults
+  pman.app_input->ProcessPackages = linear_solver_example::ProcessPackages;
+
+  // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up
+  auto manager_status = pman.ParthenonInitEnv(argc, argv);
+  if (manager_status == ParthenonStatus::complete) {
+    pman.ParthenonFinalize();
+    return 0;
+  }
+  if (manager_status == ParthenonStatus::error) {
+    pman.ParthenonFinalize();
+    return 1;
+  }
+  // Now that ParthenonInit has been called and setup succeeded, the code can now
+  // make use of MPI and Kokkos
+
+  pman.ParthenonInitPackagesAndMesh();
+
+  // This needs to be scoped so that the driver object is destructed before Finalize
+  bool success = true;
+  {
+    // Initialize the driver
+    linear_solver_example::LinearSolverDriver driver(
+        pman.pinput.get(), pman.app_input.get(), pman.pmesh.get());
+
+    // This line actually runs the simulation
+    auto driver_status = driver.Execute();
+    if (driver_status != parthenon::DriverStatus::complete) success = false;
+    // Go through all the solutions that registered a solution quality and
+    // make sure they are all doing better than some thresholds
+    for (auto &[label, err] : driver.final_rms_error) {
+      auto res = driver.final_rms_residual[label];
+      if (res > 1.e-10 || res != res) success = false;
+      if (err > 1.e-12 || err != err) success = false;
+    }
+  }
+  // call MPI_Finalize and Kokkos::finalize if necessary
+  pman.ParthenonFinalize();
+
+  // MPI and Kokkos can no longer be used
+  return static_cast<int>(!success);
+}
diff --git a/example/linear_solvers/parthenon_app_inputs.cpp b/example/linear_solvers/parthenon_app_inputs.cpp
new file mode 100644
index 0000000000000..3d2e9e65bcfcf
--- /dev/null
+++ b/example/linear_solvers/parthenon_app_inputs.cpp
@@ -0,0 +1,43 @@
+//========================================================================================
+// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#include <math.h>
+#include <sstream>
+#include <string>
+
+#include <parthenon/package.hpp>
+
+#include "config.hpp"
+#include "defs.hpp"
+#include "helmholtz_package.hpp"
+#include "poisson_cell_package.hpp"
+#include "poisson_nodal_package.hpp"
+#include "utils/error_checking.hpp"
+
+using namespace parthenon::package::prelude;
+using namespace parthenon;
+
+// *************************************************//
+// redefine some weakly linked parthenon functions *//
+// *************************************************//
+
+namespace linear_solver_example {
+
+Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
+  Packages_t packages;
+  packages.Add(poisson_cell_package::Initialize(pin.get()));
+  packages.Add(poisson_nodal_package::Initialize(pin.get()));
+  packages.Add(helmholtz_package::Initialize(pin.get()));
+  return packages;
+}
+
+} // namespace linear_solver_example
diff --git a/example/linear_solvers/parthinput.poisson b/example/linear_solvers/parthinput.poisson
new file mode 100644
index 0000000000000..9117ba3c056ff
--- /dev/null
+++ b/example/linear_solvers/parthinput.poisson
@@ -0,0 +1,85 @@
+# ========================================================================================
+#  Parthenon performance portable AMR framework
+#  Copyright(C) 2020-2023 The Parthenon collaboration
+#  Licensed under the 3-clause BSD License, see LICENSE file for details
+# ========================================================================================
+#  (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved.
+#
+#  This program was produced under U.S. Government contract 89233218CNA000001 for Los
+#  Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+#  for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+#  in the program are reserved by Triad National Security, LLC, and the U.S. Department
+#  of Energy/National Nuclear Security Administration. The Government is granted for
+#  itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+#  license in this material to reproduce, prepare derivative works, distribute copies to
+#  the public, perform publicly and display publicly, and to permit others to do so.
+# ========================================================================================
+
+<parthenon/job>
+problem_id = poisson
+
+<parthenon/mesh>
+refinement = static
+multigrid = true 
+
+nx1 = 64
+x1min = -1.0
+x1max = 1.0
+ix1_bc = outflow
+ox1_bc = outflow
+
+nx2 = 64
+x2min = -1.0
+x2max = 1.0
+ix2_bc = outflow
+ox2_bc = outflow
+
+nx3 = 1
+x3min = 0.0
+x3max = 1.0
+ix3_bc = periodic
+ox3_bc = periodic
+
+<parthenon/meshblock>
+nx1 = 32
+nx2 = 32
+nx3 = 1
+
+<parthenon/time>
+#nlim = -1
+#tlim = 1.0
+#integrator = rk2
+#ncycle_out_mesh = -10000
+
+<parthenon/output0>
+file_type = hdf5
+dt = 0.05
+variables = poisson.res_err, poisson.u, poisson.x, poisson.r, poisson.rhs, poisson.p, poisson.s, poisson.t, poisson.v, poisson.exact
+ghost_zones = true 
+
+<parthenon/static_refinement0> 
+x1min = -1.0
+x1max = -0.75
+x2min = -1.0
+x2max = -0.75
+level = 3
+
+<poisson_nodal>
+solver = BiCGSTAB # or MG
+flux_correct = true
+diagonal_alpha = 0.0
+
+x0 = 0.0
+y0 = 0.0
+z0 = 0.0
+radius = 0.5
+interior_D = 100.0
+exterior_D = 1.0
+
+<poisson_nodal/solver_params>
+precondition = true
+max_iterations = 15
+residual_tolerance = 1.e-8
+print_per_step = true
+smoother = SRJ2
+do_FAS = true
diff --git a/example/linear_solvers/poisson_cell_equation.hpp b/example/linear_solvers/poisson_cell_equation.hpp
new file mode 100644
index 0000000000000..edfa7d0bb28ee
--- /dev/null
+++ b/example/linear_solvers/poisson_cell_equation.hpp
@@ -0,0 +1,318 @@
+//========================================================================================
+// (C) (or copyright) 2023-2024. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef EXAMPLE_LINEAR_SOLVERS_POISSON_CELL_EQUATION_HPP_
+#define EXAMPLE_LINEAR_SOLVERS_POISSON_CELL_EQUATION_HPP_
+
+#include <memory>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <kokkos_abstraction.hpp>
+#include <parthenon/package.hpp>
+
+#include "poisson_cell_package.hpp"
+
+namespace poisson_cell_package {
+
+// This class implement methods for calculating A.x = y and returning the diagonal of A,
+// where A is the the matrix representing the discretized Poisson equation on the grid.
+// Here we implement the Laplace operator in terms of a flux divergence to (potentially)
+// consistently deal with coarse fine boundaries on the grid. Only the routines Ax and
+// SetDiagonal need to be defined for interfacing this with solvers. The other methods
+// are internal, but can't be marked private or protected because they launch kernels
+// on device.
+template <class var_t, class D_t>
+class PoissonEquation {
+ public:
+  bool do_flux_cor = false;
+  bool set_flux_boundary = false;
+  bool include_flux_dx = false;
+
+  using IndependentVars = parthenon::TypeList<var_t>;
+
+  PoissonEquation(parthenon::ParameterInput *pin, const std::string &label) {
+    do_flux_cor = pin->GetOrAddBoolean(label, "flux_correct", false);
+    set_flux_boundary = pin->GetOrAddBoolean(label, "set_flux_boundary", false);
+    include_flux_dx =
+        (pin->GetOrAddString(label, "boundary_prolongation", "Linear") == "Constant");
+  }
+
+  // Add tasks to calculate the result of the matrix A (which is implicitly defined by
+  // this class) being applied to x_t and store it in field out_t
+  parthenon::TaskID Ax(parthenon::TaskList &tl, parthenon::TaskID depends_on,
+                       std::shared_ptr<parthenon::MeshData<Real>> &md_mat,
+                       std::shared_ptr<parthenon::MeshData<Real>> &md_in,
+                       std::shared_ptr<parthenon::MeshData<Real>> &md_out) {
+    auto flux_res = tl.AddTask(depends_on, CalculateFluxes, md_mat, md_in);
+    if (set_flux_boundary) {
+      flux_res = tl.AddTask(flux_res, SetFluxBoundaries, md_mat, md_in, include_flux_dx);
+    }
+    if (do_flux_cor && !(md_mat->grid.type == parthenon::GridType::two_level_composite)) {
+      auto start_flxcor =
+          tl.AddTask(flux_res, parthenon::StartReceiveFluxCorrections, md_in);
+      auto send_flxcor =
+          tl.AddTask(flux_res, parthenon::LoadAndSendFluxCorrections, md_in);
+      auto recv_flxcor =
+          tl.AddTask(start_flxcor, parthenon::ReceiveFluxCorrections, md_in);
+      flux_res = tl.AddTask(recv_flxcor, parthenon::SetFluxCorrections, md_in);
+    }
+    return tl.AddTask(flux_res, FluxMultiplyMatrix, md_in, md_out);
+  }
+
+  // Calculate an approximation to the diagonal of the matrix A and store it in diag_t.
+  // For a uniform grid or when flux correction is ignored, this diagonal calculation
+  // is exact. Exactness is (probably) not required since it is just used in Jacobi
+  // iterations.
+  parthenon::TaskStatus SetDiagonal(std::shared_ptr<parthenon::MeshData<Real>> &md_mat,
+                                    std::shared_ptr<parthenon::MeshData<Real>> &md_diag) {
+    using namespace parthenon;
+    const int ndim = md_mat->GetMeshPointer()->ndim;
+    IndexRange ib = md_mat->GetBoundsI(IndexDomain::interior, te);
+    IndexRange jb = md_mat->GetBoundsJ(IndexDomain::interior, te);
+    IndexRange kb = md_mat->GetBoundsK(IndexDomain::interior, te);
+
+    auto pkg = md_mat->GetMeshPointer()->packages.Get("poisson_cell_package");
+    const auto alpha = pkg->Param<Real>("diagonal_alpha");
+
+    int nblocks = md_mat->NumBlocks();
+    std::vector<bool> include_block(nblocks, true);
+
+    auto desc_mat = parthenon::MakePackDescriptor<D_t>(md_mat.get());
+    auto desc_diag = parthenon::MakePackDescriptor<var_t>(md_diag.get());
+    auto pack_mat = desc_mat.GetPack(md_mat.get(), include_block);
+    auto pack_diag = desc_diag.GetPack(md_diag.get(), include_block);
+    using TE = parthenon::TopologicalElement;
+    parthenon::par_for(
+        "StoreDiagonal", 0, pack_mat.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+        KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          const auto &coords = pack_mat.GetCoordinates(b);
+          // Build the unigrid diagonal of the matrix
+          Real dx1 = coords.template Dxc<X1DIR>(k, j, i);
+          Real diag_elem = -(pack_mat(b, TE::F1, D_t(), k, j, i) +
+                             pack_mat(b, TE::F1, D_t(), k, j, i + 1)) /
+                               (dx1 * dx1) -
+                           alpha;
+          if (ndim > 1) {
+            Real dx2 = coords.template Dxc<X2DIR>(k, j, i);
+            diag_elem -= (pack_mat(b, TE::F2, D_t(), k, j, i) +
+                          pack_mat(b, TE::F2, D_t(), k, j + 1, i)) /
+                         (dx2 * dx2);
+          }
+          if (ndim > 2) {
+            Real dx3 = coords.template Dxc<X3DIR>(k, j, i);
+            diag_elem -= (pack_mat(b, TE::F3, D_t(), k, j, i) +
+                          pack_mat(b, TE::F3, D_t(), k + 1, j, i)) /
+                         (dx3 * dx3);
+          }
+          pack_diag(b, te, var_t(), k, j, i) = diag_elem;
+        });
+    return TaskStatus::complete;
+  }
+
+  static parthenon::TaskStatus
+  CalculateFluxes(std::shared_ptr<parthenon::MeshData<Real>> &md_mat,
+                  std::shared_ptr<parthenon::MeshData<Real>> &md) {
+    using namespace parthenon;
+    const int ndim = md->GetMeshPointer()->ndim;
+    using TE = parthenon::TopologicalElement;
+    TE te = TE::CC;
+    IndexRange ib = md->GetBoundsI(IndexDomain::interior, te);
+    IndexRange jb = md->GetBoundsJ(IndexDomain::interior, te);
+    IndexRange kb = md->GetBoundsK(IndexDomain::interior, te);
+
+    using TE = parthenon::TopologicalElement;
+
+    int nblocks = md->NumBlocks();
+    std::vector<bool> include_block(nblocks, true);
+
+    auto desc = parthenon::MakePackDescriptor<var_t>(md.get(), {}, {PDOpt::WithFluxes});
+    auto pack = desc.GetPack(md.get(), include_block);
+    auto desc_mat = parthenon::MakePackDescriptor<D_t>(md_mat.get(), {});
+    auto pack_mat = desc_mat.GetPack(md_mat.get(), include_block);
+    parthenon::par_for(
+        "CaclulateFluxes", 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+        KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          const auto &coords = pack.GetCoordinates(b);
+          Real dx1 = coords.template Dxc<X1DIR>(k, j, i);
+          pack.flux(b, X1DIR, var_t(), k, j, i) =
+              pack_mat(b, TE::F1, D_t(), k, j, i) / dx1 *
+              (pack(b, te, var_t(), k, j, i - 1) - pack(b, te, var_t(), k, j, i));
+          if (i == ib.e)
+            pack.flux(b, X1DIR, var_t(), k, j, i + 1) =
+                pack_mat(b, TE::F1, D_t(), k, j, i + 1) / dx1 *
+                (pack(b, te, var_t(), k, j, i) - pack(b, te, var_t(), k, j, i + 1));
+
+          if (ndim > 1) {
+            Real dx2 = coords.template Dxc<X2DIR>(k, j, i);
+            pack.flux(b, X2DIR, var_t(), k, j, i) =
+                pack_mat(b, TE::F2, D_t(), k, j, i) *
+                (pack(b, te, var_t(), k, j - 1, i) - pack(b, te, var_t(), k, j, i)) / dx2;
+            if (j == jb.e)
+              pack.flux(b, X2DIR, var_t(), k, j + 1, i) =
+                  pack_mat(b, TE::F2, D_t(), k, j + 1, i) *
+                  (pack(b, te, var_t(), k, j, i) - pack(b, te, var_t(), k, j + 1, i)) /
+                  dx2;
+          }
+
+          if (ndim > 2) {
+            Real dx3 = coords.template Dxc<X3DIR>(k, j, i);
+            pack.flux(b, X3DIR, var_t(), k, j, i) =
+                pack_mat(b, TE::F3, D_t(), k, j, i) *
+                (pack(b, te, var_t(), k - 1, j, i) - pack(b, te, var_t(), k, j, i)) / dx3;
+            if (k == kb.e)
+              pack.flux(b, X2DIR, var_t(), k + 1, j, i) =
+                  pack_mat(b, TE::F3, D_t(), k + 1, j, i) *
+                  (pack(b, te, var_t(), k, j, i) - pack(b, te, var_t(), k + 1, j, i)) /
+                  dx3;
+          }
+        });
+    return TaskStatus::complete;
+  }
+
+  static parthenon::TaskStatus
+  SetFluxBoundaries(std::shared_ptr<parthenon::MeshData<Real>> &md_mat,
+                    std::shared_ptr<parthenon::MeshData<Real>> &md, bool do_flux_dx) {
+    using namespace parthenon;
+    const int ndim = md->GetMeshPointer()->ndim;
+    IndexRange ib = md->GetBoundsI(IndexDomain::interior);
+    IndexRange jb = md->GetBoundsJ(IndexDomain::interior);
+    IndexRange kb = md->GetBoundsK(IndexDomain::interior);
+
+    using TE = parthenon::TopologicalElement;
+
+    int nblocks = md->NumBlocks();
+    std::vector<bool> include_block(nblocks, true);
+
+    auto desc = parthenon::MakePackDescriptor<var_t>(md.get(), {}, {PDOpt::WithFluxes});
+    auto desc_mat = parthenon::MakePackDescriptor<D_t>(md.get());
+    auto pack = desc.GetPack(md.get(), include_block);
+    auto pack_mat = desc_mat.GetPack(md_mat.get(), include_block);
+    const std::size_t scratch_size_in_bytes = 0;
+    const std::size_t scratch_level = 1;
+
+    const parthenon::Indexer3D idxers[6]{
+        parthenon::Indexer3D(kb, jb, {ib.s, ib.s}),
+        parthenon::Indexer3D(kb, jb, {ib.e + 1, ib.e + 1}),
+        parthenon::Indexer3D(kb, {jb.s, jb.s}, ib),
+        parthenon::Indexer3D(kb, {jb.e + 1, jb.e + 1}, ib),
+        parthenon::Indexer3D({kb.s, kb.s}, jb, ib),
+        parthenon::Indexer3D({kb.e + 1, kb.e + 1}, jb, ib)};
+    constexpr int x1off[6]{-1, 1, 0, 0, 0, 0};
+    constexpr int x2off[6]{0, 0, -1, 1, 0, 0};
+    constexpr int x3off[6]{0, 0, 0, 0, -1, 1};
+    constexpr TE tes[6]{TE::F1, TE::F1, TE::F2, TE::F2, TE::F3, TE::F3};
+    constexpr int dirs[6]{X1DIR, X1DIR, X2DIR, X2DIR, X3DIR, X3DIR};
+    parthenon::par_for_outer(
+        DEFAULT_OUTER_LOOP_PATTERN, "SetFluxBoundaries", DevExecSpace(),
+        scratch_size_in_bytes, scratch_level, 0, pack.GetNBlocks() - 1,
+        KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b) {
+          const auto &coords = pack.GetCoordinates(b);
+          const int gid = pack.GetGID(b);
+          const int level = pack.GetLevel(b, 0, 0, 0);
+          const Real dxs[3]{coords.template Dxc<X1DIR>(), coords.template Dxc<X2DIR>(),
+                            coords.template Dxc<X3DIR>()};
+          for (int face = 0; face < ndim * 2; ++face) {
+            const Real dx = dxs[dirs[face] - 1];
+            const auto &idxer = idxers[face];
+            const auto dir = dirs[face];
+            const auto te = tes[face];
+            // Impose the zero Dirichlet boundary condition at the actual boundary
+            if (pack.IsPhysicalBoundary(b, x3off[face], x2off[face], x1off[face])) {
+              const int koff = x3off[face] > 0 ? -1 : 0;
+              const int joff = x2off[face] > 0 ? -1 : 0;
+              const int ioff = x1off[face] > 0 ? -1 : 0;
+              const int sign = x1off[face] + x2off[face] + x3off[face];
+              parthenon::par_for_inner(
+                  DEFAULT_INNER_LOOP_PATTERN, member, 0, idxer.size() - 1,
+                  [&](const int idx) {
+                    const auto [k, j, i] = idxer(idx);
+                    pack.flux(b, dir, var_t(), k, j, i) =
+                        sign * pack_mat(b, te, D_t(), k, j, i) *
+                        pack(b, var_t(), k + koff, j + joff, i + ioff) / (0.5 * dx);
+                  });
+            }
+            // Correct for size of neighboring zone at fine-coarse boundary when using
+            // constant prolongation
+            if (do_flux_dx &&
+                pack.GetLevel(b, x3off[face], x2off[face], x1off[face]) == level - 1) {
+              parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, member, 0,
+                                       idxer.size() - 1, [&](const int idx) {
+                                         const auto [k, j, i] = idxer(idx);
+                                         pack.flux(b, dir, var_t(), k, j, i) /= 1.5;
+                                       });
+            }
+          }
+        });
+    return TaskStatus::complete;
+  }
+
+  // Calculate A in_t = out_t (in the region covered by md) for a given set of fluxes
+  // calculated with in_t (which have possibly been corrected at coarse fine boundaries)
+  static parthenon::TaskStatus
+  FluxMultiplyMatrix(std::shared_ptr<parthenon::MeshData<Real>> &md,
+                     std::shared_ptr<parthenon::MeshData<Real>> &md_out) {
+    using namespace parthenon;
+    const int ndim = md->GetMeshPointer()->ndim;
+    using TE = parthenon::TopologicalElement;
+    TE te = TE::CC;
+    IndexRange ib = md->GetBoundsI(IndexDomain::interior, te);
+    IndexRange jb = md->GetBoundsJ(IndexDomain::interior, te);
+    IndexRange kb = md->GetBoundsK(IndexDomain::interior, te);
+
+    auto pkg = md->GetMeshPointer()->packages.Get("poisson_cell_package");
+    const auto alpha = pkg->Param<Real>("diagonal_alpha");
+
+    int nblocks = md->NumBlocks();
+    std::vector<bool> include_block(nblocks, true);
+
+    static auto desc =
+        parthenon::MakePackDescriptor<var_t>(md.get(), {}, {PDOpt::WithFluxes});
+    static auto desc_out = parthenon::MakePackDescriptor<var_t>(md_out.get());
+    auto pack = desc.GetPack(md.get(), include_block);
+    auto pack_out = desc_out.GetPack(md_out.get(), include_block);
+    parthenon::par_for(
+        "FluxMultiplyMatrix", 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s,
+        ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          const auto &coords = pack.GetCoordinates(b);
+          Real dx1 = coords.template Dxc<X1DIR>(k, j, i);
+          pack_out(b, te, var_t(), k, j, i) = -alpha * pack(b, te, var_t(), k, j, i);
+          pack_out(b, te, var_t(), k, j, i) +=
+              (pack.flux(b, X1DIR, var_t(), k, j, i) -
+               pack.flux(b, X1DIR, var_t(), k, j, i + 1)) /
+              dx1;
+
+          if (ndim > 1) {
+            Real dx2 = coords.template Dxc<X2DIR>(k, j, i);
+            pack_out(b, te, var_t(), k, j, i) +=
+                (pack.flux(b, X2DIR, var_t(), k, j, i) -
+                 pack.flux(b, X2DIR, var_t(), k, j + 1, i)) /
+                dx2;
+          }
+
+          if (ndim > 2) {
+            Real dx3 = coords.template Dxc<X3DIR>(k, j, i);
+            pack_out(b, te, var_t(), k, j, i) +=
+                (pack.flux(b, X3DIR, var_t(), k, j, i) -
+                 pack.flux(b, X3DIR, var_t(), k + 1, j, i)) /
+                dx3;
+          }
+        });
+    return TaskStatus::complete;
+  }
+};
+
+} // namespace poisson_cell_package
+
+#endif // EXAMPLE_LINEAR_SOLVERS_POISSON_CELL_EQUATION_HPP_
diff --git a/example/linear_solvers/poisson_cell_package.cpp b/example/linear_solvers/poisson_cell_package.cpp
new file mode 100644
index 0000000000000..cff1416514307
--- /dev/null
+++ b/example/linear_solvers/poisson_cell_package.cpp
@@ -0,0 +1,232 @@
+//========================================================================================
+// (C) (or copyright) 2021-2024. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+#include <algorithm>
+#include <cmath>
+#include <cstdio>
+#include <limits>
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <bvals/boundary_conditions_generic.hpp>
+#include <coordinates/coordinates.hpp>
+#include <parthenon/driver.hpp>
+#include <parthenon/package.hpp>
+#include <solvers/bicgstab_solver.hpp>
+#include <solvers/cg_solver.hpp>
+#include <solvers/solver_utils.hpp>
+
+#include "defs.hpp"
+#include "kokkos_abstraction.hpp"
+#include "poisson_cell_equation.hpp"
+#include "poisson_cell_package.hpp"
+
+using namespace parthenon::package::prelude;
+using parthenon::HostArray1D;
+namespace poisson_cell_package {
+
+using namespace parthenon;
+using namespace parthenon::BoundaryFunction;
+// We need to register FixedFace boundary conditions by hand since they can't
+// be chosen in the parameter input file. FixedFace boundary conditions assume
+// Dirichlet booundary conditions on the face of the domain and linearly extrapolate
+// into the ghosts to ensure the linear reconstruction on the block face obeys the
+// chosen boundary condition. Just setting the ghost zones of CC variables to a fixed
+// value results in poor MG convergence because the effective BC at the face
+// changes with MG level.
+
+// Build type that selects only variables within the poisson_cell namespace. Internal
+// solver variables have the namespace of input variables prepended, so they will also be
+// selected by this type.
+struct any_poisson_cell : public parthenon::variable_names::base_t<true> {
+  template <class... Ts>
+  KOKKOS_INLINE_FUNCTION any_poisson_cell(Ts &&...args)
+      : base_t<true>(std::forward<Ts>(args)...) {}
+  static std::string name() { return "poisson_cell[.].*"; }
+};
+
+template <CoordinateDirection DIR, BCSide SIDE>
+auto GetBC() {
+  return [](std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) -> void {
+    using namespace parthenon;
+    using namespace parthenon::BoundaryFunction;
+    GenericBC<DIR, SIDE, BCType::FixedFace, any_poisson_cell>(rc, coarse, 0.0);
+  };
+}
+
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
+  auto pkg = std::make_shared<StateDescriptor>("poisson_cell_package");
+
+  // Set boundary conditions for Poisson variables
+  using BF = parthenon::BoundaryFace;
+  pkg->UserBoundaryFunctions[BF::inner_x1].push_back(GetBC<X1DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::inner_x2].push_back(GetBC<X2DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::inner_x3].push_back(GetBC<X3DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::outer_x1].push_back(GetBC<X1DIR, BCSide::Outer>());
+  pkg->UserBoundaryFunctions[BF::outer_x2].push_back(GetBC<X2DIR, BCSide::Outer>());
+  pkg->UserBoundaryFunctions[BF::outer_x3].push_back(GetBC<X3DIR, BCSide::Outer>());
+
+  Real diagonal_alpha = pin->GetOrAddReal("poisson_cell", "diagonal_alpha", 0.0);
+  pkg->AddParam<>("diagonal_alpha", diagonal_alpha);
+
+  std::string solver = pin->GetOrAddString("poisson_cell", "solver", "MG");
+  pkg->AddParam<>("solver", solver);
+
+  bool use_exact_rhs = pin->GetOrAddBoolean("poisson_cell", "use_exact_rhs", false);
+  pkg->AddParam<>("use_exact_rhs", use_exact_rhs);
+
+  std::string prolong =
+      pin->GetOrAddString("poisson_cell", "boundary_prolongation", "Linear");
+
+  using PoissEq = poisson_cell_package::PoissonEquation<u, D>;
+  PoissEq eq(pin, "poisson_cell");
+  pkg->AddParam<>("poisson_cell_equation", eq, parthenon::Params::Mutability::Mutable);
+
+  std::shared_ptr<parthenon::solvers::SolverBase> psolver;
+  using prolongator_t = parthenon::solvers::ProlongationBlockInteriorDefault;
+  using preconditioner_t = parthenon::solvers::MGSolver<PoissEq, prolongator_t>;
+  if (solver == "MG") {
+    psolver = std::make_shared<parthenon::solvers::MGSolver<PoissEq, prolongator_t>>(
+        "base", u_label, rhs_label, pin, "poisson_cell/solver_params",
+        PoissEq(pin, "poisson_cell"));
+  } else if (solver == "CG") {
+    psolver = std::make_shared<parthenon::solvers::CGSolver<PoissEq, preconditioner_t>>(
+        "base", u_label, rhs_label, pin, "poisson_cell/solver_params",
+        PoissEq(pin, "poisson_cell"));
+  } else if (solver == "BiCGSTAB") {
+    psolver =
+        std::make_shared<parthenon::solvers::BiCGSTABSolver<PoissEq, preconditioner_t>>(
+            "base", u_label, rhs_label, pin, "poisson_cell/solver_params",
+            PoissEq(pin, "poisson_cell"));
+  } else {
+    PARTHENON_FAIL("Unknown solver type.");
+  }
+  pkg->AddParam<>("solver_pointer", psolver);
+
+  using namespace parthenon::refinement_ops;
+  auto mD = Metadata(
+      {Metadata::Independent, Metadata::OneCopy, Metadata::Face, Metadata::GMGRestrict});
+  mD.RegisterRefinementOps<ProlongateSharedLinear, RestrictAverage>();
+
+  // Holds the discretized version of D in \nabla \cdot D(\vec{x}) \nabla u = rhs. D = 1
+  // for the standard Poisson equation.
+  pkg->AddField(D::name(), mD);
+
+  std::vector<MetadataFlag> flags{Metadata::Cell,        Metadata::Independent,
+                                  Metadata::FillGhost,   Metadata::WithFluxes,
+                                  Metadata::GMGRestrict, Metadata::GMGProlongate};
+  auto mflux_comm = Metadata(flags);
+  if (prolong == "Linear") {
+    mflux_comm.RegisterRefinementOps<ProlongateSharedLinear, RestrictAverage>();
+  } else if (prolong == "Constant") {
+    mflux_comm.RegisterRefinementOps<ProlongatePiecewiseConstant, RestrictAverage>();
+  } else {
+    PARTHENON_FAIL("Unknown prolongation method for Poisson boundaries.");
+  }
+  // u is the solution vector that starts with an initial guess and then gets updated
+  // by the solver
+  pkg->AddField(u::name(), mflux_comm);
+
+  auto m_no_ghost = Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy});
+  // rhs is the field that contains the desired rhs side
+  pkg->AddField(rhs::name(), m_no_ghost);
+
+  // Auxillary field for storing the exact solution when it is known
+  pkg->AddField(exact::name(), m_no_ghost);
+
+  return pkg;
+}
+
+parthenon::TaskStatus
+SetVector(parthenon::ParameterInput *pin, bool use_exponential,
+          std::shared_ptr<parthenon::MeshData<parthenon::Real>> md) {
+  using namespace parthenon;
+  Real x0 = pin->GetOrAddReal("poisson_cell", "x0", 0.0);
+  Real y0 = pin->GetOrAddReal("poisson_cell", "y0", 0.0);
+  Real z0 = pin->GetOrAddReal("poisson_cell", "z0", 0.0);
+  Real radius0 = pin->GetOrAddReal("poisson_cell", "radius", 0.1);
+  const int ndim = md->GetMeshPointer()->ndim;
+
+  auto desc = MakePackDescriptor<u>(md.get());
+  auto pack = desc.GetPack(md.get());
+
+  using TE = parthenon::TopologicalElement;
+  auto ib = md->GetBoundsI(IndexDomain::entire, TE::CC);
+  auto jb = md->GetBoundsJ(IndexDomain::entire, TE::CC);
+  auto kb = md->GetBoundsK(IndexDomain::entire, TE::CC);
+
+  parthenon::par_for(
+      "PoissonCell::Ax", 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        const auto &coords = pack.GetCoordinates(b);
+        Real x1 = coords.X<1, TE::CC>(i);
+        Real x2 = coords.X<2, TE::CC>(j);
+        Real x3 = coords.X<2, TE::CC>(k);
+        Real rad = (x1 - x0) * (x1 - x0);
+        if (ndim > 1) rad += (x2 - y0) * (x2 - y0);
+        if (ndim > 2) rad += (x3 - z0) * (x3 - z0);
+        rad = std::sqrt(rad);
+
+        pack(b, TE::CC, u(), k, j, i) = rad < radius0 ? 1.0 : 0.0;
+        if (use_exponential) pack(b, TE::CC, u(), k, j, i) = -exp(-10.0 * rad * rad);
+      });
+  return TaskStatus::complete;
+}
+
+KOKKOS_FUNCTION
+bool InsideLRegion(int ndim, Real x, Real y, Real z) {
+  bool inside1 = (x < -0.25) && (x > -0.75);
+  if (ndim > 1) inside1 = inside1 && (y < 0.5) && (y > -0.5);
+  if (ndim > 2) inside1 = inside1 && (z < 0.25) && (z > -0.25);
+
+  bool inside2 = (x < 0.5) && (x > -0.75);
+  if (ndim > 1) inside2 = inside2 && (y < -0.25) && (y > -0.75);
+  if (ndim > 2) inside2 = inside2 && (z < 0.25) && (z > -0.25);
+
+  return inside1 || inside2;
+}
+
+parthenon::TaskStatus SetD(parthenon::ParameterInput *pin,
+                           std::shared_ptr<parthenon::MeshData<parthenon::Real>> md) {
+  using namespace parthenon;
+  Real interior_D = pin->GetOrAddReal("poisson_cell", "interior_D", 1.0);
+  Real exterior_D = pin->GetOrAddReal("poisson_cell", "exterior_D", 1.0);
+  const int ndim = md->GetMeshPointer()->ndim;
+
+  auto desc = MakePackDescriptor<D>(md.get());
+  auto pack = desc.GetPack(md.get());
+
+  using TE = parthenon::TopologicalElement;
+  for (auto te : {TE::F1, TE::F2, TE::F3}) {
+    auto ib = md->GetBoundsI(IndexDomain::entire, te);
+    auto jb = md->GetBoundsJ(IndexDomain::entire, te);
+    auto kb = md->GetBoundsK(IndexDomain::entire, te);
+
+    parthenon::par_for(
+        "PoissonCell::Ax", 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+        KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          const auto &coords = pack.GetCoordinates(b);
+          Real x1 = coords.X<1>(te, k, j, i);
+          Real x2 = coords.X<2>(te, k, j, i);
+          Real x3 = coords.X<3>(te, k, j, i);
+
+          pack(b, te, D(), k, j, i) =
+              InsideLRegion(ndim, x1, x2, x3) ? interior_D : exterior_D;
+        });
+  }
+  return TaskStatus::complete;
+}
+
+} // namespace poisson_cell_package
diff --git a/example/linear_solvers/poisson_cell_package.hpp b/example/linear_solvers/poisson_cell_package.hpp
new file mode 100644
index 0000000000000..6bdb5bc87d82d
--- /dev/null
+++ b/example/linear_solvers/poisson_cell_package.hpp
@@ -0,0 +1,55 @@
+//========================================================================================
+// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef EXAMPLE_LINEAR_SOLVERS_POISSON_CELL_PACKAGE_HPP_
+#define EXAMPLE_LINEAR_SOLVERS_POISSON_CELL_PACKAGE_HPP_
+
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <kokkos_abstraction.hpp>
+#include <parthenon/package.hpp>
+
+#include "linear_solver_driver.hpp"
+#include "variable_type.hpp"
+
+namespace poisson_cell_package {
+using namespace parthenon::package::prelude;
+
+VARIABLE(poisson, D);
+VARIABLE(poisson, u);
+VARIABLE(poisson, rhs);
+VARIABLE(poisson, exact);
+
+// Meshdata container labels
+inline const std::string u_label = "cell_u";
+inline const std::string rhs_label = "cell_rhs";
+inline const std::string exact_label = "cell_exact";
+
+// This just provides a convenient short hand for TE::CC and will make it
+// easier for testing solves with different topological elements in the
+// future (although other types of fields require significantly different
+// condition boundary implementations)
+constexpr parthenon::TopologicalElement te = parthenon::TopologicalElement::CC;
+
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
+parthenon::TaskStatus SetVector(parthenon::ParameterInput *pin, bool use_exponential,
+                                std::shared_ptr<parthenon::MeshData<parthenon::Real>> md);
+parthenon::TaskStatus SetD(parthenon::ParameterInput *pin,
+                           std::shared_ptr<parthenon::MeshData<parthenon::Real>> md);
+void AddTaskRegion(parthenon::TaskCollection &tc,
+                   linear_solver_example::LinearSolverDriver *driver);
+} // namespace poisson_cell_package
+
+#endif // EXAMPLE_LINEAR_SOLVERS_POISSON_CELL_PACKAGE_HPP_
diff --git a/example/linear_solvers/poisson_nodal_equation.hpp b/example/linear_solvers/poisson_nodal_equation.hpp
new file mode 100644
index 0000000000000..f99d97c4008ed
--- /dev/null
+++ b/example/linear_solvers/poisson_nodal_equation.hpp
@@ -0,0 +1,155 @@
+//========================================================================================
+// (C) (or copyright) 2023-2024. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef EXAMPLE_LINEAR_SOLVERS_POISSON_NODAL_EQUATION_HPP_
+#define EXAMPLE_LINEAR_SOLVERS_POISSON_NODAL_EQUATION_HPP_
+
+#include <memory>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <kokkos_abstraction.hpp>
+#include <parthenon/package.hpp>
+
+namespace poisson_nodal_package {
+using namespace parthenon::package::prelude;
+// This class implement methods for calculating A.x = y and returning the diagonal of A,
+// where A is the the matrix representing the discretized Poisson equation on the grid.
+// Here we implement the Laplace operator in terms of a flux divergence to (potentially)
+// consistently deal with coarse fine boundaries on the grid. Only the routines Ax and
+// SetDiagonal need to be defined for interfacing this with solvers. The other methods
+// are internal, but can't be marked private or protected because they launch kernels
+// on device.
+template <class var_t>
+class PoissonEquation {
+ public:
+  using IndependentVars = parthenon::TypeList<var_t>;
+
+  PoissonEquation(parthenon::ParameterInput *pin, const std::string &label) {}
+
+  parthenon::TaskID Ax(parthenon::TaskList &tl, parthenon::TaskID depends_on,
+                       std::shared_ptr<parthenon::MeshData<Real>> & /*md_mat*/,
+                       std::shared_ptr<parthenon::MeshData<Real>> &md_in,
+                       std::shared_ptr<parthenon::MeshData<Real>> &md_out) {
+    return tl.AddTask(depends_on, AxImpl, md_in, md_out);
+  }
+
+  static parthenon::TaskStatus
+  AxImpl(std::shared_ptr<parthenon::MeshData<Real>> &md_in,
+         std::shared_ptr<parthenon::MeshData<Real>> &md_out) {
+    using namespace parthenon;
+    auto pkg = md_in->GetMeshPointer()->packages.Get("poisson_nodal_package");
+    const auto alpha = pkg->Param<Real>("diagonal_alpha");
+
+    constexpr auto te = TopologicalElement::NN;
+    const int ndim = md_in->GetMeshPointer()->ndim;
+    IndexRange ib = md_in->GetBoundsI(IndexDomain::interior, te);
+    IndexRange jb = md_in->GetBoundsJ(IndexDomain::interior, te);
+    IndexRange kb = md_in->GetBoundsK(IndexDomain::interior, te);
+
+    auto desc = parthenon::MakePackDescriptor<var_t>(md_in.get());
+    auto pack_in = desc.GetPack(md_in.get());
+    auto pack_out = desc.GetPack(md_out.get());
+
+    parthenon::par_for(
+        "PoissonNodal::Ax", 0, pack_in.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s,
+        ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          const auto &coords = pack_in.GetCoordinates(b);
+          const Real dx1 = coords.template Dxc<X1DIR>(k, j, i);
+          const Real dx2 = coords.template Dxc<X2DIR>(k, j, i);
+          const Real dx3 = coords.template Dxc<X3DIR>(k, j, i);
+
+          Real Ax = alpha * pack_in(b, te, var_t(), k, j, i);
+          Ax += pack_in(b, te, var_t(), k, j, i) * 2.0 / (dx1 * dx1);
+          Ax -= (pack_in(b, te, var_t(), k, j, i + 1) +
+                 pack_in(b, te, var_t(), k, j, i - 1)) /
+                (dx1 * dx1);
+          if (ndim > 1) {
+            Ax += pack_in(b, te, var_t(), k, j, i) * 2.0 / (dx2 * dx2);
+            Ax -= (pack_in(b, te, var_t(), k, j + 1, i) +
+                   pack_in(b, te, var_t(), k, j - 1, i)) /
+                  (dx2 * dx2);
+          }
+          if (ndim > 2) {
+            Ax += pack_in(b, te, var_t(), k, j, i) * 2.0 / (dx3 * dx3);
+            Ax -= (pack_in(b, te, var_t(), k + 1, j, i) +
+                   pack_in(b, te, var_t(), k - 1, j, i)) /
+                  (dx3 * dx3);
+          }
+          pack_out(b, te, var_t(), k, j, i) = Ax;
+        });
+    return TaskStatus::complete;
+  }
+
+  static parthenon::TaskStatus SetBoundary(std::shared_ptr<parthenon::MeshData<Real>> &md,
+                                           bool coarse) {
+    using namespace parthenon;
+
+    constexpr auto te = TopologicalElement::NN;
+    const int ndim = md->GetMeshPointer()->ndim;
+    CellLevel cl = coarse ? CellLevel::coarse : CellLevel::same;
+    IndexRange ib = md->GetBoundsI(cl, IndexDomain::interior, te);
+    IndexRange jb = md->GetBoundsJ(cl, IndexDomain::interior, te);
+    IndexRange kb = md->GetBoundsK(cl, IndexDomain::interior, te);
+
+    std::set<PDOpt> opts{};
+    if (coarse) opts.emplace(PDOpt::Coarse);
+    auto desc = parthenon::MakePackDescriptor<var_t>(md.get(), {}, opts);
+    auto pack = desc.GetPack(md.get(), GetBlockSelector::OnPhysicalBoundary());
+
+    parthenon::par_for(
+        "PoissonNodal::SetBoundary", 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e,
+        ib.s, ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          const int oi = TopologicalOffsetI(te) * ((ib.e == i) - (ib.s == i));
+          const int oj = TopologicalOffsetJ(te) * ((jb.e == j) - (jb.s == j));
+          const int ok = TopologicalOffsetK(te) * ((kb.e == k) - (kb.s == k));
+          if (pack.IsPhysicalBoundary(b, ok, oj, oi)) pack(b, te, var_t(), k, j, i) = 0.0;
+        });
+    return TaskStatus::complete;
+  }
+
+  parthenon::TaskStatus
+  SetDiagonal(std::shared_ptr<parthenon::MeshData<Real>> & /*md_mat*/,
+              std::shared_ptr<parthenon::MeshData<Real>> &md_diag) {
+    using namespace parthenon;
+    const int ndim = md_diag->GetMeshPointer()->ndim;
+    constexpr auto te = TopologicalElement::NN;
+    IndexRange ib = md_diag->GetBoundsI(IndexDomain::interior, te);
+    IndexRange jb = md_diag->GetBoundsJ(IndexDomain::interior, te);
+    IndexRange kb = md_diag->GetBoundsK(IndexDomain::interior, te);
+
+    auto pkg = md_diag->GetMeshPointer()->packages.Get("poisson_nodal_package");
+    const auto alpha = pkg->Param<Real>("diagonal_alpha");
+
+    auto desc = parthenon::MakePackDescriptor<var_t>(md_diag.get());
+    auto pack_diag = desc.GetPack(md_diag.get());
+    parthenon::par_for(
+        "StoreDiagonal", 0, pack_diag.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s,
+        ib.e, KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+          const auto &coords = pack_diag.GetCoordinates(b);
+          // Build the diagonal of the matrix
+          Real dx1 = coords.template Dxc<X1DIR>(k, j, i);
+          Real dx2 = coords.template Dxc<X2DIR>(k, j, i);
+          Real dx3 = coords.template Dxc<X3DIR>(k, j, i);
+          pack_diag(b, te, var_t(), k, j, i) = alpha + 2.0 / (dx1 * dx1) +
+                                               (ndim > 1) * 2.0 / (dx2 * dx2) +
+                                               (ndim > 2) * 2.0 / (dx3 * dx3);
+        });
+    return TaskStatus::complete;
+  }
+};
+
+} // namespace poisson_nodal_package
+
+#endif // EXAMPLE_LINEAR_SOLVERS_POISSON_NODAL_EQUATION_HPP_
diff --git a/example/linear_solvers/poisson_nodal_package.cpp b/example/linear_solvers/poisson_nodal_package.cpp
new file mode 100644
index 0000000000000..811b38502748b
--- /dev/null
+++ b/example/linear_solvers/poisson_nodal_package.cpp
@@ -0,0 +1,186 @@
+//========================================================================================
+// (C) (or copyright) 2021-2024. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+#include <algorithm>
+#include <cmath>
+#include <cstdio>
+#include <limits>
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <bvals/boundary_conditions_generic.hpp>
+#include <coordinates/coordinates.hpp>
+#include <parthenon/driver.hpp>
+#include <parthenon/package.hpp>
+#include <solvers/bicgstab_solver.hpp>
+#include <solvers/cg_solver.hpp>
+#include <solvers/solver_utils.hpp>
+
+#include "defs.hpp"
+#include "kokkos_abstraction.hpp"
+#include "linear_solver_driver.hpp"
+#include "poisson_nodal_equation.hpp"
+#include "poisson_nodal_package.hpp"
+
+using namespace parthenon::package::prelude;
+using parthenon::HostArray1D;
+namespace poisson_nodal_package {
+
+using namespace parthenon;
+using namespace parthenon::BoundaryFunction;
+// We need to register FixedFace boundary conditions by hand since they can't
+// be chosen in the parameter input file. FixedFace boundary conditions assume
+// Dirichlet booundary conditions on the face of the domain and linearly extrapolate
+// into the ghosts to ensure the linear reconstruction on the block face obeys the
+// chosen boundary condition. Just setting the ghost zones of CC variables to a fixed
+// value results in poor MG convergence because the effective BC at the face
+// changes with MG level.
+
+// Build type that selects only variables within the poisson_nodal namespace. Internal
+// solver variables have the namespace of input variables prepended, so they will also be
+// selected by this type.
+struct any_poisson_nodal : public parthenon::variable_names::base_t<true> {
+  template <class... Ts>
+  KOKKOS_INLINE_FUNCTION any_poisson_nodal(Ts &&...args)
+      : base_t<true>(std::forward<Ts>(args)...) {}
+  static std::string name() { return "poisson_nodal[.].*"; }
+};
+
+template <CoordinateDirection DIR, BCSide SIDE>
+auto GetBC() {
+  return [](std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) -> void {
+    using namespace parthenon;
+    using namespace parthenon::BoundaryFunction;
+    GenericBC<DIR, SIDE, BCType::FixedFace, any_poisson_nodal>(rc, coarse, 0.0);
+  };
+}
+
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
+  auto pkg = std::make_shared<StateDescriptor>("poisson_nodal_package");
+
+  // Set boundary conditions for Poisson variables
+  using BF = parthenon::BoundaryFace;
+  pkg->UserBoundaryFunctions[BF::inner_x1].push_back(GetBC<X1DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::inner_x2].push_back(GetBC<X2DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::inner_x3].push_back(GetBC<X3DIR, BCSide::Inner>());
+  pkg->UserBoundaryFunctions[BF::outer_x1].push_back(GetBC<X1DIR, BCSide::Outer>());
+  pkg->UserBoundaryFunctions[BF::outer_x2].push_back(GetBC<X2DIR, BCSide::Outer>());
+  pkg->UserBoundaryFunctions[BF::outer_x3].push_back(GetBC<X3DIR, BCSide::Outer>());
+
+  Real diagonal_alpha = pin->GetOrAddReal("poisson_nodal", "diagonal_alpha", 0.0);
+  pkg->AddParam<>("diagonal_alpha", diagonal_alpha);
+
+  std::string solver = pin->GetOrAddString("poisson_nodal", "solver", "MG");
+  pkg->AddParam<>("solver", solver);
+
+  bool use_exact_rhs = pin->GetOrAddBoolean("poisson_nodal", "use_exact_rhs", false);
+  pkg->AddParam<>("use_exact_rhs", use_exact_rhs);
+
+  std::string prolong =
+      pin->GetOrAddString("poisson_nodal", "boundary_prolongation", "Linear");
+
+  using PoissEq = poisson_nodal_package::PoissonEquation<u>;
+  PoissEq eq(pin, "poisson_nodal");
+  pkg->AddParam<>("poisson_nodal_equation", eq, parthenon::Params::Mutability::Mutable);
+
+  std::shared_ptr<parthenon::solvers::SolverBase> psolver;
+  using prolongator_t = parthenon::solvers::ProlongationBlockInteriorDefault;
+  using preconditioner_t = parthenon::solvers::MGSolver<PoissEq, prolongator_t>;
+
+  const std::string base_label = "base";
+  const std::string u_label = "nodal_u";
+  const std::string rhs_label = "nodal_rhs";
+  if (solver == "MG") {
+    psolver = std::make_shared<parthenon::solvers::MGSolver<PoissEq, prolongator_t>>(
+        base_label, u_label, rhs_label, pin, "poisson_nodal/solver_params",
+        PoissEq(pin, "poisson_nodal"));
+  } else if (solver == "CG") {
+    psolver = std::make_shared<parthenon::solvers::CGSolver<PoissEq, preconditioner_t>>(
+        base_label, u_label, rhs_label, pin, "poisson_nodal/solver_params",
+        PoissEq(pin, "poisson_nodal"));
+  } else if (solver == "BiCGSTAB") {
+    psolver =
+        std::make_shared<parthenon::solvers::BiCGSTABSolver<PoissEq, preconditioner_t>>(
+            base_label, u_label, rhs_label, pin, "poisson_nodal/solver_params",
+            PoissEq(pin, "poisson_nodal"));
+  } else {
+    PARTHENON_FAIL("Unknown solver type.");
+  }
+  pkg->AddParam<>("solver_pointer", psolver);
+
+  using namespace parthenon::refinement_ops;
+
+  std::vector<MetadataFlag> flags{Metadata::Node,        Metadata::Independent,
+                                  Metadata::FillGhost,   Metadata::WithFluxes,
+                                  Metadata::GMGRestrict, Metadata::GMGProlongate};
+  auto mflux_comm = Metadata(flags);
+  if (prolong == "Linear") {
+    mflux_comm.RegisterRefinementOps<ProlongateSharedLinear, RestrictAverage>();
+  } else if (prolong == "Constant") {
+    mflux_comm.RegisterRefinementOps<ProlongatePiecewiseConstant, RestrictAverage>();
+  } else {
+    PARTHENON_FAIL("Unknown prolongation method for Poisson boundaries.");
+  }
+  // u is the solution vector that starts with an initial guess and then gets updated
+  // by the solver
+  pkg->AddField(u::name(), mflux_comm);
+
+  auto m_no_ghost = Metadata({Metadata::Node, Metadata::Derived, Metadata::OneCopy});
+  // rhs is the field that contains the desired rhs side
+  pkg->AddField(rhs::name(), m_no_ghost);
+
+  // Auxillary field for storing the exact solution when it is known
+  pkg->AddField(exact::name(), m_no_ghost);
+
+  return pkg;
+}
+
+parthenon::TaskStatus
+SetVector(parthenon::ParameterInput *pin, bool use_exponential,
+          std::shared_ptr<parthenon::MeshData<parthenon::Real>> md) {
+  using namespace parthenon;
+  Real x0 = pin->GetOrAddReal("poisson_nodal", "x0", 0.0);
+  Real y0 = pin->GetOrAddReal("poisson_nodal", "y0", 0.0);
+  Real z0 = pin->GetOrAddReal("poisson_nodal", "z0", 0.0);
+  Real radius0 = pin->GetOrAddReal("poisson_nodal", "radius", 0.1);
+  const int ndim = md->GetMeshPointer()->ndim;
+
+  auto desc = MakePackDescriptor<u>(md.get());
+  auto pack = desc.GetPack(md.get());
+
+  using TE = parthenon::TopologicalElement;
+  auto ib = md->GetBoundsI(IndexDomain::entire, TE::NN);
+  auto jb = md->GetBoundsJ(IndexDomain::entire, TE::NN);
+  auto kb = md->GetBoundsK(IndexDomain::entire, TE::NN);
+
+  parthenon::par_for(
+      "PoissonNodal::Ax", 0, pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        const auto &coords = pack.GetCoordinates(b);
+        Real x1 = coords.X<1, TE::NN>(i);
+        Real x2 = coords.X<2, TE::NN>(j);
+        Real x3 = coords.X<2, TE::NN>(k);
+        Real rad = (x1 - x0) * (x1 - x0);
+        if (ndim > 1) rad += (x2 - y0) * (x2 - y0);
+        if (ndim > 2) rad += (x3 - z0) * (x3 - z0);
+        rad = std::sqrt(rad);
+
+        pack(b, TE::NN, u(), k, j, i) = rad < radius0 ? 1.0 : 0.0;
+        if (use_exponential) pack(b, TE::NN, u(), k, j, i) = -exp(-10.0 * rad * rad);
+      });
+  return TaskStatus::complete;
+}
+
+} // namespace poisson_nodal_package
diff --git a/example/linear_solvers/poisson_nodal_package.hpp b/example/linear_solvers/poisson_nodal_package.hpp
new file mode 100644
index 0000000000000..a37396798caf4
--- /dev/null
+++ b/example/linear_solvers/poisson_nodal_package.hpp
@@ -0,0 +1,38 @@
+//========================================================================================
+// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef EXAMPLE_LINEAR_SOLVERS_POISSON_NODAL_PACKAGE_HPP_
+#define EXAMPLE_LINEAR_SOLVERS_POISSON_NODAL_PACKAGE_HPP_
+
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <kokkos_abstraction.hpp>
+#include <parthenon/package.hpp>
+
+#include "variable_type.hpp"
+
+namespace poisson_nodal_package {
+using namespace parthenon::package::prelude;
+
+VARIABLE(poisson_nodal, u);
+VARIABLE(poisson_nodal, rhs);
+VARIABLE(poisson_nodal, exact);
+
+std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
+parthenon::TaskStatus SetVector(parthenon::ParameterInput *pin, bool use_exponential,
+                                std::shared_ptr<parthenon::MeshData<parthenon::Real>> md);
+} // namespace poisson_nodal_package
+
+#endif // EXAMPLE_LINEAR_SOLVERS_POISSON_NODAL_PACKAGE_HPP_
diff --git a/example/linear_solvers/variable_type.hpp b/example/linear_solvers/variable_type.hpp
new file mode 100644
index 0000000000000..1cd97a2b71d50
--- /dev/null
+++ b/example/linear_solvers/variable_type.hpp
@@ -0,0 +1,27 @@
+//========================================================================================
+// (C) (or copyright) 2021-2024. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef EXAMPLE_LINEAR_SOLVERS_VARIABLE_TYPE_HPP_
+#define EXAMPLE_LINEAR_SOLVERS_VARIABLE_TYPE_HPP_
+
+#include <string>
+#include <utility>
+
+#define VARIABLE(ns, varname)                                                            \
+  struct varname : public parthenon::variable_names::base_t<false> {                     \
+    template <class... Ts>                                                               \
+    KOKKOS_INLINE_FUNCTION varname(Ts &&...args)                                         \
+        : parthenon::variable_names::base_t<false>(std::forward<Ts>(args)...) {}         \
+    static std::string name() { return #ns "." #varname; }                               \
+  }
+
+#endif // EXAMPLE_LINEAR_SOLVERS_VARIABLE_TYPE_HPP_
diff --git a/src/basic_types.hpp b/src/basic_types.hpp
index 6f930f3b19f1f..38ac10f98bc62 100644
--- a/src/basic_types.hpp
+++ b/src/basic_types.hpp
@@ -23,6 +23,7 @@
 #include <Kokkos_Core.hpp>
 
 #include "config.hpp"
+#include "utils/error_checking.hpp"
 
 namespace parthenon {
 
@@ -166,7 +167,7 @@ enum class TopologicalElement : std::size_t {
   E3 = 8,
   NN = 9
 };
-enum class TopologicalType { Cell, Face, Edge, Node };
+enum class TopologicalType : std::size_t { Cell = 0, Face = 3, Edge = 6, Node = 9 };
 
 KOKKOS_FORCEINLINE_FUNCTION
 constexpr TopologicalType GetTopologicalType(TopologicalElement el) {
@@ -183,6 +184,13 @@ constexpr TopologicalType GetTopologicalType(TopologicalElement el) {
   }
 }
 
+KOKKOS_FORCEINLINE_FUNCTION
+constexpr std::size_t GetNumberOfElements(TopologicalType tt) {
+  using TT = TopologicalType;
+  if (tt == TT::Cell || tt == TT::Node) return 1;
+  return 3;
+}
+
 inline std::vector<TopologicalElement> GetTopologicalElements(TopologicalType tt) {
   using TE = TopologicalElement;
   using TT = TopologicalType;
@@ -212,16 +220,30 @@ TopologicalElement GetTopologicalElementInDir(const TopologicalType tt,
 using TE = TopologicalElement;
 // Returns one if the I coordinate of el is offset from the zone center coordinates,
 // and zero otherwise
-inline constexpr int TopologicalOffsetI(TE el) {
+KOKKOS_FORCEINLINE_FUNCTION
+constexpr int TopologicalOffsetI(TE el) {
   return (el == TE::F1 || el == TE::E2 || el == TE::E3 || el == TE::NN);
 }
-inline constexpr int TopologicalOffsetJ(TE el) {
+KOKKOS_FORCEINLINE_FUNCTION
+constexpr int TopologicalOffsetJ(TE el) {
   return (el == TE::F2 || el == TE::E3 || el == TE::E1 || el == TE::NN);
 }
-inline constexpr int TopologicalOffsetK(TE el) {
+KOKKOS_FORCEINLINE_FUNCTION
+constexpr int TopologicalOffsetK(TE el) {
   return (el == TE::F3 || el == TE::E2 || el == TE::E1 || el == TE::NN);
 }
 
+KOKKOS_FORCEINLINE_FUNCTION
+constexpr int TopologicalOffset(CoordinateDirection dir, TE el) {
+  if (X1DIR == dir)
+    return TopologicalOffsetI(el);
+  else if (X2DIR == dir)
+    return TopologicalOffsetJ(el);
+  else if (X3DIR == dir)
+    return TopologicalOffsetK(el);
+  return 0;
+}
+
 // Returns wether or not topological element containee is a boundary of
 // topological element container
 inline constexpr bool IsSubmanifold(TopologicalElement containee,
diff --git a/src/coordinates/uniform_coordinates.hpp b/src/coordinates/uniform_coordinates.hpp
index c3b63fac326e3..6dddf3b0edf05 100644
--- a/src/coordinates/uniform_coordinates.hpp
+++ b/src/coordinates/uniform_coordinates.hpp
@@ -191,6 +191,30 @@ class UniformCoordinates {
     return 0.0;
   }
 
+  template <int dir>
+  KOKKOS_FORCEINLINE_FUNCTION Real X(TopologicalElement el, const int k, const int j,
+                                     const int i) const {
+    static_assert(dir > 0 && dir < 4);
+    using TE = TopologicalElement;
+    if (el == TE::CC)
+      return X<dir, TE::CC>(k, j, i);
+    else if (el == TE::F1)
+      return X<dir, TE::F1>(k, j, i);
+    else if (el == TE::F2)
+      return X<dir, TE::F2>(k, j, i);
+    else if (el == TE::F3)
+      return X<dir, TE::F3>(k, j, i);
+    else if (el == TE::E1)
+      return X<dir, TE::E1>(k, j, i);
+    else if (el == TE::E2)
+      return X<dir, TE::E2>(k, j, i);
+    else if (el == TE::E3)
+      return X<dir, TE::E3>(k, j, i);
+    else if (el == TE::NN)
+      return X<dir, TE::NN>(k, j, i);
+    return 0.0;
+  }
+
   template <int dir, TopologicalElement el>
   KOKKOS_FORCEINLINE_FUNCTION Real Scale(const int k, const int j, const int i) const {
     if constexpr (dir > 0 && dir < 4) return 1.0;
diff --git a/src/interface/mesh_data.hpp b/src/interface/mesh_data.hpp
index 0a81d91548058..2e4cc843720f9 100644
--- a/src/interface/mesh_data.hpp
+++ b/src/interface/mesh_data.hpp
@@ -238,6 +238,13 @@ class MeshData {
     return IndexRange{-1, -2};
   }
 
+  template <class... Ts>
+  IndexShape GetCellBounds(Ts &&...args) const {
+    if (block_data_.size() > 0)
+      return block_data_[0]->GetBlockPointer()->GetCellBounds(std::forward<Ts>(args)...);
+    return IndexShape();
+  }
+
   template <class... Args>
   void Add(Args &&...args) {
     for (const auto &pbd : block_data_) {
diff --git a/src/mesh/domain.hpp b/src/mesh/domain.hpp
index 65a869b7a39ee..f752b7f3d7628 100644
--- a/src/mesh/domain.hpp
+++ b/src/mesh/domain.hpp
@@ -64,9 +64,20 @@ namespace parthenon {
 // containing the ghost zones above
 // and below the interior
 //
+// For non-cell centered topoplogical types, the
+// inner_* and outer_* domains include the outer
+// layer of elements that are shared between
+// blocks or sit on the boundary surface of the
+// domain. unshared_interior runs over only the
+// positions that are not shared with another
+// block or with the physical boundary of the
+// domain.
 enum class IndexDomain {
   entire,
   interior,
+  unshared_interior,
+  inner,
+  outer,
   inner_x1,
   outer_x1,
   inner_x2,
@@ -80,6 +91,32 @@ bool DomainTouchesOuterGhosts(const IndexDomain domain) {
          (domain == IndexDomain::outer_x2) || (domain == IndexDomain::outer_x3);
 }
 
+// inner_x1 and outer_x1 imply IndexDomain::entire for directions not
+// equal to X1DIR, encode that here
+KOKKOS_FORCEINLINE_FUNCTION
+IndexDomain DecayIndexDomain(CoordinateDirection dir, IndexDomain domain) noexcept {
+  if (domain == IndexDomain::inner_x1) {
+    if (dir == X1DIR) return IndexDomain::inner;
+    return IndexDomain::entire;
+  } else if (domain == IndexDomain::outer_x1) {
+    if (dir == X1DIR) return IndexDomain::outer;
+    return IndexDomain::entire;
+  } else if (domain == IndexDomain::inner_x2) {
+    if (dir == X2DIR) return IndexDomain::inner;
+    return IndexDomain::entire;
+  } else if (domain == IndexDomain::outer_x2) {
+    if (dir == X2DIR) return IndexDomain::outer;
+    return IndexDomain::entire;
+  } else if (domain == IndexDomain::inner_x3) {
+    if (dir == X3DIR) return IndexDomain::inner;
+    return IndexDomain::entire;
+  } else if (domain == IndexDomain::outer_x3) {
+    if (dir == X3DIR) return IndexDomain::outer;
+    return IndexDomain::entire;
+  }
+  return domain;
+}
+
 //! \class IndexVolume
 //  \brief Defines the dimensions of a shape of indices
 //
@@ -164,6 +201,40 @@ class IndexShape {
     return out;
   }
 
+  KOKKOS_INLINE_FUNCTION int GetStartIdx(CoordinateDirection dir, IndexDomain domain,
+                                         TE el = TE::CC) const noexcept {
+    const int idx = dir - 1;
+    domain = DecayIndexDomain(dir, domain);
+    switch (domain) {
+    case IndexDomain::interior:
+      return x_[idx].s;
+    case IndexDomain::unshared_interior:
+      return x_[idx].s + TopologicalOffset(dir, el);
+    case IndexDomain::outer:
+      return entire_ncells_[idx] == 1 ? 0 : x_[idx].e + 1;
+    default:
+      return 0;
+    }
+  }
+
+  KOKKOS_INLINE_FUNCTION int GetEndIdx(CoordinateDirection dir, IndexDomain domain,
+                                       TE el = TE::CC) const noexcept {
+    const int idx = dir - 1;
+    domain = DecayIndexDomain(dir, domain);
+    switch (domain) {
+    case IndexDomain::interior:
+      return entire_ncells_[idx] == 1 ? 0 : x_[idx].e + TopologicalOffset(dir, el);
+    case IndexDomain::unshared_interior:
+      return entire_ncells_[idx] == 1 ? 0 : x_[idx].e;
+    case IndexDomain::inner:
+      return x_[idx].s == 0 ? 0 : x_[idx].s - 1 + TopologicalOffset(dir, el);
+    default:
+      return entire_ncells_[idx] == 1
+                 ? 0
+                 : entire_ncells_[idx] - 1 + TopologicalOffset(dir, el);
+    }
+  }
+
   KOKKOS_INLINE_FUNCTION const IndexRange GetBoundsI(const IndexDomain &domain,
                                                      TE el = TE::CC) const noexcept {
     return (domain == IndexDomain::interior && el == TE::CC)
@@ -187,74 +258,32 @@ class IndexShape {
 
   KOKKOS_INLINE_FUNCTION int is(const IndexDomain &domain,
                                 TE el = TE::CC) const noexcept {
-    switch (domain) {
-    case IndexDomain::interior:
-      return x_[0].s;
-    case IndexDomain::outer_x1:
-      return entire_ncells_[0] == 1 ? 0 : x_[0].e + 1 + TopologicalOffsetI(el);
-    default:
-      return 0;
-    }
+    return GetStartIdx(X1DIR, domain, el);
   }
 
   KOKKOS_INLINE_FUNCTION int js(const IndexDomain &domain,
                                 TE el = TE::CC) const noexcept {
-    switch (domain) {
-    case IndexDomain::interior:
-      return x_[1].s;
-    case IndexDomain::outer_x2:
-      return entire_ncells_[1] == 1 ? 0 : x_[1].e + 1 + TopologicalOffsetJ(el);
-    default:
-      return 0;
-    }
+    return GetStartIdx(X2DIR, domain, el);
   }
 
   KOKKOS_INLINE_FUNCTION int ks(const IndexDomain &domain,
                                 TE el = TE::CC) const noexcept {
-    switch (domain) {
-    case IndexDomain::interior:
-      return x_[2].s;
-    case IndexDomain::outer_x3:
-      return entire_ncells_[2] == 1 ? 0 : x_[2].e + 1 + TopologicalOffsetK(el);
-    default:
-      return 0;
-    }
+    return GetStartIdx(X3DIR, domain, el);
   }
 
   KOKKOS_INLINE_FUNCTION int ie(const IndexDomain &domain,
                                 TE el = TE::CC) const noexcept {
-    switch (domain) {
-    case IndexDomain::interior:
-      return entire_ncells_[0] == 1 ? 0 : x_[0].e + TopologicalOffsetI(el);
-    case IndexDomain::inner_x1:
-      return x_[0].s == 0 ? 0 : x_[0].s - 1;
-    default:
-      return entire_ncells_[0] == 1 ? 0 : entire_ncells_[0] - 1 + TopologicalOffsetI(el);
-    }
+    return GetEndIdx(X1DIR, domain, el);
   }
 
   KOKKOS_INLINE_FUNCTION int je(const IndexDomain &domain,
                                 TE el = TE::CC) const noexcept {
-    switch (domain) {
-    case IndexDomain::interior:
-      return entire_ncells_[1] == 1 ? 0 : x_[1].e + TopologicalOffsetJ(el);
-    case IndexDomain::inner_x2:
-      return x_[1].s == 0 ? 0 : x_[1].s - 1;
-    default:
-      return entire_ncells_[1] == 1 ? 0 : entire_ncells_[1] - 1 + TopologicalOffsetJ(el);
-    }
+    return GetEndIdx(X2DIR, domain, el);
   }
 
   KOKKOS_INLINE_FUNCTION int ke(const IndexDomain &domain,
                                 TE el = TE::CC) const noexcept {
-    switch (domain) {
-    case IndexDomain::interior:
-      return entire_ncells_[2] == 1 ? 0 : x_[2].e + TopologicalOffsetK(el);
-    case IndexDomain::inner_x3:
-      return x_[2].s == 0 ? 0 : x_[2].s - 1;
-    default:
-      return entire_ncells_[2] == 1 ? 0 : entire_ncells_[2] - 1 + TopologicalOffsetK(el);
-    }
+    return GetEndIdx(X3DIR, domain, el);
   }
 
   KOKKOS_INLINE_FUNCTION int ncellsi(const IndexDomain &domain,
diff --git a/src/mesh/mesh-gmg.cpp b/src/mesh/mesh-gmg.cpp
index e51ad8e5daec9..4d14a61dc7e7d 100644
--- a/src/mesh/mesh-gmg.cpp
+++ b/src/mesh/mesh-gmg.cpp
@@ -58,6 +58,9 @@ void Mesh::SetMeshBlockNeighbors(
     const auto &loc = pmb->loc;
     auto neighbors = forest.FindNeighbors(loc, grid_id);
 
+    // Set this blocks ownership
+    pmb->ownership = DetermineOwnership(loc, neighbors, newly_refined);
+
     // Build NeighborBlocks for unique neighbors
     for (const auto &nloc : neighbors) {
       auto gid = forest.GetGid(nloc.global_loc);
diff --git a/src/mesh/meshblock.hpp b/src/mesh/meshblock.hpp
index ddb35d16926d2..a6f7ffd401267 100644
--- a/src/mesh/meshblock.hpp
+++ b/src/mesh/meshblock.hpp
@@ -184,6 +184,8 @@ class MeshBlock : public std::enable_shared_from_this<MeshBlock> {
   std::vector<NeighborBlock> gmg_finer_neighbors;
   std::vector<NeighborBlock> gmg_leaf_neighbors;
 
+  block_ownership_t ownership;
+
   BoundaryFlag boundary_flag[6];
 
   bool IsPhysicalBoundary(BoundaryFace bf) const {
diff --git a/src/pack/sparse_pack.hpp b/src/pack/sparse_pack.hpp
index d1d6dc73dfd6d..98008f39de10c 100644
--- a/src/pack/sparse_pack.hpp
+++ b/src/pack/sparse_pack.hpp
@@ -158,13 +158,18 @@ class SparsePack : public SparsePackBase {
            physical_bnd_flag;
   }
 
-  KOKKOS_INLINE_FUNCTION int GetGID(const int b) const { return block_props_(b, 27); }
+  KOKKOS_INLINE_FUNCTION bool IsOwned(const int b, const int off3, const int off2,
+                                      const int off1) const {
+    return block_props_(b, (off1 + 1) + 3 * ((off2 + 1) + 3 * (off3 + 1)) + 27) == 1;
+  }
+
+  KOKKOS_INLINE_FUNCTION int GetGID(const int b) const { return block_props_(b, 54); }
 
   int GetLevelHost(const int b, const int off3, const int off2, const int off1) const {
     return block_props_h_(b, (off1 + 1) + 3 * ((off2 + 1) + 3 * (off3 + 1)));
   }
 
-  int GetGIDHost(const int b) const { return block_props_h_(b, 27); }
+  int GetGIDHost(const int b) const { return block_props_h_(b, 54); }
 
   // Number of components of a variable on a block
   template <typename T>
@@ -223,6 +228,10 @@ class SparsePack : public SparsePackBase {
     return (... && ContainsHost(b, Args()));
   }
 
+  KOKKOS_INLINE_FUNCTION auto &GetTopologicalType(const int b, const int idx) const {
+    return pack_(0, b, idx).topological_type;
+  }
+
   // Informational
   auto LabelHost(int b, int idx) const { return pack_h_(0, b, idx).label(); }
   template <typename... Vars>
@@ -230,9 +239,18 @@ class SparsePack : public SparsePackBase {
 
   // operator() overloads
   using TE = TopologicalElement;
+  KOKKOS_INLINE_FUNCTION auto &operator()(const int b, const int el_idx,
+                                          const int idx) const {
+    PARTHENON_DEBUG_REQUIRE(
+        el_idx < GetNumberOfElements(pack_(0, b, idx).topological_type),
+        "Asking for an element index that doesn't exist for this TT.");
+    return pack_(el_idx, b, idx);
+  }
+
   KOKKOS_INLINE_FUNCTION auto &operator()(const int b, const TE el, const int idx) const {
     return pack_(static_cast<int>(el) % 3, b, idx);
   }
+
   KOKKOS_INLINE_FUNCTION auto &operator()(const int b, const int idx) const {
     PARTHENON_DEBUG_REQUIRE(pack_(0, b, idx).topological_type == TopologicalType::Cell,
                             "Suppressed topological element index assumes that this is a "
diff --git a/src/pack/sparse_pack_base.cpp b/src/pack/sparse_pack_base.cpp
index c98de31535daa..83aaeca7201e8 100644
--- a/src/pack/sparse_pack_base.cpp
+++ b/src/pack/sparse_pack_base.cpp
@@ -158,7 +158,7 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc,
 
   // This array stores refinement levels of current block and all neighboring blocks.
   const Indexer3D bp_idxer({-1, 1}, {-1, 1}, {-1, 1});
-  pack.block_props_ = block_props_t("block_props", nblocks, bp_idxer.size() + 1);
+  pack.block_props_ = block_props_t("block_props", nblocks, 2 * bp_idxer.size() + 1);
   pack.block_props_h_ = Kokkos::create_mirror_view(pack.block_props_);
 
   pack.coords_ = coords_t(ViewOfViewAlloc("coords"), desc.flat ? max_size : nblocks);
@@ -208,6 +208,8 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc,
       // Currently not storing neighbor gids
     }
 
+    // Include information about whether or not surface elements of this block
+    // are on physical boundaries of the domain
     for (int oxb = -1; oxb <= 1; ++oxb) {
       for (int oxa = -1; oxa <= 1; ++oxa) {
         if (pmb->IsPhysicalBoundary(inner_x1))
@@ -231,6 +233,12 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc,
       }
     }
 
+    // Include information about this blocks possible ownership of shared elements
+    for (int idx = 0; idx < bp_idxer.size(); ++idx) {
+      const auto [ok, oj, oi] = bp_idxer(idx);
+      pack.block_props_h_(blidx, bp_idxer.size() + idx) = pmb->ownership(oi, oj, ok);
+    }
+
     for (int i = 0; i < nvar; ++i) {
       pack.bounds_h_(0, blidx, i) = idx;
       for (const auto &[var_name, uid] : desc.var_groups[i]) {
diff --git a/src/solvers/mg_solver.hpp b/src/solvers/mg_solver.hpp
index 1de93a4b8bb40..d85504fa8c206 100644
--- a/src/solvers/mg_solver.hpp
+++ b/src/solvers/mg_solver.hpp
@@ -229,11 +229,6 @@ class MGSolver : public SolverBase, MGSolverCounter {
                     std::shared_ptr<MeshData<Real>> &md_xnew, double weight) {
     using namespace parthenon;
     const int ndim = md_rhs->GetMeshPointer()->ndim;
-    using TE = parthenon::TopologicalElement;
-    TE te = TE::CC;
-    IndexRange ib = md_rhs->GetBoundsI(IndexDomain::interior, te);
-    IndexRange jb = md_rhs->GetBoundsJ(IndexDomain::interior, te);
-    IndexRange kb = md_rhs->GetBoundsK(IndexDomain::interior, te);
 
     int nblocks = md_rhs->NumBlocks();
     std::vector<bool> include_block(nblocks, true);
@@ -254,18 +249,15 @@ class MGSolver : public SolverBase, MGSolverCounter {
     const int scratch_level = 0;
     parthenon::par_for_outer(
         DEFAULT_OUTER_LOOP_PATTERN, "Jacobi", DevExecSpace(), scratch_size, scratch_level,
-        0, pack_rhs.GetNBlocks() - 1, kb.s, kb.e,
-        KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) {
-          const int nvars = pack_rhs.GetUpperBound(b) - pack_rhs.GetLowerBound(b) + 1;
-          for (int c = 0; c < nvars; ++c) {
-            Real *Ax = &pack_Ax(b, te, c, k, jb.s, ib.s);
-            Real *diag = &pack_diag(b, te, c, k, jb.s, ib.s);
-            Real *prhs = &pack_rhs(b, te, c, k, jb.s, ib.s);
-            Real *xo = &pack_xold(b, te, c, k, jb.s, ib.s);
-            Real *xn = &pack_xnew(b, te, c, k, jb.s, ib.s);
-            // Use ptr arithmetic to get the number of points we need to go over
-            // (including ghost zones) to get from (k, jb.s, ib.s) to (k, jb.e, ib.e)
-            const int npoints = &pack_Ax(b, te, c, k, jb.e, ib.e) - Ax + 1;
+        0, pack_rhs.GetNBlocks() - 1,
+        KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b) {
+          LoopOverBlockVarsAndTEs(b, pack_rhs, [&](TopologicalElement te, int c) {
+            Real *Ax = pack_Ax(b, te, c).data();
+            Real *diag = pack_diag(b, te, c).data();
+            Real *prhs = pack_rhs(b, te, c).data();
+            Real *xo = pack_xold(b, te, c).data();
+            Real *xn = pack_xnew(b, te, c).data();
+            const int npoints = pack_Ax(b, te, c).size();
             parthenon::par_for_inner(
                 DEFAULT_INNER_LOOP_PATTERN, member, 0, npoints - 1, [&](const int idx) {
                   const Real off_diag = Ax[idx] - diag[idx] * xo[idx];
@@ -273,7 +265,7 @@ class MGSolver : public SolverBase, MGSolverCounter {
                   xn[idx] =
                       weight * robust::ratio(val, diag[idx]) + (1.0 - weight) * xo[idx];
                 });
-          }
+          });
         });
     return TaskStatus::complete;
   }
diff --git a/src/solvers/solver_utils.hpp b/src/solvers/solver_utils.hpp
index d24f8e3c5026d..2683a770384fb 100644
--- a/src/solvers/solver_utils.hpp
+++ b/src/solvers/solver_utils.hpp
@@ -14,23 +14,54 @@
 #define SOLVERS_SOLVER_UTILS_HPP_
 
 #include <algorithm>
+#include <cstdio>
 #include <limits>
 #include <memory>
+#include <set>
 #include <string>
 #include <utility>
 #include <vector>
 
 #include "kokkos_abstraction.hpp"
 
-#define PARTHENON_INTERNALSOLVERVARIABLE(base, varname)                                  \
-  struct varname : public parthenon::variable_names::base_t<false> {                     \
-    template <class... Ts>                                                               \
-    KOKKOS_INLINE_FUNCTION varname(Ts &&...args)                                         \
-        : parthenon::variable_names::base_t<false>(std::forward<Ts>(args)...) {}         \
-    static std::string name() { return base::name() + "." #varname; }                    \
+namespace parthenon {
+
+template <class pack_t, class func_t>
+KOKKOS_FORCEINLINE_FUNCTION void LoopOverBlockVarsAndTEs(const int b, pack_t &pack,
+                                                         func_t func) {
+  const int nvars = pack.GetUpperBound(b) - pack.GetLowerBound(b) + 1;
+  for (int c = 0; c < nvars; ++c) {
+    const auto tt = pack.GetTopologicalType(b, c);
+    const auto nel = GetNumberOfElements(tt);
+    for (int el = 0; el < nel; ++el) {
+      const auto te = GetTopologicalElementInDir(tt, el);
+      func(te, c);
+    }
   }
+}
 
-namespace parthenon {
+template <class TL_t>
+TaskStatus PrintFields(const std::shared_ptr<MeshData<Real>> &md_a, std::string label) {
+  using TE = parthenon::TopologicalElement;
+  IndexRange ib = md_a->GetBoundsI(IndexDomain::interior, TE::NN);
+  IndexRange jb = md_a->GetBoundsJ(IndexDomain::interior, TE::NN);
+  IndexRange kb = md_a->GetBoundsK(IndexDomain::interior, TE::NN);
+
+  printf("%s\n", label.c_str());
+  static auto desc = parthenon::MakePackDescriptorFromTypeList<TL_t>(md_a.get());
+  auto pack_a = desc.GetPack(md_a.get());
+  parthenon::par_for(
+      DEFAULT_LOOP_PATTERN, "print", parthenon::DevExecSpace(), 0,
+      pack_a.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
+        printf("[%i](%i, %i, %i) ", b, k, j, i);
+        LoopOverBlockVarsAndTEs(b, pack_a, [&](TopologicalElement te, int c) {
+          printf("[%i, %i] = %e ", c, static_cast<int>(te), pack_a(b, te, c, k, j, i));
+        });
+        printf("\n");
+      });
+  return TaskStatus::complete;
+}
 
 namespace solvers {
 
@@ -184,31 +215,22 @@ TaskStatus CopyData(const std::shared_ptr<MeshData<Real>> &md) {
 template <class TL, bool only_fine_on_composite = true>
 TaskStatus CopyData(const std::shared_ptr<MeshData<Real>> &md_in,
                     const std::shared_ptr<MeshData<Real>> &md_out) {
-  using TE = parthenon::TopologicalElement;
-  TE te = TE::CC;
-  IndexRange ib = md_in->GetBoundsI(IndexDomain::entire, te);
-  IndexRange jb = md_in->GetBoundsJ(IndexDomain::entire, te);
-  IndexRange kb = md_in->GetBoundsK(IndexDomain::entire, te);
-
   static auto desc = parthenon::MakePackDescriptorFromTypeList<TL>(md_in.get());
   auto pack_in = desc.GetPack(md_in.get(), only_fine_on_composite);
   auto pack_out = desc.GetPack(md_out.get(), only_fine_on_composite);
   const int scratch_size = 0;
   const int scratch_level = 0;
-  // Warning: This inner loop strategy only works because we are using IndexDomain::entire
-  const int npoints_inner = (kb.e - kb.s + 1) * (jb.e - jb.s + 1) * (ib.e - ib.s + 1);
   parthenon::par_for_outer(
       DEFAULT_OUTER_LOOP_PATTERN, "CopyData", DevExecSpace(), scratch_size, scratch_level,
       0, pack_in.GetNBlocks() - 1,
       KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b) {
-        const int nvars = pack_in.GetUpperBound(b) - pack_in.GetLowerBound(b) + 1;
-        for (int c = 0; c < nvars; ++c) {
-          Real *in = &pack_in(b, te, c, kb.s, jb.s, ib.s);
-          Real *out = &pack_out(b, te, c, kb.s, jb.s, ib.s);
-          parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, member, 0,
-                                   npoints_inner - 1,
+        LoopOverBlockVarsAndTEs(b, pack_in, [&](TopologicalElement te, int c) {
+          const int npts = pack_in(b, te, c).size();
+          Real const *const in = pack_in(b, te, c).data();
+          Real *out = pack_out(b, te, c).data();
+          parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, member, 0, npts - 1,
                                    [&](const int idx) { out[idx] = in[idx]; });
-        }
+        });
       });
   return TaskStatus::complete;
 }
@@ -216,8 +238,6 @@ TaskStatus CopyData(const std::shared_ptr<MeshData<Real>> &md_in,
 template <class TL, bool only_fine_on_composite = true>
 TaskStatus SetToZero(const std::shared_ptr<MeshData<Real>> &md) {
   int nblocks = md->NumBlocks();
-  using TE = parthenon::TopologicalElement;
-  TE te = TE::CC;
   static auto desc = [&] {
     if constexpr (isTypeList<TL>::value) {
       return parthenon::MakePackDescriptorFromTypeList<TL>(md.get());
@@ -233,17 +253,12 @@ TaskStatus SetToZero(const std::shared_ptr<MeshData<Real>> &md) {
       DEFAULT_OUTER_LOOP_PATTERN, "SetFieldsToZero", DevExecSpace(),
       scratch_size_in_bytes, scratch_level, 0, pack.GetNBlocks() - 1,
       KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b) {
-        auto cb = GetIndexShape(pack(b, te, 0), ng);
-        const auto &coords = pack.GetCoordinates(b);
-        IndexRange ib = cb.GetBoundsI(IndexDomain::interior, te);
-        IndexRange jb = cb.GetBoundsJ(IndexDomain::interior, te);
-        IndexRange kb = cb.GetBoundsK(IndexDomain::interior, te);
-        const int nvars = pack.GetUpperBound(b) - pack.GetLowerBound(b) + 1;
-        for (int c = 0; c < nvars; ++c) {
-          parthenon::par_for_inner(
-              DEFAULT_INNER_LOOP_PATTERN, member, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
-              [&](int k, int j, int i) { pack(b, te, c, k, j, i) = 0.0; });
-        }
+        LoopOverBlockVarsAndTEs(b, pack, [&](TopologicalElement te, int c) {
+          const int npts = pack(b, te, c).size();
+          Real *out = pack(b, te, c).data();
+          parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, member, 0, npts - 1,
+                                   [&](const int idx) { out[idx] = 0.0; });
+        });
       });
   return TaskStatus::complete;
 }
@@ -307,12 +322,6 @@ TaskStatus AddFieldsAndStoreInteriorSelect(const std::shared_ptr<MeshData<Real>>
                                            const std::shared_ptr<MeshData<Real>> &md_out,
                                            Real wa = 1.0, Real wb = 1.0,
                                            bool only_interior_blocks = false) {
-  using TE = parthenon::TopologicalElement;
-  TE te = TE::CC;
-  IndexRange ib = md_a->GetBoundsI(IndexDomain::entire, te);
-  IndexRange jb = md_a->GetBoundsJ(IndexDomain::entire, te);
-  IndexRange kb = md_a->GetBoundsK(IndexDomain::entire, te);
-
   int nblocks = md_a->NumBlocks();
   std::vector<bool> include_block(nblocks, true);
   if (only_interior_blocks) {
@@ -327,21 +336,19 @@ TaskStatus AddFieldsAndStoreInteriorSelect(const std::shared_ptr<MeshData<Real>>
   auto pack_out = desc.GetPack(md_out.get(), include_block, only_fine_on_composite);
   const int scratch_size = 0;
   const int scratch_level = 0;
-  // Warning: This inner loop strategy only works because we are using IndexDomain::entire
-  const int npoints_inner = (kb.e - kb.s + 1) * (jb.e - jb.s + 1) * (ib.e - ib.s + 1);
   parthenon::par_for_outer(
       DEFAULT_OUTER_LOOP_PATTERN, "AddFieldsAndStore", DevExecSpace(), scratch_size,
       scratch_level, 0, pack_a.GetNBlocks() - 1,
       KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b) {
-        const int nvars = pack_a.GetUpperBound(b) - pack_a.GetLowerBound(b) + 1;
-        for (int c = 0; c < nvars; ++c) {
-          Real *avar = &pack_a(b, te, c, kb.s, jb.s, ib.s);
-          Real *bvar = &pack_b(b, te, c, kb.s, jb.s, ib.s);
-          Real *out = &pack_out(b, te, c, kb.s, jb.s, ib.s);
+        LoopOverBlockVarsAndTEs(b, pack_a, [&](TopologicalElement te, int c) {
+          const int npts = pack_a(b, te, c).size();
+          Real const *const avar = pack_a(b, te, c).data();
+          Real const *const bvar = pack_b(b, te, c).data();
+          Real *out = pack_out(b, te, c).data();
           parthenon::par_for_inner(
-              DEFAULT_INNER_LOOP_PATTERN, member, 0, npoints_inner - 1,
+              DEFAULT_INNER_LOOP_PATTERN, member, 0, npts - 1,
               [&](const int idx) { out[idx] = wa * avar[idx] + wb * bvar[idx]; });
-        }
+        });
       });
   return TaskStatus::complete;
 }
@@ -504,10 +511,13 @@ TaskStatus DotProductLocal(const std::shared_ptr<MeshData<Real>> &md_a,
                            const std::shared_ptr<MeshData<Real>> &md_b,
                            AllReduce<Real> *adotb) {
   using TE = parthenon::TopologicalElement;
-  TE te = TE::CC;
-  IndexRange ib = md_a->GetBoundsI(IndexDomain::interior, te);
-  IndexRange jb = md_a->GetBoundsJ(IndexDomain::interior, te);
-  IndexRange kb = md_a->GetBoundsK(IndexDomain::interior, te);
+  // We iterate over the nodal index range since this encompasses all possible
+  // active cells for every topological type, then mask out elements that aren't
+  // owned/required by a given block
+  IndexRange ib = md_a->GetBoundsI(IndexDomain::interior, TE::NN);
+  IndexRange jb = md_a->GetBoundsJ(IndexDomain::interior, TE::NN);
+  IndexRange kb = md_a->GetBoundsK(IndexDomain::interior, TE::NN);
+  const int ndim = md_a->GetMeshPointer()->ndim;
 
   static auto desc = parthenon::MakePackDescriptorFromTypeList<TL>(md_a.get());
   auto pack_a = desc.GetPack(md_a.get());
@@ -517,12 +527,21 @@ TaskStatus DotProductLocal(const std::shared_ptr<MeshData<Real>> &md_a,
       parthenon::loop_pattern_mdrange_tag, "DotProduct", DevExecSpace(), 0,
       pack_a.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
       KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lsum) {
-        const int nvars = pack_a.GetUpperBound(b) - pack_a.GetLowerBound(b) + 1;
         // TODO(LFR): If this becomes a bottleneck, exploit hierarchical parallelism and
         //            pull the loop over vars outside of the innermost loop to promote
         //            vectorization.
-        for (int c = 0; c < nvars; ++c)
-          lsum += pack_a(b, te, c, k, j, i) * pack_b(b, te, c, k, j, i);
+        LoopOverBlockVarsAndTEs(b, pack_a, [&](TopologicalElement te, int c) {
+          const int oi = TopologicalOffsetI(te) * ((ib.e == i) - (ib.s == i));
+          const int oj = TopologicalOffsetJ(te) * ((jb.e == j) - (jb.s == j));
+          const int ok = TopologicalOffsetK(te) * ((kb.e == k) - (kb.s == k));
+          // Mask cell centered directions if loop bounds take them into ghosts
+          const int maski = (TopologicalOffsetI(te) || ndim < 1) ? 1 : (i < ib.e);
+          const int maskj = (TopologicalOffsetJ(te) || ndim < 2) ? 1 : (j < jb.e);
+          const int maskk = (TopologicalOffsetK(te) || ndim < 3) ? 1 : (k < kb.e);
+          if (pack_a.IsOwned(b, ok, oj, oi) *
+              (!pack_a.IsPhysicalBoundary(b, ok, oj, oi)) * maski * maskj * maskk)
+            lsum += pack_a(b, te, c, k, j, i) * pack_b(b, te, c, k, j, i);
+        });
       },
       Kokkos::Sum<Real>(gsum));
   adotb->val += gsum;
@@ -551,6 +570,50 @@ TaskID DotProduct(TaskID dependency_in, TaskList &tl, AllReduce<Real> *adotb,
   return finish_global_adotb;
 }
 
+template <class TL>
+TaskStatus ConstantBC(std::shared_ptr<MeshData<Real>> &md, bool coarse, Real val) {
+  using TE = TopologicalElement;
+  const int ndim = md->GetMeshPointer()->ndim;
+
+  std::set<PDOpt> opts{};
+  if (coarse) opts.emplace(PDOpt::Coarse);
+  auto desc = parthenon::MakePackDescriptorFromTypeList<TL>(
+      md.get(), std::vector<MetadataFlag>{}, opts);
+  auto pack = desc.GetPack(md.get(), GetBlockSelector::OnPhysicalBoundary());
+  const auto cellbounds = md->GetCellBounds(coarse ? CellLevel::coarse : CellLevel::same);
+
+  const int scratch_size = 0;
+  const int scratch_level = 0;
+  Indexer3D offset_idxer({-ndim > 2, ndim > 2}, {-ndim > 1, ndim > 1},
+                         {-ndim > 0, ndim > 0});
+  parthenon::par_for_outer(
+      DEFAULT_OUTER_LOOP_PATTERN, "DoBCs", DevExecSpace(), scratch_size, scratch_level, 0,
+      pack.GetNBlocks() - 1, 0, offset_idxer.size() - 1,
+      KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int oidx) {
+        IndexDomain idomains[3]{IndexDomain::inner, IndexDomain::unshared_interior,
+                                IndexDomain::outer};
+        const auto offset_tup = offset_idxer(oidx);
+        const int ok = std::get<0>(offset_tup);
+        const int oj = std::get<1>(offset_tup);
+        const int oi = std::get<2>(offset_tup);
+        if (pack.IsPhysicalBoundary(b, ok, oj, oi)) {
+          LoopOverBlockVarsAndTEs(b, pack, [&](TopologicalElement te, int c) {
+            const auto ib = cellbounds.GetBoundsI(idomains[oi + 1], te);
+            const auto jb = cellbounds.GetBoundsJ(idomains[oj + 1], te);
+            const auto kb = cellbounds.GetBoundsK(idomains[ok + 1], te);
+            Indexer3D idxer({kb.s, kb.e}, {jb.s, jb.e}, {ib.s, ib.e});
+            parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, member, 0,
+                                     idxer.size() - 1, [&](int idx) {
+                                       const auto [k, j, i] = idxer(idx);
+                                       pack(b, te, c, k, j, i) = val;
+                                     });
+          });
+        }
+      });
+
+  return TaskStatus::complete;
+}
+
 } // namespace utils
 
 } // namespace solvers