From 0cf4ec61c2ce7f3334b11c41399feaf2c7679611 Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 29 Apr 2026 14:54:33 +0200 Subject: [PATCH] Introducing possibility to use dedicated linear solver for multisegment wells --- .../flow/BlackoilModelParameters.cpp | 12 ++ .../flow/BlackoilModelParameters.hpp | 15 ++ .../wells/MultisegmentWellEquations.cpp | 168 +++++++++++++++--- .../wells/MultisegmentWellEquations.hpp | 18 ++ .../wells/MultisegmentWellGeneric.hpp | 5 + opm/simulators/wells/WellInterfaceGeneric.hpp | 2 + 6 files changed, 199 insertions(+), 21 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelParameters.cpp b/opm/simulators/flow/BlackoilModelParameters.cpp index 05bbc7bdb1e..2e83491187c 100644 --- a/opm/simulators/flow/BlackoilModelParameters.cpp +++ b/opm/simulators/flow/BlackoilModelParameters.cpp @@ -86,6 +86,9 @@ BlackoilModelParameters::BlackoilModelParameters() local_well_solver_control_switching_ = Parameters::Get(); use_implicit_ipr_ = Parameters::Get(); check_group_constraints_inner_well_iterations_ = Parameters::Get(); + local_well_eq_linear_solver_ = Parameters::Get(); + local_well_eq_reuse_preconditioner_ = Parameters::Get(); + local_well_eq_use_legacy_global_operator_ = Parameters::Get(); nonlinear_solver_ = Parameters::Get(); const auto approach = Parameters::Get(); if (approach == "jacobi") { @@ -256,6 +259,15 @@ void BlackoilModelParameters::registerParameters() ("Compute implict IPR for stability checks and stable solution search"); Parameters::Register ("Allow checking of group constraints during inner well iterations"); + Parameters::Register + ("Linear solver preset or JSON file path for multisegment well equations. " + "Default is 'umfpack'. Use a filename ending in '.json' to provide a full solver configuration."); + Parameters::Register + ("If true, reuse the multisegment-well flexible solver and only update preconditioner " + "when matrix structure is unchanged. Default false recreates the solver each iteration."); + Parameters::Register + ("If true, use the legacy UMFPACK-based multisegment-well implementation for the global ISTL operator path. " + "Default true keeps previous behavior for apply/recover/extract while local well solves use the configured flexible solver."); Parameters::Register ("Maximum outer iterations in network solver before relaxing tolerance"); Parameters::Register diff --git a/opm/simulators/flow/BlackoilModelParameters.hpp b/opm/simulators/flow/BlackoilModelParameters.hpp index 41697958de0..86f1e7cd6ad 100644 --- a/opm/simulators/flow/BlackoilModelParameters.hpp +++ b/opm/simulators/flow/BlackoilModelParameters.hpp @@ -143,6 +143,9 @@ struct UseAverageDensityMsWells { static constexpr bool value = false; }; struct LocalWellSolveControlSwitching { static constexpr bool value = true; }; struct UseImplicitIpr { static constexpr bool value = true; }; struct CheckGroupConstraintsInnerWellIterations { static constexpr bool value = true; }; +struct LocalWellEqLinearSolver { static constexpr auto value = "umfpack"; }; +struct LocalWellEqReusePreconditioner { static constexpr bool value = false; }; +struct LocalWellEqUseLegacyGlobalOperator { static constexpr bool value = true; }; // Network solver parameters struct NetworkMaxStrictOuterIterations { static constexpr int value = 10; }; @@ -332,6 +335,18 @@ struct BlackoilModelParameters /// Whether to allow checking/changing to group controls during inner well iterations bool check_group_constraints_inner_well_iterations_; + /// Well-equation linear solver specification. Can be a solver name (e.g., "umfpack") + /// or a path to a JSON solver configuration file. + std::string local_well_eq_linear_solver_; + + /// If true, reuse the existing flexible solver and only update its preconditioner + /// when the matrix sparsity structure is unchanged. + bool local_well_eq_reuse_preconditioner_; + + /// If true, use the legacy UMFPACK-based implementation for multisegment well + /// contributions that are applied through the global ISTL operator. + bool local_well_eq_use_legacy_global_operator_; + /// Maximum number of iterations in the network solver before relaxing tolerance int network_max_strict_outer_iterations_; diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp index def1efa40f2..cbdb6d320fa 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.cpp +++ b/opm/simulators/wells/MultisegmentWellEquations.cpp @@ -22,7 +22,13 @@ #include #include +#include +#include +#include + +#include #include +#include #include #include @@ -52,9 +58,26 @@ namespace Opm { template MultisegmentWellEquations:: MultisegmentWellEquations(const MultisegmentWellGeneric& well, const ParallelWellInfo& parallel_well_info) - : well_(well) + : local_solver_reuse_preconditioner_(well.modelParameters().local_well_eq_reuse_preconditioner_) + , use_legacy_global_operator_(well.modelParameters().local_well_eq_use_legacy_global_operator_) + , well_(well) , parallelB_(duneB_, parallel_well_info) { + setupLocalSolverConfig(); +} + +template +void MultisegmentWellEquations::setupLocalSolverConfig() +{ + FlowLinearSolverParameters linear_params; + linear_params.reset(); + linear_params.linsolver_ = well_.modelParameters().local_well_eq_linear_solver_; + local_solver_prm_ = setupPropertyTree(linear_params, + /*linearSolverMaxIterSet=*/false, + /*linearSolverReductionSet=*/false, + /*tpsaSetup=*/false); + local_solver_prm_.put("solver", local_solver_prm_.get("solver", "umfpack")); + local_solver_prm_.put("verbosity", local_solver_prm_.get("verbosity", 0)); } template @@ -135,6 +158,7 @@ init(const int numPerfs, // Store the global index of well perforated cells cells_ = cells; + matrix_structure_changed_ = true; } template @@ -146,7 +170,58 @@ void MultisegmentWellEquations::clear() resWell_ = 0.0; if constexpr (std::is_same_v) { duneDSolver_.reset(); + flexibleDsolver_.reset(); + duneDOperator_.reset(); + } +} + +template +typename MultisegmentWellEquations::BVectorWell +MultisegmentWellEquations::solveLocalDSystem(const BVectorWell& rhs) const +{ + if (!flexibleDsolver_) { + OPM_THROW(std::logic_error, + "MultisegmentWell local linear solver was used before createSolver() initialized it"); + } + + BVectorWell x(rhs.size()); + x = Scalar{0.0}; + auto rhs_copy = rhs; + Dune::InverseOperatorResult result; + flexibleDsolver_->apply(x, rhs_copy, result); + if (!result.converged) { + OPM_THROW(std::runtime_error, + "MultisegmentWell local linear solver failed to converge for well equation system"); + } + return x; +} + +template +typename MultisegmentWellEquations::BVectorWell +MultisegmentWellEquations::solveLegacyGlobalDSystem(const BVectorWell& rhs) const +{ +#if HAVE_SUITESPARSE_UMFPACK + if (!duneDSolver_) { + OPM_THROW(std::logic_error, + "MultisegmentWell legacy UMFPACK solver was used before createSolver() initialized it"); } + return mswellhelpers::applyUMFPack(*duneDSolver_, rhs); +#else + OPM_THROW(std::runtime_error, + "MultisegmentWell legacy global-operator support requires SuiteSparse/UMFPACK. " + "Disable LocalWellEqUseLegacyGlobalOperator to use the flexible solver path instead."); +#endif +} + +template +typename MultisegmentWellEquations::BVectorWell +MultisegmentWellEquations::solveGlobalOperatorDSystem(const BVectorWell& rhs) const +{ + if (use_legacy_global_operator_) { + return solveLegacyGlobalDSystem(rhs); + } + + return solveLocalDSystem(rhs); } template @@ -162,7 +237,7 @@ apply(const BVector& x, BVector& Ax) const // the single process to complete the computation. // invDBx = duneD^-1 * Bx_ if constexpr (std::is_same_v) { - const BVectorWell invDBx = mswellhelpers::applyUMFPack(*duneDSolver_, Bx); + const BVectorWell invDBx = solveGlobalOperatorDSystem(Bx); // Ax.size() == 0 indicates that there are no active perforations on this process. // Then, Ax does not need to be updated by the following calculation. if (Ax.size() > 0) { @@ -184,7 +259,7 @@ apply(BVector& r) const // because the other processes would remain idle while waiting for // the single process to complete the computation. // invDrw_ = duneD^-1 * resWell_ - const BVectorWell invDrw = mswellhelpers::applyUMFPack(*duneDSolver_, resWell_); + const BVectorWell invDrw = solveGlobalOperatorDSystem(resWell_); // r = r - duneC_^T * invDrw duneC_.mmtv(invDrw, r); } @@ -194,20 +269,40 @@ apply(BVector& r) const template void MultisegmentWellEquations::createSolver() { -#if HAVE_SUITESPARSE_UMFPACK if constexpr (std::is_same_v) { - OPM_THROW(std::runtime_error, "MultisegmentWell support requires UMFPACK, " - "and UMFPACK does not support float"); + OPM_THROW(std::runtime_error, "MultisegmentWell local well equations are only supported for double precision"); } else { - if (duneDSolver_) { - return; + const bool can_reuse = local_solver_reuse_preconditioner_ + && !matrix_structure_changed_ + && static_cast(flexibleDsolver_); + + if (can_reuse) { + flexibleDsolver_->preconditioner().update(); + } else { + duneDOperator_ = std::make_shared(duneD_); + flexibleDsolver_ = std::make_shared>( + *duneDOperator_, + local_solver_prm_, + std::function(), + 0); } - duneDSolver_ = std::make_shared>(duneD_, 0); - } + + if (use_legacy_global_operator_) { +#if HAVE_SUITESPARSE_UMFPACK + if (!duneDSolver_) { + duneDSolver_ = std::make_shared>(duneD_, 0); + } #else - OPM_THROW(std::runtime_error, "MultisegmentWell support requires UMFPACK. " - "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile."); + OPM_THROW(std::runtime_error, + "MultisegmentWell legacy global-operator support requires SuiteSparse/UMFPACK. " + "Disable LocalWellEqUseLegacyGlobalOperator to use the flexible solver path instead."); #endif + } else { + duneDSolver_.reset(); + } + + matrix_structure_changed_ = false; + } } template @@ -218,7 +313,7 @@ MultisegmentWellEquations::solve() const // It is ok to do this on each process instead of only on one, // because the other processes would remain idle while waiting for // the single process to complete the computation. - return mswellhelpers::applyUMFPack(*duneDSolver_, resWell_); + return solveLocalDSystem(resWell_); } else { return {}; @@ -233,7 +328,7 @@ MultisegmentWellEquations::solve(const BV // It is ok to do this on each process instead of only on one, // because the other processes would remain idle while waiting for // the single process to complete the computation. - return mswellhelpers::applyUMFPack(*duneDSolver_, rhs); + return solveLocalDSystem(rhs); } else { return {}; @@ -253,7 +348,7 @@ recoverSolutionWell(const BVector& x, BVectorWell& xw) const // It is ok to do this on each process instead of only on one, // because the other processes would remain idle while waiting for // the single process to complete the computation. - xw = mswellhelpers::applyUMFPack(*duneDSolver_, resWell); + xw = solveGlobalOperatorDSystem(resWell); } } @@ -333,10 +428,8 @@ void MultisegmentWellEquations:: extract(SparseMatrixAdapter& jacobian) const { if constexpr (std::is_same_v) { - const auto invDuneD = mswellhelpers::invertWithUMFPack(duneD_.M(), - numWellEq, - *duneDSolver_); - + const auto addExtractContribution = [&](const auto& invDuneD) + { // We need to change matrix A as follows // A -= C^T D^-1 B // D is a (nseg x nseg) block matrix with (4 x 4) blocks. @@ -347,12 +440,12 @@ extract(SparseMatrixAdapter& jacobian) const // i.e. the columns of B/C have no more than one nonzero. for (std::size_t rowC = 0; rowC < duneC_.N(); ++rowC) { for (auto colC = duneC_[rowC].begin(), - endC = duneC_[rowC].end(); colC != endC; ++colC) { + endC = duneC_[rowC].end(); colC != endC; ++colC) { // map the well perforated cell index to global cell index const auto row_index = cells_[colC.index()]; for (std::size_t rowB = 0; rowB < duneB_.N(); ++rowB) { for (auto colB = duneB_[rowB].begin(), - endB = duneB_[rowB].end(); colB != endB; ++colB) { + endB = duneB_[rowB].end(); colB != endB; ++colB) { const auto col_index = cells_[colB.index()]; OffDiagMatrixBlockWellType tmp1; detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type()); @@ -363,6 +456,39 @@ extract(SparseMatrixAdapter& jacobian) const } } } + }; + + if (use_legacy_global_operator_) { +#if HAVE_SUITESPARSE_UMFPACK + const auto invDuneD = mswellhelpers::invertWithUMFPack(duneD_.M(), + numWellEq, + *duneDSolver_); + addExtractContribution(invDuneD); +#else + OPM_THROW(std::runtime_error, + "MultisegmentWell legacy global-operator support requires SuiteSparse/UMFPACK. " + "Disable LocalWellEqUseLegacyGlobalOperator to use the flexible solver path instead."); +#endif + } else { + std::vector> invDuneD( + duneD_.N(), std::vector(duneD_.N(), DiagMatrixBlockWellType(0.0))); + const int nrows = static_cast(duneD_.N()); + for (int row = 0; row < nrows; ++row) { + BVectorWell e(nrows); + e = Scalar{0.0}; + for (int localEq = 0; localEq < numWellEq; ++localEq) { + e[row][localEq] = 1.0; + const auto col = solveLocalDSystem(e); + for (int col_idx = 0; col_idx < nrows; ++col_idx) { + for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + invDuneD[col_idx][row][eq_idx][localEq] = col[col_idx][eq_idx]; + } + } + e[row][localEq] = 0.0; + } + } + addExtractContribution(invDuneD); + } } } diff --git a/opm/simulators/wells/MultisegmentWellEquations.hpp b/opm/simulators/wells/MultisegmentWellEquations.hpp index 866e4b5dbe0..b786f516631 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.hpp +++ b/opm/simulators/wells/MultisegmentWellEquations.hpp @@ -23,17 +23,21 @@ #define OPM_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED #include +#include #include #include #include #include #include #include +#include #include #include namespace Dune { +template class MatrixAdapter; +template class FlexibleSolver; template class UMFPack; } @@ -152,6 +156,15 @@ class MultisegmentWellEquations EmptyType>; // TODO: c++20: add no_unique_address mutable UMFPackSolver duneDSolver_; + using DiagOperator = Dune::MatrixAdapter; + using FlexibleDiagSolver = std::shared_ptr>; + mutable std::shared_ptr duneDOperator_; + mutable FlexibleDiagSolver flexibleDsolver_; + mutable bool matrix_structure_changed_{true}; + mutable Opm::PropertyTree local_solver_prm_; + bool local_solver_reuse_preconditioner_{false}; + bool use_legacy_global_operator_{true}; + // residuals of the well equations BVectorWell resWell_; @@ -162,6 +175,11 @@ class MultisegmentWellEquations // Wrapper for the parallel application of B for distributed wells mswellhelpers::ParallellMSWellB parallelB_; + + BVectorWell solveLocalDSystem(const BVectorWell& rhs) const; + BVectorWell solveLegacyGlobalDSystem(const BVectorWell& rhs) const; + BVectorWell solveGlobalOperatorDSystem(const BVectorWell& rhs) const; + void setupLocalSolverConfig(); }; } diff --git a/opm/simulators/wells/MultisegmentWellGeneric.hpp b/opm/simulators/wells/MultisegmentWellGeneric.hpp index 06fa7dc1a18..631396d3aa2 100644 --- a/opm/simulators/wells/MultisegmentWellGeneric.hpp +++ b/opm/simulators/wells/MultisegmentWellGeneric.hpp @@ -48,6 +48,11 @@ class MultisegmentWellGeneric /// number of segments for this well int numberOfSegments() const; + const typename WellInterfaceGeneric::ModelParameters& modelParameters() const + { + return baseif_.modelParameters(); + } + protected: explicit MultisegmentWellGeneric(WellInterfaceGeneric& baseif); diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 54835a459b2..31dfcdf83a6 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -127,6 +127,8 @@ class WellInterfaceGeneric { int currentStep() const { return this->current_step_; } + const ModelParameters& modelParameters() const { return this->param_; } + int pvtRegionIdx() const { return pvtRegionIdx_; } const GuideRate* guideRate() const { return guide_rate_; }