diff --git a/example/diffusion/diffusion_driver.cpp b/example/diffusion/diffusion_driver.cpp
index 0ea3c389ddd65..e9b29a632a154 100644
--- a/example/diffusion/diffusion_driver.cpp
+++ b/example/diffusion/diffusion_driver.cpp
@@ -58,6 +58,9 @@ TaskCollection DiffusionDriver::MakeTaskCollection() {
   auto pkg = pmesh->packages.Get("diffusion_package");
   auto psolver =
       pkg->Param<std::shared_ptr<parthenon::solvers::SolverBase>>("solver_pointer");
+  const auto alpha = pkg->Param<Real>("diagonal_alpha");
+  auto peqs = pkg->Param<std::shared_ptr<diffusion_package::DiffusionEquation<u, D>>>(
+      "diffusion_equation");
 
   auto partitions = pmesh->GetDefaultBlockPartitions();
   const int num_partitions = partitions.size();
@@ -69,21 +72,33 @@ TaskCollection DiffusionDriver::MakeTaskCollection() {
     // SetDiffusionCoefficient
     auto set_d = tl.AddTask(none, TF(SetDiffusionCoefficient), md, tm.dt);
 
-    auto &md_u = pmesh->mesh_data.Add("u", md, {u::name()});
-    auto &md_rhs = pmesh->mesh_data.Add("rhs", md, {u::name()});
+    auto &md_deltau = psolver->AddSolutionMeshData(pmesh, md, /*shallow=*/false);
+    auto &md_rhs = psolver->AddRHSMeshData(pmesh, md);
 
-    // SetRHS
-    auto set_rhs = tl.AddTask(set_d, (SetRHS), md, md_rhs);
+    // Set the rhs
+    // We are solving for Δu using  (alpha - dt ∇ D ∇) Δu = (dt ∇ D ∇) u_old. The
+    // diffusion_equation class defines the operator A = alpha - dt ∇ D ∇, so that
+    // we have A Δu = rhs with rhs = (alpha - A) u_old.
+    auto comm =
+        AddBoundaryExchangeTasks<BoundaryType::any>(none, tl, md, pmesh->multilevel);
+    auto Au = peqs->Ax(tl, comm | set_d, md, md, md_rhs);
+    auto set_rhs =
+        tl.AddTask(Au, solvers::utils::AddFieldsAndStore<parthenon::TypeList<u>>, md,
+                   md_rhs, md_rhs, alpha, -1.0);
 
     // Set initial solution guess to zero
-    auto zero_u = tl.AddTask(set_rhs, TF(solvers::utils::SetToZero<u>), md_u);
+    auto zero_u = tl.AddTask(set_rhs, TF(solvers::utils::SetToZero<u>), md_deltau);
+    psolver->initial_guess_is_zero = true;
     auto setup = psolver->AddSetupTasks(tl, zero_u, i, pmesh);
     auto solve = psolver->AddTasks(tl, setup, i, pmesh);
-    auto copy_back =
-        tl.AddTask(solve, TF(solvers::utils::CopyData<parthenon::TypeList<u>>), md_u, md);
 
-    auto new_dt = tl.AddTask(
-        copy_back, parthenon::Update::EstimateTimestep<MeshData<Real>>, md.get());
+    // Update to u = u_0 + Δu
+    auto update_u =
+        tl.AddTask(solve, solvers::utils::AddFieldsAndStore<parthenon::TypeList<u>>, md,
+                   md_deltau, md, 1.0, 1.0);
+
+    // Update the timestep
+    tl.AddTask(update_u, parthenon::Update::EstimateTimestep<MeshData<Real>>, md.get());
   }
   return tc;
 }
diff --git a/example/diffusion/diffusion_driver.hpp b/example/diffusion/diffusion_driver.hpp
index b23fddf73d36d..f0af76a530d80 100644
--- a/example/diffusion/diffusion_driver.hpp
+++ b/example/diffusion/diffusion_driver.hpp
@@ -20,6 +20,7 @@
 #include <kokkos_abstraction.hpp>
 #include <parthenon/driver.hpp>
 #include <parthenon/package.hpp>
+#include <solvers/solver_base.hpp>
 
 namespace diffusion_example {
 using namespace parthenon::driver::prelude;
@@ -33,9 +34,33 @@ class DiffusionDriver : public EvolutionDriver {
   }
   // This next function essentially defines the driver.
   TaskCollection MakeTaskCollection();
-  TaskListStatus Step();
+  TaskListStatus Step() override;
 
   // DriverStatus Execute() override;
+  void OutputDownstreamCycleDiagnostics() override {
+    auto pkg = pmesh->packages.Get("diffusion_package");
+    auto solver_type = pkg->Param<std::string>("solver");
+    auto psolver =
+        pkg->Param<std::shared_ptr<parthenon::solvers::SolverBase>>("solver_pointer");
+    int v_cycles = psolver->GetFinalIterations();
+    if (solver_type == "BiCGSTAB") v_cycles *= 2;
+    std::cout << " v-cycles=" << v_cycles;
+  }
+  
+  void PostExecute(DriverStatus status) override {
+    EvolutionDriver::PostExecute(status);
+    if (parthenon::Globals::my_rank == 0) {
+      auto pkg = pmesh->packages.Get("diffusion_package");
+      if (pkg->Param<bool>("report_timings")) {
+        printf("\nTiming data\n-----------\n");
+        auto psolver =
+            pkg->Param<std::shared_ptr<parthenon::solvers::SolverBase>>("solver_pointer");
+        std::cout << "Solver breakdown: \n" << psolver->solver_timings;
+        psolver->solver_timings.clear();
+      }
+      
+    }
+  }
 
  private:
   LowStorageIntegrator integrator;
diff --git a/example/diffusion/diffusion_equation.hpp b/example/diffusion/diffusion_equation.hpp
index ae7dc1c01575e..418d4ec6838db 100644
--- a/example/diffusion/diffusion_equation.hpp
+++ b/example/diffusion/diffusion_equation.hpp
@@ -23,6 +23,7 @@
 #include <parthenon/package.hpp>
 
 #include "diffusion_package.hpp"
+#include "raw_memory_indexer.hpp"
 
 namespace diffusion_package {
 // Calculate A in_t = out_t (in the region covered by md) for a given set of fluxes
@@ -49,29 +50,45 @@ FluxMultiplyMatrix(std::shared_ptr<parthenon::MeshData<Real>> &md,
   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 idxer =
+      RawMemoryIndexer::IJ(IndexDomain::interior, 0, md.get(), TE::CC, TE::CC);
+  const int max_inner_raw = idxer.GetMaxNinnerRaw();
+  auto scratch_size_in_bytes =
+      0 * parthenon::ScratchPad1D<Real>::shmem_size(max_inner_raw);
+  const std::size_t scratch_level = 1;
+  const int ioff = (ndim > 0);
+  const int joff = (ndim > 1);
+  const int koff = (ndim > 2);
+  parthenon::par_for_outer(
+      DEFAULT_OUTER_LOOP_PATTERN, "FluxMultiplyMatrix", DevExecSpace(),
+      scratch_size_in_bytes, scratch_level, 0, pack.GetNBlocks() - 1, 0,
+      idxer.GetNouter() - 1,
+      KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int idx_out) {
+        const auto [ks, js, is] = idxer.GetStartIndices(idx_out);
+        const auto ninner = idxer.GetNinnerRaw(idx_out);
+        const auto raw_idx_start = idxer.GetStartingRawFlatIdx(idx_out);
+
         const auto &coords = pack.GetCoordinates(b);
-        Real dx1 = coords.template Dxc<X1DIR>(k, j, i);
-        pack_out(b, var_t(), k, j, i) = -alpha * pack(b, var_t(), k, j, i);
-        pack_out(b, 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, 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, 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;
-        }
+        const Real invdx1 = ioff / coords.template Dxc<X1DIR>(ks, js, is);
+        const Real invdx2 = joff / coords.template Dxc<X2DIR>(ks, js, is);
+        const Real invdx3 = koff / coords.template Dxc<X3DIR>(ks, js, is);
+
+        const Real *const flx1_up = &pack.flux(b, X1DIR, var_t(), ks, js, is + 1);
+        const Real *const flx1_lo = &pack.flux(b, X1DIR, var_t(), ks, js, is);
+        const Real *const flx2_up = &pack.flux(b, X2DIR, var_t(), ks, js + joff, is);
+        const Real *const flx2_lo = &pack.flux(b, X2DIR, var_t(), ks, js, is);
+        const Real *const flx3_up = &pack.flux(b, X3DIR, var_t(), ks + koff, js, is);
+        const Real *const flx3_lo = &pack.flux(b, X3DIR, var_t(), ks, js, is);
+        const Real *const in = &pack(b, var_t(), ks, js, is);
+        Real *out = &pack_out(b, var_t(), ks, js, is);
+        parthenon::par_for_inner(
+            DEFAULT_INNER_LOOP_PATTERN, member, 0, ninner - 1, [&](const int idx) {
+              const Real dfx = (flx1_up[idx] - flx1_lo[idx]) * invdx1;
+              const Real dfy = (flx2_up[idx] - flx2_lo[idx]) * invdx2;
+              const Real dfz = (flx3_up[idx] - flx3_lo[idx]) * invdx3;
+              out[idx] = alpha * in[idx] + dfx + dfy + dfz;
+            });
       });
   return TaskStatus::complete;
 }
@@ -151,19 +168,19 @@ class DiffusionEquation {
           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) -
+          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) +
+            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) +
+            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);
           }
@@ -183,48 +200,48 @@ class DiffusionEquation {
 
     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 pack = desc.GetPack(md.get());
     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, var_t(), k, j, i - 1) - pack(b, 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, var_t(), k, j, i) - pack(b, var_t(), k, j, i + 1));
+    auto pack_mat = desc_mat.GetPack(md_mat.get());
 
-          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, var_t(), k, j - 1, i) - pack(b, 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, var_t(), k, j, i) - pack(b, var_t(), k, j + 1, i)) / dx2;
-          }
+    for (int dim = 0; dim < ndim; ++dim) {
+      const auto te = dim == 0 ? TE::F1 : (dim == 1 ? TE::F2 : TE::F3);
+      const int ioff = dim == 0;
+      const int joff = dim == 1;
+      const int koff = dim == 2;
+      const int dir = dim + 1;
+      // The diffusion coefficient arrays are cell mem aligned
+      const auto idxer =
+          RawMemoryIndexer::IJ(IndexDomain::interior, 0, md.get(), te, TE::CC);
+      const std::size_t scratch_level = 1;
+      auto scratch_size_in_bytes =
+          0 * parthenon::ScratchPad1D<Real>::shmem_size(idxer.GetMaxNinnerRaw());
+      parthenon::par_for_outer(
+          DEFAULT_OUTER_LOOP_PATTERN, "CalculateFluxes", DevExecSpace(),
+          scratch_size_in_bytes, scratch_level, 0, pack.GetNBlocks() - 1, 0,
+          idxer.GetNouter() - 1,
+          KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int idx_out) {
+            const auto [ks, js, is] = idxer.GetStartIndices(idx_out);
+            const auto ninner = idxer.GetNinnerRaw(idx_out);
+            const auto raw_idx_start = idxer.GetStartingRawFlatIdx(idx_out);
+
+            const auto &coords = pack.GetCoordinates(b);
+            const Real inv_dx = ioff / coords.template Dxc<X1DIR>(ks, js, is) +
+                                joff / coords.template Dxc<X2DIR>(ks, js, is) +
+                                koff / coords.template Dxc<X3DIR>(ks, js, is);
+
+            Real *flx = &(pack.flux(b, dir, var_t(), ks, js, is));
+            const Real *const D = &pack_mat(b, te, D_t(), ks, js, is);
+            const Real *const vup = &pack(b, var_t(), ks, js, is);
+            const Real *const vlo = &pack(b, var_t(), ks - koff, js - joff, is - ioff);
+            parthenon::par_for_inner(
+                DEFAULT_INNER_LOOP_PATTERN, member, 0, ninner - 1, [&](const int idx) {
+                  flx[idx] = -D[idx] * (vup[idx] - vlo[idx]) * inv_dx;
+                });
+          });
+    }
 
-          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, var_t(), k - 1, j, i) - pack(b, 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, var_t(), k, j, i) - pack(b, var_t(), k + 1, j, i)) / dx3;
-          }
-        });
     return TaskStatus::complete;
   }
 
@@ -305,6 +322,60 @@ class DiffusionEquation {
         });
     return TaskStatus::complete;
   }
+
+  static parthenon::TaskStatus SetBoundary(std::shared_ptr<parthenon::MeshData<Real>> &md,
+                                           bool coarse) {
+    using namespace parthenon;
+    const int ndim = md->GetMeshPointer()->ndim;
+
+    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());
+    if (pack.GetNBlocks() > 0) {
+      CellLevel cl = coarse ? CellLevel::coarse : CellLevel::same;
+      IndexRange ib = md->GetBoundsI(cl, IndexDomain::interior);
+      IndexRange jb = md->GetBoundsJ(cl, IndexDomain::interior);
+      IndexRange kb = md->GetBoundsK(cl, IndexDomain::interior);
+
+      const int scratch_size_in_bytes = 0;
+      const std::size_t scratch_level = 1;
+      parthenon::par_for_outer(
+          DEFAULT_OUTER_LOOP_PATTERN, "SetBoundaries", DevExecSpace(),
+          scratch_size_in_bytes, scratch_level, 0, pack.GetNBlocks() - 1, -(ndim > 2),
+          (ndim > 2), -(ndim > 1), (ndim > 1), -1, 1,
+          KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, int ok, int oj,
+                        int oi) {
+            const int tot_offset = std::abs(ok) + std::abs(oj) + std::abs(oi);
+
+            auto get_lower = [](int offset, auto bound) {
+              if (offset != 0) {
+                return offset > 0 ? bound.e : bound.s;
+              } else {
+                return bound.s;
+              }
+            };
+            auto get_upper = [](int offset, auto bound) {
+              if (offset != 0) {
+                return offset > 0 ? bound.e : bound.s;
+              } else {
+                return bound.e;
+              }
+            };
+
+            if (tot_offset == 1 && pack.IsPhysicalBoundary(b, ok, oj, oi)) {
+              parthenon::par_for_inner(
+                  DEFAULT_INNER_LOOP_PATTERN, member, get_lower(ok, kb),
+                  get_upper(ok, kb), get_lower(oj, jb), get_upper(oj, jb),
+                  get_lower(oi, ib), get_upper(oi, ib),
+                  [&](const int k, const int j, const int i) {
+                    pack(b, var_t(), k + ok, j + oj, i + oi) = -pack(b, var_t(), k, j, i);
+                  });
+            }
+          });
+    }
+    return TaskStatus::complete;
+  }
 };
 
 } // namespace diffusion_package
diff --git a/example/diffusion/diffusion_package.cpp b/example/diffusion/diffusion_package.cpp
index 60b726d9809b9..5eeebbd457b38 100644
--- a/example/diffusion/diffusion_package.cpp
+++ b/example/diffusion/diffusion_package.cpp
@@ -76,25 +76,35 @@ struct any_diffusion : public parthenon::variable_names::base_t<true> {
 };
 
 template <CoordinateDirection DIR, BCSide SIDE>
-auto GetBC() {
-  return [](std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) -> void {
+auto GetBC(Real val = 0.0) {
+  return [val](std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse) -> void {
     using namespace parthenon;
     using namespace parthenon::BoundaryFunction;
-    GenericBC<DIR, SIDE, BCType::FixedFace, any_diffusion>(rc, coarse, 0.0);
+    GenericBC<DIR, SIDE, BCType::FixedFace, any_diffusion>(rc, coarse, val);
   };
 }
 
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto pkg = std::make_shared<StateDescriptor>("diffusion_package");
 
+  auto u_bounds =
+      pin->GetOrAddVector<Real>("diffusion", "boundary_u", {0.0}, "Boundary us.");
+  if (u_bounds.size() == 1) u_bounds = std::vector<Real>(6, u_bounds[0]);
+
   // Set boundary conditions for Diffusion 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>());
+  pkg->UserBoundaryFunctions[BF::inner_x1].push_back(
+      GetBC<X1DIR, BCSide::Inner>(u_bounds[0]));
+  pkg->UserBoundaryFunctions[BF::inner_x2].push_back(
+      GetBC<X2DIR, BCSide::Inner>(u_bounds[2]));
+  pkg->UserBoundaryFunctions[BF::inner_x3].push_back(
+      GetBC<X3DIR, BCSide::Inner>(u_bounds[4]));
+  pkg->UserBoundaryFunctions[BF::outer_x1].push_back(
+      GetBC<X1DIR, BCSide::Outer>(u_bounds[1]));
+  pkg->UserBoundaryFunctions[BF::outer_x2].push_back(
+      GetBC<X2DIR, BCSide::Outer>(u_bounds[3]));
+  pkg->UserBoundaryFunctions[BF::outer_x3].push_back(
+      GetBC<X3DIR, BCSide::Outer>(u_bounds[5]));
 
   // probably should stay 1.0
   Real diagonal_alpha = pin->GetOrAddReal("diffusion", "diagonal_alpha", 1.0);
@@ -102,8 +112,6 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
   Real cfl = pin->GetOrAddReal("diffusion", "cfl", 1.0);
   pkg->AddParam<>("cfl", cfl);
-  pkg->AddParam<>("dt", 1.0,
-                  parthenon::Params::Mutability::Mutable); // hold timestep for controls
 
   Real t0 = pin->GetOrAddReal("diffusion", "t0", 0.001);
   pkg->AddParam<>("t0", t0);
@@ -117,13 +125,29 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   std::string prolong =
       pin->GetOrAddString("diffusion", "boundary_prolongation", "Linear");
 
+  pkg->AddParam<>("diffusion_coefficient", DiffusionCoefficient(pin));
+
   using PoissEq = diffusion_package::DiffusionEquation<u, D>;
-  PoissEq eq(pin, "diffusion");
-  pkg->AddParam<>("diffusion_equation", eq, parthenon::Params::Mutability::Mutable);
+  pkg->AddParam<>("diffusion_equation", std::make_shared<PoissEq>(pin, "diffusion"));
+  
+  bool report_timings =
+      pin->GetOrAddBoolean("diffusion", "report_timings", false,
+                           "Report different timings of the diffusion solver "
+                           "at the end of the calculation.");
+  pkg->AddParam<>("report_timings", report_timings);
+  if (parthenon::Globals::my_rank == 0 && report_timings) {
+    parthenon::Task::enable_timing = true;
+    parthenon::Task::enable_timing_chunks = false;
+  }
 
   std::shared_ptr<parthenon::solvers::SolverBase> psolver;
-  using prolongator_t = parthenon::solvers::ProlongationBlockInteriorDefault;
-  using preconditioner_t = parthenon::solvers::MGSolver<PoissEq, prolongator_t>;
+
+  using prolongator_t = parthenon::solvers::ProlongationBlockInteriorZeroDirichlet;
+  using restrictor_t = parthenon::solvers::RestrictionCombined;
+
+  using preconditioner_t =
+      parthenon::solvers::MGSolver<PoissEq, prolongator_t, restrictor_t>;
+
   if (solver == "MG") {
     psolver = std::make_shared<parthenon::solvers::MGSolver<PoissEq, prolongator_t>>(
         "base", "u", "rhs", pin, "diffusion/solver_params", PoissEq(pin, "diffusion"));
@@ -142,7 +166,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
   using namespace parthenon::refinement_ops;
   auto mD = Metadata({Metadata::Independent, Metadata::OneCopy, Metadata::Face,
-                      Metadata::GMGRestrict, Metadata::FillGhost});
+                      Metadata::GMGRestrict, Metadata::CellMemAligned});
   mD.RegisterRefinementOps<ProlongateSharedLinear, RestrictAverage>();
 
   // Holds the discretized version of D in \nabla \cdot D(\vec{x}) \nabla u = rhs. D = 1
@@ -176,7 +200,7 @@ Real EstimateTimestep(MeshData<Real> *md) {
   std::shared_ptr<StateDescriptor> pkg =
       md->GetMeshPointer()->packages.Get("diffusion_package");
   const auto &cfl = pkg->Param<Real>("cfl");
-  const auto &old_dt = pkg->Param<Real>("dt");
+  const auto profile_D = pkg->Param<DiffusionCoefficient>("diffusion_coefficient");
 
   auto desc = parthenon::MakePackDescriptor<diffusion_package::D>(md);
   auto pack = desc.GetPack(md);
@@ -187,8 +211,6 @@ Real EstimateTimestep(MeshData<Real> *md) {
 
   const int ndim = md->GetMeshPointer()->ndim;
 
-  constexpr static Real ONE_FOURTH = 0.25;
-
   Real min_dt;
   parthenon::par_reduce(
       parthenon::loop_pattern_mdrange_tag, PARTHENON_AUTO_LABEL, DevExecSpace(), 0,
@@ -199,29 +221,18 @@ Real EstimateTimestep(MeshData<Real> *md) {
         const Real dx2 = coords.Dxc<X2DIR>(k, j, i);
         const Real dx3 = coords.Dxc<X3DIR>(k, j, i);
 
-        // average face coefficients and divide by 2
-        const Real d1 = (pack(b, TE::F1, diffusion_package::D(), k, j, i + 1) +
-                         pack(b, TE::F1, diffusion_package::D(), k, j, i)) *
-                        ONE_FOURTH;
-        lmin_dt = std::min(lmin_dt, dx1 * dx1 / d1);
-        if (ndim > 1) {
-          const Real d2 = (pack(b, TE::F2, diffusion_package::D(), k, j + 1, i) +
-                           pack(b, TE::F2, diffusion_package::D(), k, j, i)) *
-                          ONE_FOURTH;
-          lmin_dt = std::min(lmin_dt, dx2 * dx2 / d2);
-        }
-        if (ndim > 2) {
-          const Real d3 = (pack(b, TE::F3, diffusion_package::D(), k + 1, j, i) +
-                           pack(b, TE::F3, diffusion_package::D(), k, j, i)) *
-                          ONE_FOURTH;
-          lmin_dt = std::min(lmin_dt, dx3 * dx3 / d3);
-        }
+        const Real x1 = coords.Xc<X1DIR>(k, j, i);
+        const Real x2 = coords.Xc<X2DIR>(k, j, i);
+        const Real x3 = coords.Xc<X3DIR>(k, j, i);
+
+        const Real D = profile_D(x1, x2, x3);
+        Real dtc = dx1 * dx1 / D;
+        if (ndim > 1) dtc = std::min(dtc, dx2 * dx2 / D);
+        if (ndim > 2) dtc = std::min(dtc, dx3 * dx3 / D);
+        lmin_dt = std::min(lmin_dt, dtc);
       },
       Kokkos::Min<Real>(min_dt));
-
-  const Real new_dt = cfl * min_dt * old_dt; // need to scale by dt as D := D * dt
-  pkg->UpdateParam<Real>("dt", new_dt);
-  return new_dt;
+  return cfl * min_dt;
 } // EstimateTimestep
 
 parthenon::TaskStatus SetRHS(std::shared_ptr<MeshData<Real>> md,
@@ -230,10 +241,7 @@ parthenon::TaskStatus SetRHS(std::shared_ptr<MeshData<Real>> md,
   const auto alpha = pkg->Param<Real>("diagonal_alpha");
   auto desc = parthenon::MakePackDescriptor<diffusion_package::u>(md.get());
   auto pack = desc.GetPack(md.get());
-
-  // holds rhs
-  auto desc_rhs = parthenon::MakePackDescriptor<diffusion_package::u>(md_rhs.get());
-  auto pack_rhs = desc_rhs.GetPack(md_rhs.get());
+  auto pack_rhs = desc.GetPack(md_rhs.get());
 
   IndexRange ib = md->GetBoundsI(IndexDomain::interior);
   IndexRange jb = md->GetBoundsJ(IndexDomain::interior);
@@ -243,7 +251,7 @@ parthenon::TaskStatus SetRHS(std::shared_ptr<MeshData<Real>> md,
       "SetRHS", 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_rhs(b, diffusion_package::u(), k, j, i) = // rhs
-            -alpha * pack(b, diffusion_package::u(), k, j, i);
+            alpha * pack(b, diffusion_package::u(), k, j, i);
       });
   return TaskStatus::complete;
 } // SetRHS
@@ -257,7 +265,7 @@ parthenon::TaskStatus SetDiffusionCoefficient(std::shared_ptr<MeshData<Real>> md
       parthenon::MakePackDescriptor<diffusion_package::u, diffusion_package::D>(md.get());
   auto pack = desc.GetPack(md.get());
 
-  const bool constant_coeff = pkg->Param<bool>("constant_coefficient");
+  const auto profile_D = pkg->Param<DiffusionCoefficient>("diffusion_coefficient");
 
   for (auto te : {TE::F1, TE::F2, TE::F3}) {
     IndexRange ib = md->GetBoundsI(IndexDomain::interior, te);
@@ -274,12 +282,15 @@ parthenon::TaskStatus SetDiffusionCoefficient(std::shared_ptr<MeshData<Real>> md
           const Real u = 0.5 * (pack(b, TE::CC, diffusion_package::u(), k - offset_x3,
                                      j - offset_x2, i - offset_x1) +
                                 pack(b, TE::CC, diffusion_package::u(), k, j, i));
-          if (constant_coeff) {
-            pack(b, te, diffusion_package::D(), k, j, i) = 1.0 * dt;
-          } else {
-            pack(b, te, diffusion_package::D(), k, j, i) =
-                10.0 * std::sqrt(std::fabs(u)) * dt;
-          }
+
+          const auto &coords = pack.GetCoordinates(b);
+          const Real x1 =
+              offset_x1 ? coords.Xc<1>(k, j, i) : coords.X<1, TE::F1>(k, j, i);
+          const Real x2 =
+              offset_x2 ? coords.Xc<2>(k, j, i) : coords.X<2, TE::F2>(k, j, i);
+          const Real x3 =
+              offset_x3 ? coords.Xc<3>(k, j, i) : coords.X<3, TE::F3>(k, j, i);
+          pack(b, te, diffusion_package::D(), k, j, i) = dt * profile_D(x1, x2, x3);
         });
   }
   return TaskStatus::complete;
diff --git a/example/diffusion/diffusion_package.hpp b/example/diffusion/diffusion_package.hpp
index c9e224a18d1f1..5a0a42e339ca8 100644
--- a/example/diffusion/diffusion_package.hpp
+++ b/example/diffusion/diffusion_package.hpp
@@ -35,6 +35,36 @@ using namespace parthenon::package::prelude;
 VARIABLE(diffusion, D);
 VARIABLE(diffusion, u);
 
+struct DiffusionCoefficient {
+  Real Dright{1.0};
+  Real Dleft{1.0};
+  Real amplitude{1.0};
+  Real wavelength{1.0};
+  DiffusionCoefficient(parthenon::ParameterInput *pin) {
+    const bool constant_coeff =
+        pin->GetOrAddBoolean("diffusion", "constant_coefficient", true);
+    Dleft = pin->GetOrAddReal("diffusion", "Dleft", 1.e3,
+                              "Value of diffusion coefficient to the left.");
+    Dright = pin->GetOrAddReal("diffusion", "Dright", 1.e8,
+                               "Value of diffusion coefficient to the right.");
+    amplitude = pin->GetOrAddReal("diffusion", "amplitude", 0.15,
+                                  "Spatial amplitude of boundary.");
+    wavelength = pin->GetOrAddReal("diffusion", "wavelength", 1.0,
+                                   "Spatial wavelength of boundary.");
+    if (constant_coeff) {
+      Dleft = 1.0;
+      Dright = 1.0;
+    }
+  }
+
+  KOKKOS_FORCEINLINE_FUNCTION
+  Real operator()(Real x, Real y, Real z) const {
+    const Real xcrit = amplitude * sin(2.0 * M_PI * y / wavelength);
+    if (x >= xcrit) return Dright;
+    return Dleft;
+  }
+};
+
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);
 TaskStatus SetRHS(std::shared_ptr<MeshData<Real>> md,
                   std::shared_ptr<MeshData<Real>> md_rhs);
diff --git a/example/diffusion/parthenon_app_inputs.cpp b/example/diffusion/parthenon_app_inputs.cpp
index 8664cad3321ce..5d6d431d533a2 100644
--- a/example/diffusion/parthenon_app_inputs.cpp
+++ b/example/diffusion/parthenon_app_inputs.cpp
@@ -34,6 +34,7 @@ void ProblemGenerator(Mesh *pm, ParameterInput *pin, MeshData<Real> *md) {
   auto pmb = md->GetBlockData(0)->GetBlockPointer();
   const int ndim = md->GetMeshPointer()->ndim;
 
+  Real scale = pin->GetOrAddReal("diffusion", "scale", 1.0);
   Real x0 = pin->GetOrAddReal("diffusion", "x0", 0.0);
   Real y0 = pin->GetOrAddReal("diffusion", "y0", 0.0);
   Real z0 = pin->GetOrAddReal("diffusion", "z0", 0.0);
@@ -42,8 +43,7 @@ void ProblemGenerator(Mesh *pm, ParameterInput *pin, MeshData<Real> *md) {
   const bool constant_coeff =
       pin->GetOrAddBoolean("diffusion", "constant_coefficient", true);
 
-  auto desc =
-      parthenon::MakePackDescriptor<diffusion_package::u, diffusion_package::D>(md);
+  auto desc = parthenon::MakePackDescriptor<diffusion_package::u>(md);
   auto pack = desc.GetPack(md);
 
   using TE = parthenon::TopologicalElement;
@@ -51,16 +51,17 @@ void ProblemGenerator(Mesh *pm, ParameterInput *pin, MeshData<Real> *md) {
   auto ib = cellbounds.GetBoundsI(IndexDomain::entire);
   auto jb = cellbounds.GetBoundsJ(IndexDomain::entire);
   auto kb = cellbounds.GetBoundsK(IndexDomain::entire);
+
   pmb->par_for(
       "Diffusion::ProblemGenerator", 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.Xc<1>(i);
         Real x2 = coords.Xc<2>(j);
-        Real x3 = coords.Xc<2>(k);
+        Real x3 = coords.Xc<3>(k);
         Real x1f = coords.X<1, TE::F1>(k, j, i);
         Real x2f = coords.X<2, TE::F2>(k, j, i);
-        Real x3f = coords.X<2, TE::F3>(k, j, i);
+        Real x3f = coords.X<3, TE::F3>(k, j, i);
         Real dx1 = coords.Dxc<1>(k, j, i);
         Real dx2 = coords.Dxc<2>(k, j, i);
         Real dx3 = coords.Dxc<3>(k, j, i);
@@ -75,17 +76,7 @@ void ProblemGenerator(Mesh *pm, ParameterInput *pin, MeshData<Real> *md) {
           return std::exp(exponent);
         };
         const Real val = profile(x1, x2, x3);
-        pack(b, diffusion_package::u(), k, j, i) = val;
-
-        if (constant_coeff) {
-          pack(b, TE::F1, diffusion_package::D(), k, j, i) = 1.0 * dt;
-          pack(b, TE::F2, diffusion_package::D(), k, j, i) = 1.0 * dt;
-          pack(b, TE::F3, diffusion_package::D(), k, j, i) = 1.0 * dt;
-        } else {
-          pack(b, TE::F1, diffusion_package::D(), k, j, i) = profile(x1f, x2, x3) * dt;
-          pack(b, TE::F2, diffusion_package::D(), k, j, i) = profile(x1, x2f, x3) * dt;
-          pack(b, TE::F3, diffusion_package::D(), k, j, i) = profile(x1, x2, x3f) * dt;
-        }
+        pack(b, diffusion_package::u(), k, j, i) = scale * val;
       });
 }
 
diff --git a/example/diffusion/parthinput.diffusion b/example/diffusion/parthinput.diffusion
index cd849285a783b..0898cf2ed6a98 100644
--- a/example/diffusion/parthinput.diffusion
+++ b/example/diffusion/parthinput.diffusion
@@ -21,66 +21,86 @@ problem_id = diffusion
 <parthenon/mesh>
 refinement = static
 multigrid = true 
+base_block_coarsenings = 3
 
 nx1 = 64
-x1min = -1.0
-x1max = 1.0
+x1min = -0.5
+x1max = 0.5
 ix1_bc = outflow 
 ox1_bc = outflow
 
 nx2 = 64
-x2min = -1.0
-x2max = 1.0
-ix2_bc = outflow
-ox2_bc = outflow
+x2min = -0.5
+x2max = 0.5
+ix2_bc = periodic
+ox2_bc = periodic
 
 nx3 = 1
-x3min = 0.0
-x3max = 1.0
+x3min = -0.5
+x3max = 0.5
 ix3_bc = periodic
 ox3_bc = periodic
 
 <parthenon/meshblock>
-nx1 = 32
-nx2 = 32
+nx1 = 64
+nx2 = 64
 nx3 = 1
 
 <parthenon/time>
-nlim = -1
+nlim = 50
 tlim = 0.02
 integrator = rk1
 #ncycle_out_mesh = -10000
 
 <parthenon/output0>
 file_type = hdf5
-dt = 0.0001
+dt = 0.1
 variables = diffusion.u
 ghost_zones = false
 
 <parthenon/static_refinement0> 
-x1min = -0.15
-x1max = +0.15
-x2min = -0.15
-x2max = +0.15
-level = 3
+x1min = -0.07
+x1max = +0.07
+x2min = -0.07
+x2max = +0.07
+x3min = -0.07
+x3max = +0.07
+level = 7
 
 <diffusion>
-#solver = BiCGSTAB # or MG
-solver = MG
+solver = BiCGSTAB
 flux_correct = true
-diagonal_alpha = 1.0
-cfl = 10.0
-set_flux_boundary = true
+cfl = 1e12
+set_flux_boundary = false
+constant_coefficient = false
+boundary_prolongation = Constant
+
+# Set up the operator
+diagonal_alpha = 1e6
+wavelength = 1.0
+amplitude = 0.15
+Dleft = 1.e3
+Dright = 1.e8
+
+# Set up boundary conditions
+boundary_u = 1.e-5, 1.0, 1.e-5, 1.e-5, 1.e-5, 1.e-5
 
+# Set up constant initial conditions
+scale = 1.e-5
 x0 = 0.0
 y0 = 0.0
 z0 = 0.0
-t0 = 0.001
+t0 = 1e200
 
 <diffusion/solver_params>
 precondition = true
-max_iterations = 20
-residual_tolerance = 1.e-9
-print_per_step = true
-smoother = SRJ2
+max_iterations = 200
+max_coarsenings = 10000000
+absolute_residual_tolerance = 1.e-16
+relative_residual_tolerance = 1.e-8
+print_per_step = false
+presmoother = SRJ1
+postsmoother = SRJ3
 do_FAS = true
+block_interior_prolongation = Kwak
+volume_weight = true
diff --git a/example/diffusion/raw_memory_indexer.hpp b/example/diffusion/raw_memory_indexer.hpp
new file mode 100644
index 0000000000000..99672377bd7db
--- /dev/null
+++ b/example/diffusion/raw_memory_indexer.hpp
@@ -0,0 +1,105 @@
+//========================================================================================
+// (C) (or copyright) 2026. 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_DIFFUSION_RAW_MEMORY_INDEXER_HPP_
+#define EXAMPLE_DIFFUSION_RAW_MEMORY_INDEXER_HPP_
+
+#include <parthenon/package.hpp>
+
+using namespace parthenon::package::prelude;
+
+struct RawMemoryIndexer {
+  using ID = parthenon::IndexDomain;
+  using TE = parthenon::TopologicalElement;
+  RawMemoryIndexer(int inner_length, ID domain, int halo, parthenon::MeshData<Real> *md,
+                   TE domain_te, TE memory_te = TE::CC)
+      : inner_length(inner_length), idxer_entire(md->GetBoundsK(ID::entire, memory_te),
+                                                 md->GetBoundsJ(ID::entire, memory_te),
+                                                 md->GetBoundsI(ID::entire, memory_te)) {
+    auto pmesh = md->GetMeshPointer();
+    const int ndim = pmesh->ndim;
+    ib = md->GetBoundsI(domain, domain_te);
+    jb = md->GetBoundsJ(domain, domain_te);
+    kb = md->GetBoundsK(domain, domain_te);
+    if (ndim > 0) {
+      ib.s -= halo;
+      ib.e += halo;
+    }
+    if (ndim > 1) {
+      jb.s -= halo;
+      jb.e += halo;
+    }
+    if (ndim > 2) {
+      kb.s -= halo;
+      kb.e += halo;
+    }
+    idxer = parthenon::Indexer3D({kb.s, kb.e}, {jb.s, jb.e}, {ib.s, ib.e});
+    PARTHENON_REQUIRE(memory_te == TE::CC || memory_te == TE::NN,
+                      "Only two kinds of memory layouts for topological elements.");
+  }
+
+  static RawMemoryIndexer IJ(ID domain, parthenon::MeshData<Real> *md, TE domain_te,
+                             TE memory_te = TE::CC) {
+    return RawMemoryIndexer::IJ(domain, 0, md, domain_te, memory_te);
+  }
+
+  static RawMemoryIndexer IJ(ID domain, int halo, parthenon::MeshData<Real> *md,
+                             TE domain_te, TE memory_te = TE::CC) {
+    RawMemoryIndexer idxer(0, domain, halo, md, domain_te, memory_te);
+    const int ni = idxer.ib.e - idxer.ib.s + 1;
+    const int nj = idxer.jb.e - idxer.jb.s + 1;
+    idxer.inner_length = ni * nj;
+    return idxer;
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  auto GetStartIndices(int outer_idx) const { return idxer(outer_idx * inner_length); }
+
+  KOKKOS_INLINE_FUNCTION
+  int GetNinnerRaw(int outer_idx) const {
+    auto [ks, js, is] = idxer(outer_idx * inner_length);
+    auto [ke, je, ie] = idxer(
+        std::min((outer_idx + 1) * inner_length - 1, static_cast<int>(idxer.size()) - 1));
+    return idxer_entire.GetFlatIdx(ke, je, ie) - idxer_entire.GetFlatIdx(ks, js, is) + 1;
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  int GetNouter() const {
+    return idxer.size() / inner_length + (idxer.size() % inner_length > 0);
+  }
+
+  int GetMaxNinnerRaw() const {
+    int max_ninner_raw{0};
+    for (int i = 0; i < GetNouter(); ++i) {
+      max_ninner_raw = std::max(max_ninner_raw, GetNinnerRaw(i));
+    }
+    return max_ninner_raw;
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  int GetStartingRawFlatIdx(int outer_idx) const {
+    auto [ks, js, is] = idxer(outer_idx * inner_length);
+    return idxer_entire.GetFlatIdx(ks, js, is);
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  auto GetCurrentIndices(int starting_raw_flat_idx, int inner_idx) const {
+    return idxer_entire(starting_raw_flat_idx + inner_idx);
+  }
+
+  int inner_length;
+  parthenon::IndexRange ib, jb, kb;
+  parthenon::Indexer3D idxer_entire;
+  parthenon::Indexer3D idxer;
+};
+
+#endif // EXAMPLE_DIFFUSION_RAW_MEMORY_INDEXER_HPP_
diff --git a/src/bvals/comms/boundary_communication.cpp b/src/bvals/comms/boundary_communication.cpp
index 076ab8085ca4e..77fe828f161a4 100644
--- a/src/bvals/comms/boundary_communication.cpp
+++ b/src/bvals/comms/boundary_communication.cpp
@@ -59,7 +59,7 @@ TaskStatus SendBoundBufsWithRestrictOption(std::shared_ptr<MeshData<Real>> &md,
     InitializeBufferCache<bound_type>(md, &(pmesh->boundary_comm_map), &cache, SendKey,
                                       true);
 
-  auto [rebuild, nbound, other_communication_unfinished] =
+  auto [rebuild, nbound, other_communication_unfinished, any_sparse] =
       CheckSendBufferCacheForRebuild<bound_type, true>(md);
 
   if (nbound == 0) {
@@ -101,56 +101,87 @@ TaskStatus SendBoundBufsWithRestrictOption(std::shared_ptr<MeshData<Real>> &md,
   auto &sending_nonzero_flags = cache.sending_non_zero_flags;
   auto &sending_nonzero_flags_h = cache.sending_non_zero_flags_h;
 
-  Kokkos::parallel_for(
-      PARTHENON_AUTO_LABEL,
-      Kokkos::TeamPolicy<>(parthenon::DevExecSpace(), nbound, Kokkos::AUTO),
-      KOKKOS_LAMBDA(parthenon::team_mbr_t team_member) {
-        const int b = team_member.league_rank();
+  if (any_sparse) {
+    Kokkos::parallel_for(
+        PARTHENON_AUTO_LABEL,
+        Kokkos::TeamPolicy<>(parthenon::DevExecSpace(), nbound, Kokkos::AUTO),
+        KOKKOS_LAMBDA(parthenon::team_mbr_t team_member) {
+          const int b = team_member.league_rank();
+
+          if (!bnd_info(b).allocated || bnd_info(b).same_to_same) {
+            Kokkos::single(Kokkos::PerTeam(team_member),
+                           [&]() { sending_nonzero_flags(b) = false; });
+            return;
+          }
+          Real threshold = bnd_info(b).var.allocation_threshold;
+          bool non_zero[3]{false, false, false};
+          int idx_offset = 0;
+          for (int it = 0; it < bnd_info(b).ntopological_elements; ++it) {
+            auto &idxer = bnd_info(b).idxer[it];
+            const int iel = static_cast<int>(bnd_info(b).topo_idx[it]) % 3;
+            const int Ni = idxer.template EndIdx<5>() - idxer.template StartIdx<5>() + 1;
+            Kokkos::parallel_reduce(
+                Kokkos::TeamThreadRange<>(team_member, idxer.size() / Ni),
+                [&](const int idx, bool &lnon_zero) {
+                  const auto [t, u, v, k, j, i] = idxer(idx * Ni);
+                  Real const *const var = &bnd_info(b).var(iel, t, u, v, k, j, i);
+                  Real *buf = &bnd_info(b).buf(idx * Ni + idx_offset);
 
-        if (!bnd_info(b).allocated || bnd_info(b).same_to_same) {
-          Kokkos::single(Kokkos::PerTeam(team_member),
-                         [&]() { sending_nonzero_flags(b) = false; });
-          return;
-        }
-        Real threshold = bnd_info(b).var.allocation_threshold;
-        bool non_zero[3]{false, false, false};
-        int idx_offset = 0;
-        for (int it = 0; it < bnd_info(b).ntopological_elements; ++it) {
-          auto &idxer = bnd_info(b).idxer[it];
-          const int iel = static_cast<int>(bnd_info(b).topo_idx[it]) % 3;
-          const int Ni = idxer.template EndIdx<5>() - idxer.template StartIdx<5>() + 1;
-          Kokkos::parallel_reduce(
-              Kokkos::TeamThreadRange<>(team_member, idxer.size() / Ni),
-              [&](const int idx, bool &lnon_zero) {
-                const auto [t, u, v, k, j, i] = idxer(idx * Ni);
-                Real *var = &bnd_info(b).var(iel, t, u, v, k, j, i);
-                Real *buf = &bnd_info(b).buf(idx * Ni + idx_offset);
-
-                Kokkos::parallel_for(Kokkos::ThreadVectorRange<>(team_member, Ni),
-                                     [&](int m) { buf[m] = var[m]; });
-
-                bool mnon_zero = false;
-                Kokkos::parallel_reduce(
-                    Kokkos::ThreadVectorRange<>(team_member, Ni),
-                    [&](int m, bool &llnon_zero) {
-                      llnon_zero = llnon_zero || (std::abs(buf[m]) >= threshold);
-                    },
-                    Kokkos::LOr<bool, parthenon::DevMemSpace>(mnon_zero));
-
-                lnon_zero = lnon_zero || mnon_zero;
-                if (bound_type == BoundaryType::flxcor_send) lnon_zero = true;
-              },
-              Kokkos::LOr<bool, parthenon::DevMemSpace>(non_zero[iel]));
-          idx_offset += idxer.size();
-        }
-        Kokkos::single(Kokkos::PerTeam(team_member), [&]() {
-          sending_nonzero_flags(b) = non_zero[0] || non_zero[1] || non_zero[2];
+                  Kokkos::parallel_for(Kokkos::ThreadVectorRange<>(team_member, Ni),
+                                       [&](int m) { buf[m] = var[m]; });
+
+                  bool mnon_zero = false;
+                  Kokkos::parallel_reduce(
+                      Kokkos::ThreadVectorRange<>(team_member, Ni),
+                      [&](int m, bool &llnon_zero) {
+                        llnon_zero = llnon_zero || (std::abs(buf[m]) >= threshold);
+                      },
+                      Kokkos::LOr<bool, parthenon::DevMemSpace>(mnon_zero));
+
+                  lnon_zero = lnon_zero || mnon_zero;
+                  if (bound_type == BoundaryType::flxcor_send) lnon_zero = true;
+                },
+                Kokkos::LOr<bool, parthenon::DevMemSpace>(non_zero[iel]));
+            idx_offset += idxer.size();
+          }
+          Kokkos::single(Kokkos::PerTeam(team_member), [&]() {
+            sending_nonzero_flags(b) = non_zero[0] || non_zero[1] || non_zero[2];
+          });
         });
-      });
+    if (Globals::sparse_config.enabled)
+      Kokkos::deep_copy(sending_nonzero_flags_h, sending_nonzero_flags);
+  } else {
+    Kokkos::parallel_for(
+        PARTHENON_AUTO_LABEL,
+        Kokkos::TeamPolicy<>(parthenon::DevExecSpace(), nbound, Kokkos::AUTO),
+        KOKKOS_LAMBDA(parthenon::team_mbr_t team_member) {
+          const int b = team_member.league_rank();
+
+          if (bnd_info(b).same_to_same) return;
+
+          int idx_offset = 0;
+          for (int it = 0; it < bnd_info(b).ntopological_elements; ++it) {
+            auto &idxer = bnd_info(b).idxer[it];
+            const int iel = static_cast<int>(bnd_info(b).topo_idx[it]) % 3;
+            const int Ni = idxer.template EndIdx<5>() - idxer.template StartIdx<5>() + 1;
+            Kokkos::parallel_for(
+                Kokkos::TeamThreadRange<>(team_member, idxer.size() / Ni),
+                [&](const int idx) {
+                  const auto [t, u, v, k, j, i] = idxer(idx * Ni);
+                  Real const *const var = &bnd_info(b).var(iel, t, u, v, k, j, i);
+                  Real *buf = &bnd_info(b).buf(idx * Ni + idx_offset);
+                  Kokkos::parallel_for(Kokkos::ThreadVectorRange<>(team_member, Ni),
+                                       [&](int m) { buf[m] = var[m]; });
+                });
+            idx_offset += idxer.size();
+          }
+        });
+    for (int ibuf = 0; ibuf < cache.buf_vec.size(); ++ibuf)
+      sending_nonzero_flags_h(ibuf) = true;
+  }
 
   // Send buffers
-  if (Globals::sparse_config.enabled)
-    Kokkos::deep_copy(sending_nonzero_flags_h, sending_nonzero_flags);
+
 #ifdef MPI_PARALLEL
   if (bound_type == BoundaryType::any || bound_type == BoundaryType::nonlocal)
     Kokkos::fence();
@@ -307,7 +338,18 @@ TaskStatus SetBounds(std::shared_ptr<MeshData<Real>> &md) {
           Real fac = ftemp; // Can't capture structured bindings
           const int iel = static_cast<int>(tel) % 3;
           const int Ni = idxer.template EndIdx<5>() - idxer.template StartIdx<5>() + 1;
-          if (bnd_info(b).buf_allocated && bnd_info(b).allocated) {
+          if (bnd_info(b).buf_allocated && bnd_info(b).allocated &&
+              tel == TopologicalElement::CC && lcoord_trans.IsIdentity()) {
+            Kokkos::parallel_for(
+                Kokkos::TeamThreadRange<>(team_member, idxer.size() / Ni),
+                [&](const int idx) {
+                  Real const *const buf = &bnd_info(b).buf(idx * Ni + idx_offset);
+                  const auto [t, u, v, k, j, i] = idxer(idx * Ni);
+                  Real *v_ = &var(iel, t, u, v, k, j, i);
+                  Kokkos::parallel_for(Kokkos::ThreadVectorRange<>(team_member, Ni),
+                                       [&](int m) { v_[m] = buf[m]; });
+                });
+          } else if (bnd_info(b).buf_allocated && bnd_info(b).allocated) {
             Kokkos::parallel_for(
                 Kokkos::TeamThreadRange<>(team_member, idxer.size() / Ni),
                 [&](const int idx) {
diff --git a/src/bvals/comms/bvals_utils.hpp b/src/bvals/comms/bvals_utils.hpp
index 131a7ee58338f..bc64a87d429d1 100644
--- a/src/bvals/comms/bvals_utils.hpp
+++ b/src/bvals/comms/bvals_utils.hpp
@@ -162,14 +162,17 @@ inline auto CheckSendBufferCacheForRebuild(std::shared_ptr<MeshData<Real>> md) {
   using namespace loops::shorthands;
   BvarsSubCache_t &cache = md->GetBvarsCache().GetSubCache(BOUND_TYPE, SENDER);
 
-  bool rebuild = false;
-  bool other_communication_unfinished = false;
-  int nbound = 0;
+  bool rebuild{false};
+  bool other_communication_unfinished{false};
+  int nbound{0};
+  bool any_sparse{false};
   ForEachBoundary<BOUND_TYPE>(md, [&](auto pmb, sp_mbd_t rc, const nb_t &nb,
                                       const sp_cv_t v) {
     const std::size_t ibuf = cache.idx_vec[nbound];
     auto &buf = *(cache.buf_vec[ibuf]);
 
+    any_sparse = any_sparse || v->IsSparse();
+
     if (!buf.IsAvailableForWrite()) other_communication_unfinished = true;
 
     if (v->IsAllocated()) {
@@ -186,7 +189,7 @@ inline auto CheckSendBufferCacheForRebuild(std::shared_ptr<MeshData<Real>> md) {
     }
     ++nbound;
   });
-  return std::make_tuple(rebuild, nbound, other_communication_unfinished);
+  return std::make_tuple(rebuild, nbound, other_communication_unfinished, any_sparse);
 }
 
 template <BoundaryType BOUND_TYPE, bool SENDER>
diff --git a/src/mesh/forest/logical_coordinate_transformation.cpp b/src/mesh/forest/logical_coordinate_transformation.cpp
index 34d4bec27a6f1..ee9fcd1be7e1f 100644
--- a/src/mesh/forest/logical_coordinate_transformation.cpp
+++ b/src/mesh/forest/logical_coordinate_transformation.cpp
@@ -114,6 +114,7 @@ ComposeTransformations(const LogicalCoordinateTransformation &first,
   for (int dir : {0, 1, 2})
     out.dir_connection_inverse[out.dir_connection[dir]] = dir;
   out.use_offset = first.use_offset && second.use_offset;
+  out.is_identity = first.is_identity && second.is_identity;
   return out;
 }
 
diff --git a/src/mesh/forest/logical_coordinate_transformation.hpp b/src/mesh/forest/logical_coordinate_transformation.hpp
index 0e0d0bc293cf6..b31795dc910ed 100644
--- a/src/mesh/forest/logical_coordinate_transformation.hpp
+++ b/src/mesh/forest/logical_coordinate_transformation.hpp
@@ -37,15 +37,19 @@ struct LogicalCoordinateTransformation {
   KOKKOS_INLINE_FUNCTION
   LogicalCoordinateTransformation()
       : dir_connection{0, 1, 2}, dir_connection_inverse{0, 1, 2},
-        dir_flip{false, false, false}, offset{0, 0, 0} {};
+        dir_flip{false, false, false}, offset{0, 0, 0}, is_identity{true} {};
 
   void SetDirection(CoordinateDirection origin, CoordinateDirection neighbor,
                     bool reversed = false) {
     dir_connection[origin - 1] = neighbor - 1;
     dir_connection_inverse[neighbor - 1] = origin - 1;
     dir_flip[origin - 1] = reversed;
+    is_identity = is_identity && (origin == neighbor) && !reversed;
   }
 
+  KOKKOS_INLINE_FUNCTION
+  bool IsIdentity() const { return is_identity; }
+
   LogicalLocation Transform(const LogicalLocation &loc_in,
                             std::int64_t destination) const;
   LogicalLocation InverseTransform(const LogicalLocation &loc_in,
@@ -98,6 +102,7 @@ struct LogicalCoordinateTransformation {
   }
 
   bool use_offset = false;
+  bool is_identity;
   std::array<int, 3> offset;
   std::array<int, 3> dir_connection, dir_connection_inverse;
   std::array<bool, 3> dir_flip;
diff --git a/src/mesh/meshblock.hpp b/src/mesh/meshblock.hpp
index facdf5415fae6..21c23006f805e 100644
--- a/src/mesh/meshblock.hpp
+++ b/src/mesh/meshblock.hpp
@@ -207,8 +207,7 @@ class MeshBlock : public std::enable_shared_from_this<MeshBlock> {
   BoundaryFlag boundary_flag[6];
 
   bool IsPhysicalBoundary(BoundaryFace bf) const {
-    // TODO(LFR): Should we only return true if this is set to user?
-    return boundary_flag[bf] != BoundaryFlag::block;
+    return boundary_flag[bf] == BoundaryFlag::user;
   }
 
   bool IsPhysicalBoundary() const {
diff --git a/src/solvers/solver_base.hpp b/src/solvers/solver_base.hpp
index af2b824eb2f8f..29c182b2d02fd 100644
--- a/src/solvers/solver_base.hpp
+++ b/src/solvers/solver_base.hpp
@@ -81,6 +81,27 @@ class SolverBase {
 
   const std::vector<std::string> &GetFieldLabels() const { return sol_fields; }
 
+  // We do not include a helper function for getting MeshData on the base container
+  // since it may contain an arbitrary set of fields
+
+  template <class T>
+  auto &AddRHSMeshData(Mesh *pmesh, T &&other, bool shallow = false) {
+    if (shallow)
+      return pmesh->mesh_data.AddShallow(GetRHSContainerLabel(), std::forward<T>(other),
+                                         GetFieldLabels());
+    return pmesh->mesh_data.Add(GetRHSContainerLabel(), std::forward<T>(other),
+                                GetFieldLabels());
+  }
+
+  template <class T>
+  auto &AddSolutionMeshData(Mesh *pmesh, T &&other, bool shallow = false) {
+    if (shallow)
+      return pmesh->mesh_data.AddShallow(GetSolutionContainerLabel(),
+                                         std::forward<T>(other), GetFieldLabels());
+    return pmesh->mesh_data.Add(GetSolutionContainerLabel(), std::forward<T>(other),
+                                GetFieldLabels());
+  }
+
   bool initial_guess_is_zero{false};
 
   static inline TimingAccumulatorDictionary solver_timings;
@@ -99,8 +120,8 @@ class SolverBase {
   std::string container_rhs;
 
   Real initial_residual{-1.0};
-  Real final_residual;
-  int final_iteration;
+  Real final_residual{-1.0};
+  int final_iteration{-1};
 };
 
 } // namespace solvers