diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4d9af1b0..73117b75 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -147,6 +147,7 @@ set(HEADER_FILES
     ${INC_DIR}/nektar_interface/bounding_box_intersection.hpp
     ${INC_DIR}/nektar_interface/cell_id_translation.hpp
     ${INC_DIR}/nektar_interface/coordinate_mapping.hpp
+    ${INC_DIR}/nektar_interface/solver_base/eqnsys_base.hpp
     ${INC_DIR}/nektar_interface/expansion_looping/basis_evaluate_base.hpp
     ${INC_DIR}/nektar_interface/expansion_looping/expansion_looping.hpp
     ${INC_DIR}/nektar_interface/expansion_looping/geom_to_expansion_builder.hpp
@@ -176,6 +177,7 @@ set(HEADER_FILES
     ${INC_DIR}/nektar_interface/geometry_transport/remote_geom_3d.hpp
     ${INC_DIR}/nektar_interface/geometry_transport/shape_mapping.hpp
     ${INC_DIR}/nektar_interface/parameter_store.hpp
+    ${INC_DIR}/nektar_interface/solver_base/partsys_base.hpp
     ${INC_DIR}/nektar_interface/particle_boundary_conditions.hpp
     ${INC_DIR}/nektar_interface/particle_cell_mapping/coarse_lookup_map.hpp
     ${INC_DIR}/nektar_interface/particle_cell_mapping/coarse_mappers_base.hpp
@@ -209,6 +211,7 @@ set(HEADER_FILES
     ${INC_DIR}/nektar_interface/particle_interface.hpp
     ${INC_DIR}/nektar_interface/particle_mesh_interface.hpp
     ${INC_DIR}/nektar_interface/special_functions.hpp
+    ${INC_DIR}/nektar_interface/solver_base/time_evolved_eqnsys_base.hpp
     ${INC_DIR}/nektar_interface/utilities.hpp
     ${INC_DIR}/nektar_interface/utility_mesh_cartesian.hpp
     ${INC_DIR}/nektar_interface/utility_mesh_plotting.hpp
diff --git a/include/nektar_interface/solver_base/eqnsys_base.hpp b/include/nektar_interface/solver_base/eqnsys_base.hpp
new file mode 100644
index 00000000..e1e9b222
--- /dev/null
+++ b/include/nektar_interface/solver_base/eqnsys_base.hpp
@@ -0,0 +1,142 @@
+#ifndef __EQNSYS_BASE_H_
+#define __EQNSYS_BASE_H_
+
+#include <SolverUtils/EquationSystem.h>
+
+#include <type_traits>
+
+#include "nektar_interface/solver_base/partsys_base.hpp"
+#include "nektar_interface/utilities.hpp"
+
+namespace LU = Nektar::LibUtilities;
+namespace SD = Nektar::SpatialDomains;
+namespace SU = Nektar::SolverUtils;
+
+namespace NESO::Solvers {
+
+// Fwd declare NEKEQNSYS, PARTSYS as classes
+class NEKEQNSYS;
+class PARTSYS;
+
+/**
+ * @brief Base class for Nektar++ equation systems, coupled to a NESO-Particles
+ * particle system derived from NESO::Solvers::PartSysBase.
+ */
+template <typename NEKEQNSYS, typename PARTSYS>
+class EqnSysBase : public NEKEQNSYS {
+  // Template param must derive from Nektar's EquationSystem base class
+  static_assert(std::is_base_of<SU::EquationSystem, NEKEQNSYS>(),
+                "Template arg to EqnSysBase must derive from "
+                "Nektar::SolverUtils::EquationSystem");
+  // Particle system must derive from PartSysBase
+  static_assert(
+      std::is_base_of<PartSysBase, PARTSYS>(),
+      "PARTSYS template arg to EqnSysBase must derive from PartSysBase");
+
+protected:
+  EqnSysBase(const LU::SessionReaderSharedPtr &session,
+             const SD::MeshGraphSharedPtr &graph)
+      : NEKEQNSYS(session, graph), field_to_index(session->GetVariables()),
+        required_fld_names() {
+
+    // If number of particles / number per cell was set in config; construct the
+    // particle system
+    int num_parts_per_cell, num_parts_tot;
+    session->LoadParameter(PartSysBase::NUM_PARTS_TOT_STR, num_parts_tot, -1);
+    session->LoadParameter(PartSysBase::NUM_PARTS_PER_CELL_STR,
+                           num_parts_per_cell, -1);
+    this->particles_enabled = num_parts_tot > 0 || num_parts_per_cell > 0;
+    if (this->particles_enabled) {
+      this->particle_sys = std::make_shared<PARTSYS>(session, graph);
+    }
+  }
+
+  /// Field name => index mapper
+  NESO::NektarFieldIndexMap field_to_index;
+
+  /// Particle system
+  std::shared_ptr<PARTSYS> particle_sys;
+
+  /// Flag identifying whether particles were enabled in the config file
+  bool particles_enabled;
+
+  /// List of field names required by the solver
+  std::vector<std::string> required_fld_names;
+
+  /// Placeholder for subclasses to override; called in v_InitObject()
+  virtual void load_params(){};
+
+  /**
+   * @brief Check that all required fields are defined. All fields must have the
+   * same number of quad points for now.
+   */
+  void validate_fields() {
+    int npts_exp = NEKEQNSYS::GetNpoints();
+    for (auto &fld_name : this->required_fld_names) {
+      int idx = this->field_to_index.get_idx(fld_name);
+      // Check field exists
+
+      std::string err_msg = "Required field [" + fld_name + "] is not defined.";
+      NESOASSERT(idx >= 0, err_msg.c_str());
+
+      // Check fields all have the same number of quad points
+      int npts = this->m_fields[idx]->GetNpoints();
+      err_msg = "Expecting " + std::to_string(npts_exp) +
+                " quad points, but field '" + fld_name + "' has " +
+                std::to_string(npts) +
+                ". Check NUMMODES is the same for all required fields.";
+      NESOASSERT(npts == npts_exp, err_msg.c_str());
+    }
+  }
+
+  /** @brief Write particle params to stdout on task 0. Ensures they appear just
+   * after fluid params are written by nektar.
+   *
+   * */
+  virtual void v_DoInitialise() override {
+    if (this->m_session->GetComm()->TreatAsRankZero() &&
+        this->particles_enabled) {
+      particle_sys->add_params_report();
+    }
+    NEKEQNSYS::v_DoInitialise();
+  }
+
+  /**
+   * @brief Free particle system memory after solver loop has finished.
+   * Prevent further overrides to guarantee that subclasses do the same.
+   */
+  virtual void v_DoSolve() override final {
+    NEKEQNSYS::v_DoSolve();
+    if (this->particle_sys) {
+      this->particle_sys->free();
+    }
+  }
+
+  /**
+   * @brief Initialise the equation system, then check required fields are set
+   * and load parameters.
+   */
+  virtual void v_InitObject(bool create_fields) override {
+    NEKEQNSYS::v_InitObject(create_fields);
+
+    // Ensure that the session file defines all required variables
+    validate_fields();
+
+    // Load parameters
+    load_params();
+  }
+
+  /**
+   * @brief Utility function to zero a Nektar++ array of arrays
+   *
+   * @param arr Array-of-arrays to zero
+   */
+  void zero_array_of_arrays(Array<OneD, Array<OneD, NekDouble>> &arr) {
+    for (auto ii = 0; ii < arr.size(); ii++) {
+      Vmath::Zero(arr[ii].size(), arr[ii], 1);
+    }
+  }
+};
+
+} // namespace NESO::Solvers
+#endif
\ No newline at end of file
diff --git a/include/nektar_interface/solver_base/partsys_base.hpp b/include/nektar_interface/solver_base/partsys_base.hpp
new file mode 100644
index 00000000..8b30416b
--- /dev/null
+++ b/include/nektar_interface/solver_base/partsys_base.hpp
@@ -0,0 +1,241 @@
+#ifndef __PARTSYS_BASE_H_
+#define __PARTSYS_BASE_H_
+
+#include <SolverUtils/EquationSystem.h>
+#include <mpi.h>
+#include <nektar_interface/geometry_transport/halo_extension.hpp>
+#include <nektar_interface/particle_interface.hpp>
+#include <neso_particles.hpp>
+#include <type_traits>
+
+namespace LU = Nektar::LibUtilities;
+namespace SD = Nektar::SpatialDomains;
+
+namespace NESO::Particles {
+
+/// Struct used to set common options for particle systems
+struct PartSysOptions {
+  int extend_halos_offset = 0;
+};
+
+class PartSysBase {
+
+public:
+  // Some parameter names used in solver config files
+  inline static const std::string NUM_PARTS_TOT_STR = "num_particles_total";
+  inline static const std::string NUM_PARTS_PER_CELL_STR =
+      "num_particles_per_cell";
+  inline static const std::string PART_OUTPUT_FREQ_STR = "particle_output_freq";
+
+  /// Total number of particles in simulation
+  int64_t num_parts_tot;
+
+  /// NESO-Particles ParticleGroup
+  ParticleGroupSharedPtr particle_group;
+  /// Compute target
+  SYCLTargetSharedPtr sycl_target;
+
+  /**
+   * @brief Write particle parameter values to stdout for any parameter.
+   *  @see also report_param()
+   */
+  inline void add_params_report() {
+    std::cout << "Particle settings:" << std::endl;
+    for (auto const &[param_lbl, param_str_val] : this->param_vals_to_report) {
+      std::cout << "  " << param_lbl << ": " << param_str_val << std::endl;
+    }
+    std::cout << "============================================================="
+                 "=========="
+              << std::endl
+              << std::endl;
+  }
+
+  /**
+   * @brief Clear up memory related to the particle system
+   */
+  inline void free() {
+    if (this->h5part_exists) {
+      this->h5part->close();
+    }
+    this->particle_group->free();
+    this->sycl_target->free();
+    this->particle_mesh_interface->free();
+  };
+
+  /**
+   *  @brief Write particle properties to an output file.
+   *
+   *  @param step Time step number.
+   */
+  inline void write(const int step) {
+    if (this->h5part_exists) {
+      if (this->sycl_target->comm_pair.rank_parent == 0) {
+        nprint("Writing particle properties at step", step);
+      }
+      this->h5part->write();
+    } else {
+      if (this->sycl_target->comm_pair.rank_parent == 0) {
+        nprint("Ignoring call to write particle data because an output file "
+               "wasn't set up. init_output() not called?");
+      }
+    }
+  }
+
+protected:
+  /**
+   * Protected constructor to prohibit direct instantiation.
+   *
+   *  @param session Nektar++ session to use for parameters and simulation
+   * specification.
+   *  @param graph Nektar++ MeshGraph on which particles exist.
+   *  @param comm (optional) MPI communicator to use - default MPI_COMM_WORLD.
+   *
+   */
+  PartSysBase(const LU::SessionReaderSharedPtr session,
+              const SD::MeshGraphSharedPtr graph, ParticleSpec particle_spec,
+              MPI_Comm comm = MPI_COMM_WORLD,
+              PartSysOptions options = PartSysOptions())
+      : session(session), graph(graph), comm(comm), h5part_exists(false),
+        ndim(graph->GetSpaceDimension()) {
+
+    read_params();
+
+    // Store options
+    this->options = options;
+
+    // Create interface between particles and nektar++
+    this->particle_mesh_interface =
+        std::make_shared<ParticleMeshInterface>(graph, 0, this->comm);
+    extend_halos_fixed_offset(this->options.extend_halos_offset,
+                              this->particle_mesh_interface);
+    this->sycl_target =
+        std::make_shared<SYCLTarget>(0, particle_mesh_interface->get_comm());
+    this->nektar_graph_local_mapper = std::make_shared<NektarGraphLocalMapper>(
+        this->sycl_target, this->particle_mesh_interface);
+    this->domain = std::make_shared<Domain>(this->particle_mesh_interface,
+                                            this->nektar_graph_local_mapper);
+
+    // Create ParticleGroup
+    this->particle_group =
+        std::make_shared<ParticleGroup>(domain, particle_spec, sycl_target);
+
+    // Set up map between cell indices
+    this->cell_id_translation = std::make_shared<CellIDTranslation>(
+        this->sycl_target, this->particle_group->cell_id_dat,
+        this->particle_mesh_interface);
+  }
+
+  /// Object used to map to/from nektar geometry ids to 0,N-1
+  std::shared_ptr<CellIDTranslation> cell_id_translation;
+  /// MPI communicator
+  MPI_Comm comm;
+  /// NESO-Particles domain.
+  DomainSharedPtr domain;
+  /// Pointer to Nektar Meshgraph object
+  SD::MeshGraphSharedPtr graph;
+  /// HDF5 output file
+  std::shared_ptr<H5Part> h5part;
+  /// HDF5 output file flag
+  bool h5part_exists;
+  /// Number of spatial dimensions being used
+  const int ndim;
+  /// Mapping instance to map particles into nektar++ elements.
+  std::shared_ptr<NektarGraphLocalMapper> nektar_graph_local_mapper;
+  /// Options struct
+  PartSysOptions options;
+  /// HMesh instance that allows particles to move over nektar++ meshes.
+  ParticleMeshInterfaceSharedPtr particle_mesh_interface;
+  /// Pointer to Session object
+  LU::SessionReaderSharedPtr session;
+
+  /**
+   * @brief Set up per-step particle output
+   *
+   *  @param fname Output filename. Default is 'particle_trajectory.h5part'.
+   *  @param args Remaining arguments (variable length) should be sym instances
+   *  indicating which ParticleDats are to be written.
+   */
+  template <typename... T>
+  inline void init_output(std::string fname, T... args) {
+    if (this->h5part_exists) {
+      if (this->sycl_target->comm_pair.rank_parent == 0) {
+        nprint("Ignoring (duplicate?) call to init_output().");
+      }
+    } else {
+      // Create H5Part instance
+      this->h5part =
+          std::make_shared<H5Part>(fname, this->particle_group, args...);
+      this->h5part_exists = true;
+    }
+  }
+
+  /**
+   * @brief Store particle param values in a map.Values are reported later via
+   * add_params_report()
+   *
+   * @param label Label to attach to the parameter
+   * @param value Value of the parameter that was set
+   */
+  template <typename T> void report_param(std::string label, T val) {
+    // form stringstream and store string value in private map
+    std::stringstream ss;
+    ss << val;
+    param_vals_to_report[label] = ss.str();
+  }
+
+  /**
+   * @brief Read some parameters associated with all particle systems.
+   */
+  inline void read_params() {
+
+    // Read total number of particles / number per cell from config
+    int num_parts_per_cell, num_parts_tot;
+    this->session->LoadParameter(NUM_PARTS_TOT_STR, num_parts_tot, -1);
+    this->session->LoadParameter(NUM_PARTS_PER_CELL_STR, num_parts_per_cell,
+                                 -1);
+
+    if (num_parts_tot > 0) {
+      this->num_parts_tot = num_parts_tot;
+      if (num_parts_per_cell > 0) {
+        nprint("Ignoring value of '" + NUM_PARTS_PER_CELL_STR +
+               "' because  "
+               "'" +
+               NUM_PARTS_TOT_STR + "' was specified.");
+      }
+    } else {
+      if (num_parts_per_cell > 0) {
+        // Determine the global number of elements
+        const int num_elements_local = this->graph->GetNumElements();
+        int num_elements_global;
+        MPICHK(MPI_Allreduce(&num_elements_local, &num_elements_global, 1,
+                             MPI_INT, MPI_SUM, this->comm));
+
+        // compute the global number of particles
+        this->num_parts_tot =
+            ((int64_t)num_elements_global) * num_parts_per_cell;
+
+        report_param("Number of particles per cell/element",
+                     num_parts_per_cell);
+      } else {
+        nprint("Particles disabled (Neither '" + NUM_PARTS_TOT_STR +
+               "' or "
+               "'" +
+               NUM_PARTS_PER_CELL_STR + "' are set)");
+      }
+    }
+    report_param("Total number of particles", this->num_parts_tot);
+
+    // Output frequency
+    int particle_output_freq;
+    this->session->LoadParameter(PART_OUTPUT_FREQ_STR, particle_output_freq, 0);
+    report_param("Output frequency (steps)", particle_output_freq);
+  }
+
+private:
+  /// Map containing parameter name,value pairs to be written to stdout when the
+  /// nektar equation system is initialised. Populated with report_param().
+  std::map<std::string, std::string> param_vals_to_report;
+};
+
+} // namespace NESO::Particles
+#endif
\ No newline at end of file
diff --git a/include/nektar_interface/solver_base/time_evolved_eqnsys_base.hpp b/include/nektar_interface/solver_base/time_evolved_eqnsys_base.hpp
new file mode 100644
index 00000000..ddc0d572
--- /dev/null
+++ b/include/nektar_interface/solver_base/time_evolved_eqnsys_base.hpp
@@ -0,0 +1,64 @@
+#ifndef __TIME_EVOLVD_EQNSYS_BASE_H_
+#define __TIME_EVOLVD_EQNSYS_BASE_H_
+
+#include <SolverUtils/UnsteadySystem.h>
+
+#include "eqnsys_base.hpp"
+
+#include <type_traits>
+
+namespace SU = Nektar::SolverUtils;
+
+namespace NESO::Solvers {
+
+/**
+ * @brief Base class for time-evolving Nektar++ equation systems, based on
+ * Nektar::SolverUtils::UnsteadySystem, coupled to a NESO-Particles
+ * particle system derived from NESO::Solvers::PartSysBase.
+ */
+template <typename NEKEQNSYS, typename PARTSYS>
+class TimeEvoEqnSysBase : public EqnSysBase<NEKEQNSYS, PARTSYS> {
+  /// Template param must derive from Nektar's UnsteadySystem class
+  static_assert(std::is_base_of<SU::UnsteadySystem, NEKEQNSYS>(),
+                "Template arg to TimeEvoEqnSysBase must derive from "
+                "Nektar::SolverUtils::UnsteadySystem");
+
+protected:
+  TimeEvoEqnSysBase(const LU::SessionReaderSharedPtr &session,
+                    const SD::MeshGraphSharedPtr &graph)
+      : EqnSysBase<NEKEQNSYS, PARTSYS>(session, graph),
+        int_fld_names(std::vector<std::string>()) {}
+
+  /// Names of fields that will be time integrated
+  std::vector<std::string> int_fld_names;
+
+  /**
+   * @brief Load common parameters associated with all time evolved equation
+   * systems.
+   */
+  virtual void load_params() override {
+    EqnSysBase<NEKEQNSYS, PARTSYS>::load_params();
+
+    // No additional params yet
+  };
+
+  /** @brief Check that the names of fields identified as time-evolving are
+   * valid.
+   *
+   */
+  virtual void v_InitObject(bool create_fields) override {
+    EqnSysBase<NEKEQNSYS, PARTSYS>::v_InitObject(create_fields);
+    // Tell UnsteadySystem to only integrate a subset of fields in time
+    this->m_intVariables.resize(this->int_fld_names.size());
+    for (auto ii = 0; ii < this->int_fld_names.size(); ii++) {
+      int var_idx = this->field_to_index.get_idx(this->int_fld_names[ii]);
+      ASSERTL0(var_idx >= 0,
+               "Setting time integration vars - GetIntFieldNames() "
+               "returned an invalid field name.");
+      this->m_intVariables[ii] = var_idx;
+    }
+  }
+};
+
+} // namespace NESO::Solvers
+#endif
\ No newline at end of file
diff --git a/include/nektar_interface/utilities.hpp b/include/nektar_interface/utilities.hpp
index d663e14b..f0ac0e7b 100644
--- a/include/nektar_interface/utilities.hpp
+++ b/include/nektar_interface/utilities.hpp
@@ -63,8 +63,19 @@ class NektarFieldIndexMap {
    *
    * @param field_name Name of field to get index for.
    * @returns Non-negative integer if field exists.
+   * @throws std::out_of_range if field doesn't exist.
    */
   int at(std::string field_name) { return this->field_to_index.at(field_name); }
+
+  /**
+   * Unlike std::map, [] operator can be used for read access, but not for
+   * write.
+   *
+   * @param field_name Name of field to get index for.
+   * @returns Non-negative integer if field exists.
+   * @throws std::out_of_range if field doesn't exist.
+   */
+  int operator[](std::string field_name) { return at(field_name); }
 };
 
 /**
diff --git a/solvers/SimpleSOL/EquationSystems/SOLSystem.cpp b/solvers/SimpleSOL/EquationSystems/SOLSystem.cpp
index 9f95757a..849586c8 100644
--- a/solvers/SimpleSOL/EquationSystems/SOLSystem.cpp
+++ b/solvers/SimpleSOL/EquationSystems/SOLSystem.cpp
@@ -10,84 +10,59 @@ std::string SOLSystem::class_name =
 
 SOLSystem::SOLSystem(const LU::SessionReaderSharedPtr &session,
                      const SD::MeshGraphSharedPtr &graph)
-    : UnsteadySystem(session, graph),
-      m_field_to_index(session->GetVariables()) {
+    : TimeEvoEqnSysBase<SU::UnsteadySystem, NeutralParticleSystem>(session,
+                                                                   graph) {
 
   // m_spacedim isn't set at this point, for some reason; use mesh dim instead
   NESOASSERT(graph->GetSpaceDimension() == 1 || graph->GetSpaceDimension() == 2,
              "Unsupported mush dimension for SOLSystem - must be 1 or 2.");
   if (graph->GetSpaceDimension() == 2) {
-    m_required_flds = {"rho", "rhou", "rhov", "E"};
+    this->required_fld_names = {"rho", "rhou", "rhov", "E"};
   } else {
-    m_required_flds = {"rho", "rhou", "E"};
-  }
-  m_int_fld_names = std::vector<std::string>(m_required_flds);
-}
-
-/**
- * Check all required fields are defined
- */
-void SOLSystem::validate_field_list() {
-  for (auto &fld_name : m_required_flds) {
-    ASSERTL0(m_field_to_index.get_idx(fld_name) >= 0,
-             "Required field [" + fld_name + "] is not defined.");
+    this->required_fld_names = {"rho", "rhou", "E"};
   }
+  int_fld_names = std::vector<std::string>(this->required_fld_names);
 }
 
 /**
  * @brief Initialization object for SOLSystem class.
  */
 void SOLSystem::v_InitObject(bool DeclareField) {
-  validate_field_list();
-  UnsteadySystem::v_InitObject(DeclareField);
-
-  // Tell UnsteadySystem to only integrate a subset of fields in time
-  // (Ignore fields that don't have a time derivative)
-  m_intVariables.resize(m_int_fld_names.size());
-  for (auto ii = 0; ii < m_int_fld_names.size(); ii++) {
-    int var_idx = m_field_to_index.get_idx(m_int_fld_names[ii]);
-    ASSERTL0(var_idx >= 0, "Setting time integration vars - GetIntFieldNames() "
-                           "returned an invalid field name.");
-    m_intVariables[ii] = var_idx;
-  }
+  TimeEvoEqnSysBase<SU::UnsteadySystem, NeutralParticleSystem>::v_InitObject(
+      DeclareField);
 
   for (int i = 0; i < m_fields.size(); i++) {
     // Use BwdTrans to make sure initial condition is in solution space
     m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), m_fields[i]->UpdatePhys());
   }
 
-  m_var_converter = MemoryManager<VariableConverter>::AllocateSharedPtr(
+  this->var_converter = MemoryManager<VariableConverter>::AllocateSharedPtr(
       m_session, m_spacedim);
 
   ASSERTL0(m_session->DefinesSolverInfo("UPWINDTYPE"),
            "No UPWINDTYPE defined in session.");
 
-  // Store velocity field indices for the Riemann solver.
-  m_vec_locs = Array<OneD, Array<OneD, NekDouble>>(1);
-  m_vec_locs[0] = Array<OneD, NekDouble>(m_spacedim);
+  // Store velocity field indices in the format required by the Riemann solver.
+  this->vel_fld_indices = Array<OneD, Array<OneD, NekDouble>>(1);
+  this->vel_fld_indices[0] = Array<OneD, NekDouble>(m_spacedim);
   for (int i = 0; i < m_spacedim; ++i) {
-    m_vec_locs[0][i] = 1 + i;
+    this->vel_fld_indices[0][i] = 1 + i;
   }
 
   // Loading parameters from session file
-  m_session->LoadParameter("Gamma", m_gamma, 1.4);
+  m_session->LoadParameter("Gamma", this->gamma, 1.4);
 
   // Setting up advection and diffusion operators
   init_advection();
 
   // Set up forcing/source term objects.
-  m_forcing = SU::Forcing::Load(m_session, shared_from_this(), m_fields,
-                                m_fields.size());
+  this->fluid_src_terms = SU::Forcing::Load(m_session, shared_from_this(),
+                                            m_fields, m_fields.size());
 
   m_ode.DefineOdeRhs(&SOLSystem::explicit_time_int, this);
   m_ode.DefineProjection(&SOLSystem::do_ode_projection, this);
 }
 
-/**
- * @brief Destructor for SOLSystem class.
- */
-SOLSystem::~SOLSystem() {}
-
 /**
  * @brief Initialisation, including creation of advection object.
  */
@@ -99,9 +74,9 @@ void SOLSystem::init_advection() {
   std::string adv_type, riemann_type;
   m_session->LoadSolverInfo("AdvectionType", adv_type, "WeakDG");
 
-  m_adv = SU::GetAdvectionFactory().CreateInstance(adv_type, adv_type);
+  this->adv_obj = SU::GetAdvectionFactory().CreateInstance(adv_type, adv_type);
 
-  m_adv->SetFluxVector(&SOLSystem::get_flux_vector, this);
+  this->adv_obj->SetFluxVector(&SOLSystem::get_flux_vector, this);
 
   // Setting up Riemann solver for advection operator
   m_session->LoadSolverInfo("UpwindType", riemann_type, "Average");
@@ -116,8 +91,8 @@ void SOLSystem::init_advection() {
   riemann_solver->SetVector("N", &SOLSystem::get_trace_norms, this);
 
   // Concluding initialisation of advection / diffusion operators
-  m_adv->SetRiemannSolver(riemann_solver);
-  m_adv->InitObject(m_session, m_fields);
+  this->adv_obj->SetRiemannSolver(riemann_solver);
+  this->adv_obj->InitObject(m_session, m_fields);
 }
 
 /**
@@ -154,7 +129,7 @@ void SOLSystem::explicit_time_int(
   }
 
   // Add forcing terms
-  for (auto &x : m_forcing) {
+  for (auto &x : this->fluid_src_terms) {
     x->Apply(m_fields, in_arr, out_arr, time);
   }
 }
@@ -178,11 +153,11 @@ void SOLSystem::do_advection(
     const Array<OneD, const Array<OneD, NekDouble>> &fwd,
     const Array<OneD, const Array<OneD, NekDouble>> &bwd) {
   // Only fields up to and including the energy need to be advected
-  int num_fields_to_advect = m_field_to_index.get_idx("E") + 1;
+  int num_fields_to_advect = this->field_to_index.get_idx("E") + 1;
 
   Array<OneD, Array<OneD, NekDouble>> adv_vel(m_spacedim);
-  m_adv->Advect(num_fields_to_advect, m_fields, adv_vel, in_arr, out_arr, time,
-                fwd, bwd);
+  this->adv_obj->Advect(num_fields_to_advect, m_fields, adv_vel, in_arr,
+                        out_arr, time, fwd, bwd);
 }
 
 /**
@@ -194,8 +169,8 @@ void SOLSystem::do_advection(
 void SOLSystem::get_flux_vector(
     const Array<OneD, const Array<OneD, NekDouble>> &fields_vals,
     TensorOfArray3D<NekDouble> &flux) {
-  const auto rho_idx = m_field_to_index.get_idx("rho");
-  const auto E_idx = m_field_to_index.get_idx("E");
+  const auto rho_idx = this->field_to_index.get_idx("rho");
+  const auto E_idx = this->field_to_index.get_idx("E");
   // Energy is the last field of relevance, regardless of mesh dimension
   const auto num_vars = E_idx + 1;
   const auto num_pts = fields_vals[0].size();
@@ -222,7 +197,7 @@ void SOLSystem::get_flux_vector(
       vel_vals_pt[dim] = field_vals_pt[dim + 1] * oneOrho;
     }
 
-    NekDouble pressure = m_var_converter->GetPressure(field_vals_pt.data());
+    NekDouble pressure = this->var_converter->GetPressure(field_vals_pt.data());
     NekDouble e_plus_P = field_vals_pt[E_idx] + pressure;
     for (auto dim = 0; dim < m_spacedim; ++dim) {
       // Flux vector for the velocity fields
diff --git a/solvers/SimpleSOL/EquationSystems/SOLSystem.hpp b/solvers/SimpleSOL/EquationSystems/SOLSystem.hpp
index 430b7605..1c111f42 100644
--- a/solvers/SimpleSOL/EquationSystems/SOLSystem.hpp
+++ b/solvers/SimpleSOL/EquationSystems/SOLSystem.hpp
@@ -1,8 +1,9 @@
 #ifndef __SIMPLESOL_SOLSYSTEM_H_
 #define __SIMPLESOL_SOLSYSTEM_H_
 
+#include "../ParticleSystems/neutral_particles.hpp"
+#include "nektar_interface/solver_base/time_evolved_eqnsys_base.hpp"
 #include "nektar_interface/utilities.hpp"
-
 #include <CompressibleFlowSolver/Misc/VariableConverter.h>
 #include <LocalRegions/Expansion2D.h>
 #include <LocalRegions/Expansion3D.h>
@@ -22,7 +23,8 @@ namespace SU = Nektar::SolverUtils;
 
 namespace NESO::Solvers {
 
-class SOLSystem : public SU::UnsteadySystem {
+class SOLSystem
+    : public TimeEvoEqnSysBase<SU::UnsteadySystem, NeutralParticleSystem> {
 public:
   friend class MemoryManager<SOLSystem>;
 
@@ -39,24 +41,19 @@ class SOLSystem : public SU::UnsteadySystem {
   /// Name of class.
   static std::string class_name;
 
-  virtual ~SOLSystem();
-
 protected:
   SOLSystem(const LU::SessionReaderSharedPtr &session,
             const SD::MeshGraphSharedPtr &graph);
 
-  SU::AdvectionSharedPtr m_adv;
-  NektarFieldIndexMap m_field_to_index;
-  // Forcing terms
-  std::vector<SU::ForcingSharedPtr> m_forcing;
-  NekDouble m_gamma;
-  /// Names of fields that will be time integrated
-  std::vector<std::string> m_int_fld_names;
-  /// Names of fields required by the solver
-  std::vector<std::string> m_required_flds;
-  // Auxiliary object to convert variables
-  VariableConverterSharedPtr m_var_converter;
-  Array<OneD, Array<OneD, NekDouble>> m_vec_locs;
+  /// Advection object
+  SU::AdvectionSharedPtr adv_obj;
+  /// Fluid forcing / source terms
+  std::vector<SU::ForcingSharedPtr> fluid_src_terms;
+  /// Gamma value, read from config
+  NekDouble gamma;
+  /// Auxiliary object to convert variables
+  VariableConverterSharedPtr var_converter;
+  Array<OneD, Array<OneD, NekDouble>> vel_fld_indices;
 
   void do_advection(const Array<OneD, const Array<OneD, NekDouble>> &in_arr,
                     Array<OneD, Array<OneD, NekDouble>> &out_arr,
@@ -78,7 +75,7 @@ class SOLSystem : public SU::UnsteadySystem {
   get_flux_vector(const Array<OneD, const Array<OneD, NekDouble>> &physfield,
                   TensorOfArray3D<NekDouble> &flux);
 
-  NekDouble get_gamma() { return m_gamma; }
+  NekDouble get_gamma() { return this->gamma; }
 
   const Array<OneD, const Array<OneD, NekDouble>> &get_trace_norms() {
     return m_traceNormals;
@@ -89,13 +86,12 @@ class SOLSystem : public SU::UnsteadySystem {
    * (velocity field indices, in this case)
    */
   const Array<OneD, const Array<OneD, NekDouble>> &get_vec_locs() {
-    return m_vec_locs;
+    return this->vel_fld_indices;
   }
 
   void init_advection();
 
   virtual void v_InitObject(bool DeclareField) override;
-  void validate_field_list();
 };
 
 } // namespace NESO::Solvers
diff --git a/solvers/SimpleSOL/EquationSystems/SOLWithParticlesSystem.cpp b/solvers/SimpleSOL/EquationSystems/SOLWithParticlesSystem.cpp
index 9da36aea..bd06905f 100644
--- a/solvers/SimpleSOL/EquationSystems/SOLWithParticlesSystem.cpp
+++ b/solvers/SimpleSOL/EquationSystems/SOLWithParticlesSystem.cpp
@@ -11,17 +11,16 @@ std::string SOLWithParticlesSystem::class_name =
 SOLWithParticlesSystem::SOLWithParticlesSystem(
     const LU::SessionReaderSharedPtr &session,
     const SD::MeshGraphSharedPtr &graph)
-    : SOLSystem(session, graph), m_field_to_index(session->GetVariables()) {
+    : SOLSystem(session, graph) {
 
-  m_particle_sys = std::make_shared<NeutralParticleSystem>(session, graph);
-  m_required_flds.push_back("E_src");
-  m_required_flds.push_back("rho_src");
-  m_required_flds.push_back("rhou_src");
-  m_required_flds.push_back("rhov_src");
-  m_required_flds.push_back("T");
+  this->required_fld_names.push_back("E_src");
+  this->required_fld_names.push_back("rho_src");
+  this->required_fld_names.push_back("rhou_src");
+  this->required_fld_names.push_back("rhov_src");
+  this->required_fld_names.push_back("T");
 
   // mass recording diagnostic creation
-  m_diag_mass_recording_enabled =
+  this->mass_recording_enabled =
       session->DefinesParameter("mass_recording_step");
 }
 
@@ -31,13 +30,13 @@ SOLWithParticlesSystem::SOLWithParticlesSystem(
  */
 void SOLWithParticlesSystem::update_temperature() {
   // Need values of rho,rhou,[rhov],[rhow],E
-  int num_fields_for_T_calc = m_field_to_index.get_idx("E") + 1;
+  int num_fields_for_T_calc = this->field_to_index.get_idx("E") + 1;
   Array<OneD, Array<OneD, NekDouble>> phys_vals(num_fields_for_T_calc);
   for (int i = 0; i < num_fields_for_T_calc; ++i) {
     phys_vals[i] = m_fields[i]->GetPhys();
   }
-  auto temperature = m_fields[m_field_to_index.get_idx("T")];
-  m_var_converter->GetTemperature(phys_vals, temperature->UpdatePhys());
+  auto temperature = m_fields[this->field_to_index.get_idx("T")];
+  this->var_converter->GetTemperature(phys_vals, temperature->UpdatePhys());
   // Update coeffs - may not be needed?
   temperature->FwdTrans(temperature->GetPhys(), temperature->UpdateCoeffs());
 }
@@ -61,51 +60,46 @@ void SOLWithParticlesSystem::v_InitObject(bool DeclareField) {
     idx++;
   }
 
-  m_particle_sys->setup_project(
+  particle_sys->setup_project(
       m_discont_fields["rho_src"], m_discont_fields["rhou_src"],
       m_discont_fields["rhov_src"], m_discont_fields["E_src"]);
 
-  m_particle_sys->setup_evaluate_n(m_discont_fields["rho"]);
-  m_particle_sys->setup_evaluate_T(m_discont_fields["T"]);
+  particle_sys->setup_evaluate_n(m_discont_fields["rho"]);
+  particle_sys->setup_evaluate_T(m_discont_fields["T"]);
 
   m_diag_mass_recording = std::make_shared<MassRecording<MR::DisContField>>(
-      m_session, m_particle_sys, m_discont_fields["rho"]);
+      m_session, this->particle_sys, m_discont_fields["rho"]);
 }
 
-/**
- * @brief Destructor for SOLWithParticlesSystem class.
- */
-SOLWithParticlesSystem::~SOLWithParticlesSystem() { m_particle_sys->free(); }
-
 bool SOLWithParticlesSystem::v_PostIntegrate(int step) {
   // Writes a step of the particle trajectory.
   if (m_num_write_particle_steps > 0) {
     if ((step % m_num_write_particle_steps) == 0) {
-      m_particle_sys->write(step);
-      m_particle_sys->write_source_fields();
+      particle_sys->write(step);
+      particle_sys->write_source_fields();
     }
   }
 
-  if (m_diag_mass_recording_enabled) {
+  if (this->mass_recording_enabled) {
     m_diag_mass_recording->compute(step);
   }
 
-  m_solver_callback_handler.call_post_integrate(this);
+  this->solver_callback_handler.call_post_integrate(this);
   return SOLSystem::v_PostIntegrate(step);
 }
 
 bool SOLWithParticlesSystem::v_PreIntegrate(int step) {
-  m_solver_callback_handler.call_pre_integrate(this);
+  this->solver_callback_handler.call_pre_integrate(this);
 
-  if (m_diag_mass_recording_enabled) {
+  if (this->mass_recording_enabled) {
     m_diag_mass_recording->compute_initial_fluid_mass();
   }
   // Update Temperature field
   update_temperature();
   // Integrate the particle system to the requested time.
-  m_particle_sys->integrate(m_time + m_timestep, m_part_timestep);
+  particle_sys->integrate(m_time + m_timestep, m_part_timestep);
   // Project onto the source fields
-  m_particle_sys->project_source_terms();
+  particle_sys->project_source_terms();
 
   return SOLSystem::v_PreIntegrate(step);
 }
diff --git a/solvers/SimpleSOL/EquationSystems/SOLWithParticlesSystem.hpp b/solvers/SimpleSOL/EquationSystems/SOLWithParticlesSystem.hpp
index 054481ad..007fbe7e 100644
--- a/solvers/SimpleSOL/EquationSystems/SOLWithParticlesSystem.hpp
+++ b/solvers/SimpleSOL/EquationSystems/SOLWithParticlesSystem.hpp
@@ -19,7 +19,7 @@ class SOLWithParticlesSystem : public SOLSystem {
   static std::string class_name;
 
   /// Callback handler to call user defined callbacks.
-  SolverCallbackHandler<SOLWithParticlesSystem> m_solver_callback_handler;
+  SolverCallbackHandler<SOLWithParticlesSystem> solver_callback_handler;
 
   // Object that allows optional recording of stats related to mass conservation
   std::shared_ptr<MassRecording<MR::DisContField>> m_diag_mass_recording;
@@ -38,15 +38,9 @@ class SOLWithParticlesSystem : public SOLSystem {
   SOLWithParticlesSystem(const LU::SessionReaderSharedPtr &session,
                          const SD::MeshGraphSharedPtr &graph);
 
-  virtual ~SOLWithParticlesSystem();
-
 protected:
   // Flag to toggle mass conservation checking
-  bool m_diag_mass_recording_enabled;
-  // Map of field name to field index
-  NESO::NektarFieldIndexMap m_field_to_index;
-  // Particles system object
-  std::shared_ptr<NeutralParticleSystem> m_particle_sys;
+  bool mass_recording_enabled;
   // Number of particle timesteps per fluid timestep.
   int m_num_part_substeps;
   // Number of time steps between particle trajectory step writes.
diff --git a/solvers/SimpleSOL/ParticleSystems/neutral_particles.hpp b/solvers/SimpleSOL/ParticleSystems/neutral_particles.hpp
index 21f3af76..21642a3d 100644
--- a/solvers/SimpleSOL/ParticleSystems/neutral_particles.hpp
+++ b/solvers/SimpleSOL/ParticleSystems/neutral_particles.hpp
@@ -4,9 +4,9 @@
 #include <nektar_interface/function_evaluation.hpp>
 #include <nektar_interface/function_projection.hpp>
 #include <nektar_interface/particle_interface.hpp>
+#include <nektar_interface/solver_base/partsys_base.hpp>
 #include <nektar_interface/utilities.hpp>
 #include <neso_particles.hpp>
-
 #include <particle_utility/particle_initialisation_line.hpp>
 #include <particle_utility/position_distribution.hpp>
 
@@ -44,14 +44,21 @@ inline double expint_barry_approx(const double x) {
   return std::exp(-x) / (G + (1 - G) * std::exp(-(x / (1 - G)))) * logfactor;
 }
 
-class NeutralParticleSystem {
+class NeutralParticleSystem : public PartSysBase {
+  inline static ParticleSpec particle_spec{
+      ParticleProp(Sym<REAL>("POSITION"), 2, true),
+      ParticleProp(Sym<INT>("CELL_ID"), 1, true),
+      ParticleProp(Sym<INT>("PARTICLE_ID"), 2),
+      ParticleProp(Sym<REAL>("COMPUTATIONAL_WEIGHT"), 1),
+      ParticleProp(Sym<REAL>("SOURCE_DENSITY"), 1),
+      ParticleProp(Sym<REAL>("SOURCE_ENERGY"), 1),
+      ParticleProp(Sym<REAL>("SOURCE_MOMENTUM"), 2),
+      ParticleProp(Sym<REAL>("ELECTRON_DENSITY"), 1),
+      ParticleProp(Sym<REAL>("ELECTRON_TEMPERATURE"), 1),
+      ParticleProp(Sym<REAL>("MASS"), 1),
+      ParticleProp(Sym<REAL>("VELOCITY"), 3)};
+
 protected:
-  LU::SessionReaderSharedPtr session;
-  SD::MeshGraphSharedPtr graph;
-  MPI_Comm comm;
-  const double tol;
-  const int ndim = 2;
-  bool h5part_exists;
   double simulation_time;
 
   /**
@@ -120,13 +127,7 @@ class NeutralParticleSystem {
   /// Disable (implicit) copies.
   NeutralParticleSystem &operator=(NeutralParticleSystem const &a) = delete;
 
-  ~NeutralParticleSystem() {
-    std::cout << "NeutralParticleSystem dtor @" << this << std::endl;
-  }
-  /// Global number of particles in the simulation.
-  int64_t num_particles;
-  /// Average number of particles per cell (element) in the simulation.
-  int64_t num_particles_per_cell;
+  ~NeutralParticleSystem() {}
   /// Total number of particles added on this MPI rank.
   uint64_t total_num_particles_added;
   /// Mass of particles
@@ -137,22 +138,10 @@ class NeutralParticleSystem {
   double particle_number_density;
   // PARTICLE_ID value used to flag particles for removal from the simulation
   const int particle_remove_key = -1;
-  /// HMesh instance that allows particles to move over nektar++ meshes.
-  ParticleMeshInterfaceSharedPtr particle_mesh_interface;
-  /// Compute target.
-  SYCLTargetSharedPtr sycl_target;
-  /// Mapping instance to map particles into nektar++ elements.
-  std::shared_ptr<NektarGraphLocalMapper> nektar_graph_local_mapper;
-  /// NESO-Particles domain.
-  DomainSharedPtr domain;
-  /// NESO-Particles ParticleGroup containing charged particles.
-  ParticleGroupSharedPtr particle_group;
   /// Method to apply particle boundary conditions.
   std::shared_ptr<NektarCartesianPeriodic> periodic_bc;
   /// Method to map to/from nektar geometry ids to 0,N-1 used by NESO-Particles
   std::shared_ptr<CellIDTranslation> cell_id_translation;
-  /// Trajectory writer for particles.
-  std::shared_ptr<H5Part> h5part;
 
   // Factors to convert nektar units to units required by ionisation calc
   double t_to_SI;
@@ -171,41 +160,10 @@ class NeutralParticleSystem {
   NeutralParticleSystem(LU::SessionReaderSharedPtr session,
                         SD::MeshGraphSharedPtr graph,
                         MPI_Comm comm = MPI_COMM_WORLD)
-      : session(session), graph(graph), comm(comm), tol(1.0e-8),
-        h5part_exists(false), simulation_time(0.0) {
+      : PartSysBase(session, graph, particle_spec, comm), simulation_time(0.0) {
     this->total_num_particles_added = 0;
     this->debug_write_fields_count = 0;
 
-    // Read the number of requested particles per cell.
-    int tmp_int;
-    this->session->LoadParameter("num_particles_per_cell", tmp_int);
-    this->num_particles_per_cell = tmp_int;
-
-    // Reduce the global number of elements
-    const int num_elements_local = this->graph->GetNumElements();
-    int num_elements_global;
-    MPICHK(MPI_Allreduce(&num_elements_local, &num_elements_global, 1, MPI_INT,
-                         MPI_SUM, this->comm));
-
-    // compute the global number of particles
-    this->num_particles =
-        ((int64_t)num_elements_global) * this->num_particles_per_cell;
-
-    this->session->LoadParameter("num_particles_total", tmp_int);
-    if (tmp_int > -1) {
-      this->num_particles = tmp_int;
-    }
-
-    // Create interface between particles and nektar++
-    this->particle_mesh_interface =
-        std::make_shared<ParticleMeshInterface>(graph, 0, this->comm);
-    this->sycl_target =
-        std::make_shared<SYCLTarget>(0, particle_mesh_interface->get_comm());
-    this->nektar_graph_local_mapper = std::make_shared<NektarGraphLocalMapper>(
-        this->sycl_target, this->particle_mesh_interface);
-    this->domain = std::make_shared<Domain>(this->particle_mesh_interface,
-                                            this->nektar_graph_local_mapper);
-
     // Load scaling parameters from session
     double Rs, pInf, rhoInf, uInf;
     get_from_session(this->session, "GasConstant", Rs, 1.0);
@@ -241,23 +199,6 @@ class NeutralParticleSystem {
     double L_to_SI = 1;
     this->t_to_SI = L_to_SI / vel_to_SI;
 
-    // Create ParticleGroup
-    ParticleSpec particle_spec{
-        ParticleProp(Sym<REAL>("POSITION"), 2, true),
-        ParticleProp(Sym<INT>("CELL_ID"), 1, true),
-        ParticleProp(Sym<INT>("PARTICLE_ID"), 2),
-        ParticleProp(Sym<REAL>("COMPUTATIONAL_WEIGHT"), 1),
-        ParticleProp(Sym<REAL>("SOURCE_DENSITY"), 1),
-        ParticleProp(Sym<REAL>("SOURCE_ENERGY"), 1),
-        ParticleProp(Sym<REAL>("SOURCE_MOMENTUM"), 2),
-        ParticleProp(Sym<REAL>("ELECTRON_DENSITY"), 1),
-        ParticleProp(Sym<REAL>("ELECTRON_TEMPERATURE"), 1),
-        ParticleProp(Sym<REAL>("MASS"), 1),
-        ParticleProp(Sym<REAL>("VELOCITY"), 3)};
-
-    this->particle_group = std::make_shared<ParticleGroup>(
-        this->domain, particle_spec, this->sycl_target);
-
     this->particle_remover =
         std::make_shared<ParticleRemover>(this->sycl_target);
 
@@ -300,14 +241,14 @@ class NeutralParticleSystem {
     if (this->particle_number_density < 0.0) {
       this->particle_weight = 1.0;
       this->particle_number_density =
-          this->num_particles / particle_region_volume;
+          this->num_parts_tot / particle_region_volume;
     } else {
       const double number_physical_particles =
           this->particle_number_density * particle_region_volume;
       this->particle_weight =
-          (this->num_particles == 0)
+          (this->num_parts_tot == 0)
               ? 0.0
-              : number_physical_particles / this->num_particles;
+              : number_physical_particles / this->num_parts_tot;
     }
 
     // get seed from file
@@ -462,7 +403,7 @@ class NeutralParticleSystem {
 
     const int num_particles_per_line =
         add_proportion *
-        (((double)this->num_particles / ((double)total_lines)));
+        (((double)this->num_parts_tot / ((double)total_lines)));
 
     const long rank = this->sycl_target->comm_pair.rank_parent;
 
@@ -623,18 +564,6 @@ class NeutralParticleSystem {
         profile_elapsed(t0, profile_timestamp()));
   }
 
-  /**
-   *  Free the object before MPI_Finalize is called.
-   */
-  inline void free() {
-    if (this->h5part_exists) {
-      this->h5part->close();
-    }
-    this->particle_group->free();
-    this->particle_mesh_interface->free();
-    this->sycl_target->free();
-  };
-
   /**
    *  Integrate the particle system forward in time to the requested time using
    *  at most the requested time step.
diff --git a/solvers/SimpleSOL/SimpleSOL.cpp b/solvers/SimpleSOL/SimpleSOL.cpp
index c764bc6e..ee8ccfe8 100644
--- a/solvers/SimpleSOL/SimpleSOL.cpp
+++ b/solvers/SimpleSOL/SimpleSOL.cpp
@@ -7,43 +7,14 @@
 //
 ///////////////////////////////////////////////////////////////////////////////
 
-#include <LibUtilities/BasicUtils/SessionReader.h>
-#include <SolverUtils/Driver.h>
-
 #include "SimpleSOL.hpp"
-
-namespace LU = Nektar::LibUtilities;
-namespace SD = Nektar::SpatialDomains;
-namespace SU = Nektar::SolverUtils;
+#include "solvers/solver_runner.hpp"
 
 namespace NESO::Solvers {
 int run_SimpleSOL(int argc, char *argv[]) {
-  try {
-    // Create session reader.
-    auto session = LU::SessionReader::CreateInstance(argc, argv);
-
-    // Read the mesh and create a MeshGraph object.
-    auto graph = SD::MeshGraph::Read(session);
-
-    // Create driver.
-    std::string driverName;
-    session->LoadSolverInfo("Driver", driverName, "Standard");
-    auto drv =
-        SU::GetDriverFactory().CreateInstance(driverName, session, graph);
-
-    // Execute driver
-    drv->Execute();
-
-    // Finalise session
-    session->Finalise();
-  } catch (const std::runtime_error &e) {
-    std::cerr << "Error: " << e.what() << std::endl;
-    return 1;
-  } catch (const std::string &eStr) {
-    std::cerr << "Error: " << eStr << std::endl;
-    return 2;
-  }
-
+  SolverRunner solver_runner(argc, argv);
+  solver_runner.execute();
+  solver_runner.finalise();
   return 0;
 }
 
diff --git a/solvers/SimpleSOL/SimpleSOL.hpp b/solvers/SimpleSOL/SimpleSOL.hpp
index 028e3cfd..b32a296a 100644
--- a/solvers/SimpleSOL/SimpleSOL.hpp
+++ b/solvers/SimpleSOL/SimpleSOL.hpp
@@ -6,8 +6,8 @@
 // Description: Header for the SimpleSOL solver.
 //
 ///////////////////////////////////////////////////////////////////////////////
-namespace NESO {
-namespace Solvers {
+namespace NESO::Solvers {
+
 int run_SimpleSOL(int argc, char *argv[]);
-}
-} // namespace NESO
\ No newline at end of file
+
+} // namespace NESO::Solvers
\ No newline at end of file
diff --git a/test/integration/solvers/SimpleSOL/test_SimpleSOL.cpp b/test/integration/solvers/SimpleSOL/test_SimpleSOL.cpp
index 2ff84ab8..4e7a3851 100644
--- a/test/integration/solvers/SimpleSOL/test_SimpleSOL.cpp
+++ b/test/integration/solvers/SimpleSOL/test_SimpleSOL.cpp
@@ -41,9 +41,9 @@ TEST_F(SimpleSOLTest, 2DWithParticles) {
     auto equation_system = std::dynamic_pointer_cast<SOLWithParticlesSystem>(
         solver_runner.driver->GetEqu()[0]);
 
-    equation_system->m_solver_callback_handler.register_pre_integrate(
+    equation_system->solver_callback_handler.register_pre_integrate(
         callback_pre);
-    equation_system->m_solver_callback_handler.register_post_integrate(
+    equation_system->solver_callback_handler.register_post_integrate(
         callback_post);
 
     solver_runner.execute();