diff --git a/src/executioners/steady_executioner.cpp b/src/executioners/steady_executioner.cpp
index b5994d430..f8ac5d293 100644
--- a/src/executioners/steady_executioner.cpp
+++ b/src/executioners/steady_executioner.cpp
@@ -13,7 +13,7 @@ SteadyExecutioner::Solve() const
 {
   // Advance time step.
   _problem->_preprocessors.Solve();
-  _problem->GetOperator()->Solve(*(_problem->_f));
+  _problem->GetOperator()->Solve();
   _problem->_postprocessors.Solve();
 
   // Output data
diff --git a/src/executioners/transient_executioner.cpp b/src/executioners/transient_executioner.cpp
index ddd7fa067..a947395e6 100644
--- a/src/executioners/transient_executioner.cpp
+++ b/src/executioners/transient_executioner.cpp
@@ -27,7 +27,7 @@ TransientExecutioner::Step(double dt, int it) const
 
   // Advance time step.
   _problem->_preprocessors.Solve(_t);
-  _problem->_ode_solver->Step(*(_problem->_f), _t, dt);
+  _problem->GetOperator()->Step(_t, dt);
   _problem->_postprocessors.Solve(_t);
 
   // Output data
diff --git a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp
index d4bcbdd0e..238ca0c36 100644
--- a/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp
+++ b/src/formulations/ComplexMaxwell/complex_maxwell_formulation.cpp
@@ -218,8 +218,8 @@ ComplexMaxwellOperator::Solve(mfem::Vector & X)
 
   auto * jac_z = jac.As<mfem::ComplexHypreParMatrix>()->GetSystemMatrix();
 
-  _problem._jacobian_solver->SetOperator(*jac_z);
-  _problem._jacobian_solver->Mult(rhs, u);
+  _jacobian_solver->SetOperator(*jac_z);
+  _jacobian_solver->Mult(rhs, u);
   sqlf.RecoverFEMSolution(u, lf, *_u);
 
   _problem._gridfunctions.GetRef(_trial_var_names.at(0)) = _u->real();
diff --git a/src/formulations/Dual/dual_formulation.cpp b/src/formulations/Dual/dual_formulation.cpp
index e85002533..56a21370c 100644
--- a/src/formulations/Dual/dual_formulation.cpp
+++ b/src/formulations/Dual/dual_formulation.cpp
@@ -59,20 +59,16 @@ DualFormulation::DualFormulation(std::string alpha_coef_name,
 }
 
 void
-DualFormulation::ConstructJacobianPreconditioner()
+DualFormulation::ConstructJacobianSolver()
 {
   auto precond =
-      std::make_shared<mfem::HypreAMS>(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0));
+      std::make_unique<mfem::HypreAMS>(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0));
 
   precond->SetSingularProblem();
   precond->SetPrintLevel(-1);
 
-  GetProblem()->_jacobian_preconditioner = precond;
-}
+  GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond));
 
-void
-DualFormulation::ConstructJacobianSolver()
-{
   ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG, {._max_iteration = 1000});
 }
 
@@ -238,12 +234,12 @@ DualOperator::ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vecto
   logger.info("{} ImplicitSolve: {} seconds", typeid(this).name(), sw);
 }
 
-void
-DualOperator::UpdateOffsets()
+int
+DualOperator::GetSolutionVectorSize() const
 {
   // Blocks for solution vector are smaller than the operator size for DualOperator,
   // as curl is stored separately. Block operator only has the HCurl TrueVSize;
-  return UpdateOffsetsWithSize(_trial_variables.size() - 1);
+  return static_cast<int>(GetTrialVariablesSize() - 1);
 }
 
 } // namespace hephaestus
diff --git a/src/formulations/Dual/dual_formulation.hpp b/src/formulations/Dual/dual_formulation.hpp
index 58d55005d..5a09f8923 100644
--- a/src/formulations/Dual/dual_formulation.hpp
+++ b/src/formulations/Dual/dual_formulation.hpp
@@ -16,8 +16,6 @@ class DualFormulation : public TimeDomainEMFormulation
 
   ~DualFormulation() override = default;
 
-  void ConstructJacobianPreconditioner() override;
-
   void ConstructJacobianSolver() override;
 
   void ConstructOperator() override;
@@ -74,7 +72,7 @@ class DualOperator : public TimeDomainEquationSystemProblemOperator
   mfem::ParGridFunction * _dv{nullptr}; // HDiv vector field
 
 protected:
-  void UpdateOffsets() override;
+  int GetSolutionVectorSize() const override;
 
   std::unique_ptr<mfem::ParDiscreteLinearOperator> _curl;
 };
diff --git a/src/formulations/HCurl/hcurl_formulation.cpp b/src/formulations/HCurl/hcurl_formulation.cpp
index 435e932f7..bdb5acee0 100644
--- a/src/formulations/HCurl/hcurl_formulation.cpp
+++ b/src/formulations/HCurl/hcurl_formulation.cpp
@@ -67,20 +67,16 @@ HCurlFormulation::ConstructOperator()
 }
 
 void
-HCurlFormulation::ConstructJacobianPreconditioner()
+HCurlFormulation::ConstructJacobianSolver()
 {
   auto precond =
-      std::make_shared<mfem::HypreAMS>(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0));
+      std::make_unique<mfem::HypreAMS>(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0));
 
   precond->SetSingularProblem();
   precond->SetPrintLevel(-1);
 
-  GetProblem()->_jacobian_preconditioner = precond;
-}
+  GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond));
 
-void
-HCurlFormulation::ConstructJacobianSolver()
-{
   ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG);
 }
 
diff --git a/src/formulations/HCurl/hcurl_formulation.hpp b/src/formulations/HCurl/hcurl_formulation.hpp
index b1c4b8ba4..dc9e491a1 100644
--- a/src/formulations/HCurl/hcurl_formulation.hpp
+++ b/src/formulations/HCurl/hcurl_formulation.hpp
@@ -18,8 +18,6 @@ class HCurlFormulation : public TimeDomainEMFormulation
 
   void ConstructOperator() override;
 
-  void ConstructJacobianPreconditioner() override;
-
   void ConstructJacobianSolver() override;
 
   void RegisterGridFunctions() override;
diff --git a/src/formulations/Magnetostatic/statics_formulation.cpp b/src/formulations/Magnetostatic/statics_formulation.cpp
index 05e0f0d08..25ddd21e5 100644
--- a/src/formulations/Magnetostatic/statics_formulation.cpp
+++ b/src/formulations/Magnetostatic/statics_formulation.cpp
@@ -38,20 +38,16 @@ StaticsFormulation::StaticsFormulation(std::string alpha_coef_name, std::string
 }
 
 void
-StaticsFormulation::ConstructJacobianPreconditioner()
+StaticsFormulation::ConstructJacobianSolver()
 {
-  std::shared_ptr<mfem::HypreAMS> precond{std::make_shared<mfem::HypreAMS>(
-      GetProblem()->_gridfunctions.Get(_h_curl_var_name)->ParFESpace())};
+  auto precond = std::make_unique<mfem::HypreAMS>(
+      GetProblem()->_gridfunctions.Get(_h_curl_var_name)->ParFESpace());
 
   precond->SetSingularProblem();
   precond->SetPrintLevel(-1);
 
-  GetProblem()->_jacobian_preconditioner = precond;
-}
+  GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond));
 
-void
-StaticsFormulation::ConstructJacobianSolver()
-{
   ConstructJacobianSolverWithOptions(SolverType::HYPRE_FGMRES,
                                      {._max_iteration = 100, ._k_dim = 10});
 }
@@ -152,8 +148,8 @@ StaticsOperator::Solve(mfem::Vector & X)
 
   // Define and apply a parallel FGMRES solver for AX=B with the AMS
   // preconditioner from hypre.
-  _problem._jacobian_solver->SetOperator(curl_mu_inv_curl);
-  _problem._jacobian_solver->Mult(rhs_tdofs, sol_tdofs);
+  _jacobian_solver->SetOperator(curl_mu_inv_curl);
+  _jacobian_solver->Mult(rhs_tdofs, sol_tdofs);
   blf.RecoverFEMSolution(sol_tdofs, lf, gf);
 
   logger.info("{} Solve: {} seconds", typeid(this).name(), sw);
diff --git a/src/formulations/Magnetostatic/statics_formulation.hpp b/src/formulations/Magnetostatic/statics_formulation.hpp
index cf45d6f8a..b0a901695 100644
--- a/src/formulations/Magnetostatic/statics_formulation.hpp
+++ b/src/formulations/Magnetostatic/statics_formulation.hpp
@@ -14,8 +14,6 @@ class StaticsFormulation : public SteadyStateEMFormulation
 
   ~StaticsFormulation() override = default;
 
-  void ConstructJacobianPreconditioner() override;
-
   void ConstructJacobianSolver() override;
 
   void ConstructOperator() override;
diff --git a/src/problem_builders/problem_builder_base.cpp b/src/problem_builders/problem_builder_base.cpp
index 49c40900b..d5e4ef8dd 100644
--- a/src/problem_builders/problem_builder_base.cpp
+++ b/src/problem_builders/problem_builder_base.cpp
@@ -3,13 +3,6 @@
 namespace hephaestus
 {
 
-Problem::~Problem()
-{
-  // Ensure that all owned memory is properly freed!
-  _f.reset();
-  _ode_solver.reset();
-}
-
 void
 Problem::Update()
 {
@@ -97,15 +90,15 @@ ProblemBuilder::SetSolverOptions(hephaestus::InputParameters & solver_options)
 }
 
 void
-ProblemBuilder::SetJacobianPreconditioner(std::shared_ptr<mfem::Solver> preconditioner)
+ProblemBuilder::SetJacobianPreconditioner(std::unique_ptr<mfem::Solver> preconditioner)
 {
-  GetProblem()->_jacobian_preconditioner = preconditioner;
+  GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(preconditioner));
 }
 
 void
-ProblemBuilder::SetJacobianSolver(std::shared_ptr<mfem::Solver> jacobian_solver)
+ProblemBuilder::SetJacobianSolver(std::unique_ptr<mfem::Solver> jacobian_solver)
 {
-  GetProblem()->_jacobian_solver = jacobian_solver;
+  GetProblem()->GetOperator()->SetJacobianSolver(std::move(jacobian_solver));
 }
 
 void
@@ -223,26 +216,16 @@ ProblemBuilder::AddSource(std::string source_name, std::shared_ptr<hephaestus::S
 }
 
 void
-ProblemBuilder::ConstructJacobianPreconditioner()
+ProblemBuilder::ConstructJacobianSolver()
 {
-  auto precond = std::make_shared<mfem::HypreBoomerAMG>();
+  auto precond = std::make_unique<mfem::HypreBoomerAMG>();
   precond->SetPrintLevel(GetGlobalPrintLevel());
 
-  GetProblem()->_jacobian_preconditioner = precond;
-}
+  GetProblem()->GetOperator()->SetJacobianPreconditioner(std::move(precond));
 
-void
-ProblemBuilder::ConstructJacobianSolver()
-{
   ConstructJacobianSolverWithOptions(SolverType::HYPRE_GMRES);
 }
 
-void
-ProblemBuilder::ConstructBlockVector()
-{
-  GetProblem()->_f = std::make_unique<mfem::BlockVector>();
-}
-
 void
 ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams default_params)
 {
@@ -258,14 +241,13 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams
       solver_options.GetOptionalParam<int>("PrintLevel", default_params._print_level);
   const auto k_dim = solver_options.GetOptionalParam<unsigned int>("KDim", default_params._k_dim);
 
-  auto preconditioner =
-      std::dynamic_pointer_cast<mfem::HypreSolver>(GetProblem()->_jacobian_preconditioner);
+  auto preconditioner = GetProblem()->GetOperator()->JacobianPreconditioner<mfem::HypreSolver>();
 
   switch (type)
   {
     case SolverType::HYPRE_PCG:
     {
-      auto solver = std::make_shared<mfem::HyprePCG>(GetProblem()->_comm);
+      auto solver = std::make_unique<mfem::HyprePCG>(GetProblem()->_comm);
 
       solver->SetTol(tolerance);
       solver->SetAbsTol(abs_tolerance);
@@ -275,12 +257,12 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams
       if (preconditioner)
         solver->SetPreconditioner(*preconditioner);
 
-      GetProblem()->_jacobian_solver = solver;
+      GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver));
       break;
     }
     case SolverType::HYPRE_GMRES:
     {
-      auto solver = std::make_shared<mfem::HypreGMRES>(GetProblem()->_comm);
+      auto solver = std::make_unique<mfem::HypreGMRES>(GetProblem()->_comm);
 
       solver->SetTol(tolerance);
       solver->SetAbsTol(abs_tolerance);
@@ -291,12 +273,12 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams
       if (preconditioner)
         solver->SetPreconditioner(*preconditioner);
 
-      GetProblem()->_jacobian_solver = solver;
+      GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver));
       break;
     }
     case SolverType::HYPRE_FGMRES:
     {
-      auto solver = std::make_shared<mfem::HypreFGMRES>(GetProblem()->_comm);
+      auto solver = std::make_unique<mfem::HypreFGMRES>(GetProblem()->_comm);
 
       solver->SetTol(tolerance);
       solver->SetMaxIter(max_iter);
@@ -306,25 +288,25 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams
       if (preconditioner)
         solver->SetPreconditioner(*preconditioner);
 
-      GetProblem()->_jacobian_solver = solver;
+      GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver));
       break;
     }
     case SolverType::HYPRE_AMG:
     {
-      auto solver = std::make_shared<mfem::HypreBoomerAMG>();
+      auto solver = std::make_unique<mfem::HypreBoomerAMG>();
 
       solver->SetTol(tolerance);
       solver->SetMaxIter(max_iter);
       solver->SetPrintLevel(print_level);
 
-      GetProblem()->_jacobian_solver = solver;
+      GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver));
       break;
     }
     case SolverType::SUPER_LU:
     {
-      auto solver = std::make_shared<hephaestus::SuperLUSolver>(GetProblem()->_comm);
+      auto solver = std::make_unique<hephaestus::SuperLUSolver>(GetProblem()->_comm);
 
-      GetProblem()->_jacobian_solver = solver;
+      GetProblem()->GetOperator()->SetJacobianSolver(std::move(solver));
       break;
     }
     default:
@@ -338,14 +320,14 @@ ProblemBuilder::ConstructJacobianSolverWithOptions(SolverType type, SolverParams
 void
 ProblemBuilder::ConstructNonlinearSolver()
 {
-  auto nl_solver = std::make_shared<mfem::NewtonSolver>(GetProblem()->_comm);
+  auto nl_solver = std::make_unique<mfem::NewtonSolver>(GetProblem()->_comm);
 
   // Defaults to one iteration, without further nonlinear iterations
   nl_solver->SetRelTol(0.0);
   nl_solver->SetAbsTol(0.0);
   nl_solver->SetMaxIter(1);
 
-  GetProblem()->_nonlinear_solver = nl_solver;
+  GetProblem()->GetOperator()->SetNonlinearSolver(std::move(nl_solver));
 }
 
 void
@@ -390,13 +372,10 @@ ProblemBuilder::FinalizeProblem(bool build_operator)
     ConstructOperator();
   }
 
-  ConstructBlockVector();
-
   InitializeAuxSolvers();
   InitializeSources();
   InitializeOperator();
 
-  ConstructJacobianPreconditioner();
   ConstructJacobianSolver();
   ConstructNonlinearSolver();
 
diff --git a/src/problem_builders/problem_builder_base.hpp b/src/problem_builders/problem_builder_base.hpp
index 55d881def..39d85ac62 100644
--- a/src/problem_builders/problem_builder_base.hpp
+++ b/src/problem_builders/problem_builder_base.hpp
@@ -12,14 +12,14 @@ namespace hephaestus
 {
 
 /// Forwards declaration.
-class ProblemOperatorInterface;
+class ProblemOperatorBase;
 
 /// Base problem class.
 class Problem
 {
 public:
   Problem() = default;
-  virtual ~Problem();
+  virtual ~Problem() = default;
 
   std::shared_ptr<mfem::ParMesh> _pmesh{nullptr};
   hephaestus::BCMap _bc_map;
@@ -30,13 +30,6 @@ class Problem
   hephaestus::Outputs _outputs;
   hephaestus::InputParameters _solver_options;
 
-  std::unique_ptr<mfem::ODESolver> _ode_solver{nullptr};
-  std::unique_ptr<mfem::BlockVector> _f{nullptr};
-
-  std::shared_ptr<mfem::Solver> _jacobian_preconditioner{nullptr};
-  std::shared_ptr<mfem::Solver> _jacobian_solver{nullptr};
-  std::shared_ptr<mfem::NewtonSolver> _nonlinear_solver{nullptr};
-
   hephaestus::FECollections _fecs;
   hephaestus::FESpaces _fespaces;
   hephaestus::GridFunctions _gridfunctions;
@@ -46,7 +39,7 @@ class Problem
   int _num_procs;
 
   /// Returns a pointer to the operator. See derived classes.
-  [[nodiscard]] virtual hephaestus::ProblemOperatorInterface * GetOperator() const = 0;
+  [[nodiscard]] virtual hephaestus::ProblemOperatorBase * GetOperator() const = 0;
 
   /// Virtual method to construct the operator. Call for default problems.
   virtual void ConstructOperator() = 0;
@@ -80,8 +73,8 @@ class ProblemBuilder
   void SetSources(hephaestus::Sources & sources);
   void SetOutputs(hephaestus::Outputs & outputs);
   void SetSolverOptions(hephaestus::InputParameters & solver_options);
-  void SetJacobianPreconditioner(std::shared_ptr<mfem::Solver> preconditioner);
-  void SetJacobianSolver(std::shared_ptr<mfem::Solver> solver);
+  void SetJacobianPreconditioner(std::unique_ptr<mfem::Solver> preconditioner);
+  void SetJacobianSolver(std::unique_ptr<mfem::Solver> solver);
   void SetCoefficients(hephaestus::Coefficients & coefficients);
 
   void AddFESpace(std::string fespace_name,
@@ -100,10 +93,8 @@ class ProblemBuilder
   virtual void RegisterAuxSolvers() = 0;
   virtual void RegisterCoefficients() = 0;
 
-  virtual void ConstructJacobianPreconditioner();
   virtual void ConstructJacobianSolver();
   virtual void ConstructNonlinearSolver();
-  virtual void ConstructBlockVector();
 
   virtual void ConstructOperator() = 0;
   virtual void ConstructTimestepper() = 0;
diff --git a/src/problem_builders/time_domain_problem_builder.cpp b/src/problem_builders/time_domain_problem_builder.cpp
index b1a62099e..916857a8f 100644
--- a/src/problem_builders/time_domain_problem_builder.cpp
+++ b/src/problem_builders/time_domain_problem_builder.cpp
@@ -3,17 +3,6 @@
 namespace hephaestus
 {
 
-void
-TimeDomainProblem::Update()
-{
-  // 1. Call superclass. Updates the offsets and block vector.
-  Problem::Update();
-
-  // 2. The dimensions of the problem operator have now changed. We must call
-  // the ode solver's Init method to resize its internal vector.
-  _ode_solver->Init(*GetOperator());
-}
-
 std::vector<mfem::ParGridFunction *>
 TimeDomainProblemBuilder::RegisterTimeDerivatives(std::vector<std::string> gridfunction_names,
                                                   hephaestus::GridFunctions & gridfunctions)
@@ -59,8 +48,10 @@ TimeDomainProblemBuilder::InitializeOperator()
 void
 TimeDomainProblemBuilder::ConstructTimestepper()
 {
-  GetProblem()->_ode_solver = std::make_unique<mfem::BackwardEulerSolver>();
-  GetProblem()->_ode_solver->Init(*(GetProblem()->GetOperator()));
+  auto ode_solver = std::make_unique<mfem::BackwardEulerSolver>();
+  ode_solver->Init(*GetProblem()->GetOperator());
+
+  GetProblem()->GetOperator()->SetODESolver(std::move(ode_solver));
 }
 
 } // namespace hephaestus
diff --git a/src/problem_builders/time_domain_problem_builder.hpp b/src/problem_builders/time_domain_problem_builder.hpp
index 50c5e3caf..231b20586 100644
--- a/src/problem_builders/time_domain_problem_builder.hpp
+++ b/src/problem_builders/time_domain_problem_builder.hpp
@@ -31,8 +31,6 @@ class TimeDomainProblem : public Problem
     _problem_operator = std::make_unique<hephaestus::TimeDomainProblemOperator>(*this);
   }
 
-  void Update() override;
-
 private:
   std::unique_ptr<hephaestus::TimeDomainProblemOperator> _problem_operator{nullptr};
 };
diff --git a/src/problem_operators/equation_system_problem_operator.hpp b/src/problem_operators/equation_system_problem_operator.hpp
index 5e4ec50e6..5f583b4e8 100644
--- a/src/problem_operators/equation_system_problem_operator.hpp
+++ b/src/problem_operators/equation_system_problem_operator.hpp
@@ -1,6 +1,6 @@
 #pragma once
 #include "problem_operator.hpp"
-#include "problem_operator_interface.hpp"
+#include "problem_operator_base.hpp"
 #include "equation_system_interface.hpp"
 
 namespace hephaestus
diff --git a/src/problem_operators/problem_operator.cpp b/src/problem_operators/problem_operator.cpp
deleted file mode 100644
index 5e9efa441..000000000
--- a/src/problem_operators/problem_operator.cpp
+++ /dev/null
@@ -1,12 +0,0 @@
-#include "problem_operator.hpp"
-
-namespace hephaestus
-{
-
-void
-ProblemOperator::UpdateOperatorWidthAndHeight()
-{
-  width = height = _true_offsets.Last();
-}
-
-} // namespace hephaestus
\ No newline at end of file
diff --git a/src/problem_operators/problem_operator.hpp b/src/problem_operators/problem_operator.hpp
index 1aba00690..12d151272 100644
--- a/src/problem_operators/problem_operator.hpp
+++ b/src/problem_operators/problem_operator.hpp
@@ -2,21 +2,25 @@
 #include "../common/pfem_extras.hpp"
 #include "hephaestus_solvers.hpp"
 #include "problem_builder_base.hpp"
-#include "problem_operator_interface.hpp"
+#include "problem_operator_base.hpp"
 
 namespace hephaestus
 {
 /// Steady-state problem operator with no equation system.
-class ProblemOperator : public mfem::Operator, public ProblemOperatorInterface
+class ProblemOperator : public mfem::Operator, public ProblemOperatorBase
 {
 public:
-  ProblemOperator(hephaestus::Problem & problem) : ProblemOperatorInterface(problem) {}
+  ProblemOperator(hephaestus::Problem & problem) : ProblemOperatorBase(problem) {}
   ~ProblemOperator() override = default;
 
+  virtual void Solve() { Solve(*_block_vector); }
+
   virtual void Solve(mfem::Vector & X) {}
   void Mult(const mfem::Vector & x, mfem::Vector & y) const override {}
 
-  void UpdateOperatorWidthAndHeight() final;
+protected:
+  int & Width() final { return mfem::Operator::width; }
+  int & Height() final { return mfem::Operator::height; }
 };
 
 } // namespace hephaestus
\ No newline at end of file
diff --git a/src/problem_operators/problem_operator_base.cpp b/src/problem_operators/problem_operator_base.cpp
new file mode 100644
index 000000000..55aeeb746
--- /dev/null
+++ b/src/problem_operators/problem_operator_base.cpp
@@ -0,0 +1,101 @@
+#include "problem_operator_base.hpp"
+
+namespace hephaestus
+{
+
+ProblemOperatorBase::ProblemOperatorBase(hephaestus::Problem & problem) : _problem{problem}
+{
+  _block_vector = std::make_unique<mfem::BlockVector>();
+}
+
+void
+ProblemOperatorBase::SetTrialVariables()
+{
+  SetTrialVariableNames();
+
+  _trial_variables = _problem._gridfunctions.Get(_trial_var_names);
+}
+
+int
+ProblemOperatorBase::GetSolutionVectorSize() const
+{
+  return static_cast<int>(GetTrialVariablesSize());
+}
+
+int
+ProblemOperatorBase::GetTrialVariablesSize() const
+{
+  return static_cast<int>(_trial_variables.size());
+}
+
+void
+ProblemOperatorBase::UpdateOffsets()
+{
+  if (GetSolutionVectorSize() > GetTrialVariablesSize())
+  {
+    MFEM_ABORT("Solution vector size (" << GetSolutionVectorSize()
+                                        << ") cannot exceed the number of trial variables ("
+                                        << GetTrialVariablesSize() << ").");
+  }
+
+  _block_true_offsets.SetSize(GetSolutionVectorSize() + 1);
+  _true_offsets.SetSize(GetTrialVariablesSize() + 1);
+
+  _block_true_offsets[0] = _true_offsets[0] = 0;
+  for (int i = 0; i < GetTrialVariablesSize(); ++i)
+  {
+    mfem::ParFiniteElementSpace * fespace = _trial_variables.at(i)->ParFESpace();
+
+    if (i < GetSolutionVectorSize())
+      _block_true_offsets[i + 1] = fespace->GetTrueVSize();
+
+    _true_offsets[i + 1] = fespace->GetVSize();
+  }
+
+  // Partial sum over values to calculate offsets.
+  _block_true_offsets.PartialSum();
+  _true_offsets.PartialSum();
+
+  // Update block vectors.
+  _true_x.Update(_block_true_offsets);
+  _true_rhs.Update(_block_true_offsets);
+
+  // Set width and height member variables.
+  Width() = Height() = _true_offsets.Last();
+}
+
+void
+ProblemOperatorBase::UpdateBlockVector(mfem::BlockVector & X)
+{
+  X.Update(_true_offsets);
+
+  for (int i = 0; i < GetTrialVariablesSize(); ++i)
+  {
+    mfem::ParGridFunction * trial_var = _trial_variables.at(i);
+
+    trial_var->MakeRef(trial_var->ParFESpace(), X, _true_offsets[i]);
+    *trial_var = 0.0;
+  }
+}
+
+void
+ProblemOperatorBase::Init()
+{
+  SetTrialVariables();
+
+  UpdateOffsets();
+  UpdateBlockVector(*_block_vector);
+}
+
+void
+ProblemOperatorBase::Update()
+{
+  logger.debug("Update called for ProblemOperatorBase.");
+
+  // Recalculate the offsets from gridfunction trial variables.
+  UpdateOffsets();
+
+  UpdateBlockVector(*_block_vector);
+}
+
+}
\ No newline at end of file
diff --git a/src/problem_operators/problem_operator_base.hpp b/src/problem_operators/problem_operator_base.hpp
new file mode 100644
index 000000000..323c373ab
--- /dev/null
+++ b/src/problem_operators/problem_operator_base.hpp
@@ -0,0 +1,107 @@
+#pragma once
+#include "mfem.hpp"
+#include "problem_builder_base.hpp"
+#include <vector>
+#include <string>
+
+namespace hephaestus
+{
+class ProblemOperatorBase
+{
+public:
+  ProblemOperatorBase() = delete;
+
+  virtual ~ProblemOperatorBase() = default;
+
+  /// Initialize the problem operator.
+  virtual void Init();
+
+  /// Update the problem operator after a mesh change.
+  virtual void Update();
+
+  /// Set the Jacobian preconditioner.
+  void SetJacobianPreconditioner(std::unique_ptr<mfem::Solver> preconditioner)
+  {
+    _jacobian_preconditioner = std::move(preconditioner);
+  }
+
+  /// Set the Jacobian solver.
+  void SetJacobianSolver(std::unique_ptr<mfem::Solver> solver)
+  {
+    _jacobian_solver = std::move(solver);
+  }
+
+  /// Set the nonlinear solver.
+  void SetNonlinearSolver(std::unique_ptr<mfem::NewtonSolver> nl_solver)
+  {
+    _nonlinear_solver = std::move(nl_solver);
+  }
+
+  /// Accessor for Jacobian preconditioner.
+  template <class TSolver>
+  TSolver * JacobianPreconditioner() const
+  {
+    return static_cast<TSolver *>(_jacobian_preconditioner.get());
+  }
+
+protected:
+  /// Use of protected constructor to only allow construction by derived classes.
+  /// All problem operator classes are built on-top of this class and it should not
+  /// be possible to use directly.
+  explicit ProblemOperatorBase(hephaestus::Problem & problem);
+
+  /// Set trial variables names. Override in derived classes.
+  virtual void SetTrialVariableNames() {}
+
+  /// Set gridfunction of trial variables from the trial variable names.
+  virtual void SetTrialVariables();
+
+  /// Returns a reference to the operator's width.
+  virtual int & Width() = 0;
+
+  /// Returns a reference to the operator's height.
+  virtual int & Height() = 0;
+
+  /// Returns the solution vector size. By default this will be the same as the
+  /// number of trial variables.
+  virtual int GetSolutionVectorSize() const;
+
+  /// Returns the number of trial variables.
+  int GetTrialVariablesSize() const;
+
+  // Reference to the current problem.
+  hephaestus::Problem & _problem;
+
+  /// Vector of names of state gridfunctions used in formulation, ordered by
+  /// appearance in block vector during solve.
+  std::vector<std::string> _trial_var_names;
+
+  /// Vector of trial variables set by SetTrialVariables.
+  std::vector<mfem::ParGridFunction *> _trial_variables;
+
+  /// Arrays of offsets.
+  mfem::Array<int> _true_offsets, _block_true_offsets;
+
+  /// Block vectors.
+  mfem::BlockVector _true_x, _true_rhs;
+
+  /// Store the Jacobian preconditioner.
+  std::unique_ptr<mfem::Solver> _jacobian_preconditioner{nullptr};
+
+  /// Store the Jacobian solver.
+  std::unique_ptr<mfem::Solver> _jacobian_solver{nullptr};
+
+  /// Store the non-linear solver.
+  std::unique_ptr<mfem::NewtonSolver> _nonlinear_solver{nullptr};
+
+  /// Store the block vector.
+  std::unique_ptr<mfem::BlockVector> _block_vector{nullptr};
+
+private:
+  /// Update the block vectors and offsets after a mesh change.
+  void UpdateOffsets();
+
+  /// Update a block vector. Should be called after the offsets have been updated.
+  void UpdateBlockVector(mfem::BlockVector & X);
+};
+}
\ No newline at end of file
diff --git a/src/problem_operators/problem_operator_interface.cpp b/src/problem_operators/problem_operator_interface.cpp
deleted file mode 100644
index 5fd4d291b..000000000
--- a/src/problem_operators/problem_operator_interface.cpp
+++ /dev/null
@@ -1,90 +0,0 @@
-#include "problem_operator_interface.hpp"
-
-namespace hephaestus
-{
-
-void
-ProblemOperatorInterface::SetTrialVariables()
-{
-  SetTrialVariableNames();
-
-  _trial_variables = _problem._gridfunctions.Get(_trial_var_names);
-}
-
-void
-ProblemOperatorInterface::UpdateOffsets()
-{
-  UpdateOffsetsWithSize(_trial_variables.size());
-}
-
-void
-ProblemOperatorInterface::UpdateOffsetsWithSize(const size_t block_vector_size)
-{
-  if (block_vector_size > _trial_variables.size())
-  {
-    MFEM_ABORT("Solution vector size (" << block_vector_size
-                                        << ") cannot exceed the number of trial variables ("
-                                        << _trial_variables.size() << ").");
-  }
-
-  _block_true_offsets.SetSize(block_vector_size + 1);
-  _true_offsets.SetSize(_trial_variables.size() + 1);
-
-  _block_true_offsets[0] = _true_offsets[0] = 0;
-  for (size_t i = 0; i < _trial_variables.size(); ++i)
-  {
-    mfem::ParFiniteElementSpace * fespace = _trial_variables.at(i)->ParFESpace();
-
-    if (i < block_vector_size)
-      _block_true_offsets[i + 1] = fespace->GetTrueVSize();
-
-    _true_offsets[i + 1] = fespace->GetVSize();
-  }
-
-  // Partial sum over values to calculate offsets.
-  _block_true_offsets.PartialSum();
-  _true_offsets.PartialSum();
-
-  // Update block vectors.
-  _true_x.Update(_block_true_offsets);
-  _true_rhs.Update(_block_true_offsets);
-
-  // Set width and height member variables.
-  UpdateOperatorWidthAndHeight();
-}
-
-void
-ProblemOperatorInterface::UpdateBlockVector(mfem::BlockVector & X)
-{
-  X.Update(_true_offsets);
-
-  for (size_t i = 0; i < _trial_variables.size(); ++i)
-  {
-    mfem::ParGridFunction * trial_var = _trial_variables.at(i);
-
-    trial_var->MakeRef(trial_var->ParFESpace(), X, _true_offsets[i]);
-    *trial_var = 0.0;
-  }
-}
-
-void
-ProblemOperatorInterface::Init()
-{
-  SetTrialVariables();
-
-  UpdateOffsets();
-  UpdateBlockVector(*_problem._f);
-}
-
-void
-ProblemOperatorInterface::Update()
-{
-  logger.debug("Update called for ProblemOperatorInterface.");
-
-  // Recalculate the offsets from gridfunction trial variables.
-  UpdateOffsets();
-
-  UpdateBlockVector(*_problem._f);
-}
-
-}
\ No newline at end of file
diff --git a/src/problem_operators/problem_operator_interface.hpp b/src/problem_operators/problem_operator_interface.hpp
deleted file mode 100644
index af91e7759..000000000
--- a/src/problem_operators/problem_operator_interface.hpp
+++ /dev/null
@@ -1,41 +0,0 @@
-#pragma once
-#include "problem_builder_base.hpp"
-
-namespace hephaestus
-{
-/// Interface inherited by ProblemOperator and TimeDomainProblemOperator. Removes duplicated code in both classes.
-class ProblemOperatorInterface
-{
-public:
-  ProblemOperatorInterface(hephaestus::Problem & problem) : _problem(problem) {}
-  virtual ~ProblemOperatorInterface() = default;
-
-  virtual void UpdateOperatorWidthAndHeight() {}
-
-  virtual void Init();
-
-  virtual void Update();
-
-  mfem::Array<int> _true_offsets, _block_true_offsets;
-
-  mfem::BlockVector _true_x, _true_rhs;
-  mfem::OperatorHandle _equation_system_operator;
-
-protected:
-  virtual void SetTrialVariableNames() {}
-  virtual void SetTrialVariables();
-  virtual void UpdateOffsets();
-
-  virtual void UpdateBlockVector(mfem::BlockVector & X);
-
-  void UpdateOffsetsWithSize(size_t soln_vector_size);
-
-  // Reference to the current problem.
-  hephaestus::Problem & _problem;
-
-  // Vector of names of state gridfunctions used in formulation, ordered by appearance in block
-  // vector during solve.
-  std::vector<std::string> _trial_var_names;
-  std::vector<mfem::ParGridFunction *> _trial_variables;
-};
-}
\ No newline at end of file
diff --git a/src/problem_operators/time_domain_equation_system_problem_operator.cpp b/src/problem_operators/time_domain_equation_system_problem_operator.cpp
index 94b6383f9..ccc210c42 100644
--- a/src/problem_operators/time_domain_equation_system_problem_operator.cpp
+++ b/src/problem_operators/time_domain_equation_system_problem_operator.cpp
@@ -55,21 +55,20 @@ TimeDomainEquationSystemProblemOperator::Update()
     // Rebuild the jacobian preconditioner.
     auto first_pfespace = GetEquationSystem()->_test_pfespaces.at(0);
 
-    auto precond = std::make_shared<mfem::HypreAMS>(first_pfespace);
+    auto precond = std::make_unique<mfem::HypreAMS>(first_pfespace);
 
     precond->SetSingularProblem();
     precond->SetPrintLevel(-1);
 
-    _problem._jacobian_preconditioner = precond;
+    _jacobian_preconditioner = std::move(precond);
 
     // Set new preconditioner.
-    std::static_pointer_cast<mfem::HyprePCG>(_problem._jacobian_solver)
-        ->SetPreconditioner(
-            *std::static_pointer_cast<mfem::HypreSolver>(_problem._jacobian_preconditioner));
+    static_cast<mfem::HyprePCG *>(_jacobian_solver.get())
+        ->SetPreconditioner(*static_cast<mfem::HypreSolver *>(_jacobian_preconditioner.get()));
 
     // Set Jacobian matrix.
     auto * matrix = GetEquationSystem()->JacobianOperatorHandle().As<mfem::HypreParMatrix>();
-    _problem._jacobian_solver->SetOperator(*matrix);
+    _jacobian_solver->SetOperator(*matrix);
   }
 }
 
@@ -89,9 +88,9 @@ TimeDomainEquationSystemProblemOperator::ImplicitSolve(const double dt,
   _problem._coefficients.SetTime(GetTime());
   BuildEquationSystemOperator(dt);
 
-  _problem._nonlinear_solver->SetSolver(*_problem._jacobian_solver);
-  _problem._nonlinear_solver->SetOperator(*GetEquationSystem());
-  _problem._nonlinear_solver->Mult(_true_rhs, _true_x);
+  _nonlinear_solver->SetSolver(*_jacobian_solver);
+  _nonlinear_solver->SetOperator(*GetEquationSystem());
+  _nonlinear_solver->Mult(_true_rhs, _true_x);
 
   GetEquationSystem()->RecoverFEMSolution(_true_x, _problem._gridfunctions);
 }
diff --git a/src/problem_operators/time_domain_equation_system_problem_operator.hpp b/src/problem_operators/time_domain_equation_system_problem_operator.hpp
index 3aa518753..5e83e6cbd 100644
--- a/src/problem_operators/time_domain_equation_system_problem_operator.hpp
+++ b/src/problem_operators/time_domain_equation_system_problem_operator.hpp
@@ -1,7 +1,7 @@
 #pragma once
 #include "../common/pfem_extras.hpp"
 #include "time_domain_problem_operator.hpp"
-#include "problem_operator_interface.hpp"
+#include "problem_operator_base.hpp"
 #include "equation_system_interface.hpp"
 
 namespace hephaestus
diff --git a/src/problem_operators/time_domain_problem_operator.cpp b/src/problem_operators/time_domain_problem_operator.cpp
index ecde0de8c..a4ffa7211 100644
--- a/src/problem_operators/time_domain_problem_operator.cpp
+++ b/src/problem_operators/time_domain_problem_operator.cpp
@@ -21,9 +21,16 @@ GetTimeDerivativeNames(std::vector<std::string> gridfunction_names)
 }
 
 void
-TimeDomainProblemOperator::UpdateOperatorWidthAndHeight()
+TimeDomainProblemOperator::Update()
 {
-  width = height = _true_offsets.Last();
+  ProblemOperatorBase::Update();
+
+  // The dimensions of the problem operator have now changed. We must call the
+  // ODE solver's Init method to resize its internal vector.
+  if (_ode_solver)
+  {
+    _ode_solver->Init(*this);
+  }
 }
 
 } // namespace hephaestus
\ No newline at end of file
diff --git a/src/problem_operators/time_domain_problem_operator.hpp b/src/problem_operators/time_domain_problem_operator.hpp
index 095bd4fc6..d36f17365 100644
--- a/src/problem_operators/time_domain_problem_operator.hpp
+++ b/src/problem_operators/time_domain_problem_operator.hpp
@@ -2,7 +2,7 @@
 #include "../common/pfem_extras.hpp"
 #include "hephaestus_solvers.hpp"
 #include "problem_builder_base.hpp"
-#include "problem_operator_interface.hpp"
+#include "problem_operator_base.hpp"
 
 namespace hephaestus
 {
@@ -13,16 +13,35 @@ std::vector<std::string> GetTimeDerivativeNames(std::vector<std::string> gridfun
 
 /// Problem operator for time-dependent problems with no equation system. The user will need to subclass this since the solve is not
 /// implemented.
-class TimeDomainProblemOperator : public mfem::TimeDependentOperator,
-                                  public ProblemOperatorInterface
+class TimeDomainProblemOperator : public mfem::TimeDependentOperator, public ProblemOperatorBase
 {
 public:
-  TimeDomainProblemOperator(hephaestus::Problem & problem) : ProblemOperatorInterface(problem) {}
+  TimeDomainProblemOperator(hephaestus::Problem & problem) : ProblemOperatorBase(problem) {}
   ~TimeDomainProblemOperator() override = default;
 
-  void UpdateOperatorWidthAndHeight() final;
-
   void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override {}
+
+  /// Set the ODE solver.
+  void SetODESolver(std::unique_ptr<mfem::ODESolver> ode_solver)
+  {
+    _ode_solver = std::move(ode_solver);
+  }
+
+  /// Wrapper around the ODE solver's Step method using the block vector.
+  virtual void Step(mfem::real_t & t, mfem::real_t & dt)
+  {
+    _ode_solver->Step(*_block_vector, t, dt);
+  }
+
+  void Update() override;
+
+protected:
+  int & Width() final { return mfem::TimeDependentOperator::width; }
+  int & Height() final { return mfem::TimeDependentOperator::height; }
+
+private:
+  /// Store the ODE solver.
+  std::unique_ptr<mfem::ODESolver> _ode_solver{nullptr};
 };
 
 } // namespace hephaestus
\ No newline at end of file