Merge branch 'LSCG' into 'master'

Add LeastSquareCG-Solver from Eigen See merge request ogs/ogs!5150
ufz · Nov 29, 2024 · ffc324f · ffc324f
2 parents 8da24ed + 51c29bd
commit ffc324f
Show file tree

Hide file tree

Showing 13 changed files with 195 additions and 10 deletions.
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
@@ -388,6 +388,9 @@ std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
         case EigenOption::PreconType::DIAGONAL:
             return createIterativeSolver<
                 Solver, Eigen::DiagonalPreconditioner<double>>();
+        case EigenOption::PreconType::LeastSquareDIAGONAL:
+            return createIterativeSolver<
+                Solver, Eigen::LeastSquareDiagonalPreconditioner<double>>();
         case EigenOption::PreconType::ILUT:
             // TODO for this preconditioner further options can be passed.
             // see
@@ -429,6 +432,11 @@ std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
         {
             return createIterativeSolver<EigenCGSolver>(precon_type);
         }
+        case EigenOption::SolverType::LeastSquareCG:
+        {
+            return createIterativeSolver<Eigen::LeastSquaresConjugateGradient>(
+                precon_type);
+        }
         case EigenOption::SolverType::GMRES:
         {
 #ifdef USE_EIGEN_UNSUPPORTED
@@ -488,6 +496,7 @@ 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:
@@ -498,12 +507,19 @@ EigenLinearSolver::EigenLinearSolver(std::string const& /*solver_name*/,
         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(

diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.h b/MathLib/LinAlg/Eigen/EigenLinearSolver.h
@@ -71,9 +71,13 @@ class EigenLinearSolver final
                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();

diff --git a/MathLib/LinAlg/Eigen/EigenOption.cpp b/MathLib/LinAlg/Eigen/EigenOption.cpp
@@ -38,6 +38,10 @@ EigenOption::SolverType EigenOption::getSolverType(
     {
         return SolverType::CG;
     }
+    if (solver_name == "LeastSquareCG")
+    {
+        return SolverType::LeastSquareCG;
+    }
     if (solver_name == "BiCGSTAB")
     {
         return SolverType::BiCGSTAB;
@@ -81,6 +85,10 @@ EigenOption::PreconType EigenOption::getPreconType(
     {
         return PreconType::DIAGONAL;
     }
+    if (precon_name == "LeastSquareDIAGONAL")
+    {
+        return PreconType::LeastSquareDIAGONAL;
+    }
     if (precon_name == "ILUT")
     {
         return PreconType::ILUT;
@@ -95,6 +103,8 @@ std::string EigenOption::getSolverName(SolverType const solver_type)
     {
         case SolverType::CG:
             return "CG";
+        case SolverType::LeastSquareCG:
+            return "LeastSquareCG";
         case SolverType::BiCGSTAB:
             return "BiCGSTAB";
         case SolverType::BiCGSTABL:
@@ -121,6 +131,8 @@ std::string EigenOption::getPreconName(PreconType const precon_type)
             return "NONE";
         case PreconType::DIAGONAL:
             return "DIAGONAL";
+        case PreconType::LeastSquareDIAGONAL:
+            return "LeastSquareDIAGONAL";
         case PreconType::ILUT:
             return "ILUT";
     }

diff --git a/MathLib/LinAlg/Eigen/EigenOption.h b/MathLib/LinAlg/Eigen/EigenOption.h
@@ -21,6 +21,7 @@ struct EigenOption final
     enum class SolverType : short
     {
         CG,
+        LeastSquareCG,
         BiCGSTAB,
         BiCGSTABL,
         IDRS,
@@ -35,6 +36,7 @@ struct EigenOption final
     {
         NONE,
         DIAGONAL,
+        LeastSquareDIAGONAL,
         ILUT
     };
 

diff --git a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h b/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h
@@ -51,6 +51,9 @@ class EigenLisLinearSolver final
 
     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_;

diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp b/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)

diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.h b/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 @@ class PETScLinearSolver
     /// 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
 };
 

diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/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());
 

diff --git a/ProcessLib/HeatTransportBHE/Tests.cmake b/ProcessLib/HeatTransportBHE/Tests.cmake
@@ -34,12 +34,26 @@ AddTest(
     WRAPPER time
     TESTER vtkdiff
     REQUIREMENTS NOT OGS_USE_MPI
-    RUNTIME 20
+    RUNTIME 3
     DIFF_DATA
     sandwich_ts_10_t_600.000000.vtu sandwich_algebraic_bc_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 0 1e-6
     sandwich_ts_10_t_600.000000.vtu sandwich_algebraic_bc_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 5e-9
 )
 
+AddTest(
+    NAME HeatTransportBHE_1U_3D_bhe_sandwich_algebraicBC_LSCG
+    PATH Parabolic/T/3D_BHE_Sandwich
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS sandwich_algebraicBC_LSCG.xml
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    RUNTIME 3
+    DIFF_DATA
+    sandwich_ts_10_t_600.000000.vtu sandwich_algebraic_bc_LSCG_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 0 1e-6
+    sandwich_ts_10_t_600.000000.vtu sandwich_algebraic_bc_LSCG_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 5e-9
+)
+
 AddTest(
     NAME HeatTransportBHE_1U_3D_bhe_sandwich_fixed_power
     PATH Parabolic/T/3D_BHE_Sandwich
@@ -57,7 +71,7 @@ AddTest(
 AddTest(
     NAME HeatTransportBHE_1U_3D_bhe_sandwich_fixed_power_algebraicBC
     PATH Parabolic/T/3D_BHE_Sandwich
-    RUNTIME 100
+    RUNTIME 3
     EXECUTABLE ogs
     EXECUTABLE_ARGS sandwich_fixed_power_algebraicBC.xml
     WRAPPER time
@@ -68,6 +82,20 @@ AddTest(
     sandwich_fixed_power_ts_10_t_600.000000.vtu sandwich_fixed_power_algebraic_bc_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 1e-6
 )
 
+AddTest(
+    NAME HeatTransportBHE_1U_3D_bhe_sandwich_fixed_power_algebraicBC_LSCG
+    PATH Parabolic/T/3D_BHE_Sandwich
+    RUNTIME 3
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS sandwich_fixed_power_algebraicBC_LSCG.xml
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    sandwich_fixed_power_ts_10_t_600.000000.vtu sandwich_fixed_power_algebraic_bc_LSCG_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 0 5e-3
+    sandwich_fixed_power_ts_10_t_600.000000.vtu sandwich_fixed_power_algebraic_bc_LSCG_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 5e-6
+)
+
 AddTest(
     NAME HeatTransportBHE_1U_3D_beier_sandbox
     PATH Parabolic/T/3D_Beier_sandbox
@@ -96,6 +124,20 @@ AddTest(
     beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_newton_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 1e-13
 )
 
+AddTest(
+    NAME HeatTransportBHE_1U_3D_MassLumping
+    PATH Parabolic/T/3D_Beier_sandbox
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS beier_sandbox_MassLumping.xml
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    RUNTIME 20
+    DIFF_DATA
+    beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_mass_lumping_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 0 1e-6
+    beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_mass_lumping_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 1e-4
+)
+
 AddTest(
     NAME HeatTransportBHE_1U_3D_beier_sandbox_binary_curve
     PATH Parabolic/T/3D_Beier_sandbox
@@ -118,24 +160,24 @@ AddTest(
     WRAPPER time
     TESTER vtkdiff
     REQUIREMENTS NOT OGS_USE_MPI
-    RUNTIME 20
+    RUNTIME 3
     DIFF_DATA
     beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_algebraic_bc_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 0 5e-7
     beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_algebraic_bc_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 5e-10
 )
 
 AddTest(
-    NAME HeatTransportBHE_1U_3D_MassLumping
+    NAME HeatTransportBHE_1U_3D_beier_sandbox_algebraicBC_LSCG
     PATH Parabolic/T/3D_Beier_sandbox
     EXECUTABLE ogs
-    EXECUTABLE_ARGS beier_sandbox_MassLumping.xml
+    EXECUTABLE_ARGS beier_sandbox_algebraicBC_LSCG.xml
     WRAPPER time
     TESTER vtkdiff
     REQUIREMENTS NOT OGS_USE_MPI
-    RUNTIME 20
+    RUNTIME 3
     DIFF_DATA
-    beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_mass_lumping_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 0 1e-6
-    beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_mass_lumping_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 1e-4
+    beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_algebraic_bc_LSCG_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 0 5e-7
+    beier_sandbox_ts_10_t_600.000000.vtu beier_sandbox_algebraic_bc_LSCG_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 5e-10
 )
 
 AddTest(
@@ -155,7 +197,7 @@ AddTest(
 AddTest(
     NAME HeatTransportBHE_1U_beier_sandbox_fixed_power_constant_flow_algebraicBC
     PATH Parabolic/T/3D_Beier_sandbox
-    RUNTIME 220
+    RUNTIME 3
     EXECUTABLE ogs
     EXECUTABLE_ARGS fixed_power_constant_flow_algebraicBC.xml
     WRAPPER time
@@ -166,6 +208,20 @@ AddTest(
     fixed_power_constant_flow_ts_10_t_600.000000.vtu fixed_power_constant_flow_algebraic_bc_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 5e-9
 )
 
+AddTest(
+    NAME HeatTransportBHE_1U_beier_sandbox_fixed_power_constant_flow_algebraicBC_LSCG
+    PATH Parabolic/T/3D_Beier_sandbox
+    RUNTIME 3
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS fixed_power_constant_flow_algebraicBC_LSCG.xml
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    fixed_power_constant_flow_ts_10_t_600.000000.vtu fixed_power_constant_flow_algebraic_bc_LSCG_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 0 1e-4
+    fixed_power_constant_flow_ts_10_t_600.000000.vtu fixed_power_constant_flow_algebraic_bc_LSCG_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 5e-9
+)
+
 AddTest(
     NAME HeatTransportBHE_coaxial_pipe_3D_deep_BHE_CXA
     PATH Parabolic/T/3D_deep_BHE

diff --git a/Tests/Data/Parabolic/T/3D_BHE_Sandwich/sandwich_algebraicBC_LSCG.xml b/Tests/Data/Parabolic/T/3D_BHE_Sandwich/sandwich_algebraicBC_LSCG.xml
@@ -0,0 +1,15 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProjectDiff base_file="sandwich.prj">
+    <add sel="/*/processes/process">
+        <use_algebraic_bc>true</use_algebraic_bc>
+    </add>
+    <add sel="/*/processes/process">
+        <weighting_factor>100</weighting_factor>
+    </add>
+    <add sel="/*/processes/process">
+        <linear>true</linear>
+    </add>
+    <replace sel="/*/time_loop/output/prefix/text()">sandwich_algebraic_bc_LSCG</replace>
+    <replace sel="/*/linear_solvers/linear_solver/eigen/solver_type/text()">LeastSquareCG</replace>
+    <replace sel="/*/linear_solvers/linear_solver/eigen/precon_type/text()">LeastSquareDIAGONAL</replace>
+</OpenGeoSysProjectDiff>
diff --git a/Tests/Data/Parabolic/T/3D_BHE_Sandwich/sandwich_fixed_power_algebraicBC_LSCG.xml b/Tests/Data/Parabolic/T/3D_BHE_Sandwich/sandwich_fixed_power_algebraicBC_LSCG.xml
@@ -0,0 +1,15 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProjectDiff base_file="sandwich_fixed_power.prj">
+    <add sel="/*/processes/process">
+        <use_algebraic_bc>true</use_algebraic_bc>
+    </add>
+    <add sel="/*/processes/process">
+        <weighting_factor>100</weighting_factor>
+    </add>
+    <add sel="/*/processes/process">
+        <linear>true</linear>
+    </add>
+    <replace sel="/*/time_loop/output/prefix/text()">sandwich_fixed_power_algebraic_bc_LSCG</replace>
+    <replace sel="/*/linear_solvers/linear_solver/eigen/solver_type/text()">LeastSquareCG</replace>
+    <replace sel="/*/linear_solvers/linear_solver/eigen/precon_type/text()">LeastSquareDIAGONAL</replace>
+</OpenGeoSysProjectDiff>
diff --git a/Tests/Data/Parabolic/T/3D_Beier_sandbox/beier_sandbox_algebraicBC_LSCG.xml b/Tests/Data/Parabolic/T/3D_Beier_sandbox/beier_sandbox_algebraicBC_LSCG.xml
@@ -0,0 +1,15 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProjectDiff base_file="beier_sandbox.prj">
+    <add sel="/*/processes/process">
+        <use_algebraic_bc>true</use_algebraic_bc>
+    </add>
+    <add sel="/*/processes/process">
+        <weighting_factor>100</weighting_factor>
+    </add>
+    <add sel="/*/processes/process">
+        <linear>true</linear>
+    </add>
+    <replace sel="/*/time_loop/output/prefix/text()">beier_sandbox_algebraic_bc_LSCG</replace>
+    <replace sel="/*/linear_solvers/linear_solver/eigen/solver_type/text()">LeastSquareCG</replace>
+    <replace sel="/*/linear_solvers/linear_solver/eigen/precon_type/text()">LeastSquareDIAGONAL</replace>
+</OpenGeoSysProjectDiff>