diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp index 1f0c394355fe60bba38875674a9e2e4f6b07b775..2dca2852594c06020d50c4390a5820eb60c77710 100644 --- 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 index 1955f1fd8495e3609af9556119b066bc3e725887..2f58fc83738f5137132abb4778f6b7cb7e615b19 100644 --- a/MathLib/LinAlg/Eigen/EigenLinearSolver.h +++ b/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(); diff --git a/MathLib/LinAlg/Eigen/EigenOption.cpp b/MathLib/LinAlg/Eigen/EigenOption.cpp index b9df0ff226ace71eb98ea7d1eb4a5b82b7092eaa..07b11082bffb88ec04189ebedbc4577f79d287c1 100644 --- 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 index 54815c948b92818850031565d9b59dc20e1c1741..217d69eae31e94ca97a8514c7eb2efb17c2428cf 100644 --- 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 index 23763277d0a9241931e0537cb5330f1eb52f5654..3543d65533131c26542f07f35a48713e88a9368e 100644 --- a/MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h +++ b/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_; diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp index 7310848f3ea60f118df33340b5dca8ea1ed13060..1fc14a8bd588b019e9d3dbf442787d5759dd2154 100644 --- 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 index 0dfcbc359560f5945d8996c3b1ac0ff1a04ac58c..59c1ad0d8d2ba9d16a41aa6a339b7ee9d7512287 100644 --- 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 @@ 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 }; diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp index eb533380fd16d64f63b39752765ace7c16bba02c..516ffa207d8dced39eeda68dc48909cd6a7ffb3a 100644 --- 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 index 269f8d3089d311116151f077352c64231fc671fd..ddaa11ce26d82c80d1c9cfb4cddff222fd306dd8 100644 --- 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 new file mode 100644 index 0000000000000000000000000000000000000000..7aae83c25dd8cd11ff32bd32b88b3f070c8a9b20 --- /dev/null +++ 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 new file mode 100644 index 0000000000000000000000000000000000000000..50024b03ab9ca6c09871ea0cac7dd710af59a062 --- /dev/null +++ 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 new file mode 100644 index 0000000000000000000000000000000000000000..9c9ee6554bb0b061c5715547b009dd6fce7f21d4 --- /dev/null +++ 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> diff --git a/Tests/Data/Parabolic/T/3D_Beier_sandbox/fixed_power_constant_flow_algebraicBC_LSCG.xml b/Tests/Data/Parabolic/T/3D_Beier_sandbox/fixed_power_constant_flow_algebraicBC_LSCG.xml new file mode 100644 index 0000000000000000000000000000000000000000..8cc4e749e8a3d3d7c8a7eb969c423318ee09473c --- /dev/null +++ b/Tests/Data/Parabolic/T/3D_Beier_sandbox/fixed_power_constant_flow_algebraicBC_LSCG.xml @@ -0,0 +1,15 @@ +<?xml version="1.0" encoding="ISO-8859-1"?> +<OpenGeoSysProjectDiff base_file="fixed_power_constant_flow.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()">fixed_power_constant_flow_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>