diff --git a/CHANGELOG.md b/CHANGELOG.md
index 13430ffd2d2de..5eefff1f09379 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1315]](https://github.com/parthenon-hpc-lab/parthenon/pull/1315) Add user specifiable BCs to solvers
 - [[PR 1192]](https://github.com/parthenon-hpc-lab/parthenon/pull/1192) Coalesced buffer communication
 - [[PR 1314]](https://github.com/parthenon-hpc-lab/parthenon/pull/1314) Add option of user specified BCs in AddBoundaryExchangeTasks
 - [[PR 1244]](https://github.com/parthenon-hpc-lab/parthenon/pull/1244) Add TaskCollection timeout capability
diff --git a/src/solvers/bicgstab_solver.hpp b/src/solvers/bicgstab_solver.hpp
index 946c95457e7a8..50ca039c8396c 100644
--- a/src/solvers/bicgstab_solver.hpp
+++ b/src/solvers/bicgstab_solver.hpp
@@ -19,6 +19,7 @@
 #include <utility>
 #include <vector>
 
+#include "bvals/comms/bvals_in_one.hpp"
 #include "interface/mesh_data.hpp"
 #include "interface/meshblock_data.hpp"
 #include "interface/state_descriptor.hpp"
@@ -27,6 +28,7 @@
 #include "solvers/solver_base.hpp"
 #include "solvers/solver_utils.hpp"
 #include "tasks/tasks.hpp"
+#include "utils/reductions.hpp"
 #include "utils/type_list.hpp"
 
 namespace parthenon {
@@ -63,42 +65,34 @@ struct BiCGSTABParams {
   }
 };
 
-// The equations class must include a template method
+struct BiCGSTABSolverCounter {
+  static inline std::size_t id{0};
+};
+
+// The equations_t class must include a template method
 //
 //   template <class x_t, class y_t, class TL_t>
 //   TaskID Ax(TL_t &tl, TaskID depends_on, std::shared_ptr<MeshData<Real>> &md)
 //
 // that takes a field associated with x_t and applies
 // the matrix A to it and stores the result in y_t.
-template <class equations, class preconditioner_t = MGSolver<equations>>
-class BiCGSTABSolver : public SolverBase {
-  using FieldTL = typename equations::IndependentVars;
-
-  std::vector<std::string> sol_fields;
-  // Name of user defined container that should contain information required to
-  // calculate the matrix part of the matrix vector product
-  std::string container_base;
-  // User defined container in which the solution will reside, only needs to contain
-  // sol_fields
-  // TODO(LFR): Also allow for an initial guess to come in here
-  std::string container_u;
-  // User defined container containing the rhs vector, only needs to contain sol_fields
-  std::string container_rhs;
+template <class equations_t, class preconditioner_t = MGSolver<equations_t>>
+class BiCGSTABSolver : public SolverBase, BiCGSTABSolverCounter {
+  using FieldTL = typename equations_t::IndependentVars;
+
   // Internal containers for solver which create deep copies of sol_fields
   std::string container_rhat0, container_v, container_h, container_s;
   std::string container_t, container_r, container_p, container_x, container_diag;
-
-  static inline std::size_t id{0};
+  BValOnMDFunc_t BCFunc;
 
  public:
   BiCGSTABSolver(const std::string &container_base, const std::string &container_u,
                  const std::string &container_rhs, ParameterInput *pin,
-                 const std::string &input_block, equations eq_in = equations())
+                 const std::string &input_block, equations_t eq_in = equations_t())
       : preconditioner(container_base, container_u, container_rhs, pin, input_block,
                        eq_in),
-        container_base(container_base), container_u(container_u),
-        container_rhs(container_rhs), params_(pin, input_block), iter_counter(0),
-        eqs_(eq_in) {
+        SolverBase(container_base, container_u, container_rhs), params_(pin, input_block),
+        iter_counter(0), eqs_(eq_in) {
     FieldTL::IterateTypes(
         [this](auto t) { this->sol_fields.push_back(decltype(t)::name()); });
     std::string solver_id = "bicgstab" + std::to_string(id++);
@@ -111,6 +105,17 @@ class BiCGSTABSolver : public SolverBase {
     container_p = solver_id + "_p";
     container_x = solver_id + "_x";
     container_diag = solver_id + "_diag";
+    if constexpr (has_SetBoundary<equations_t>::value) {
+      BCFunc = equations_t::SetBoundary;
+    } else {
+      BCFunc = ApplyBoundaryConditionsOnCoarseOrFineMD;
+    }
+  }
+
+  TaskID Ax(TaskList &tl, TaskID dependence, std::shared_ptr<MeshData<Real>> &md_mat,
+            std::shared_ptr<MeshData<Real>> &md_in,
+            std::shared_ptr<MeshData<Real>> &md_out) {
+    return eqs_.Ax(tl, dependence, md_mat, md_in, md_out);
   }
 
   TaskID AddSetupTasks(TaskList &tl, TaskID dependence, int partition, Mesh *pmesh) {
@@ -120,7 +125,7 @@ class BiCGSTABSolver : public SolverBase {
       auto partitions = pmesh->GetDefaultBlockPartitions();
       auto &md = pmesh->mesh_data.Add(container_base, partitions[partition]);
       auto &md_diag = pmesh->mesh_data.Add(container_diag, md, sol_fields);
-      return tl.AddTask(dependence, &equations::SetDiagonal, &eqs_, md, md_diag);
+      return tl.AddTask(dependence, &equations_t::SetDiagonal, &eqs_, md, md_diag);
     } else {
       return dependence;
     }
@@ -217,8 +222,8 @@ class BiCGSTABSolver : public SolverBase {
     }
 
     // 2. v <- A u
-    auto comm =
-        AddBoundaryExchangeTasks<BoundaryType::any>(precon1, itl, md_u, multilevel);
+    auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(precon1, itl, md_u,
+                                                            multilevel, BCFunc);
     auto get_v = eqs_.Ax(itl, comm, md_base, md_u, md_v);
 
     // 3. rhat0v <- (rhat0, v)
@@ -271,8 +276,8 @@ class BiCGSTABSolver : public SolverBase {
     }
 
     // 7. t <- A u
-    auto pre_t_comm =
-        AddBoundaryExchangeTasks<BoundaryType::any>(precon2, itl, md_u, multilevel);
+    auto pre_t_comm = AddBoundaryExchangeTasks<BoundaryType::any>(precon2, itl, md_u,
+                                                                  multilevel, BCFunc);
     auto get_t = eqs_.Ax(itl, pre_t_comm, md_base, md_u, md_t);
 
     // 8. omega <- (t,s) / (t,t)
@@ -365,7 +370,7 @@ class BiCGSTABSolver : public SolverBase {
   int iter_counter;
   AllReduce<Real> rtr, pAp, rhat0v, rhat0r, ts, tt, residual, rhs2;
   Real rhat0r_old;
-  equations eqs_;
+  equations_t eqs_;
   std::string container_;
 };
 
diff --git a/src/solvers/cg_solver.hpp b/src/solvers/cg_solver.hpp
index 898b16c885b70..9cde01ff10460 100644
--- a/src/solvers/cg_solver.hpp
+++ b/src/solvers/cg_solver.hpp
@@ -17,9 +17,11 @@
 #include <limits>
 #include <memory>
 #include <string>
+#include <type_traits>
 #include <utility>
 #include <vector>
 
+#include "bvals/comms/bvals_in_one.hpp"
 #include "interface/mesh_data.hpp"
 #include "interface/meshblock_data.hpp"
 #include "interface/state_descriptor.hpp"
@@ -28,6 +30,7 @@
 #include "solvers/solver_base.hpp"
 #include "solvers/solver_utils.hpp"
 #include "tasks/tasks.hpp"
+#include "utils/reductions.hpp"
 #include "utils/type_list.hpp"
 
 namespace parthenon {
@@ -54,41 +57,33 @@ struct CGParams {
   }
 };
 
-// The equations class must include a template method
+struct CGSolverCounter {
+  static inline std::size_t id{0};
+};
+
+// The equations_t class must include a template method
 //
 //   template <class x_t, class y_t, class TL_t>
 //   TaskID Ax(TL_t &tl, TaskID depends_on, std::shared_ptr<MeshData<Real>> &md)
 //
 // that takes a field associated with x_t and applies
 // the matrix A to it and stores the result in y_t.
-template <class equations, class preconditioner_t = MGSolver<equations>>
-class CGSolver : public SolverBase {
-  using FieldTL = typename equations::IndependentVars;
-
-  std::vector<std::string> sol_fields;
-  // Name of user defined container that should contain information required to
-  // calculate the matrix part of the matrix vector product
-  std::string container_base;
-  // User defined container in which the solution will reside, only needs to contain
-  // sol_fields
-  // TODO(LFR): Also allow for an initial guess to come in here
-  std::string container_u;
-  // User defined container containing the rhs vector, only needs to contain sol_fields
-  std::string container_rhs;
+template <class equations_t, class preconditioner_t = MGSolver<equations_t>>
+class CGSolver : public SolverBase, CGSolverCounter {
+  using FieldTL = typename equations_t::IndependentVars;
+
   // Internal containers for solver which create deep copies of sol_fields
   std::string container_x, container_r, container_v, container_p;
-
-  static inline std::size_t id{0};
+  BValOnMDFunc_t BCFunc;
 
  public:
   CGSolver(const std::string &container_base, const std::string &container_u,
            const std::string &container_rhs, ParameterInput *pin,
-           const std::string &input_block, const equations &eq_in = equations())
-      : preconditioner(container_base, container_u, container_rhs, pin, input_block,
+           const std::string &input_block, const equations_t &eq_in = equations_t())
+      : SolverBase(container_base, container_u, container_rhs),
+        preconditioner(container_base, container_u, container_rhs, pin, input_block,
                        eq_in),
-        container_base(container_base), container_u(container_u),
-        container_rhs(container_rhs), params_(pin, input_block), iter_counter(0),
-        eqs_(eq_in) {
+        params_(pin, input_block), iter_counter(0), eqs_(eq_in) {
     FieldTL::IterateTypes(
         [this](auto t) { this->sol_fields.push_back(decltype(t)::name()); });
     std::string solver_id = "cg" + std::to_string(id++);
@@ -96,6 +91,17 @@ class CGSolver : public SolverBase {
     container_r = solver_id + "_r";
     container_v = solver_id + "_v";
     container_p = solver_id + "_p";
+    if constexpr (has_SetBoundary<equations_t>::value) {
+      BCFunc = equations_t::SetBoundary;
+    } else {
+      BCFunc = ApplyBoundaryConditionsOnCoarseOrFineMD;
+    }
+  }
+
+  TaskID Ax(TaskList &tl, TaskID dependence, std::shared_ptr<MeshData<Real>> &md_mat,
+            std::shared_ptr<MeshData<Real>> &md_in,
+            std::shared_ptr<MeshData<Real>> &md_out) {
+    return eqs_.Ax(tl, dependence, md_mat, md_in, md_out);
   }
 
   TaskID AddSetupTasks(TaskList &tl, TaskID dependence, int partition, Mesh *pmesh) {
@@ -154,7 +160,7 @@ class CGSolver : public SolverBase {
                            : *res_tol;
             printf("# [0] v-cycle\n# [1] rms-residual (tol = %e) \n# [2] rms-error\n",
                    tol);
-            printf("0 %e\n", std::sqrt(solver->rhs2.val / pm->GetTotalCells()));
+            printf("\t0 %e\n", std::sqrt(solver->rhs2.val / pm->GetTotalCells()));
             return TaskStatus::complete;
           },
           this, params_.residual_tolerance, params_.relative_residual, pmesh);
@@ -199,8 +205,8 @@ class CGSolver : public SolverBase {
         this, md_u, md_p);
 
     // 4. v <- A p
-    auto comm =
-        AddBoundaryExchangeTasks<BoundaryType::any>(correct_p, itl, md_p, multilevel);
+    auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(correct_p, itl, md_p,
+                                                            multilevel, BCFunc);
     auto get_v = eqs_.Ax(itl, comm, md_base, md_p, md_v);
 
     // 5. alpha <- r dot u / p dot v (calculate denominator)
@@ -234,7 +240,7 @@ class CGSolver : public SolverBase {
         [&](CGSolver *solver, Mesh *pmesh) {
           Real rms_res = std::sqrt(solver->residual.val / pmesh->GetTotalCells());
           if (Globals::my_rank == 0 && solver->params_.print_per_step)
-            printf("\t%i %e\n", solver->iter_counter, rms_res);
+            printf("\t%i %e\n", solver->iter_counter + 1, rms_res);
           return TaskStatus::complete;
         },
         this, pmesh);
@@ -273,7 +279,7 @@ class CGSolver : public SolverBase {
   int iter_counter;
   AllReduce<Real> ru, pAp, residual, rhs2;
   Real ru_old;
-  equations eqs_;
+  equations_t eqs_;
 };
 
 } // namespace solvers
diff --git a/src/solvers/mg_solver.hpp b/src/solvers/mg_solver.hpp
index 6662b6cdffb5b..79772f1bb5852 100644
--- a/src/solvers/mg_solver.hpp
+++ b/src/solvers/mg_solver.hpp
@@ -21,6 +21,7 @@
 #include <utility>
 #include <vector>
 
+#include "bvals/comms/bvals_in_one.hpp"
 #include "interface/mesh_data.hpp"
 #include "interface/meshblock_data.hpp"
 #include "interface/state_descriptor.hpp"
@@ -29,6 +30,7 @@
 #include "solvers/solver_base.hpp"
 #include "solvers/solver_utils.hpp"
 #include "tasks/tasks.hpp"
+#include "utils/reductions.hpp"
 #include "utils/robust.hpp"
 #include "utils/type_list.hpp"
 
@@ -60,6 +62,9 @@ struct MGParams {
   }
 };
 
+struct MGSolverCounter {
+  static inline std::size_t id{0};
+};
 // The equations_t class must include a template method
 //
 //   template <class x_t, class y_t, class TL_t>
@@ -75,25 +80,13 @@ struct MGParams {
 // That stores the (possibly approximate) diagonal of matrix A in the field
 // associated with the type diag_t. This is used for Jacobi iteration.
 template <class equations_t, class prolongator_t = ProlongationBlockInteriorDefault>
-class MGSolver : public SolverBase {
-  static inline std::size_t id{0};
-
+class MGSolver : public SolverBase, MGSolverCounter {
  public:
   using FieldTL = typename equations_t::IndependentVars;
 
-  std::vector<std::string> sol_fields;
-
-  // Name of user defined container that should contain information required to
-  // calculate the matrix part of the matrix vector product
-  std::string container_base;
-  // User defined container in which the solution will reside, only needs to contain
-  // sol_fields
-  // TODO(LFR): Also allow for an initial guess to come in here
-  std::string container_u;
-  // User defined container containing the rhs vector, only needs to contain sol_fields
-  std::string container_rhs;
   // Internal containers for solver which create deep copies of sol_fields
   std::string container_res_err, container_temp, container_u0, container_diag;
+  BValOnMDFunc_t BCFunc;
 
   MGSolver(const std::string &container_base, const std::string &container_u,
            const std::string &container_rhs, ParameterInput *pin,
@@ -104,9 +97,8 @@ class MGSolver : public SolverBase {
   MGSolver(const std::string &container_base, const std::string &container_u,
            const std::string &container_rhs, MGParams params_in,
            equations_t eq_in = equations_t(), prolongator_t prol_in = prolongator_t())
-      : container_base(container_base), container_u(container_u),
-        container_rhs(container_rhs), params_(params_in), iter_counter(0), eqs_(eq_in),
-        prolongator_(prol_in) {
+      : SolverBase(container_base, container_u, container_rhs), params_(params_in),
+        iter_counter(0), eqs_(eq_in), prolongator_(prol_in) {
     FieldTL::IterateTypes(
         [this](auto t) { this->sol_fields.push_back(decltype(t)::name()); });
     std::string solver_id = "mg" + std::to_string(id++);
@@ -114,6 +106,17 @@ class MGSolver : public SolverBase {
     container_temp = solver_id + "_temp";
     container_u0 = solver_id + "_u0";
     container_diag = solver_id + "_diag";
+    if constexpr (has_SetBoundary<equations_t>::value) {
+      BCFunc = equations_t::SetBoundary;
+    } else {
+      BCFunc = ApplyBoundaryConditionsOnCoarseOrFineMD;
+    }
+  }
+
+  TaskID Ax(TaskList &tl, TaskID dependence, std::shared_ptr<MeshData<Real>> &md_mat,
+            std::shared_ptr<MeshData<Real>> &md_in,
+            std::shared_ptr<MeshData<Real>> &md_out) {
+    return eqs_.Ax(tl, dependence, md_mat, md_in, md_out);
   }
 
   TaskID AddTasks(TaskList &tl, TaskID dependence, const int partition, Mesh *pmesh) {
@@ -140,7 +143,8 @@ class MGSolver : public SolverBase {
     auto &md_res_err = pmesh->mesh_data.Add(container_res_err, md, sol_fields);
     auto &md_rhs = pmesh->mesh_data.Add(container_rhs, md, sol_fields);
     auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(mg_finest, itl, md_u,
-                                                            pmesh->multilevel);
+                                                            pmesh->multilevel, BCFunc);
+
     auto calc_pointwise_res = eqs_.Ax(itl, comm, md, md_u, md_res_err);
     calc_pointwise_res =
         itl.AddTask(calc_pointwise_res, TF(AddFieldsAndStoreInteriorSelect<FieldTL>),
@@ -288,8 +292,8 @@ class MGSolver : public SolverBase {
     auto &md_rhs = pmesh->mesh_data.Add(container_rhs, partitions[partition], sol_fields);
     auto &md_diag = pmesh->mesh_data.Add(container_diag, md_base, sol_fields);
 
-    auto comm =
-        AddBoundaryExchangeTasks<comm_boundary>(depends_on, tl, md_in, multilevel);
+    auto comm = AddBoundaryExchangeTasks<comm_boundary>(depends_on, tl, md_in, multilevel,
+                                                        BCFunc);
     auto mat_mult = eqs_.Ax(tl, comm, md_base, md_in, md_out);
     return tl.AddTask(mat_mult, TF(&MGSolver::Jacobi), this, md_rhs, md_out, md_diag,
                       md_in, md_out, omega);
@@ -436,7 +440,7 @@ class MGSolver : public SolverBase {
         // calling Ax. That being said, at least in one case commenting this line out
         // didn't seem to impact the solution.
         set_from_finer = AddBoundaryExchangeTasks<BoundaryType::gmg_same>(
-            set_from_finer, tl, md_u, multilevel);
+            set_from_finer, tl, md_u, multilevel, BCFunc);
         set_from_finer =
             tl.AddTask(set_from_finer, BTF(CopyData<FieldTL, true>), md_u, md_u0);
         // This should set the rhs only in blocks that correspond to interior nodes, the
@@ -460,7 +464,7 @@ class MGSolver : public SolverBase {
     if (level > min_level) {
       // 3. Communicate same level boundaries so that u is up to date everywhere
       auto comm_u = AddBoundaryExchangeTasks<BoundaryType::gmg_same>(pre_smooth, tl, md_u,
-                                                                     multilevel);
+                                                                     multilevel, BCFunc);
 
       // 4. Caclulate residual and store in communication field
       auto residual = eqs_.Ax(tl, comm_u, md, md_u, md_temp);
@@ -515,7 +519,7 @@ class MGSolver : public SolverBase {
       // This is required to make sure boundaries of res_err are up to date before
       // prolongation
       auto boundary = AddBoundaryExchangeTasks<BoundaryType::gmg_same>(
-          copy_over, tl, md_res_err, multilevel);
+          copy_over, tl, md_res_err, multilevel, BCFunc);
       last_task = tl.AddTask(
           boundary, BTF(SendBoundBufs<BoundaryType::gmg_prolongate_send>), md_res_err);
     }
diff --git a/src/solvers/solver_base.hpp b/src/solvers/solver_base.hpp
index c5fc1653a2e54..c3bef1e86b59e 100644
--- a/src/solvers/solver_base.hpp
+++ b/src/solvers/solver_base.hpp
@@ -22,15 +22,28 @@
 #include <vector>
 
 #include "interface/mesh_data.hpp"
-#include "interface/meshblock_data.hpp"
-#include "tasks/tasks.hpp"
-
 namespace parthenon {
 
 namespace solvers {
 
+// Used for checking if a given equations class has a SetBoundary function
+template <typename T, typename = void>
+struct has_SetBoundary : std::false_type {};
+
+template <typename T>
+struct has_SetBoundary<
+    T, std::void_t<decltype(std::declval<T>().SetBoundary(
+           std::declval<std::shared_ptr<MeshData<Real>> &>(), std::declval<bool>()))>>
+    : std::true_type {};
+
+// Solver base class
 class SolverBase {
  public:
+  SolverBase(const std::string &container_base, const std::string &container_u,
+             const std::string &container_rhs)
+      : container_base(container_base), container_u(container_u),
+        container_rhs(container_rhs) {}
+
   virtual ~SolverBase() {}
 
   virtual TaskID AddSetupTasks(TaskList &tl, TaskID dependence, int partition,
@@ -38,10 +51,34 @@ class SolverBase {
   virtual TaskID AddTasks(TaskList &tl, TaskID dependence, int partition,
                           Mesh *pmesh) = 0;
 
+  // Provide access to the underlying matrix operator for convenience
+  virtual TaskID Ax(TaskList &tl, TaskID dependence,
+                    std::shared_ptr<MeshData<Real>> &md_mat,
+                    std::shared_ptr<MeshData<Real>> &md_in,
+                    std::shared_ptr<MeshData<Real>> &md_out) = 0;
+
   Real GetFinalResidual() const { return final_residual; }
   int GetFinalIterations() const { return final_iteration; }
 
+  const std::string &GetBaseContainerLabel() const { return container_base; }
+  const std::string &GetRHSContainerLabel() const { return container_rhs; }
+  const std::string &GetSolutionContainerLabel() const { return container_u; }
+
+  const std::vector<std::string> &GetFieldLabels() const { return sol_fields; }
+
  protected:
+  // Labels of all fields included in the vector
+  std::vector<std::string> sol_fields;
+  // Name of user defined container that should contain information required to
+  // calculate the matrix part of the matrix vector product
+  std::string container_base;
+  // User defined container in which the solution will reside, only needs to contain
+  // sol_fields
+  // TODO(LFR): Also allow for an initial guess to come in here
+  std::string container_u;
+  // User defined container containing the rhs vector, only needs to contain sol_fields
+  std::string container_rhs;
+
   Real final_residual;
   int final_iteration;
 };
diff --git a/src/solvers/solver_utils.hpp b/src/solvers/solver_utils.hpp
index d24f8e3c5026d..d5a7dcceeaaa7 100644
--- a/src/solvers/solver_utils.hpp
+++ b/src/solvers/solver_utils.hpp
@@ -21,6 +21,7 @@
 #include <vector>
 
 #include "kokkos_abstraction.hpp"
+#include "utils/reductions.hpp"
 
 #define PARTHENON_INTERNALSOLVERVARIABLE(base, varname)                                  \
   struct varname : public parthenon::variable_names::base_t<false> {                     \
diff --git a/src/solvers/tridiag_solver.hpp b/src/solvers/tridiag_solver.hpp
index 1c675e171c723..8e168845febde 100644
--- a/src/solvers/tridiag_solver.hpp
+++ b/src/solvers/tridiag_solver.hpp
@@ -34,6 +34,10 @@ namespace parthenon {
 
 namespace solvers {
 
+struct TridiagSolverCounter {
+  static inline std::size_t id{0};
+};
+
 // The equations class must include a template method
 //
 //   template <class x_t, class y_t, class TL_t>
@@ -42,32 +46,20 @@ namespace solvers {
 // that takes a field associated with x_t and applies
 // the matrix A to it and stores the result in y_t.
 template <class equations>
-class TridiagSolver : public SolverBase {
+class TridiagSolver : public SolverBase, TridiagSolverCounter {
   using FieldTL = typename equations::IndependentVars;
 
-  std::vector<std::string> sol_fields;
-  // Name of user defined container that should contain information required to
-  // calculate the matrix part of the matrix vector product
-  std::string container_base;
-  // User defined container in which the solution will reside, only needs to contain
-  // sol_fields
-  // TODO(LFR): Also allow for an initial guess to come in here
-  std::string container_u;
-  // User defined container containing the rhs vector, only needs to contain sol_fields
-  std::string container_rhs;
   // Internal containers for solver which create deep copies of sol_fields
   std::string container_100, container_010, container_001, container_100_out,
       container_010_out, container_001_out, container_Aup, container_Adi, container_Alo,
       container_r;
 
-  static inline std::size_t id{0};
-
  public:
   TridiagSolver(const std::string &container_base, const std::string &container_u,
                 const std::string &container_rhs, ParameterInput *pin,
                 const std::string &input_block, const equations &eq_in = equations())
-      : container_base(container_base), container_u(container_u),
-        container_rhs(container_rhs), iter_counter(0), eqs_(eq_in),
+      : SolverBase(container_base, container_u, container_rhs), iter_counter(0),
+        eqs_(eq_in),
         print_solution_(pin->GetOrAddBoolean(input_block, "print_solution", false)) {
     FieldTL::IterateTypes(
         [this](auto t) { this->sol_fields.push_back(decltype(t)::name()); });
@@ -303,6 +295,12 @@ class TridiagSolver : public SolverBase {
   Real GetSquaredResidualSum() const { return 0.0; }
   int GetCurrentIterations() const { return 1; }
 
+  TaskID Ax(TaskList &tl, TaskID dependence, std::shared_ptr<MeshData<Real>> &md_mat,
+            std::shared_ptr<MeshData<Real>> &md_in,
+            std::shared_ptr<MeshData<Real>> &md_out) {
+    return eqs_.Ax(tl, dependence, md_mat, md_in, md_out);
+  }
+
  protected:
   int iter_counter;
   Real ru_old;