Iterative method and solver updates

framework/math/linear_solver/linear_eigen_solver.h

@@ -0,0 +1,40 @@
+// SPDX-FileCopyrightText: 2025 The OpenSn Authors <https://open-sn.github.io/opensn/>
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include "framework/math/linear_solver/linear_solver.h"
+#include "framework/math/linear_solver/linear_solver_context.h"
+
+namespace opensn
+{
+
+class LinearEigenSolver : public LinearSolver
+{
+public:
+  enum class IterativeMethod : int
+  {
+    POWER
+  };
+
+  LinearEigenSolver(IterativeMethod method, std::shared_ptr<LinearEigenContext> ctx)
+    : LinearSolver(ctx), method_(method)
+  {
+  }
+
+  std::string GetIterativeMethodName()
+  {
+    switch (method_)
+    {
+      case IterativeMethod::POWER:
+        return "POWER ITERATION";
+      default:
+        return "UNKNOWN ITERATIVE METHOD";
+    }
+  }
+
+protected:
+  IterativeMethod method_;
+};
+
+} // namespace opensn

framework/math/linear_solver/linear_solver.h

@@ -14,36 +14,23 @@ struct LinearSolverContext;
 class LinearSolver
 {
 public:
-  enum class IterativeMethod : int
-  {
-    NONE = 0,
-    CLASSIC_RICHARDSON = 1, ///< Classic Richardson (source iteration)
-    PETSC_RICHARDSON = 3,   ///< PETSc Richardson iteration
-    PETSC_GMRES = 2,        ///< PETSc GMRES iterative algorithm
-    PETSC_BICGSTAB = 4,     ///< PETSc BiCGStab iterative algorithm
-  };
-
-  LinearSolver(IterativeMethod iterative_method, std::shared_ptr<LinearSolverContext> context_ptr)
-    : iterative_method_(iterative_method), context_ptr_(context_ptr)
+  explicit LinearSolver(std::shared_ptr<LinearSolverContext> context_ptr)
+    : context_ptr_(context_ptr)
   {
   }
 
-  virtual ~LinearSolver() {}
-
-  std::shared_ptr<LinearSolverContext> GetContext() { return context_ptr_; }
+  virtual ~LinearSolver() = default;
 
   /// Set up the linaer solver
   virtual void Setup() {}
 
   /// Solve the system
   virtual void Solve() = 0;
 
+  std::shared_ptr<LinearSolverContext> GetContext() { return context_ptr_; }
+
 protected:
-  const IterativeMethod iterative_method_;
   std::shared_ptr<LinearSolverContext> context_ptr_;
-
-public:
-  static std::string IterativeMethodName(IterativeMethod iterative_method);
 };
 
 } // namespace opensn

framework/math/linear_solver/linear_solver_context.h

@@ -17,15 +17,23 @@ enum class ResidualScaleType
 };
 
 struct LinearSolverContext
+{
+  virtual int MatrixAction(Mat& matrix, Vec& vector, Vec& action) { return 0; }
+
+  virtual ~LinearSolverContext() = default;
+};
+
+struct LinearSystemContext : public LinearSolverContext
 {
   double rhs_norm = 0.0;
   double rhs_preconditioned_norm = 0.0;
   double custom_residual_scale = 1.0;
   ResidualScaleType residual_scale_type = ResidualScaleType::NONE;
+};
 
-  virtual int MatrixAction(Mat& matrix, Vec& vector, Vec& action) { return 0; }
-
-  virtual ~LinearSolverContext() = default;
+struct LinearEigenContext : public LinearSolverContext
+{
+  double eigenvalue = 0.0;
 };
 
 } // namespace opensn

framework/math/linear_solver/linear_system_solver.h

@@ -0,0 +1,52 @@
+// SPDX-FileCopyrightText: 2025 The OpenSn Authors <https://open-sn.github.io/opensn/>
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include "framework/math/linear_solver/linear_solver.h"
+#include "framework/math/linear_solver/linear_solver_context.h"
+
+namespace opensn
+{
+
+class LinearSystemSolver : public LinearSolver
+{
+public:
+  enum class IterativeMethod : int
+  {
+    NONE,
+    CLASSIC_RICHARDSON,
+    PETSC_RICHARDSON,
+    PETSC_GMRES,
+    PETSC_BICGSTAB
+  };
+
+  LinearSystemSolver(IterativeMethod method, std::shared_ptr<LinearSystemContext> ctx)
+    : LinearSolver(ctx), method_(method)
+  {
+  }
+
+  std::string GetIterativeMethodName()
+  {
+    switch (method_)
+    {
+      case IterativeMethod::NONE:
+        return "NONE";
+      case IterativeMethod::CLASSIC_RICHARDSON:
+        return "CLASSIC_RICHARDSON";
+      case IterativeMethod::PETSC_RICHARDSON:
+        return "PETSC_RICHARDSON";
+      case IterativeMethod::PETSC_GMRES:
+        return "PETSC_GMRES";
+      case IterativeMethod::PETSC_BICGSTAB:
+        return "PETSC_BICGSTAB";
+      default:
+        return "UNKNOWN ITERATIVE METHOD";
+    }
+  }
+
+protected:
+  IterativeMethod method_;
+};
+
+} // namespace opensn

framework/math/linear_solver/petsc_linear_solver.cc → framework/math/linear_solver/petsc_linear_system_solver.cc

@@ -1,15 +1,15 @@
 // SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
 // SPDX-License-Identifier: MIT
 
-#include "framework/math/linear_solver/petsc_linear_solver.h"
+#include "framework/math/linear_solver/petsc_linear_system_solver.h"
 #include "framework/runtime.h"
 
 namespace opensn
 {
 
-PETScLinearSolver::PETScLinearSolver(IterativeMethod iterative_method,
-                                     std::shared_ptr<LinearSolverContext> context_ptr)
-  : LinearSolver(iterative_method, context_ptr),
+PETScLinearSolver::PETScLinearSolver(IterativeMethod method,
+                                     std::shared_ptr<LinearSystemContext> context_ptr)
+  : LinearSystemSolver(method, context_ptr),
     A_(nullptr),
     b_(nullptr),
     x_(nullptr),
@@ -82,34 +82,25 @@ PETScLinearSolver::Setup()
     return;
 
   PreSetupCallback();
-
   KSPCreate(opensn::mpi_comm, &ksp_);
-
-  const auto petsc_iterative_method = PETScIterativeMethodName(iterative_method_);
+  const auto petsc_iterative_method = PETScIterativeMethodName();
   KSPSetType(ksp_, petsc_iterative_method.c_str());
-
   ApplyToleranceOptions();
-
-  if (iterative_method_ == IterativeMethod::PETSC_GMRES)
+  if (method_ == IterativeMethod::PETSC_GMRES)
   {
     KSPGMRESSetRestart(ksp_, tolerance_options.gmres_restart_interval);
     KSPGMRESSetBreakdownTolerance(ksp_, tolerance_options.gmres_breakdown_tolerance);
   }
-
   KSPSetInitialGuessNonzero(ksp_, PETSC_FALSE);
-
   SetOptions();
-
   SetSolverContext();
   SetConvergenceTest();
   SetMonitor();
-
   SetSystemSize();
   SetSystem();
-
   SetPreconditioner();
-
   PostSetupCallback();
+
   system_set_ = true;
 }
 
@@ -135,9 +126,9 @@ PETScLinearSolver::Solve()
 }
 
 std::string
-PETScLinearSolver::PETScIterativeMethodName(opensn::LinearSolver::IterativeMethod iterative_method)
+PETScLinearSolver::PETScIterativeMethodName()
 {
-  switch (iterative_method)
+  switch (method_)
   {
     case IterativeMethod::NONE:
       return "preonly";
@@ -148,14 +139,14 @@ PETScLinearSolver::PETScIterativeMethodName(opensn::LinearSolver::IterativeMetho
     case IterativeMethod::PETSC_BICGSTAB:
       return "bcgs";
     default:
-      throw std::runtime_error("Cannot get a PETSc option name for a non-PETSc iterative method.");
+      throw std::runtime_error("Unsupported PETSc iterative method.");
   }
 }
 
 int
 PETScLinearSolver::LinearSolverMatrixAction(Mat matrix, Vec vector, Vec action)
 {
-  LinearSolverContext* context;
+  LinearSystemContext* context;
   MatShellGetContext(matrix, &context);
 
   context->MatrixAction(matrix, vector, action);

framework/math/linear_solver/petsc_linear_solver.h → framework/math/linear_solver/petsc_linear_system_solver.h

@@ -3,7 +3,7 @@
 
 #pragma once
 
-#include "framework/math/linear_solver/linear_solver.h"
+#include "framework/math/linear_solver/linear_system_solver.h"
 #include "framework/math/linear_solver/linear_solver_context.h"
 #include <string>
 #include <utility>
@@ -13,7 +13,7 @@
 namespace opensn
 {
 
-class PETScLinearSolver : public LinearSolver
+class PETScLinearSolver : public LinearSystemSolver
 {
 public:
   struct ToleranceOptions
@@ -26,8 +26,7 @@ class PETScLinearSolver : public LinearSolver
     double gmres_breakdown_tolerance = 1.0e6;
   } tolerance_options;
 
-  PETScLinearSolver(IterativeMethod iterative_method,
-                    std::shared_ptr<LinearSolverContext> context_ptr);
+  PETScLinearSolver(IterativeMethod method, std::shared_ptr<LinearSystemContext> context_ptr);
 
   virtual ~PETScLinearSolver();
 
@@ -72,12 +71,10 @@ class PETScLinearSolver : public LinearSolver
 private:
   bool system_set_;
   bool suppress_kspsolve_;
+  std::string PETScIterativeMethodName();
 
 protected:
   static int LinearSolverMatrixAction(Mat matrix, Vec vector, Vec action);
-
-private:
-  static std::string PETScIterativeMethodName(IterativeMethod iterative_method);
 };
 
 } // namespace opensn

modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cc

@@ -311,8 +311,11 @@ DiscreteOrdinatesProblem::InitializeWGSSolvers()
       options_.verbose_inner_iterations,
       sweep_chunk);
 
-    if (groupset.iterative_method == LinearSolver::IterativeMethod::CLASSIC_RICHARDSON)
-      wgs_solvers_.push_back(std::make_shared<ClassicRichardson>(sweep_wgs_context_ptr));
+    if (groupset.iterative_method == LinearSystemSolver::IterativeMethod::CLASSIC_RICHARDSON)
+    {
+      wgs_solvers_.push_back(std::make_shared<ClassicRichardson>(
+        sweep_wgs_context_ptr, options_.verbose_inner_iterations));
+    }
     else
       wgs_solvers_.push_back(std::make_shared<WGSLinearSolver>(sweep_wgs_context_ptr));
   }

modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/sweep_wgs_context.cc

@@ -3,12 +3,12 @@
 
 #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/sweep_wgs_context.h"
 #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h"
-#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_vecops.h"
 #include "modules/linear_boltzmann_solvers/lbs_problem/preconditioning/lbs_shell_operations.h"
-#include "framework/runtime.h"
+#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_vecops.h"
 #include "framework/logging/log.h"
-#include <petscksp.h>
+#include "framework/runtime.h"
 #include "caliper/cali.h"
+#include <petscksp.h>
 #include <iomanip>
 
 namespace opensn
@@ -33,20 +33,6 @@ SweepWGSContext::SweepWGSContext(DiscreteOrdinatesProblem& do_problem,
 {
 }
 
-void
-SweepWGSContext::PreSetupCallback()
-{
-  CALI_CXX_MARK_SCOPE("SweepWGSContext::PreSetupCallback");
-
-  if (log_info)
-    log.Log() << "\n\n"
-              << "********** Solving groupset " << groupset.id << " with "
-              << LinearSolver::IterativeMethodName(groupset.iterative_method) << ".\n\n"
-              << "Quadrature number of angles: " << groupset.quadrature->abscissae.size() << "\n"
-              << "Groups " << groupset.groups.front().id << " " << groupset.groups.back().id
-              << "\n\n";
-}
-
 void
 SweepWGSContext::SetPreconditioner(KSP& solver)
 {
@@ -130,9 +116,9 @@ SweepWGSContext::PostSolveCallback()
   // Perform final sweep with converged phi and delayed psi dofs. This step is necessary for
   // Krylov methods to recover the actual solution (this includes all of the PETSc methods
   // currently used in OpenSn).
-  if (groupset.iterative_method == LinearSolver::IterativeMethod::PETSC_GMRES or
-      groupset.iterative_method == LinearSolver::IterativeMethod::PETSC_BICGSTAB or
-      (groupset.iterative_method == LinearSolver::IterativeMethod::PETSC_RICHARDSON and
+  if (groupset.iterative_method == LinearSystemSolver::IterativeMethod::PETSC_GMRES or
+      groupset.iterative_method == LinearSystemSolver::IterativeMethod::PETSC_BICGSTAB or
+      (groupset.iterative_method == LinearSystemSolver::IterativeMethod::PETSC_RICHARDSON and
        groupset.max_iterations > 1))
   {
     const auto scope = lhs_src_scope | rhs_src_scope;

modules/linear_boltzmann_solvers/discrete_ordinates_problem/iterative_methods/sweep_wgs_context.h

@@ -21,8 +21,6 @@ struct SweepWGSContext : public WGSContext
                   bool log_info,
                   std::shared_ptr<SweepChunk> sweep_chunk);
 
-  void PreSetupCallback() override;
-
   void SetPreconditioner(KSP& solver) override;
 
   std::pair<int64_t, int64_t> GetSystemSize() override;