Apply normalization of linear rectangular equation only, if solver only supports square

MathLib/LinAlg/Eigen/EigenLinearSolver.cpp

@@ -496,23 +496,30 @@ EigenLinearSolver::EigenLinearSolver(std::string const& /*solver_name*/,
                 Eigen::SparseLU<Matrix, Eigen::COLAMDOrdering<int>>;
             solver_ = std::make_unique<
                 details::EigenDirectLinearSolver<SolverType>>();
+            can_solve_rectangular_ = false;
             return;
         }
         case EigenOption::SolverType::BiCGSTAB:
         case EigenOption::SolverType::BiCGSTABL:
         case EigenOption::SolverType::CG:
-        case EigenOption::SolverType::LeastSquareCG:
         case EigenOption::SolverType::GMRES:
         case EigenOption::SolverType::IDRS:
         case EigenOption::SolverType::IDRSTABL:
             solver_ = details::createIterativeSolver(option_.solver_type,
                                                      option_.precon_type);
+            can_solve_rectangular_ = false;
+            return;
+        case EigenOption::SolverType::LeastSquareCG:
+            solver_ = details::createIterativeSolver(option_.solver_type,
+                                                     option_.precon_type);
+            can_solve_rectangular_ = true;
             return;
         case EigenOption::SolverType::PardisoLU:
         {
 #ifdef USE_MKL
             using SolverType = Eigen::PardisoLU<EigenMatrix::RawMatrixType>;
             solver_.reset(new details::EigenDirectLinearSolver<SolverType>);
+            can_solve_rectangular_ = false;
             return;
 #else
             OGS_FATAL(

MathLib/LinAlg/Eigen/EigenLinearSolver.h

@@ -71,9 +71,13 @@ public:
                MathLib::LinearSolverBehaviour const linear_solver_behaviour =
                    MathLib::LinearSolverBehaviour::RECOMPUTE);
 
+    /// Get, if the solver can handle rectangular equation systems
+    bool canSolveRectangular() const { return can_solve_rectangular_; }
+
 protected:
     EigenOption option_;
     std::unique_ptr<EigenLinearSolverBase> solver_;
+    bool can_solve_rectangular_ = false;
     void setRestart();
     void setL();
     void setS();

MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h

@@ -51,6 +51,9 @@ public:
 
     bool solve(EigenMatrix& A, EigenVector& b, EigenVector& x);
 
+    /// Get, if the solver can handle rectangular equation systems
+    bool canSolveRectangular() const { return false; }
+
 private:
     bool solve(LisMatrix& A, LisVector& b, LisVector& x);
     std::string lis_options_;

MathLib/LinAlg/PETSc/PETScLinearSolver.cpp

@@ -40,6 +40,10 @@ PETScLinearSolver::PETScLinearSolver(std::string const& prefix,
 
     KSPSetInitialGuessNonzero(solver_, PETSC_TRUE);
     KSPSetFromOptions(solver_);  // set run-time options
+
+    KSPType ksp_type;
+    KSPGetType(solver_, &ksp_type);
+    can_solve_rectangular_ = canSolverHandleRectangular(ksp_type);
 }
 
 bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector& b, PETScVector& x)

MathLib/LinAlg/PETSc/PETScLinearSolver.h

@@ -18,6 +18,8 @@
 
 #include <petscksp.h>
 
+#include <algorithm>
+#include <array>
 #include <string>
 
 #include "PETScMatrix.h"
@@ -62,10 +64,26 @@ public:
     /// Get elapsed wall clock time.
     double getElapsedTime() const { return elapsed_ctime_; }
 
+    /// Get, if the solver can handle rectangular equation systems
+    bool canSolveRectangular() const { return can_solve_rectangular_; }
+
 private:
+    bool canSolverHandleRectangular(std::string_view ksp_type)
+    {
+        // List of KSP types that can handle rectangular matrices
+        constexpr auto rectangular_solvers = std::to_array<std::string_view>(
+            {KSPGMRES, KSPFGMRES, KSPBCGS, KSPBCGSL, KSPTFQMR, KSPCGS});
+
+        return std::ranges::any_of(rectangular_solvers,
+                                   [&](const auto& solver)
+                                   { return solver == ksp_type; });
+    }
+
     KSP solver_;  ///< Solver type.
     PC pc_;       ///< Preconditioner type.
 
+    bool can_solve_rectangular_ = false;
+
     double elapsed_ctime_ = 0.0;  ///< Clock time
 };
 

NumLib/ODESolver/NonlinearSolver.cpp

@@ -170,8 +170,18 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
         sys.getRhs(*x_prev[process_id], rhs);
 
         // Normalize the linear equation system, if required
-        if (sys.requiresNormalization())
+        if (sys.requiresNormalization() &&
+            !_linear_solver.canSolveRectangular())
+        {
             sys.getAandRhsNormalized(A, rhs);
+            WARN(
+                "The equation system is rectangular, but the current linear "
+                "solver only supports square systems. "
+                "The system will be normalized, which lead to a squared "
+                "condition number and potential numerical issues. "
+                "It is recommended to use a solver that supports rectangular "
+                "equation systems for better numerical stability.");
+        }
 
         INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());