Merge branch 'CG_option' into 'master'

[ML/Eigen] set triangular matrix type for CG solver to enable multithreading

See merge request ogs/ogs!5187

Documentation/ProjectFile/prj/linear_solvers/linear_solver/eigen/t_triangular_matrix.md

+Set the triangular matrix type for the CG solver.
+Could be Lower, Upper, or LowerUpper. Default is Lower.
+OpenMP parallelization requires LowerUpper.
+See [Eigen Documentation](https://eigen.tuxfamily.org/dox/classEigen_1_1ConjugateGradient.html) for details.

MathLib/LinAlg/Eigen/EigenLinearSolver.cpp

@@ -403,10 +403,18 @@ std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
 }
 
 template <typename Mat, typename Precon>
-using EigenCGSolver = Eigen::ConjugateGradient<Mat, Eigen::Lower, Precon>;
+using EigenCGSolverL = Eigen::ConjugateGradient<Mat, Eigen::Lower, Precon>;
+
+template <typename Mat, typename Precon>
+using EigenCGSolverU = Eigen::ConjugateGradient<Mat, Eigen::Upper, Precon>;
+
+template <typename Mat, typename Precon>
+using EigenCGSolverLU =
+    Eigen::ConjugateGradient<Mat, Eigen::Lower | Eigen::Upper, Precon>;
 
 std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
-    EigenOption::SolverType solver_type, EigenOption::PreconType precon_type)
+    EigenOption::SolverType solver_type, EigenOption::PreconType precon_type,
+    EigenOption::TriangularMatrixType triangular_matrix_type)
 {
     switch (solver_type)
     {
@@ -430,7 +438,15 @@ std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
         }
         case EigenOption::SolverType::CG:
         {
-            return createIterativeSolver<EigenCGSolver>(precon_type);
+            switch (triangular_matrix_type)
+            {
+                case EigenOption::TriangularMatrixType::Upper:
+                    return createIterativeSolver<EigenCGSolverU>(precon_type);
+                case EigenOption::TriangularMatrixType::LowerUpper:
+                    return createIterativeSolver<EigenCGSolverLU>(precon_type);
+                default:
+                    return createIterativeSolver<EigenCGSolverL>(precon_type);
+            }
         }
         case EigenOption::SolverType::LeastSquareCG:
         {
@@ -505,13 +521,17 @@ EigenLinearSolver::EigenLinearSolver(std::string const& /*solver_name*/,
         case EigenOption::SolverType::GMRES:
         case EigenOption::SolverType::IDRS:
         case EigenOption::SolverType::IDRSTABL:
-            solver_ = details::createIterativeSolver(option_.solver_type,
-                                                     option_.precon_type);
+            solver_ =
+                details::createIterativeSolver(option_.solver_type,
+                                               option_.precon_type,
+                                               option_.triangular_matrix_type);
             can_solve_rectangular_ = false;
             return;
         case EigenOption::SolverType::LeastSquareCG:
-            solver_ = details::createIterativeSolver(option_.solver_type,
-                                                     option_.precon_type);
+            solver_ =
+                details::createIterativeSolver(option_.solver_type,
+                                               option_.precon_type,
+                                               option_.triangular_matrix_type);
             can_solve_rectangular_ = true;
             return;
         case EigenOption::SolverType::PardisoLU:

MathLib/LinAlg/Eigen/EigenOption.cpp

@@ -20,6 +20,7 @@ EigenOption::EigenOption()
     precon_type = PreconType::NONE;
     max_iterations = static_cast<int>(1e6);
     error_tolerance = 1.e-16;
+    triangular_matrix_type = TriangularMatrixType::Lower;
 #ifdef USE_EIGEN_UNSUPPORTED
     scaling = false;
     restart = 30;
@@ -97,6 +98,25 @@ EigenOption::PreconType EigenOption::getPreconType(
     OGS_FATAL("Unknown Eigen preconditioner type `{:s}'", precon_name);
 }
 
+EigenOption::TriangularMatrixType EigenOption::getTriangularMatrixType(
+    const std::string& triangular_matrix_name)
+{
+    if (triangular_matrix_name == "Lower")
+    {
+        return TriangularMatrixType::Lower;
+    }
+    if (triangular_matrix_name == "Upper")
+    {
+        return TriangularMatrixType::Upper;
+    }
+    if (triangular_matrix_name == "LowerUpper")
+    {
+        return TriangularMatrixType::LowerUpper;
+    }
+
+    OGS_FATAL("Unknown triangular matrix type `{:s}'", triangular_matrix_name);
+}
+
 std::string EigenOption::getSolverName(SolverType const solver_type)
 {
     switch (solver_type)
@@ -139,4 +159,19 @@ std::string EigenOption::getPreconName(PreconType const precon_type)
     return "Invalid";
 }
 
+std::string EigenOption::getTriangularMatrixName(
+    TriangularMatrixType const triangular_matrix_type)
+{
+    switch (triangular_matrix_type)
+    {
+        case TriangularMatrixType::Lower:
+            return "Lower";
+        case TriangularMatrixType::Upper:
+            return "Upper";
+        case TriangularMatrixType::LowerUpper:
+            return "LowerUpper";
+    }
+    return "Invalid";
+}
+
 }  // namespace MathLib

MathLib/LinAlg/Eigen/EigenOption.h

@@ -40,10 +40,20 @@ struct EigenOption final
         ILUT
     };
 
+    /// triangular matrix type
+    enum class TriangularMatrixType : short
+    {
+        Lower,
+        Upper,
+        LowerUpper
+    };
+
     /// Linear solver type
     SolverType solver_type;
     /// Preconditioner type
     PreconType precon_type;
+    /// Triangular Matrix Type
+    TriangularMatrixType triangular_matrix_type;
     /// Maximum iteration count
     int max_iterations;
     /// Error tolerance
@@ -82,11 +92,24 @@ struct EigenOption final
     ///      NONE is returned.
     static PreconType getPreconType(const std::string& precon_name);
 
+    /// return a triangular matrix type from the name
+    ///
+    /// @param triangular_matrix_name
+    /// @return a triangular_matrix type
+    ///      If there is no triangular matrix type matched with the given name,
+    ///      NONE is returned.
+    static TriangularMatrixType getTriangularMatrixType(
+        const std::string& triangular_matrix_name);
+
     /// return a linear solver name from the solver type
     static std::string getSolverName(SolverType const solver_type);
 
     /// return a preconditioner name from the preconditioner type
     static std::string getPreconName(PreconType const precon_type);
+
+    /// return a triangular matrix name from the preconditioner type
+    static std::string getTriangularMatrixName(
+        TriangularMatrixType const triangular_matrix_type);
 };
 
 }  // namespace MathLib

MathLib/LinAlg/Eigen/LinearSolverOptionsParser.cpp

@@ -62,6 +62,14 @@ LinearSolverOptionsParser<EigenLinearSolver>::parseNameAndOptions(
     {
         options.max_iterations = *max_iteration_step;
     }
+    if (auto triangular_matrix_type =
+            //! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__triangular_matrix}
+        config->getConfigParameterOptional<std::string>("triangular_matrix"))
+    {
+        options.triangular_matrix_type =
+            MathLib::EigenOption::getTriangularMatrixType(
+                *triangular_matrix_type);
+    }
     if (auto scaling =
             //! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__scaling}
         config->getConfigParameterOptional<bool>("scaling"))

ProcessLib/SteadyStateDiffusion/Tests.cmake

@@ -143,6 +143,24 @@ foreach(mesh_size 1e4 2e4 3e4 4e4 5e4 1e5 1e6)
     )
 endforeach()
 
+# CG TriangularMatrix test
+foreach(matrix lower upper lowerupper)
+    set(RUNTIME 3)
+    AddTest(
+        NAME SteadyStateDiffusion_cube_1x1x1_1e0_${matrix}
+        PATH Elliptic/cube_1x1x1_SteadyStateDiffusion
+        RUNTIME ${RUNTIME}
+        EXECUTABLE ogs
+        EXECUTABLE_ARGS cube_1e0_${matrix}.xml --write-prj
+        TESTER vtkdiff
+        REQUIREMENTS NOT OGS_USE_MPI
+        DIFF_DATA
+        cube_1x1x1_hex_1e0.vtu 1e0_cube_1x1x1_hex_1e0_${matrix}_ts_1_t_1.000000.vtu Linear_1_to_minus1 pressure 1e-13 1e-13
+    )
+endforeach()
+
+
+
 # Test FixedTimeStepping and fixed output times
 AddTest(
     NAME SteadyStateDiffusion_square_1x1_quad_1e1_FixedTimeStepping_FixedOutputTimes

Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e0_lower.xml

+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProjectDiff base_file="cube_1e0.prj">
+    <replace sel="/*/time_loop/output/prefix/text()" after_includes="true">1e0_{:meshname}_lower</replace>
+    <add sel="/*/linear_solvers/linear_solver/eigen" after_includes="true">
+      <triangular_matrix>Lower</triangular_matrix>
+    </add>
+</OpenGeoSysProjectDiff>

Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e0_lowerupper.xml

+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProjectDiff base_file="cube_1e0.prj">
+    <replace sel="/*/time_loop/output/prefix/text()" after_includes="true">1e0_{:meshname}_lowerupper</replace>
+    <add sel="/*/linear_solvers/linear_solver/eigen" after_includes="true">
+      <triangular_matrix>LowerUpper</triangular_matrix>
+    </add>
+</OpenGeoSysProjectDiff>

Tests/Data/Elliptic/cube_1x1x1_SteadyStateDiffusion/cube_1e0_upper.xml

+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProjectDiff base_file="cube_1e0.prj">
+    <replace sel="/*/time_loop/output/prefix/text()" after_includes="true">1e0_{:meshname}_upper</replace>
+    <add sel="/*/linear_solvers/linear_solver/eigen" after_includes="true">
+      <triangular_matrix>Upper</triangular_matrix>
+    </add>
+</OpenGeoSysProjectDiff>

web/content/docs/userguide/basics/openmp.md

@@ -23,6 +23,11 @@ Note:
 * Some linear solvers and some preconditioners from the Eigen library run with a
   single thread only, see
   [here](https://eigen.tuxfamily.org/dox/TopicMultiThreading.html).
+* For Eigen's CG solver, multi-threading can be exploited using both triangular
+  matrices (see
+  [here](https://eigen.tuxfamily.org/dox/classEigen_1_1ConjugateGradient.html).
+  This can be set in OGS using `<triangular_matrix>LowerUpper</triangular_matrix>`
+  as linear solver option.
 * When using distributed memory parallelization (MPI) it might not make sense to
   use OpenMP at the same time: the different threads might compete for the same
   CPU resources.