diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index 6ecfc9f3baa866d4831b175106c2c891b78a3900..2c218b4c432b345f8c3f5e987711ba7c0ec5468e 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -142,7 +142,7 @@ if(NOT OGS_USE_MPI)
EXECUTABLE_ARGS line_${mesh_size}_robin_right_picard.prj
WRAPPER time
TESTER vtkdiff
- ABSTOL 1e-14 RELTOL 1e-14
+ ABSTOL 4e-14 RELTOL 2e-14
DIFF_DATA
line_1_line_${mesh_size}.vtu line_${mesh_size}_robin_right_picard_pcs_0_ts_1_t_1.000000.vtu D1_left_N1_right pressure
)
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
index 3d05dfbdbb1eeed15657e26939d31d1c0b9b89cd..d6f45b29c494ba3424b0c637963d3e4b54d0e3d3 100644
--- a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
@@ -102,6 +102,53 @@ private:
T_SOLVER _solver;
};
+template <template <typename, typename> class Solver, typename Precon>
+std::unique_ptr<EigenLinearSolverBase> createIterativeSolver()
+{
+ using Slv = EigenIterativeLinearSolver<
+ Solver<EigenMatrix::RawMatrixType, Precon>>;
+ return std::unique_ptr<Slv>(new Slv);
+}
+
+template <template <typename, typename> class Solver>
+std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
+ EigenOption::PreconType precon_type)
+{
+ switch (precon_type) {
+ case EigenOption::PreconType::NONE:
+ return createIterativeSolver<Solver,
+ Eigen::IdentityPreconditioner>();
+ case EigenOption::PreconType::DIAGONAL:
+ return createIterativeSolver<
+ Solver, Eigen::DiagonalPreconditioner<double>>();
+ case EigenOption::PreconType::ILUT:
+ // TODO for this preconditioner further options can be passed.
+ // see https://eigen.tuxfamily.org/dox/classEigen_1_1IncompleteLUT.html
+ return createIterativeSolver<
+ Solver, Eigen::IncompleteLUT<double>>();
+ default:
+ OGS_FATAL("Invalid Eigen preconditioner type.");
+ }
+}
+
+template <typename Mat, typename Precon>
+using EigenCGSolver = Eigen::ConjugateGradient<Mat, Eigen::Lower, Precon>;
+
+std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
+ EigenOption::SolverType solver_type, EigenOption::PreconType precon_type)
+{
+ switch (solver_type) {
+ case EigenOption::SolverType::BiCGSTAB: {
+ return createIterativeSolver<Eigen::BiCGSTAB>(precon_type);
+ }
+ case EigenOption::SolverType::CG: {
+ return createIterativeSolver<EigenCGSolver>(precon_type);
+ }
+ default:
+ OGS_FATAL("Invalid Eigen iterative linear solver type. Aborting.");
+ }
+}
+
} // details
EigenLinearSolver::EigenLinearSolver(
@@ -115,26 +162,21 @@ EigenLinearSolver::EigenLinearSolver(
// TODO for my taste it is much too unobvious that the default solver type
// currently is SparseLU.
- switch (_option.solver_type)
- {
- case EigenOption::SolverType::SparseLU: {
- using SolverType = Eigen::SparseLU<Matrix, Eigen::COLAMDOrdering<int>>;
- _solver.reset(new details::EigenDirectLinearSolver<SolverType>);
- break;
- }
- case EigenOption::SolverType::BiCGSTAB: {
- using SolverType = Eigen::BiCGSTAB<Matrix, Eigen::DiagonalPreconditioner<double>>;
- _solver.reset(new details::EigenIterativeLinearSolver<SolverType>);
- break;
- }
- case EigenOption::SolverType::CG: {
- using SolverType = Eigen::ConjugateGradient<Matrix, Eigen::Lower, Eigen::DiagonalPreconditioner<double>>;
- _solver.reset(new details::EigenIterativeLinearSolver<SolverType>);
- break;
- }
- case EigenOption::SolverType::INVALID:
- OGS_FATAL("Invalid Eigen linear solver type. Aborting.");
+ switch (_option.solver_type) {
+ case EigenOption::SolverType::SparseLU: {
+ using SolverType =
+ Eigen::SparseLU<Matrix, Eigen::COLAMDOrdering<int>>;
+ _solver.reset(new details::EigenDirectLinearSolver<SolverType>);
+ return;
+ }
+ case EigenOption::SolverType::BiCGSTAB:
+ case EigenOption::SolverType::CG:
+ _solver = details::createIterativeSolver(_option.solver_type,
+ _option.precon_type);
+ return;
}
+
+ OGS_FATAL("Invalid Eigen linear solver type. Aborting.");
}
EigenLinearSolver::~EigenLinearSolver() = default;
@@ -147,20 +189,24 @@ void EigenLinearSolver::setOption(BaseLib::ConfigTree const& option)
if (!ptSolver)
return;
- //! \ogs_file_param{linear_solver__eigen__solver_type}
- if (auto solver_type = ptSolver->getConfigParameterOptional<std::string>("solver_type")) {
+ if (auto solver_type =
+ //! \ogs_file_param{linear_solver__eigen__solver_type}
+ ptSolver->getConfigParameterOptional<std::string>("solver_type")) {
_option.solver_type = _option.getSolverType(*solver_type);
}
- //! \ogs_file_param{linear_solver__eigen__precon_type}
- if (auto precon_type = ptSolver->getConfigParameterOptional<std::string>("precon_type")) {
+ if (auto precon_type =
+ //! \ogs_file_param{linear_solver__eigen__precon_type}
+ ptSolver->getConfigParameterOptional<std::string>("precon_type")) {
_option.precon_type = _option.getPreconType(*precon_type);
}
- //! \ogs_file_param{linear_solver__eigen__error_tolerance}
- if (auto error_tolerance = ptSolver->getConfigParameterOptional<double>("error_tolerance")) {
+ if (auto error_tolerance =
+ //! \ogs_file_param{linear_solver__eigen__error_tolerance}
+ ptSolver->getConfigParameterOptional<double>("error_tolerance")) {
_option.error_tolerance = *error_tolerance;
}
- //! \ogs_file_param{linear_solver__eigen__max_iteration_step}
- if (auto max_iteration_step = ptSolver->getConfigParameterOptional<int>("max_iteration_step")) {
+ if (auto max_iteration_step =
+ //! \ogs_file_param{linear_solver__eigen__max_iteration_step}
+ ptSolver->getConfigParameterOptional<int>("max_iteration_step")) {
_option.max_iterations = *max_iteration_step;
}
}
diff --git a/MathLib/LinAlg/Eigen/EigenOption.cpp b/MathLib/LinAlg/Eigen/EigenOption.cpp
index 36bc4a671ff8fbfbba51a5a3f0d6088f21c2b3cf..3c273491484a7404ccbfde4e5be205f7b3db49c1 100644
--- a/MathLib/LinAlg/Eigen/EigenOption.cpp
+++ b/MathLib/LinAlg/Eigen/EigenOption.cpp
@@ -8,6 +8,7 @@
*/
#include "EigenOption.h"
+#include "BaseLib/Error.h"
namespace MathLib
{
@@ -22,27 +23,26 @@ EigenOption::EigenOption()
EigenOption::SolverType EigenOption::getSolverType(const std::string &solver_name)
{
-#define RETURN_SOLVER_ENUM_IF_SAME_STRING(str, TypeName) \
- if (#TypeName==(str)) return SolverType::TypeName;
-
- RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, CG);
- RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, BiCGSTAB);
- RETURN_SOLVER_ENUM_IF_SAME_STRING(solver_name, SparseLU);
-
- return SolverType::INVALID;
-#undef RETURN_SOLVER_ENUM_IF_SAME_STRING
+ if (solver_name == "CG")
+ return SolverType::CG;
+ if (solver_name == "BiCGSTAB")
+ return SolverType::BiCGSTAB;
+ if (solver_name == "SparseLU")
+ return SolverType::SparseLU;
+
+ OGS_FATAL("Unknown Eigen solver type `%s'", solver_name.c_str());
}
EigenOption::PreconType EigenOption::getPreconType(const std::string &precon_name)
{
-#define RETURN_PRECOM_ENUM_IF_SAME_STRING(str, TypeName) \
- if (#TypeName==(str)) return PreconType::TypeName;
-
- RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, NONE);
- RETURN_PRECOM_ENUM_IF_SAME_STRING(precon_name, DIAGONAL);
-
- return PreconType::NONE;
-#undef RETURN_PRECOM_ENUM_IF_SAME_STRING
+ if (precon_name == "NONE")
+ return PreconType::NONE;
+ if (precon_name == "DIAGONAL")
+ return PreconType::DIAGONAL;
+ if (precon_name == "ILUT")
+ return PreconType::ILUT;
+
+ OGS_FATAL("Unknown Eigen preconditioner type `%s'", precon_name.c_str());
}
} //MathLib
diff --git a/MathLib/LinAlg/Eigen/EigenOption.h b/MathLib/LinAlg/Eigen/EigenOption.h
index b0e0f34688af2ef5d4b96570a9b802b3dd865c52..5b64fc41a7b42da5fb1ab976a975c32222cba4b1 100644
--- a/MathLib/LinAlg/Eigen/EigenOption.h
+++ b/MathLib/LinAlg/Eigen/EigenOption.h
@@ -21,7 +21,6 @@ struct EigenOption final
/// Solver type
enum class SolverType : short
{
- INVALID,
CG,
BiCGSTAB,
SparseLU
@@ -31,7 +30,8 @@ struct EigenOption final
enum class PreconType : short
{
NONE,
- DIAGONAL
+ DIAGONAL,
+ ILUT
};
/// Linear solver type
diff --git a/NumLib/NamedFunctionCaller.cpp b/NumLib/NamedFunctionCaller.cpp
index b7111612e0f45bf3bee11137e8ef745d16dd414c..e7b89d90119b9177a9406a0742686e949a76afa1 100644
--- a/NumLib/NamedFunctionCaller.cpp
+++ b/NumLib/NamedFunctionCaller.cpp
@@ -204,8 +204,6 @@ double NamedFunctionCaller::call(
assert(_deferred_plugs.empty() &&
"You must call applyPlugs() before this method!");
- DBUG("Preparing call of fct #%lu %s()", function_idx,
- _named_functions[function_idx].getName().c_str());
auto const& sis_sos = _map_sink_source[function_idx];
assert(sis_sos.size() ==
_named_functions[function_idx].getArgumentNames().size());
@@ -217,17 +215,12 @@ double NamedFunctionCaller::call(
if (source >= 0) {
fct_args[sink] = call(source, unbound_arguments);
- DBUG("setting %luth argument to %g", sink, fct_args[sink]);
} else {
assert(source != _uninitialized);
fct_args[sink] = unbound_arguments[-source-1];
- DBUG("setting %luth argument to %g", sink, fct_args[sink]);
}
}
- DBUG("Finished preparing call of fct #%lu %s()", function_idx,
- _named_functions[function_idx].getName().c_str());
-
return _named_functions[function_idx].call(fct_args);
}
diff --git a/Tests/Data b/Tests/Data
index 2085e1093f2b2235d82a6e8091c07fb98b6d3dfb..077539a7a9cb452c0797ec3dc3a255df1e05392d 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 2085e1093f2b2235d82a6e8091c07fb98b6d3dfb
+Subproject commit 077539a7a9cb452c0797ec3dc3a255df1e05392d