Add additional solvers and preconditioners for contrasting material systems

Add some more non-symmetric krylov solvers (BiCGSTAB)
Add a number of new preconditioners that should help for some highly contrasting multi-material systems such as ILU, l1-Gauss Seidel, and a Chebychev preconditioner.

Author: rcarson3

src/options.toml

@@ -379,19 +379,24 @@ grain_file = "grains.txt"
         abs_tol = 1e-30
 
         # Linear solver algorithm:
-        # - "CG" = Conjugate Gradient (fastest for symmetric problems)
-        # - "GMRES" = General solver (works for any problem)
-        # - "MINRES" = Minimal Residual (good for symmetric indefinite)
-        # Use CG for ExaCMech, GMRES if unsure
+        # - "CG" = Conjugate Gradient (fastest for well-behaved symmetric problems)
+        # - "GMRES" = General solver (works for any problem, reliable but uses more memory)
+        # - "MINRES" = Minimal Residual (for symmetric problems when CG fails)
+        # - "BiCGSTAB" = Memory-efficient alternative to GMRES (faster but less robust)
+        # Use CG for typical material models, GMRES if unsure or have convergence issues
         solver = "GMRES"
 
         # Preconditioner to accelerate convergence:
         # With multi-material systems you might find better convergence
         # properties by testing out different preconditioners
         # EA/PA will automatically run with JACOBI as the full matrix is
         # not constructed
-        # - "JACOBI" = diagonal scaling (simple / fast, always works but usually increases iteration count)
-        # - "AMG" = Algebraic MultiGrid (usually lowers iteration count and can be faster for large problems)
+        # - "JACOBI" = diagonal scaling (simple/fast, works everywhere but slow convergence)
+        # - "AMG" = Algebraic MultiGrid (fewer iterations but expensive setup, can fail on some problems)
+        # - "ILU" = Incomplete factorization (good middle-ground, useful for multi-material systems)
+        # - "L1GS" = advanced smoother (can help with multi-material systems with contrasting properties)
+        # - "CHEBYSHEV" = polynomial smoother (good for problems with multiple material scales)
+        # Try ILU / L1GS / CHEBYSHEV if JACOBI convergence is too slow
         preconditioner = "JACOBI"
 
         # Output verbosity (0 = quiet, 1+ = show iterations)

src/options/option_enum.cpp

@@ -113,7 +113,8 @@ LinearSolverType string_to_linear_solver_type(const std::string& str) {
         {"CG", LinearSolverType::CG},
         {"PCG", LinearSolverType::CG},
         {"GMRES", LinearSolverType::GMRES},
-        {"MINRES", LinearSolverType::MINRES}
+        {"MINRES", LinearSolverType::MINRES},
+        {"BICGSTAB", LinearSolverType::BICGSTAB}
     };
 
     return string_to_enum(str, mapping, LinearSolverType::NOTYPE, "linear solver");
@@ -141,7 +142,10 @@ NonlinearSolverType string_to_nonlinear_solver_type(const std::string& str) {
 PreconditionerType string_to_preconditioner_type(const std::string& str) {
     static const std::map<std::string, PreconditionerType> mapping = {
         {"JACOBI", PreconditionerType::JACOBI},
-        {"AMG", PreconditionerType::AMG}
+        {"AMG", PreconditionerType::AMG},
+        {"ILU", PreconditionerType::ILU},
+        {"L1GS", PreconditionerType::L1GS},
+        {"CHEBYSHEV", PreconditionerType::CHEBYSHEV},
     };
 
     return string_to_enum(str, mapping, PreconditionerType::NOTYPE, "preconditioner");

src/options/option_parser_v2.cpp

@@ -597,13 +597,17 @@ void ExaOptions::print_solver_options() const {
         case LinearSolverType::CG: std::cout << "Conjugate Gradient\n"; break;
         case LinearSolverType::GMRES: std::cout << "GMRES\n"; break;
         case LinearSolverType::MINRES: std::cout << "MINRES\n"; break;
+        case LinearSolverType::BICGSTAB: std::cout << "BiCGSTAB\n"; break;
         default: std::cout << "Unknown\n"; break;
     }
 
     std::cout << "    Preconditioner: ";
     switch (solvers.linear_solver.preconditioner) {
         case PreconditionerType::JACOBI: std::cout << "Jacobi\n"; break;
         case PreconditionerType::AMG: std::cout << "AMG\n"; break;
+        case PreconditionerType::ILU: std::cout << "ILU\n"; break;
+        case PreconditionerType::L1GS: std::cout << "L1GS\n"; break;
+        case PreconditionerType::CHEBYSHEV: std::cout << "CHEBYSHEV\n"; break;
         default: std::cout << "Unknown\n"; break;
     }
 

src/options/option_parser_v2.hpp

@@ -89,6 +89,7 @@ enum class LinearSolverType {
     CG,     /**< Conjugate Gradient solver */
     GMRES,  /**< Generalized Minimal Residual solver */
     MINRES, /**< Minimal Residual solver */
+    BICGSTAB, /**< BiCGSTAB Solver */
     NOTYPE  /**< Uninitialized or invalid linear solver type */
 };
 
@@ -106,7 +107,10 @@ enum class NonlinearSolverType {
  */
 enum class PreconditionerType { 
     JACOBI, /**< Jacobi preconditioner */
-    AMG,    /**< Algebraic multigrid preconditioner */
+    AMG,    /**< Algebraic multigrid preconditioner (Full assembly only) */
+    ILU,    /**< Incomplete LU factorization preconditioner (Full assembly only) */
+    L1GS,   /**< l1-scaled block Gauss-Seidel/SSOR preconditioner (Full assembly only) */
+    CHEBYSHEV,   /**< Chebyshev preconditioner (Full assembly only) */
     NOTYPE  /**< Uninitialized or invalid preconditioner type */
 };
 
@@ -580,7 +584,7 @@ struct LinearSolverOptions {
     /**
      * @brief Preconditioner type for linear solver acceleration
      */
-    PreconditionerType preconditioner = PreconditionerType::JACOBI;
+    PreconditionerType preconditioner = PreconditionerType::AMG;
 
     /**
      * @brief Absolute convergence tolerance for linear solver

src/system_driver.cpp

@@ -253,7 +253,7 @@ SystemDriver::SystemDriver(std::shared_ptr<SimulationState> sim_state)
       J_prec = mech_operator->GetPAPreconditioner();
    }
    else {
-      if (linear_solvers.solver_type == LinearSolverType::GMRES || linear_solvers.solver_type == LinearSolverType::CG) {
+      if (linear_solvers.preconditioner == PreconditionerType::AMG) {
          auto prec_amg = std::make_shared<mfem::HypreBoomerAMG>();
          HYPRE_Solver h_amg = (HYPRE_Solver) * prec_amg;
          HYPRE_Real st_val = 0.90;
@@ -280,8 +280,19 @@ SystemDriver::SystemDriver(std::shared_ptr<SimulationState> sim_state)
 
          prec_amg->SetPrintLevel(linear_solvers.print_level);
          J_prec = prec_amg;
-      }
-      else {
+      } else if (linear_solvers.preconditioner == PreconditionerType::ILU) {
+         auto J_hypreEuclid = std::make_shared<mfem::HypreEuclid>(fe_space->GetComm());
+         J_prec = J_hypreEuclid;
+      } else if (linear_solvers.preconditioner == PreconditionerType::L1GS) {
+         auto J_hypreSmoother = std::make_shared<mfem::HypreSmoother>();
+         J_hypreSmoother->SetType(mfem::HypreSmoother::l1GS);
+         J_hypreSmoother->SetPositiveDiagonal(true);
+         J_prec = J_hypreSmoother;
+      } else if (linear_solvers.preconditioner == PreconditionerType::CHEBYSHEV) {
+         auto J_hypreSmoother = std::make_shared<mfem::HypreSmoother>();
+         J_hypreSmoother->SetType(mfem::HypreSmoother::Chebyshev);
+         J_prec = J_hypreSmoother;
+      } else {
          auto J_hypreSmoother = std::make_shared<mfem::HypreSmoother>();
          J_hypreSmoother->SetType(mfem::HypreSmoother::l1Jacobi);
          J_hypreSmoother->SetPositiveDiagonal(true);
@@ -295,17 +306,20 @@ SystemDriver::SystemDriver(std::shared_ptr<SimulationState> sim_state)
    else if (linear_solvers.solver_type == LinearSolverType::CG) {
       J_solver = std::make_shared<mfem::CGSolver>(fe_space->GetComm());
    }
+   else if (linear_solvers.solver_type == LinearSolverType::BICGSTAB) {
+      J_solver = std::make_shared<mfem::BiCGSTABSolver>(fe_space->GetComm());
+   }
    else {
       J_solver = std::make_shared<mfem::MINRESSolver>(fe_space->GetComm());
    }
 
-      // The relative tolerance should be at this point or smaller
-      J_solver->SetRelTol(linear_solvers.rel_tol);
-      // The absolute tolerance could probably get even smaller then this
-      J_solver->SetAbsTol(linear_solvers.abs_tol);
-      J_solver->SetMaxIter(linear_solvers.max_iter);
-      J_solver->SetPrintLevel(linear_solvers.print_level);
-      J_solver->SetPreconditioner(*J_prec);
+   // The relative tolerance should be at this point or smaller
+   J_solver->SetRelTol(linear_solvers.rel_tol);
+   // The absolute tolerance could probably get even smaller then this
+   J_solver->SetAbsTol(linear_solvers.abs_tol);
+   J_solver->SetMaxIter(linear_solvers.max_iter);
+   J_solver->SetPrintLevel(linear_solvers.print_level);
+   J_solver->SetPreconditioner(*J_prec);
 
    auto nonlinear_solver = options.solvers.nonlinear_solver;
    newton_iter = nonlinear_solver.iter;