diff --git a/CMakeLists.txt b/CMakeLists.txt
index 359e33e1c8..ef3220ee6a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -377,6 +377,7 @@ set(BOUT_SOURCES
     ./src/sys/expressionparser.cxx
     ./src/sys/generator_context.cxx
     ./src/sys/hyprelib.cxx
+    ./src/sys/hypre_interface.cxx
     ./src/sys/msg_stack.cxx
     ./src/sys/options.cxx
     ./src/sys/options/optionparser.hxx
diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx
index 832c3fbcc0..da222b9d7b 100644
--- a/include/bout/hypre_interface.hxx
+++ b/include/bout/hypre_interface.hxx
@@ -7,11 +7,18 @@
 
 #include "bout/assert.hxx"
 #include "bout/bout_enum_class.hxx"
+#include "bout/bout_types.hxx"
 #include "bout/boutcomm.hxx"
+#include "bout/boutexception.hxx"
 #include "bout/caliper_wrapper.hxx"
 #include "bout/field.hxx"
+#include "bout/field2d.hxx"
 #include "bout/globalindexer.hxx"
 #include "bout/hyprelib.hxx"
+#include "bout/options.hxx"
+#include "bout/output.hxx"
+#include "bout/paralleltransform.hxx"
+#include "bout/region.hxx"
 #include "bout/utils.hxx"
 
 #include "HYPRE.h"
@@ -21,8 +28,11 @@
 #include "HYPRE_utilities.h"
 #include "_hypre_utilities.h"
 
+#include <algorithm>
 #include <cmath>
+#include <iterator>
 #include <memory>
+#include <vector>
 
 // BOUT_ENUM_CLASS does not work inside namespaces
 BOUT_ENUM_CLASS(HYPRE_SOLVER_TYPE, gmres, bicgstab, pcg);
@@ -60,6 +70,139 @@ int checkHypreError(int error) {
 // TODO: set sizes
 // TODO: set contiguous blocks at once
 
+/// Wrapper around HYPRE_Complex, that calls HypreFree when destroyed.
+struct HypreComplexArray {
+  HYPRE_Complex* data;
+
+  HypreComplexArray(int n) { HypreMalloc(data, sizeof(HYPRE_Complex) * n); }
+
+  ~HypreComplexArray() { HypreFree(data); }
+};
+
+/// Shared pointter to a HypreComplexArray. When the last copy is destroyed
+/// the HYPRE_Complex array inside will be free'd.
+using BCValuesPtr = std::shared_ptr<HypreComplexArray>;
+
+/*!
+ * This function modifies the input for the HYPRE_IJMatrixSetValues() routine to
+ * eliminate the boundary condition equations (see below for details on how the
+ * equations are adjusted).  It modifies the arrays ncols, rows, cols, and
+ * values.  It also returns a row_indexes array.  This can then be passed to the
+ * HYPRE_IJMatrixSetValues2() routine to set up the matrix in hypre.
+ *
+ * The arguments nb and bi_array indicate the boundary equations.  The routine
+ * returns info needed to adjust the right-hand-side and solution vector through
+ * the functions AdjustRightHandSideEquations and AdjustSolutionEquations.
+ *
+ * NOTE: It may make sense from an organizational standpoint to collect many of
+ * these arguments in a structure of some sort.
+ *
+ * Notation, assumptions, and other details:
+ *
+ * - Boundary equation i is assumed to have two coefficients
+ *
+ *      b_ii * u_i + b_ij * u_j = rhs_i
+ *
+ * - We also assume that each boundary equation has only one retained unknown j
+ *   on the boundary row
+ *
+ *      b_ii * u_i + b_ij * u_j = rhs_i
+ *
+ * - Any number of other equations k may couple to u_i with coefficient a_ki
+ *
+ *      a_ki * u_i + a_kj * u_j + ... = rhs_k
+ *
+ * - Each such equation k is adjusted as follows:
+ *
+ *      a_kj = a_kj - a_ki * b_ij / b_ii
+ *      a_ki = 0
+ *
+ * - Boundary equations are adjusted to be identity equations in the matrix, but
+ *   the boundary coefficients (b_ii, b_ij) are returned for use later
+ *
+ * - Right-hand-side equations are adjusted in AdjustRightHandSideEquations() as
+ *   follows: rhs_k = rhs_k - a_ki * rhs_i / b_ii
+ *
+ * - Solution unknowns are adjusted at boundaries in AdjustSolutionEquations as
+ *   follows: u_i = (rhs_i - b_ij * u_j) / b_ii
+ *
+ * - Naming conventions: Arrays starting with 'b' are boundary equation arrays
+ *   indexed by 'bnum', and arrays starting with 'a' are non-boundary arrays
+ *   (interior matrix equations) indexed by 'anum'.  When 'num' is prefixed with
+ *   a row or column number 'i', 'j', or 'k', the array holds the corresponding
+ *   local data index for that row or column (e.g., an index into the local
+ *   solution vector).  Matrix coefficients are named as above, e.g., 'bij' is
+ *   the coefficient for b_ij.
+ *
+ *   NOTE: Implementation in src/sys/hypre_interface.cxx
+ */
+struct BoundaryElimination {
+  HYPRE_Int nb;
+  HYPRE_Int* binum_array;
+  HYPRE_Int* bjnum_array;
+  HYPRE_Int* bdep_array;
+  HYPRE_Complex* bii_array;
+  HYPRE_Complex* bij_array;
+  HYPRE_Int na;
+  HYPRE_Int* aoffset_array;
+  HYPRE_Int* aknum_array;
+  HYPRE_Complex* aki_array;
+  std::vector<HYPRE_Int> reduction_order;
+  std::vector<HYPRE_Int> expansion_order;
+
+  BoundaryElimination() = delete;
+
+  BoundaryElimination(HYPRE_Int nrows, HYPRE_Int* ncols, HYPRE_BigInt* rows,
+                      HYPRE_Int** row_indexes_ptr, HYPRE_BigInt* cols,
+                      HYPRE_Complex* values,
+                      HYPRE_Int nb,         // number of boundary equations
+                      HYPRE_Int* bi_array); // row i for each boundary equation
+
+  ~BoundaryElimination() {
+    // Free arrays
+    HypreFree(binum_array);
+    HypreFree(bjnum_array);
+    HypreFree(bdep_array);
+    HypreFree(bii_array);
+    HypreFree(bij_array);
+    HypreFree(aoffset_array);
+    HypreFree(aknum_array);
+    HypreFree(aki_array);
+  }
+
+  HYPRE_Int size() const { return nb; }
+
+  /// Copy boundary-row entries from an existing full-space vector.
+  BCValuesPtr copyBoundaryRowValues(const HYPRE_Complex* values) const;
+
+  /// Evaluate the original boundary equations using the supplied full-space vector.
+  BCValuesPtr evaluateBoundaryEquations(const HYPRE_Complex* values) const;
+
+  /// Applies in-place modification of the rhs array.
+  ///
+  /// Returns boundary values needed to reconstruct the full-space solution.
+  BCValuesPtr reduceRightHandSideInPlace(HYPRE_Complex* rhs) const;
+
+  /// Apply boundary conditions to the solution.
+  /// Uses the BCValuesPtr returned from reduceRightHandSideInPlace()
+  void expandSolutionInPlace(BCValuesPtr brhs, HYPRE_Complex* solution) const;
+
+  /// Reconstruct the action of the full operator from the reduced matrix-vector product.
+  void expandMatvecResultInPlace(BCValuesPtr boundary_operator_values,
+                                 BCValuesPtr full_boundary_values,
+                                 HYPRE_Complex* result) const;
+};
+
+/// A shared pointer to a BoundaryElimination object
+using BoundaryEliminationPtr = std::shared_ptr<BoundaryElimination>;
+
+struct BoundaryEliminationState {
+  BCValuesPtr interior_values;
+  BCValuesPtr boundary_values;
+};
+
+enum class HypreVectorReadMode { standard, solution, matvec };
+
 template <class T>
 class HypreVector {
   MPI_Comm comm;
@@ -72,6 +215,8 @@ class HypreVector {
   bool have_indices{false};
   HYPRE_BigInt* I{nullptr};
   HYPRE_Complex* V{nullptr};
+  HYPRE_Complex* workV{nullptr};
+  bool cache_current{false};
   HypreLib hyprelib{};
 
 public:
@@ -88,6 +233,7 @@ public:
     checkHypreError(HYPRE_IJVectorDestroy(hypre_vector));
     HypreFree(I);
     HypreFree(V);
+    HypreFree(workV);
   }
 
   // Disable copy, at least for now: not clear that HYPRE_IJVector is
@@ -99,16 +245,18 @@ public:
       : comm(other.comm), jlower(other.jlower), jupper(other.jupper), vsize(other.vsize),
         indexConverter(other.indexConverter), location(other.location),
         initialised(other.initialised), have_indices(other.have_indices), I(other.I),
-        V(other.V) {
+        V(other.V), workV(other.workV), cache_current(other.cache_current) {
     std::swap(hypre_vector, other.hypre_vector);
     std::swap(parallel_vector, other.parallel_vector);
     other.initialised = false;
     other.have_indices = false;
     other.I = nullptr;
     other.V = nullptr;
+    other.workV = nullptr;
+    other.cache_current = false;
   }
 
-  HypreVector<T>& operator=(HypreVector<T>&& other) {
+  HypreVector<T>& operator=(HypreVector<T>&& other) noexcept {
     comm = other.comm;
     jlower = other.jlower;
     jupper = other.jupper;
@@ -123,8 +271,12 @@ public:
     other.have_indices = false;
     I = other.I;
     V = other.V;
+    workV = other.workV;
+    cache_current = other.cache_current;
     other.I = nullptr;
     other.V = nullptr;
+    other.workV = nullptr;
+    other.cache_current = false;
     return *this;
   }
 
@@ -154,7 +306,9 @@ public:
     initialised = true;
     HypreMalloc(I, vsize * sizeof(HYPRE_BigInt));
     HypreMalloc(V, vsize * sizeof(HYPRE_Complex));
+    HypreMalloc(workV, vsize * sizeof(HYPRE_Complex));
     importValuesFromField(f);
+    assemble();
   }
 
   /// Construct a vector with given index set, but don't set any values
@@ -174,22 +328,70 @@ public:
     location = CELL_LOC::centre;
     HypreMalloc(I, vsize * sizeof(HYPRE_BigInt));
     HypreMalloc(V, vsize * sizeof(HYPRE_Complex));
+    HypreMalloc(workV, vsize * sizeof(HYPRE_Complex));
   }
 
-  void assemble() {
+  void assemble(const BoundaryElimination* boundary_elimination = nullptr,
+                BoundaryEliminationState* elimination_state = nullptr) {
     CALI_CXX_MARK_FUNCTION;
-    writeCacheToHypre();
+
+    parallel_vector = nullptr;
+    checkHypreError(HYPRE_IJVectorInitialize(hypre_vector));
+    writeCacheToHypre(boundary_elimination, elimination_state);
     checkHypreError(HYPRE_IJVectorAssemble(hypre_vector));
     checkHypreError(HYPRE_IJVectorGetObject(hypre_vector,
                                             reinterpret_cast<void**>(&parallel_vector)));
+    cache_current = true;
   }
 
-  void writeCacheToHypre() {
-    checkHypreError(HYPRE_IJVectorSetValues(hypre_vector, vsize, I, V));
+  void writeCacheToHypre(const BoundaryElimination* boundary_elimination = nullptr,
+                         BoundaryEliminationState* elimination_state = nullptr) {
+    HYPRE_Complex* values = V;
+    if (boundary_elimination != nullptr) {
+      ASSERT1(elimination_state != nullptr);
+      std::copy(V, V + vsize, workV);
+      values = workV;
+      elimination_state->interior_values =
+          boundary_elimination->reduceRightHandSideInPlace(values);
+      elimination_state->boundary_values = elimination_state->interior_values;
+    } else if (elimination_state != nullptr) {
+      elimination_state->interior_values = nullptr;
+      elimination_state->boundary_values = nullptr;
+    }
+    checkHypreError(HYPRE_IJVectorSetValues(hypre_vector, vsize, I, values));
   }
 
-  void readCacheFromHypre() {
-    checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, V));
+  void readCacheFromHypre(const BoundaryElimination* boundary_elimination = nullptr,
+                          const BoundaryEliminationState* elimination_state = nullptr,
+                          HypreVectorReadMode mode = HypreVectorReadMode::standard) {
+    HYPRE_Complex* values = V;
+    if ((boundary_elimination != nullptr) and (mode != HypreVectorReadMode::standard)) {
+      values = workV;
+    }
+    checkHypreError(HYPRE_IJVectorGetValues(hypre_vector, vsize, I, values));
+    if (boundary_elimination != nullptr) {
+      switch (mode) {
+      case HypreVectorReadMode::standard:
+        break;
+      case HypreVectorReadMode::solution:
+        ASSERT1(elimination_state != nullptr);
+        ASSERT1(elimination_state->interior_values != nullptr);
+        boundary_elimination->expandSolutionInPlace(elimination_state->interior_values,
+                                                    values);
+        std::copy(values, values + vsize, V);
+        break;
+      case HypreVectorReadMode::matvec:
+        ASSERT1(elimination_state != nullptr);
+        ASSERT1(elimination_state->interior_values != nullptr);
+        ASSERT1(elimination_state->boundary_values != nullptr);
+        boundary_elimination->expandMatvecResultInPlace(
+            elimination_state->interior_values, elimination_state->boundary_values,
+            values);
+        std::copy(values, values + vsize, V);
+        break;
+      }
+    }
+    cache_current = true;
   }
 
   T toField() {
@@ -198,7 +400,9 @@ public:
     T result(indexConverter->getMesh());
     result.allocate().setLocation(location);
 
-    readCacheFromHypre();
+    if ((!cache_current) and (parallel_vector != nullptr)) {
+      readCacheFromHypre();
+    }
     // Note that this only populates boundaries to a depth of 1
     int count = 0;
     BOUT_FOR_SERIAL(i, indexConverter->getRegionAll()) {
@@ -212,6 +416,7 @@ public:
     return result;
   }
 
+  // Loads values into I and V arrays but does not assemble
   void importValuesFromField(const T& f) {
     CALI_CXX_MARK_FUNCTION;
 
@@ -226,16 +431,22 @@ public:
     }
 
     ASSERT2(vec_i == vsize);
-    // writeCacheToHypre(); // redundant assemble already performs writeCacheToHypre
-    assemble();
     have_indices = true;
+    cache_current = true;
   }
 
   HYPRE_IJVector get() { return hypre_vector; }
   const HYPRE_IJVector& get() const { return hypre_vector; }
+  const HYPRE_Complex* getValues() const { return V; }
 
-  HYPRE_ParVector getParallel() { return parallel_vector; }
-  const HYPRE_ParVector& getParallel() const { return parallel_vector; }
+  HYPRE_ParVector getParallel() {
+    ASSERT1(parallel_vector != nullptr);
+    return parallel_vector;
+  }
+  const HYPRE_ParVector& getParallel() const {
+    ASSERT1(parallel_vector != nullptr);
+    return parallel_vector;
+  }
 
   class Element {
     HypreVector<T>* vector;
@@ -274,6 +485,7 @@ public:
       ASSERT3(std::isfinite(value_));
       value = value_;
       vector->V[vec_i] = value_;
+      vector->cache_current = true;
       return *this;
     }
     Element& operator+=(BoutReal value_) {
@@ -281,6 +493,7 @@ public:
       value += value_;
       ASSERT3(std::isfinite(value));
       vector->V[vec_i] += value_;
+      vector->cache_current = true;
       return *this;
     }
     operator BoutReal() const { return value; }
@@ -303,11 +516,20 @@ public:
 
   friend void swap(HypreVector<T>& lhs, HypreVector<T>& rhs) {
     using std::swap;
+    swap(lhs.comm, rhs.comm);
+    swap(lhs.jlower, rhs.jlower);
+    swap(lhs.jupper, rhs.jupper);
+    swap(lhs.vsize, rhs.vsize);
     swap(lhs.hypre_vector, rhs.hypre_vector);
     swap(lhs.parallel_vector, rhs.parallel_vector);
     swap(lhs.indexConverter, rhs.indexConverter);
     swap(lhs.location, rhs.location);
     swap(lhs.initialised, rhs.initialised);
+    swap(lhs.have_indices, rhs.have_indices);
+    swap(lhs.I, rhs.I);
+    swap(lhs.V, rhs.V);
+    swap(lhs.workV, rhs.workV);
+    swap(lhs.cache_current, rhs.cache_current);
   }
 };
 
@@ -327,6 +549,8 @@ class HypreMatrix {
   std::vector<HYPRE_BigInt>* I;
   std::vector<std::vector<HYPRE_BigInt>>* J;
   std::vector<std::vector<HYPRE_Complex>>* V;
+  bool elimBE{false};
+  BoundaryEliminationPtr boundary_elimination;
   HypreLib hyprelib{};
 
   // todo also take care of I,J,V
@@ -350,7 +574,8 @@ public:
         index_converter(other.index_converter), location(other.location),
         initialised(other.initialised), yoffset(other.yoffset),
         parallel_transform(other.parallel_transform), assembled(other.assembled),
-        num_rows(other.num_rows), I(other.I), J(other.J), V(other.V) {
+        num_rows(other.num_rows), I(other.I), J(other.J), V(other.V),
+        elimBE(other.elimBE), boundary_elimination(other.boundary_elimination) {
     std::swap(hypre_matrix, other.hypre_matrix);
     std::swap(parallel_matrix, other.parallel_matrix);
   }
@@ -371,6 +596,8 @@ public:
     I = other.I;
     J = other.J;
     V = other.V;
+    elimBE = other.elimBE;
+    boundary_elimination = other.boundary_elimination;
     return *this;
   }
 
@@ -659,7 +886,7 @@ public:
           pw.begin(), pw.end(), std::back_inserter(positions),
           [this, ny, nz](ParallelTransform::PositionsAndWeights p) -> HYPRE_Int {
             return this->index_converter->getGlobal(
-                ind_type(p.i * ny * nz + p.j * nz + p.k, ny, nz));
+                ind_type((((p.i * ny) + p.j) * nz) + p.k, ny, nz));
           });
       std::transform(pw.begin(), pw.end(), std::back_inserter(weights),
                      [](ParallelTransform::PositionsAndWeights p) -> HYPRE_Complex {
@@ -669,19 +896,21 @@ public:
     return Element(*this, global_row, global_column, positions, weights);
   }
 
+  void setElimBE() { elimBE = true; }
+
   void assemble() {
     CALI_CXX_MARK_FUNCTION;
 
     HYPRE_BigInt num_entries = 0;
-    HYPRE_BigInt* num_cols;
+    HYPRE_Int* num_cols;
     HYPRE_BigInt* cols;
     HYPRE_BigInt* rawI;
     HYPRE_Complex* vals;
 
-    HypreMalloc(num_cols, num_rows * sizeof(HYPRE_BigInt));
+    HypreMalloc(num_cols, num_rows * sizeof(HYPRE_Int));
     for (HYPRE_BigInt i = 0; i < num_rows; ++i) {
-      num_cols[i] = (*J)[i].size();
-      num_entries += (*J)[i].size();
+      num_cols[i] = static_cast<HYPRE_Int>((*J)[i].size());
+      num_entries += num_cols[i];
     }
 
     HypreMalloc(rawI, num_rows * sizeof(HYPRE_BigInt));
@@ -697,8 +926,33 @@ public:
         entry++;
       }
     }
-    checkHypreError(
-        HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals));
+
+    // Eliminate boundary condition equations in hypre SetValues input arguments
+    if (elimBE) {
+      HYPRE_Int* bi_array;
+      HYPRE_Int* row_indexes;
+      // There must be an easier way to get nb
+      int nb = 0;
+      BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) { nb++; }
+      HypreMalloc(bi_array, nb * sizeof(HYPRE_Int));
+      nb = 0;
+      BOUT_FOR_SERIAL(i, index_converter->getRegionBndry()) {
+        bi_array[nb] = index_converter->getGlobal(i);
+        nb++;
+      }
+
+      boundary_elimination = std::make_shared<BoundaryElimination>(
+          num_rows, num_cols, rawI, &row_indexes, cols, vals, nb, bi_array);
+      HypreFree(bi_array);
+
+      checkHypreError(HYPRE_IJMatrixSetValues2(*hypre_matrix, num_rows, num_cols, rawI,
+                                               row_indexes, cols, vals));
+      HypreFree(row_indexes);
+    } else {
+      boundary_elimination = nullptr;
+      checkHypreError(
+          HYPRE_IJMatrixSetValues(*hypre_matrix, num_rows, num_cols, rawI, cols, vals));
+    }
     checkHypreError(HYPRE_IJMatrixAssemble(*hypre_matrix));
     checkHypreError(HYPRE_IJMatrixGetObject(*hypre_matrix,
                                             reinterpret_cast<void**>(&parallel_matrix)));
@@ -743,12 +997,17 @@ public:
     result.I = I; // We want the pointer to transfer so this works like a view
     result.J = J;
     result.V = V;
+    result.elimBE = elimBE;
+    result.boundary_elimination = boundary_elimination;
 
     return result;
   }
 
   HYPRE_IJMatrix get() { return *hypre_matrix; }
   const HYPRE_IJMatrix& get() const { return *hypre_matrix; }
+  const BoundaryElimination* getBoundaryElimination() const {
+    return boundary_elimination.get();
+  }
 
   HYPRE_ParCSRMatrix getParallel() { return parallel_matrix; }
   const HYPRE_ParCSRMatrix& getParallel() const { return parallel_matrix; }
@@ -758,8 +1017,33 @@ public:
   void computeAxpby(double alpha, HypreVector<T>& x, double beta, HypreVector<T>& y) {
     CALI_CXX_MARK_FUNCTION;
 
+    BoundaryEliminationState elimination_state;
+    if (boundary_elimination != nullptr) {
+      elimination_state.interior_values =
+          boundary_elimination->evaluateBoundaryEquations(x.getValues());
+      elimination_state.boundary_values =
+          std::make_shared<HypreComplexArray>(boundary_elimination->size());
+      BCValuesPtr y_boundary_values;
+      if (beta != 0.0) {
+        y_boundary_values = boundary_elimination->copyBoundaryRowValues(y.getValues());
+      }
+      for (HYPRE_Int i = 0; i < boundary_elimination->size(); ++i) {
+        const HYPRE_Complex operator_value = elimination_state.interior_values->data[i];
+        elimination_state.interior_values->data[i] = alpha * operator_value;
+        elimination_state.boundary_values->data[i] =
+            alpha * operator_value
+            + (y_boundary_values != nullptr ? beta * y_boundary_values->data[i] : 0.0);
+      }
+    }
+
+    x.assemble();
+    y.assemble();
     checkHypreError(HYPRE_ParCSRMatrixMatvec(alpha, parallel_matrix, x.getParallel(),
                                              beta, y.getParallel()));
+    y.readCacheFromHypre(boundary_elimination.get(),
+                         boundary_elimination != nullptr ? &elimination_state : nullptr,
+                         boundary_elimination != nullptr ? HypreVectorReadMode::matvec
+                                                         : HypreVectorReadMode::standard);
   }
 
   // y = A*x
@@ -767,8 +1051,7 @@ public:
   void computeAx(HypreVector<T>& x, HypreVector<T>& y) {
     CALI_CXX_MARK_FUNCTION;
 
-    checkHypreError(HYPRE_ParCSRMatrixMatvec(1.0, parallel_matrix, x.getParallel(), 0.0,
-                                             y.getParallel()));
+    computeAxpby(1.0, x, 0.0, y);
   }
 };
 
@@ -878,6 +1161,31 @@ public:
     setMaxIter(
         options["maxits"].doc("Maximum iterations for Hypre solver").withDefault(10000));
 
+    switch (solver_type) {
+    case HYPRE_SOLVER_TYPE::gmres: {
+      HYPRE_ParCSRGMRESSetKDim(solver,
+                               options["kdim"]
+                                   .doc("Set the maximum size of the Krylov space")
+                                   .withDefault(30));
+
+      if (options["skip_real_residual_check"]
+              .doc("Skip the evaluation and the check of the actual residual?")
+              .withDefault<bool>(false)) {
+        HYPRE_GMRESSetSkipRealResidualCheck(solver, 1);
+      }
+      break;
+    }
+    case HYPRE_SOLVER_TYPE::bicgstab: {
+      break;
+    }
+    case HYPRE_SOLVER_TYPE::pcg: {
+      break;
+    }
+    default: {
+      throw BoutException("Unsupported hypre_solver_type {}", toString(solver_type));
+    }
+    }
+
     HYPRE_BoomerAMGCreate(&precon);
     HYPRE_BoomerAMGSetOldDefault(precon);
 #if BOUT_HAS_CUDA
@@ -962,6 +1270,11 @@ public:
     ASSERT2(A != nullptr);
     ASSERT2(x != nullptr);
     ASSERT2(b != nullptr);
+    BoundaryEliminationState elimination_state;
+    const auto* boundary_elimination = A->getBoundaryElimination();
+    b->assemble(boundary_elimination,
+                boundary_elimination != nullptr ? &elimination_state : nullptr);
+    x->assemble();
     if (not solver_setup) {
       checkHypreError(
           solverSetup(solver, A->getParallel(), b->getParallel(), x->getParallel()));
@@ -971,6 +1284,12 @@ public:
     solve_err = checkHypreError(
         solverSolve(solver, A->getParallel(), b->getParallel(), x->getParallel()));
 
+    x->readCacheFromHypre(boundary_elimination,
+                          boundary_elimination != nullptr ? &elimination_state : nullptr,
+                          boundary_elimination != nullptr
+                              ? HypreVectorReadMode::solution
+                              : HypreVectorReadMode::standard);
+
     return solve_err;
   }
 }; // class HypreSystem
diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx
index 6a321b0adf..9138c5a647 100644
--- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx
+++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx
@@ -3,7 +3,7 @@
  *                           Using Hypre Solvers
  *
  **************************************************************************
- * Copyright 2021 - 2025 BOUT++ contributors
+ * Copyright 2021 - 2026 BOUT++ contributors
  *
  * Contact: Ben Dudson, dudson2@llnl.gov
  *
@@ -30,17 +30,31 @@
 #include "hypre3d_laplace.hxx"
 
 #include <bout/assert.hxx>
+#include <bout/bout_types.hxx>
 #include <bout/boutcomm.hxx>
+#include <bout/boutexception.hxx>
 #include <bout/caliper_wrapper.hxx>
 #include <bout/derivs.hxx>
+#include <bout/field2d.hxx>
+#include <bout/field3d.hxx>
+#include <bout/globalindexer.hxx>
 #include <bout/hypre_interface.hxx>
 #include <bout/mesh.hxx>
 #include <bout/operatorstencil.hxx>
+#include <bout/paralleltransform.hxx>
+#include <bout/region.hxx>
 #include <bout/solver.hxx>
+#include <bout/sys/range.hxx>
 #include <bout/sys/timer.hxx>
 #include <bout/utils.hxx>
 
+#include <algorithm>
 #include <cmath>
+#include <iterator>
+#include <memory>
+#include <set>
+#include <string>
+#include <vector>
 
 LaplaceHypre3d::LaplaceHypre3d(Options* opt, const CELL_LOC loc, Mesh* mesh_in, Solver*)
     : Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0), Ex(0.0), Ez(0.0),
@@ -220,15 +234,13 @@ Field3D LaplaceHypre3d::solve(const Field3D& b_in, const Field3D& x0) {
 
   rhs.importValuesFromField(b);
   solution.importValuesFromField(x0);
-  rhs.assemble();
-  solution.assemble();
 
   CALI_MARK_END("LaplaceHypre3d_solve:vectorAssemble");
 
   CALI_MARK_BEGIN("LaplaceHypre3d_solve:solve");
   // Invoke solver
   {
-    Timer timer("hypresolve");
+    const Timer timer("hypresolve");
     linearSystem.solve();
   }
 
@@ -305,7 +317,7 @@ void LaplaceHypre3d::updateMatrix3D() {
     }
 
     BoutReal C_d2f_dx2 = coords->g11()[l];
-    BoutReal C_d2f_dy2 = (coords->g22()[l] - 1.0 / coords->g_22()[l]);
+    BoutReal C_d2f_dy2 = (coords->g22()[l] - (1.0 / coords->g_22()[l]));
     BoutReal C_d2f_dz2 = coords->g33()[l];
     if (issetD) {
       C_d2f_dx2 *= D[l];
@@ -360,7 +372,7 @@ void LaplaceHypre3d::updateMatrix3D() {
   // Must add these (rather than assign) so that elements used in
   // interpolation don't overwrite each other.
   BOUT_FOR_SERIAL(l, indexer->getRegionNobndry()) {
-    BoutReal C_df_dy = (coords->G2()[l] - dJ_dy[l] / coords->J()[l]);
+    BoutReal C_df_dy = (coords->G2()[l] - (dJ_dy[l] / coords->J()[l]));
     if (issetD) {
       C_df_dy *= D[l];
     }
@@ -371,7 +383,7 @@ void LaplaceHypre3d::updateMatrix3D() {
                  / C1[l];
     }
 
-    BoutReal C_d2f_dy2 = (coords->g22()[l] - 1.0 / coords->g_22()[l]);
+    BoutReal C_d2f_dy2 = (coords->g22()[l] - (1.0 / coords->g_22()[l]));
     if (issetD) {
       C_d2f_dy2 *= D[l];
     }
@@ -411,6 +423,7 @@ void LaplaceHypre3d::updateMatrix3D() {
     operator3D.ydown(ydown)(l, l.ym().zp()) += -C_d2f_dydz;
     operator3D.ydown(ydown)(l, l.ym().zm()) += C_d2f_dydz;
   }
+  operator3D.setElimBE();
   operator3D.assemble();
 
   if (print_matrix) {
@@ -440,7 +453,7 @@ OperatorStencil<Ind3D> LaplaceHypre3d::getStencil(Mesh* localmesh,
                                                           : p.k};
                  });
 
-  OffsetInd3D zero;
+  const OffsetInd3D zero;
 
   // Add interior cells
   const std::vector<OffsetInd3D> interpolatedUpElements = {
@@ -544,7 +557,7 @@ void LaplaceHypre3d::outputVars(Options& output_options,
                                 const std::string& time_dimension) const {
   BoutReal mean_iterations = 0.0;
   BoutReal mean_amg_iterations = 0.0;
-  BoutReal rel_res_norm = linearSystem.getFinalRelResNorm();
+  const BoutReal rel_res_norm = linearSystem.getFinalRelResNorm();
 
   if (n_solves > 0) {
     // Calculate average
diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx
index d58ef6f688..a5ecaab45e 100644
--- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx
+++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx
@@ -145,7 +145,7 @@ public:
   }
 
   // Return a reference to the matrix objects representing the Laplace
-  // operator. These will be (re)construct if necessary.
+  // operator. These will be (re)constructed if necessary.
   bout::HypreMatrix<Field3D>& getMatrix3D();
   IndexerPtr<Field3D> getIndexer() { return indexer; }
 
diff --git a/src/sys/hypre_interface.cxx b/src/sys/hypre_interface.cxx
new file mode 100644
index 0000000000..101484913a
--- /dev/null
+++ b/src/sys/hypre_interface.cxx
@@ -0,0 +1,261 @@
+
+#include "bout/build_defines.hxx"
+
+#if BOUT_HAS_HYPRE
+
+#include "bout/hypre_interface.hxx"
+
+#include <numeric>
+#include <unordered_map>
+
+namespace bout {
+
+BoundaryElimination::BoundaryElimination(HYPRE_Int nrows, HYPRE_Int* ncols,
+                                         HYPRE_BigInt* rows, HYPRE_Int** row_indexes_ptr,
+                                         HYPRE_BigInt* cols, HYPRE_Complex* values,
+                                         HYPRE_Int nb, HYPRE_Int* bi_array)
+    : nb(nb) {
+  HYPRE_Int* row_indexes;
+  const auto find_local_row = [nrows, rows](HYPRE_BigInt row) -> HYPRE_Int {
+    auto row_position = std::lower_bound(rows, rows + nrows, row);
+    if ((row_position == rows + nrows) || (*row_position != row)) {
+      throw BoutException("Could not find local row {} while constructing boundary "
+                          "elimination data",
+                          row);
+    }
+    return static_cast<HYPRE_Int>(std::distance(rows, row_position));
+  };
+
+  // Create the row_indexes array
+  HypreMalloc(row_indexes, sizeof(HYPRE_Int) * nrows);
+  row_indexes[0] = 0;
+  for (HYPRE_Int i = 1; i < nrows; i++) {
+    row_indexes[i] = row_indexes[i - 1] + ncols[i - 1];
+  }
+
+  // Allocate arrays
+  HypreMalloc(binum_array, sizeof(HYPRE_Int) * nb);
+  HypreMalloc(bjnum_array, sizeof(HYPRE_Int) * nb);
+  HypreMalloc(bdep_array, sizeof(HYPRE_Int) * nb);
+  HypreMalloc(bii_array, sizeof(HYPRE_Complex) * nb);
+  HypreMalloc(bij_array, sizeof(HYPRE_Complex) * nb);
+
+  std::vector<HYPRE_Int> boundary_at_row(nrows, -1);
+  for (HYPRE_Int bnum = 0; bnum < nb; bnum++) {
+    boundary_at_row[find_local_row(bi_array[bnum])] = bnum;
+  }
+
+  struct BoundaryOccurrence {
+    HYPRE_Int rownum;
+    HYPRE_Int row_offset;
+  };
+  std::unordered_map<HYPRE_BigInt, HYPRE_Int> boundary_number_for_row;
+  boundary_number_for_row.reserve(nb);
+  for (HYPRE_Int bnum = 0; bnum < nb; bnum++) {
+    boundary_number_for_row.emplace(bi_array[bnum], bnum);
+  }
+
+  std::vector<std::vector<BoundaryOccurrence>> occurrences_by_boundary(nb);
+  for (HYPRE_Int rownum = 0; rownum < nrows; rownum++) {
+    const HYPRE_Int row_start = row_indexes[rownum];
+    for (HYPRE_Int m = 0; m < ncols[rownum]; m++) {
+      const auto boundary_position = boundary_number_for_row.find(cols[row_start + m]);
+      if (boundary_position != boundary_number_for_row.end()) {
+        occurrences_by_boundary[boundary_position->second].push_back({rownum, m});
+      }
+    }
+  }
+
+  std::vector<HYPRE_BigInt> bjrow_array(nb, -1);
+  for (HYPRE_Int bnum = 0; bnum < nb; bnum++) {
+    // Get boundary equation information and adjust boundary equations
+    HYPRE_Int i = bi_array[bnum];
+    const HYPRE_Int binum = find_local_row(i);
+    HYPRE_Int bcoeffnum = row_indexes[binum];
+    HYPRE_Complex bii{0.0}, bij{0.0};
+    HYPRE_BigInt j = -1;
+
+    for (HYPRE_Int m = 0; m < ncols[binum]; m++) {
+      if (cols[bcoeffnum + m] == i) {
+        bii = values[bcoeffnum + m];
+        values[bcoeffnum + m] = -1.0; // Identity equation (negative definite matrix)
+      } else {
+        j = cols[bcoeffnum + m];
+        bij = values[bcoeffnum + m];
+        values[bcoeffnum + m] = 0.0; // Identity equation
+      }
+    }
+    if (j < 0) {
+      throw BoutException("Boundary row {} does not contain a retained neighbour", i);
+    }
+    if (bii == 0.0) {
+      throw BoutException("Boundary row {} has zero diagonal coefficient", i);
+    }
+    ncols[binum] = 1; // Identity equation
+
+    // Update arrays
+    binum_array[bnum] = binum;
+    bjnum_array[bnum] = find_local_row(j);
+    bdep_array[bnum] = boundary_at_row[bjnum_array[bnum]];
+    bii_array[bnum] = bii;
+    bij_array[bnum] = bij;
+    bjrow_array[bnum] = j;
+  }
+
+  std::vector<HYPRE_Int> aoffsets(nb + 1, 0);
+  std::vector<HYPRE_Int> aknums;
+  std::vector<HYPRE_Complex> akis;
+  for (HYPRE_Int bnum = 0; bnum < nb; bnum++) {
+    const HYPRE_Int binum = binum_array[bnum];
+    const HYPRE_Int bjnum = bjnum_array[bnum];
+    const HYPRE_BigInt j = bjrow_array[bnum];
+
+    for (const auto& occurrence : occurrences_by_boundary[bnum]) {
+      const HYPRE_Int aknum = occurrence.rownum;
+      if (aknum == binum) {
+        continue;
+      }
+      if ((boundary_at_row[aknum] >= 0) && (aknum != bjnum)) {
+        continue;
+      }
+
+      const HYPRE_Int acoeffnum = row_indexes[aknum];
+      const HYPRE_Int aki_position = acoeffnum + occurrence.row_offset;
+      const HYPRE_Complex aki = values[aki_position];
+      if (aki == 0.0) {
+        continue;
+      }
+
+      HYPRE_Int j_position = -1;
+      for (HYPRE_Int m = 0; m < ncols[aknum]; m++) {
+        if (cols[acoeffnum + m] == j) {
+          j_position = acoeffnum + m;
+          break;
+        }
+      }
+      if (j_position < 0) {
+        continue;
+      }
+
+      values[aki_position] = 0.0; // Eliminate coupling to boundary equation
+      values[j_position] -= aki * bij_array[bnum] / bii_array[bnum];
+
+      aknums.push_back(aknum);
+      akis.push_back(aki);
+    }
+    aoffsets[bnum + 1] = static_cast<HYPRE_Int>(aknums.size());
+  }
+
+  na = static_cast<HYPRE_Int>(aknums.size());
+  HypreMalloc(aoffset_array, sizeof(HYPRE_Int) * (nb + 1));
+  std::copy(aoffsets.begin(), aoffsets.end(), aoffset_array);
+  HypreMalloc(aknum_array, sizeof(HYPRE_Int) * na);
+  std::copy(aknums.begin(), aknums.end(), aknum_array);
+  HypreMalloc(aki_array, sizeof(HYPRE_Complex) * na);
+  std::copy(akis.begin(), akis.end(), aki_array);
+
+  std::vector<HYPRE_Int> dependency_depth(nb, -1);
+  const auto get_depth = [&dependency_depth, this](const auto& self,
+                                                   HYPRE_Int bnum) -> HYPRE_Int {
+    if (dependency_depth[bnum] >= 0) {
+      return dependency_depth[bnum];
+    }
+    const HYPRE_Int dep = bdep_array[bnum];
+    if (dep < 0) {
+      dependency_depth[bnum] = 0;
+    } else {
+      dependency_depth[bnum] = 1 + self(self, dep);
+    }
+    return dependency_depth[bnum];
+  };
+  for (HYPRE_Int bnum = 0; bnum < nb; bnum++) {
+    get_depth(get_depth, bnum);
+  }
+
+  expansion_order.resize(nb);
+  std::iota(expansion_order.begin(), expansion_order.end(), 0);
+  std::sort(expansion_order.begin(), expansion_order.end(),
+            [&dependency_depth](HYPRE_Int lhs, HYPRE_Int rhs) {
+              if (dependency_depth[lhs] != dependency_depth[rhs]) {
+                return dependency_depth[lhs] < dependency_depth[rhs];
+              }
+              return lhs < rhs;
+            });
+  reduction_order = expansion_order;
+  std::reverse(reduction_order.begin(), reduction_order.end());
+
+  // Set return arguments
+  *row_indexes_ptr = row_indexes;
+}
+
+BCValuesPtr
+BoundaryElimination::copyBoundaryRowValues(const HYPRE_Complex* values) const {
+  BCValuesPtr boundary_values = std::make_shared<HypreComplexArray>(nb);
+
+  for (HYPRE_Int bnum = 0; bnum < nb; bnum++) {
+    boundary_values->data[bnum] = values[binum_array[bnum]];
+  }
+
+  return boundary_values;
+}
+
+BCValuesPtr
+BoundaryElimination::evaluateBoundaryEquations(const HYPRE_Complex* values) const {
+  BCValuesPtr boundary_values = std::make_shared<HypreComplexArray>(nb);
+
+  for (HYPRE_Int bnum = 0; bnum < nb; bnum++) {
+    boundary_values->data[bnum] = (bii_array[bnum] * values[binum_array[bnum]])
+                                  + (bij_array[bnum] * values[bjnum_array[bnum]]);
+  }
+
+  return boundary_values;
+}
+
+BCValuesPtr BoundaryElimination::reduceRightHandSideInPlace(HYPRE_Complex* rhs) const {
+
+  // Allocate array to store boundary row values
+  BCValuesPtr brhs = copyBoundaryRowValues(rhs);
+
+  for (HYPRE_Int bnum : reduction_order) {
+    for (HYPRE_Int anum = aoffset_array[bnum]; anum < aoffset_array[bnum + 1]; anum++) {
+      HYPRE_Int aknum = aknum_array[anum];
+      rhs[aknum] -= aki_array[anum] * brhs->data[bnum] / bii_array[bnum];
+    }
+    if (bdep_array[bnum] >= 0) {
+      brhs->data[bdep_array[bnum]] = rhs[bjnum_array[bnum]];
+    }
+  }
+
+  return brhs;
+}
+
+void BoundaryElimination::expandSolutionInPlace(BCValuesPtr brhs,
+                                                HYPRE_Complex* solution) const {
+
+  for (HYPRE_Int bnum : expansion_order) {
+    HYPRE_Int binum = binum_array[bnum];
+    HYPRE_Int bjnum = bjnum_array[bnum];
+    solution[binum] =
+        (brhs->data[bnum] - bij_array[bnum] * solution[bjnum]) / bii_array[bnum];
+  }
+}
+
+void BoundaryElimination::expandMatvecResultInPlace(BCValuesPtr boundary_operator_values,
+                                                    BCValuesPtr full_boundary_values,
+                                                    HYPRE_Complex* result) const {
+  for (HYPRE_Int bnum = 0; bnum < nb; bnum++) {
+    for (HYPRE_Int anum = aoffset_array[bnum]; anum < aoffset_array[bnum + 1]; anum++) {
+      HYPRE_Int aknum = aknum_array[anum];
+      result[aknum] +=
+          aki_array[anum] * boundary_operator_values->data[bnum] / bii_array[bnum];
+    }
+  }
+
+  for (HYPRE_Int bnum = 0; bnum < nb; bnum++) {
+    result[binum_array[bnum]] = full_boundary_values->data[bnum];
+  }
+}
+
+} // namespace bout
+
+#endif // BOUT_HAS_HYPRE
diff --git a/tests/integrated/test-laplace-hypre3d/test_laplace_hypre3d.py b/tests/integrated/test-laplace-hypre3d/test_laplace_hypre3d.py
index e2fa11c13a..014b91ef00 100755
--- a/tests/integrated/test-laplace-hypre3d/test_laplace_hypre3d.py
+++ b/tests/integrated/test-laplace-hypre3d/test_laplace_hypre3d.py
@@ -17,7 +17,7 @@ def test_laplace_hypre3d():
 
     success = True
     for directory, nproc in test_directories:
-        command = "test-laplace3d -d " + directory
+        command = "./test-laplace3d -d " + directory
         print("running on", nproc, "processors:", command)
         launch_safe(command, nproc=nproc)
 
diff --git a/tests/unit/include/bout/test_hypre_interface.cxx b/tests/unit/include/bout/test_hypre_interface.cxx
index ec22140647..2be393f1b6 100644
--- a/tests/unit/include/bout/test_hypre_interface.cxx
+++ b/tests/unit/include/bout/test_hypre_interface.cxx
@@ -1,15 +1,19 @@
+#include "bout/build_defines.hxx"
+
 #if BOUT_HAS_HYPRE
 
 #include "HYPRE.h"
 #include "HYPRE_IJ_mv.h"
 #include "HYPRE_parcsr_ls.h"
 
+#include <array>
 #include <cmath>
+#include <set>
+#include <vector>
 
 #include "fake_mesh_fixture.hxx"
 #include "test_extras.hxx"
 
-#include "gmock/gmock.h"
 #include "gtest/gtest.h"
 
 #include "bout/field3d.hxx"
@@ -154,13 +158,30 @@ TYPED_TEST(HypreVectorTest, Swap) {
 //////////////////////////////////////////////////
 // HypreMatrix tests
 
-class MockTransform : public ParallelTransformIdentity {
+class TestParallelTransform : public ParallelTransformIdentity {
 public:
-  explicit MockTransform(Mesh& mesh_in) : ParallelTransformIdentity(mesh_in){};
-  MOCK_METHOD(std::vector<PositionsAndWeights>, getWeightsForYUpApproximation,
-              (int i, int j, int k), (override));
-  MOCK_METHOD(std::vector<PositionsAndWeights>, getWeightsForYDownApproximation,
-              (int i, int j, int k), (override));
+  explicit TestParallelTransform(Mesh& mesh_in) : ParallelTransformIdentity(mesh_in) {}
+
+  std::vector<PositionsAndWeights> getWeightsForYUpApproximation(int i, int j,
+                                                                 int k) override {
+    last_up = {i, j, k};
+    up_calls += 1;
+    return y_up_weights;
+  }
+
+  std::vector<PositionsAndWeights> getWeightsForYDownApproximation(int i, int j,
+                                                                   int k) override {
+    last_down = {i, j, k};
+    down_calls += 1;
+    return y_down_weights;
+  }
+
+  std::vector<PositionsAndWeights> y_up_weights;
+  std::vector<PositionsAndWeights> y_down_weights;
+  std::array<int, 3> last_up{-1, -1, -1};
+  std::array<int, 3> last_down{-1, -1, -1};
+  int up_calls{0};
+  int down_calls{0};
 };
 
 template <class T>
@@ -169,7 +190,7 @@ class HypreMatrixTest : public FakeMeshFixture {
   WithQuietOutput all{output};
   T field;
   IndexerPtr<T> indexer;
-  MockTransform* pt{nullptr};
+  TestParallelTransform* pt{nullptr};
   std::vector<ParallelTransform::PositionsAndWeights> yUpWeights, yDownWeights;
   using ind_type = typename T::ind_type;
   ind_type indexA, indexB, iWU0, iWU1, iWU2, iWD0, iWD1, iWD2;
@@ -199,7 +220,8 @@ class HypreMatrixTest : public FakeMeshFixture {
     iWU1 = indexB;
     iWU2 = indexB.xp();
 
-    auto transform = bout::utils::make_unique<MockTransform>(*bout::globals::mesh);
+    auto transform =
+        bout::utils::make_unique<TestParallelTransform>(*bout::globals::mesh);
     ParallelTransform::PositionsAndWeights wUp0 = {iWU0.x(), iWU0.y(), iWU0.z(), 0.5},
                                            wUp1 = {iWU1.x(), iWU1.y(), iWU1.z(), 1.0},
                                            wUp2 = {iWU2.x(), iWU2.y(), iWU2.z(), 0.5},
@@ -208,8 +230,14 @@ class HypreMatrixTest : public FakeMeshFixture {
                                            wDown2 = {iWD2.x(), iWD2.y(), iWD2.z(), 0.5};
     yUpWeights = {wUp0, wUp1, wUp2};
     yDownWeights = {wDown0, wDown1, wDown2};
-    pt = transform.get();
-    field.getCoordinates()->setParallelTransform(std::move(transform));
+    transform->y_up_weights = yUpWeights;
+    transform->y_down_weights = yDownWeights;
+    bout::globals::mesh->getCoordinates()->setParallelTransform(std::move(transform));
+    pt = dynamic_cast<TestParallelTransform*>(
+        &bout::globals::mesh->getCoordinates()->getParallelTransform());
+    if (pt == nullptr) {
+      throw BoutException("Failed to install TestParallelTransform in HypreMatrixTest");
+    }
   }
   virtual ~HypreMatrixTest() = default;
 };
@@ -307,7 +335,9 @@ TYPED_TEST(HypreMatrixTest, SetElements) {
       HYPRE_Int ncolumns{1};
       HYPRE_Complex value;
       BOUT_OMP_SAFE(critical)
-      { HYPRE_IJMatrixGetValues(raw_matrix, 1, &ncolumns, &i_index, &j_index, &value); }
+      {
+        HYPRE_IJMatrixGetValues(raw_matrix, 1, &ncolumns, &i_index, &j_index, &value);
+      }
       if (i == j) {
         EXPECT_EQ(static_cast<BoutReal>(value),
                   static_cast<BoutReal>(this->indexer->getGlobal(i)));
@@ -382,8 +412,6 @@ auto IsHypreMatrixEqual(const HypreMatrix<T>& matrix, const HypreMatrix<T>& refe
 }
 
 TYPED_TEST(HypreMatrixTest, YUp) {
-  using namespace ::testing;
-
   HypreMatrix<TypeParam> matrix(this->indexer);
 
   if constexpr (std::is_same_v<TypeParam, FieldPerp>) {
@@ -392,15 +420,11 @@ TYPED_TEST(HypreMatrixTest, YUp) {
   }
 
   HypreMatrix<TypeParam> expected(this->indexer);
-  MockTransform* transform = this->pt;
   const BoutReal value = 42.0;
 
   if constexpr (std::is_same_v<TypeParam, Field2D>) {
     expected(this->indexA, this->indexB) = value;
   } else {
-    EXPECT_CALL(*transform, getWeightsForYUpApproximation(
-                                this->indexB.x(), this->indexA.y(), this->indexB.z()))
-        .WillOnce(Return(this->yUpWeights));
     expected(this->indexA, this->iWU0) = this->yUpWeights[0].weight * value;
     expected(this->indexA, this->iWU1) = this->yUpWeights[1].weight * value;
     expected(this->indexA, this->iWU2) = this->yUpWeights[2].weight * value;
@@ -408,6 +432,13 @@ TYPED_TEST(HypreMatrixTest, YUp) {
 
   matrix.yup()(this->indexA, this->indexB) = value;
 
+  if constexpr (std::is_same_v<TypeParam, Field3D>) {
+    EXPECT_EQ(this->pt->up_calls, 1);
+    EXPECT_EQ(this->pt->last_up[0], this->indexB.x());
+    EXPECT_EQ(this->pt->last_up[1], this->indexA.y());
+    EXPECT_EQ(this->pt->last_up[2], this->indexB.z());
+  }
+
   expected.assemble();
   matrix.assemble();
 
@@ -415,8 +446,6 @@ TYPED_TEST(HypreMatrixTest, YUp) {
 }
 
 TYPED_TEST(HypreMatrixTest, YDown) {
-  using namespace ::testing;
-
   HypreMatrix<TypeParam> matrix(this->indexer);
 
   if constexpr (std::is_same_v<TypeParam, FieldPerp>) {
@@ -425,15 +454,11 @@ TYPED_TEST(HypreMatrixTest, YDown) {
   }
 
   HypreMatrix<TypeParam> expected(this->indexer);
-  MockTransform* transform = this->pt;
   const BoutReal value = 42.0;
 
   if constexpr (std::is_same_v<TypeParam, Field2D>) {
     expected(this->indexB, this->indexA) = value;
   } else {
-    EXPECT_CALL(*transform, getWeightsForYDownApproximation(
-                                this->indexA.x(), this->indexB.y(), this->indexA.z()))
-        .WillOnce(Return(this->yDownWeights));
     expected(this->indexB, this->iWD0) = this->yDownWeights[0].weight * value;
     expected(this->indexB, this->iWD1) = this->yDownWeights[1].weight * value;
     expected(this->indexB, this->iWD2) = this->yDownWeights[2].weight * value;
@@ -441,6 +466,13 @@ TYPED_TEST(HypreMatrixTest, YDown) {
 
   matrix.ydown()(this->indexB, this->indexA) = value;
 
+  if constexpr (std::is_same_v<TypeParam, Field3D>) {
+    EXPECT_EQ(this->pt->down_calls, 1);
+    EXPECT_EQ(this->pt->last_down[0], this->indexA.x());
+    EXPECT_EQ(this->pt->last_down[1], this->indexB.y());
+    EXPECT_EQ(this->pt->last_down[2], this->indexA.z());
+  }
+
   expected.assemble();
   matrix.assemble();
 
@@ -462,4 +494,219 @@ TYPED_TEST(HypreMatrixTest, YNext0) {
   EXPECT_TRUE(IsHypreMatrixEqual(matrix, expected));
 }
 
+namespace {
+
+struct RawBoundaryEliminationSystem {
+  std::vector<HYPRE_Int> ncols;
+  std::vector<HYPRE_BigInt> rows;
+  std::vector<HYPRE_BigInt> cols;
+  std::vector<HYPRE_Complex> values;
+  std::vector<HYPRE_Int> boundary_rows;
+  HYPRE_Int* row_indexes{nullptr};
+  std::unique_ptr<bout::BoundaryElimination> elimination;
+
+  RawBoundaryEliminationSystem(std::initializer_list<HYPRE_Int> ncols_in,
+                               std::initializer_list<HYPRE_BigInt> rows_in,
+                               std::initializer_list<HYPRE_BigInt> cols_in,
+                               std::initializer_list<HYPRE_Complex> values_in,
+                               std::initializer_list<HYPRE_Int> boundary_rows_in)
+      : ncols(ncols_in), rows(rows_in), cols(cols_in), values(values_in),
+        boundary_rows(boundary_rows_in) {
+    elimination = std::make_unique<bout::BoundaryElimination>(
+        static_cast<HYPRE_Int>(rows.size()), ncols.data(), rows.data(), &row_indexes,
+        cols.data(), values.data(), static_cast<HYPRE_Int>(boundary_rows.size()),
+        boundary_rows.data());
+  }
+
+  ~RawBoundaryEliminationSystem() {
+    if (row_indexes != nullptr) {
+      HypreFree(row_indexes);
+    }
+  }
+};
+
+RawBoundaryEliminationSystem makeSingleBoundarySystem() {
+  return {{2, 3, 2},
+          {0, 1, 2},
+          {0, 1, 0, 1, 2, 1, 2},
+          {2.0, 3.0, 5.0, 7.0, 11.0, 13.0, 17.0},
+          {0}};
+}
+
+RawBoundaryEliminationSystem makeMultiCoupledBoundarySystem() {
+  return {{2, 3, 3, 2},
+          {0, 1, 2, 3},
+          {0, 1, 0, 1, 3, 0, 1, 2, 2, 3},
+          {2.0, 3.0, 5.0, 7.0, 11.0, 13.0, 17.0, 19.0, 23.0, 29.0},
+          {0}};
+}
+
+RawBoundaryEliminationSystem makeBackwardCoupledRowSystem() {
+  return {{2, 3, 3, 2},
+          {0, 1, 2, 3},
+          {0, 2, 1, 2, 3, 0, 2, 3, 1, 3},
+          {2.0, 3.0, 31.0, 37.0, 29.0, 5.0, 7.0, 11.0, 23.0, 19.0},
+          {0, 3}};
+}
+
+} // namespace
+
+TEST(BoundaryEliminationTest, TransformsSingleBoundarySystem) {
+  auto system = makeSingleBoundarySystem();
+
+  EXPECT_EQ(system.ncols[0], 1);
+  EXPECT_EQ(system.values[0], -1.0);
+  EXPECT_EQ(system.values[1], 0.0);
+  EXPECT_EQ(system.values[2], 0.0);
+  EXPECT_DOUBLE_EQ(system.values[3], -0.5);
+  EXPECT_EQ(system.values[4], 11.0);
+  EXPECT_EQ(system.values[5], 13.0);
+  EXPECT_EQ(system.values[6], 17.0);
+
+  ASSERT_NE(system.elimination, nullptr);
+  EXPECT_EQ(system.elimination->binum_array[0], 0);
+  EXPECT_EQ(system.elimination->bjnum_array[0], 1);
+  EXPECT_EQ(system.elimination->bdep_array[0], -1);
+  EXPECT_EQ(system.elimination->bii_array[0], 2.0);
+  EXPECT_EQ(system.elimination->bij_array[0], 3.0);
+  EXPECT_EQ(system.elimination->aoffset_array[0], 0);
+  EXPECT_EQ(system.elimination->aoffset_array[1], 1);
+  EXPECT_EQ(system.elimination->aknum_array[0], 1);
+  EXPECT_EQ(system.elimination->aki_array[0], 5.0);
+}
+
+TEST(BoundaryEliminationTest, ReducesRightHandSideForSingleBoundarySystem) {
+  auto system = makeSingleBoundarySystem();
+  std::array<HYPRE_Complex, 3> rhs{{19.0, 23.0, 29.0}};
+
+  auto brhs = system.elimination->reduceRightHandSideInPlace(rhs.data());
+
+  ASSERT_NE(brhs, nullptr);
+  EXPECT_EQ(brhs->data[0], 19.0);
+  EXPECT_EQ(rhs[0], 19.0);
+  EXPECT_DOUBLE_EQ(rhs[1], -24.5);
+  EXPECT_EQ(rhs[2], 29.0);
+}
+
+TEST(BoundaryEliminationTest, ExpandsSolutionForSingleBoundarySystem) {
+  auto system = makeSingleBoundarySystem();
+  auto brhs = std::make_shared<bout::HypreComplexArray>(1);
+  brhs->data[0] = 19.0;
+  std::array<HYPRE_Complex, 3> solution{{0.0, 2.0, 3.0}};
+
+  system.elimination->expandSolutionInPlace(brhs, solution.data());
+
+  EXPECT_DOUBLE_EQ(solution[0], 6.5);
+  EXPECT_EQ(solution[1], 2.0);
+  EXPECT_EQ(solution[2], 3.0);
+}
+
+TEST(BoundaryEliminationTest, ReconstructsMatvecForSingleBoundarySystem) {
+  auto system = makeSingleBoundarySystem();
+  std::array<HYPRE_Complex, 3> x{{1.0, 2.0, 3.0}};
+  std::array<HYPRE_Complex, 3> reduced_result{{-1.0, 32.0, 77.0}};
+  auto boundary_values = system.elimination->evaluateBoundaryEquations(x.data());
+
+  system.elimination->expandMatvecResultInPlace(boundary_values, boundary_values,
+                                                reduced_result.data());
+
+  EXPECT_EQ(reduced_result[0], 8.0);
+  EXPECT_EQ(reduced_result[1], 52.0);
+  EXPECT_EQ(reduced_result[2], 77.0);
+}
+
+TEST(BoundaryEliminationTest, EliminatesBoundaryCouplingsFromAllAffectedRows) {
+  auto system = makeMultiCoupledBoundarySystem();
+
+  EXPECT_EQ(system.ncols[0], 1);
+  EXPECT_EQ(system.values[0], -1.0);
+  EXPECT_EQ(system.values[1], 0.0);
+  EXPECT_EQ(system.values[2], 0.0);
+  EXPECT_DOUBLE_EQ(system.values[3], -0.5);
+  EXPECT_EQ(system.values[4], 11.0);
+  EXPECT_EQ(system.values[5], 0.0);
+  EXPECT_DOUBLE_EQ(system.values[6], -2.5);
+  EXPECT_EQ(system.values[7], 19.0);
+  EXPECT_EQ(system.values[8], 23.0);
+  EXPECT_EQ(system.values[9], 29.0);
+
+  EXPECT_EQ(system.elimination->na, 2);
+  EXPECT_EQ(system.elimination->aoffset_array[0], 0);
+  EXPECT_EQ(system.elimination->aoffset_array[1], 2);
+  EXPECT_EQ(system.elimination->aknum_array[0], 1);
+  EXPECT_EQ(system.elimination->aknum_array[1], 2);
+  EXPECT_EQ(system.elimination->aki_array[0], 5.0);
+  EXPECT_EQ(system.elimination->aki_array[1], 13.0);
+}
+
+TEST(BoundaryEliminationTest, ReducesRightHandSideForMultipleAffectedRows) {
+  auto system = makeMultiCoupledBoundarySystem();
+  std::array<HYPRE_Complex, 4> rhs{{19.0, 23.0, 31.0, 37.0}};
+
+  auto brhs = system.elimination->reduceRightHandSideInPlace(rhs.data());
+
+  ASSERT_NE(brhs, nullptr);
+  EXPECT_EQ(brhs->data[0], 19.0);
+  EXPECT_EQ(rhs[0], 19.0);
+  EXPECT_DOUBLE_EQ(rhs[1], -24.5);
+  EXPECT_DOUBLE_EQ(rhs[2], -92.5);
+  EXPECT_EQ(rhs[3], 37.0);
+}
+
+TEST(BoundaryEliminationTest, ReconstructsMatvecForMultipleAffectedRows) {
+  auto system = makeMultiCoupledBoundarySystem();
+  std::array<HYPRE_Complex, 4> x{{1.0, 2.0, 3.0, 4.0}};
+  std::array<HYPRE_Complex, 4> reduced_result{{-1.0, 32.0, 52.0, 185.0}};
+  auto boundary_values = system.elimination->evaluateBoundaryEquations(x.data());
+
+  system.elimination->expandMatvecResultInPlace(boundary_values, boundary_values,
+                                                reduced_result.data());
+
+  EXPECT_EQ(reduced_result[0], 8.0);
+  EXPECT_EQ(reduced_result[1], 52.0);
+  EXPECT_EQ(reduced_result[2], 104.0);
+  EXPECT_EQ(reduced_result[3], 185.0);
+}
+
+TEST(BoundaryEliminationTest, HandlesBackwardCoupledRows) {
+  auto system = makeBackwardCoupledRowSystem();
+
+  EXPECT_EQ(system.ncols[0], 1);
+  EXPECT_EQ(system.ncols[3], 1);
+  EXPECT_EQ(system.values[0], -1.0);
+  EXPECT_EQ(system.values[1], 0.0);
+  EXPECT_DOUBLE_EQ(system.values[2], 31.0 - (29.0 * 23.0 / 19.0));
+  EXPECT_EQ(system.values[4], 0.0);
+  EXPECT_EQ(system.values[5], 0.0);
+  EXPECT_DOUBLE_EQ(system.values[6], -0.5);
+  EXPECT_EQ(system.values[8], 0.0);
+  EXPECT_EQ(system.values[9], -1.0);
+
+  EXPECT_EQ(system.elimination->aoffset_array[0], 0);
+  EXPECT_EQ(system.elimination->aoffset_array[1], 1);
+  EXPECT_EQ(system.elimination->aoffset_array[2], 2);
+  EXPECT_EQ(system.elimination->aknum_array[0], 2);
+  EXPECT_EQ(system.elimination->aknum_array[1], 1);
+}
+
+TEST(BoundaryEliminationTest, DifferentEliminationObjectsProduceDifferentReducedSystems) {
+  auto system_a = makeSingleBoundarySystem();
+  RawBoundaryEliminationSystem system_b{{2, 3, 2},
+                                        {0, 1, 2},
+                                        {0, 1, 0, 1, 2, 1, 2},
+                                        {4.0, -1.0, 6.0, 9.0, 10.0, 13.0, 17.0},
+                                        {0}};
+  std::array<HYPRE_Complex, 3> rhs_a{{19.0, 23.0, 29.0}};
+  std::array<HYPRE_Complex, 3> rhs_b = rhs_a;
+
+  auto brhs_a = system_a.elimination->reduceRightHandSideInPlace(rhs_a.data());
+  auto brhs_b = system_b.elimination->reduceRightHandSideInPlace(rhs_b.data());
+
+  EXPECT_EQ(brhs_a->data[0], 19.0);
+  EXPECT_EQ(brhs_b->data[0], 19.0);
+  EXPECT_DOUBLE_EQ(rhs_a[1], -24.5);
+  EXPECT_DOUBLE_EQ(rhs_b[1], -5.5);
+  EXPECT_NE(rhs_a[1], rhs_b[1]);
+}
+
 #endif // BOUT_HAS_HYPRE
diff --git a/tests/unit/invert/laplace/test_laplace_hypre3d.cxx b/tests/unit/invert/laplace/test_laplace_hypre3d.cxx
index e0febdbb22..2838e09237 100644
--- a/tests/unit/invert/laplace/test_laplace_hypre3d.cxx
+++ b/tests/unit/invert/laplace/test_laplace_hypre3d.cxx
@@ -10,14 +10,18 @@
 #include "bout/invert_laplace.hxx"
 #include "gtest/gtest.h"
 
+#include "bout/bout_types.hxx"
+#include "bout/coordinates.hxx"
 #include "bout/derivs.hxx"
 #include "bout/difops.hxx"
 #include "bout/field2d.hxx"
 #include "bout/field3d.hxx"
+#include "bout/globals.hxx"
 #include "bout/griddata.hxx"
 #include "bout/hypre_interface.hxx"
 #include "bout/mesh.hxx"
 #include "bout/options.hxx"
+#include "bout/region.hxx"
 #include "bout/vecops.hxx"
 
 #include "fake_mesh_fixture.hxx"
@@ -38,7 +42,7 @@ class ForwardOperator {
     upper_y_neumann = yup_neumann;
   }
 
-  const Field3D operator()(Field3D& f) {
+  Field3D operator()(Field3D& f) {
     Field3D result = d * Laplace_perp(f, CELL_DEFAULT, "free", "RGN_NOY")
                      + (Grad(f) * Grad(c2) - DDY(c2) * DDY(f) / coords->g_22()) / c1
                      + a * f + ex * DDX(f) + ez * DDZ(f);
@@ -53,7 +57,7 @@ class ForwardOperator {
   bool inner_x_neumann, outer_x_neumann, // If false then use Dirichlet conditions
       lower_y_neumann, upper_y_neumann;
 
-  void applyBoundaries(Field3D& newF, Field3D& f) {
+  void applyBoundaries(Field3D& newF, Field3D& f) const {
     BOUT_FOR(i, f.getMesh()->getRegion3D("RGN_INNER_X")) {
       if (inner_x_neumann) {
         newF[i] = (f[i.xp()] - f[i]) / coords->dx()[i] / sqrt(coords->g_11()[i]);
@@ -103,14 +107,14 @@ class LaplaceHypre3dTest
     coef3.allocate();
 
     BOUT_FOR(i, mesh->getRegion2D("RGN_ALL")) {
-      BoutReal x = i.x() / (BoutReal)nx - 0.5;
-      BoutReal y = i.y() / (BoutReal)ny - 0.5;
+      const BoutReal x = i.x() / (BoutReal)nx - 0.5;
+      const BoutReal y = i.y() / (BoutReal)ny - 0.5;
       coef2[i] = x + y;
     }
     BOUT_FOR(i, mesh->getRegion3D("RGN_ALL")) {
-      BoutReal x = i.x() / (BoutReal)nx - 0.5;
-      BoutReal y = i.y() / (BoutReal)ny - 0.5;
-      BoutReal z = i.z() / (BoutReal)nz - 0.5;
+      const BoutReal x = i.x() / (BoutReal)nx - 0.5;
+      const BoutReal y = i.y() / (BoutReal)ny - 0.5;
+      const BoutReal z = i.z() / (BoutReal)nz - 0.5;
       f3[i] = 1e3 * exp(-0.5 * sqrt(x * x + y * y + z * z) / sigmasq);
       coef3[i] = x + y + sin(2 * 3.14159265358979323846 * z);
     }
@@ -346,6 +350,17 @@ TEST_P(LaplaceHypre3dTest, TestSolve3DGuess) {
   BOUT_FOR(i, mesh->getRegion3D("RGN_ALL")) { EXPECT_NEAR(expected[i], actual[i], tol); }
 }
 
+TEST_P(LaplaceHypre3dTest, TestSolve3DRepeated) {
+  const Field3D rhs = forward(f3);
+  const Field3D first = solver.solve(rhs);
+  const Field3D second = solver.solve(rhs);
+
+  BOUT_FOR(i, mesh->getRegion3D("RGN_ALL")) {
+    EXPECT_NEAR(f3[i], first[i], tol);
+    EXPECT_NEAR(f3[i], second[i], tol);
+  }
+}
+
 TEST_P(LaplaceHypre3dTest, TestSolvePerp) {
   FieldPerp f(1.0);
   EXPECT_THROW(solver.solve(f), BoutException);