diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp index 1f0c394355f..2dca2852594 100644 --- a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp +++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp @@ -388,6 +388,9 @@ std::unique_ptr createIterativeSolver( case EigenOption::PreconType::DIAGONAL: return createIterativeSolver< Solver, Eigen::DiagonalPreconditioner>(); + case EigenOption::PreconType::LeastSquareDIAGONAL: + return createIterativeSolver< + Solver, Eigen::LeastSquareDiagonalPreconditioner>(); case EigenOption::PreconType::ILUT: // TODO for this preconditioner further options can be passed. // see @@ -429,6 +432,11 @@ std::unique_ptr createIterativeSolver( { return createIterativeSolver(precon_type); } + case EigenOption::SolverType::LeastSquareCG: + { + return createIterativeSolver( + precon_type); + } case EigenOption::SolverType::GMRES: { #ifdef USE_EIGEN_UNSUPPORTED @@ -488,6 +496,7 @@ EigenLinearSolver::EigenLinearSolver(std::string const& /*solver_name*/, Eigen::SparseLU>; solver_ = std::make_unique< details::EigenDirectLinearSolver>(); + 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; solver_.reset(new details::EigenDirectLinearSolver); + can_solve_rectangular_ = false; return; #else OGS_FATAL( diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.h b/MathLib/LinAlg/Eigen/EigenLinearSolver.h index 1955f1fd849..2f58fc83738 100644 --- 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 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 b9df0ff226a..07b11082bff 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 54815c948b9..217d69eae31 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 23763277d0a..3543d655331 100644 --- 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 index 7310848f3ea..1fc14a8bd58 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 0dfcbc35956..59c1ad0d8d2 100644 --- a/MathLib/LinAlg/PETSc/PETScLinearSolver.h +++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.h @@ -18,6 +18,8 @@ #include +#include +#include #include #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( + {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 eb533380fd1..516ffa207d8 100644 --- a/NumLib/ODESolver/NonlinearSolver.cpp +++ b/NumLib/ODESolver/NonlinearSolver.cpp @@ -170,8 +170,18 @@ NonlinearSolverStatus NonlinearSolver::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 269f8d3089d..ddaa11ce26d 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 00000000000..7aae83c25dd --- /dev/null +++ b/Tests/Data/Parabolic/T/3D_BHE_Sandwich/sandwich_algebraicBC_LSCG.xml @@ -0,0 +1,15 @@ + + + + true + + + 100 + + + true + + sandwich_algebraic_bc_LSCG + LeastSquareCG + LeastSquareDIAGONAL + 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 00000000000..50024b03ab9 --- /dev/null +++ b/Tests/Data/Parabolic/T/3D_BHE_Sandwich/sandwich_fixed_power_algebraicBC_LSCG.xml @@ -0,0 +1,15 @@ + + + + true + + + 100 + + + true + + sandwich_fixed_power_algebraic_bc_LSCG + LeastSquareCG + LeastSquareDIAGONAL + 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 00000000000..9c9ee6554bb --- /dev/null +++ b/Tests/Data/Parabolic/T/3D_Beier_sandbox/beier_sandbox_algebraicBC_LSCG.xml @@ -0,0 +1,15 @@ + + + + true + + + 100 + + + true + + beier_sandbox_algebraic_bc_LSCG + LeastSquareCG + LeastSquareDIAGONAL + 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 00000000000..8cc4e749e8a --- /dev/null +++ b/Tests/Data/Parabolic/T/3D_Beier_sandbox/fixed_power_constant_flow_algebraicBC_LSCG.xml @@ -0,0 +1,15 @@ + + + + true + + + 100 + + + true + + fixed_power_constant_flow_algebraic_bc_LSCG + LeastSquareCG + LeastSquareDIAGONAL +