From a011244d9caf10d5422ecd6d5a63269f071e2a95 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 23 Jul 2020 10:31:37 +0200 Subject: [PATCH 1/9] refactored invDxDirect to allow for reusing the decomposition. --- opm/simulators/wells/MSWellHelpers.hpp | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index a3b0cfcc2..f92d6b732 100644 --- a/opm/simulators/wells/MSWellHelpers.hpp +++ b/opm/simulators/wells/MSWellHelpers.hpp @@ -1,6 +1,7 @@ /* Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. Copyright 2017 Statoil ASA. + Copyright 2020 Equinor ASA. This file is part of the Open Porous Media project (OPM). @@ -36,17 +37,17 @@ namespace Opm { namespace mswellhelpers { - // obtain y = D^-1 * x with a direct solver + +#if HAVE_UMFPACK + /// Applies umfpack and checks for singularity template VectorType - invDXDirect(const MatrixType& D, VectorType x) + applyUMFPack(Dune::UMFPack& linsolver, VectorType x) { -#if HAVE_UMFPACK + // The copy of x seems mandatory for calling UMFPack! VectorType y(x.size()); y = 0.; - Dune::UMFPack linsolver(D, 0); - // Object storing some statistics about the solving process Dune::InverseOperatorResult res; @@ -62,8 +63,18 @@ namespace mswellhelpers } } } - return y; + } +#endif + + // obtain y = D^-1 * x with a direct solver + template + VectorType + invDXDirect(const MatrixType& D, VectorType x) + { +#if HAVE_UMFPACK + Dune::UMFPack linsolver(D, 0); + return applyUMFPack(linsolver, x); #else // this is not thread safe OPM_THROW(std::runtime_error, "Cannot use invDXDirect() without UMFPACK. " From ce409737febb6aff5f261a3b4dfd7e48a1c27888 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Fri, 24 Jul 2020 17:58:12 +0200 Subject: [PATCH 2/9] Allow applyUMFPack to decompose if needed. --- opm/simulators/wells/MSWellHelpers.hpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index f92d6b732..a1ae2e90c 100644 --- a/opm/simulators/wells/MSWellHelpers.hpp +++ b/opm/simulators/wells/MSWellHelpers.hpp @@ -42,8 +42,13 @@ namespace mswellhelpers /// Applies umfpack and checks for singularity template VectorType - applyUMFPack(Dune::UMFPack& linsolver, VectorType x) + applyUMFPack(const MatrixType& D, std::shared_ptr >& linsolver, VectorType x) { + if (!linsolver) + { + linsolver.reset(new Dune::UMFPack(D, 0)); + } + // The copy of x seems mandatory for calling UMFPack! VectorType y(x.size()); y = 0.; @@ -52,7 +57,7 @@ namespace mswellhelpers Dune::InverseOperatorResult res; // Solve - linsolver.apply(y, x, res); + linsolver->apply(y, x, res); // Checking if there is any inf or nan in y // it will be the solution before we find a way to catch the singularity of the matrix @@ -73,8 +78,8 @@ namespace mswellhelpers invDXDirect(const MatrixType& D, VectorType x) { #if HAVE_UMFPACK - Dune::UMFPack linsolver(D, 0); - return applyUMFPack(linsolver, x); + std::shared_ptr > linsolver; + return applyUMFPack(D, linsolver, x); #else // this is not thread safe OPM_THROW(std::runtime_error, "Cannot use invDXDirect() without UMFPACK. " From 7261759065467f97e438cc71f80c29c2ff9a0391 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Fri, 24 Jul 2020 18:08:35 +0200 Subject: [PATCH 3/9] Reuse UMFPack decomposition whenever possible. We hold a shared pointer to the umpfack solver that gets reset whenever the matrix is changed. When applying the decomposition we will be recomputed if the solver pointer is null. --- opm/simulators/wells/MultisegmentWell.hpp | 6 +++++- .../wells/MultisegmentWell_impl.hpp | 21 ++++++++++++++----- 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index f2b61b95f..a09921fd2 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -234,8 +234,12 @@ namespace Opm // two off-diagonal matrices mutable OffDiagMatWell duneB_; mutable OffDiagMatWell duneC_; - // diagonal matrix for the well + // "diagonal" matrix for the well. It has offdiagonal entries for inlets and outlets. mutable DiagMatWell duneD_; + /// \brief solver for diagonal matrix + /// + /// This is a shared_ptr as MultisegmentWell is copied in computeWellPotentials... + mutable std::shared_ptr > duneDSolver_; // residuals of the well equations mutable BVectorWell resWell_; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 4da42b1dd..0189b2797 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -624,7 +624,8 @@ namespace Opm duneB_.mv(x, Bx); // invDBx = duneD^-1 * Bx_ - const BVectorWell invDBx = mswellhelpers::invDXDirect(duneD_, Bx); + // const BVectorWell invDBx = mswellhelpers::invDXDirect(duneD_, Bx); + const BVectorWell invDBx = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, Bx); // Ax = Ax - duneC_^T * invDBx duneC_.mmtv(invDBx,Ax); @@ -640,7 +641,8 @@ namespace Opm apply(BVector& r) const { // invDrw_ = duneD^-1 * resWell_ - const BVectorWell invDrw = mswellhelpers::invDXDirect(duneD_, resWell_); + // const BVectorWell invDrw = mswellhelpers::invDXDirect(duneD_, resWell_); + const BVectorWell invDrw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_); // r = r - duneC_^T * invDrw duneC_.mmtv(invDrw, r); } @@ -999,7 +1001,8 @@ namespace Opm // resWell = resWell - B * x duneB_.mmv(x, resWell); // xw = D^-1 * resWell - xw = mswellhelpers::invDXDirect(duneD_, resWell); + //xw = mswellhelpers::invDXDirect(duneD_, resWell); + xw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell); } @@ -1013,7 +1016,8 @@ namespace Opm { // We assemble the well equations, then we check the convergence, // which is why we do not put the assembleWellEq here. - const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); + //const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); + const BVectorWell dx_well = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_); updateWellState(dx_well, well_state, deferred_logger); } @@ -1923,6 +1927,7 @@ namespace Opm for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + numEq); } + duneDSolver_.reset(); } @@ -2087,6 +2092,7 @@ namespace Opm if (accelerationalPressureLossConsidered()) { handleAccelerationPressureLoss(seg, well_state); } + duneDSolver_.reset(); } @@ -2181,6 +2187,7 @@ namespace Opm duneD_[seg][seg][SPres][GTotal] -= accelerationPressureLoss.derivative(GTotal + numEq); duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + numEq); duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + numEq); + duneDSolver_.reset(); } @@ -2399,7 +2406,8 @@ namespace Opm assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, deferred_logger); - const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); + // const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); + const BVectorWell dx_well = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_); if (it > param_.strict_inner_iter_ms_wells_) relax_convergence = true; @@ -2656,6 +2664,7 @@ namespace Opm well_state.segPressDropFriction()[seg] + well_state.segPressDropAcceleration()[seg]; } + duneDSolver_.reset(); } @@ -3207,6 +3216,7 @@ namespace Opm duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq); } + duneDSolver_.reset(); } @@ -3248,6 +3258,7 @@ namespace Opm for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq); } + duneDSolver_.reset(); } From e96c4a99085a9497557dc73629f76767aa6bb365 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 18 Aug 2020 17:52:51 +0200 Subject: [PATCH 4/9] Removed unused method invDXDirect. --- opm/simulators/wells/MSWellHelpers.hpp | 17 ++--------------- 1 file changed, 2 insertions(+), 15 deletions(-) diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index a1ae2e90c..8295a481a 100644 --- a/opm/simulators/wells/MSWellHelpers.hpp +++ b/opm/simulators/wells/MSWellHelpers.hpp @@ -38,12 +38,12 @@ namespace Opm { namespace mswellhelpers { -#if HAVE_UMFPACK /// Applies umfpack and checks for singularity template VectorType applyUMFPack(const MatrixType& D, std::shared_ptr >& linsolver, VectorType x) { +#if HAVE_UMFPACK if (!linsolver) { linsolver.reset(new Dune::UMFPack(D, 0)); @@ -69,28 +69,15 @@ namespace mswellhelpers } } return y; - } -#endif - - // obtain y = D^-1 * x with a direct solver - template - VectorType - invDXDirect(const MatrixType& D, VectorType x) - { -#if HAVE_UMFPACK - std::shared_ptr > linsolver; - return applyUMFPack(D, linsolver, x); #else // this is not thread safe - OPM_THROW(std::runtime_error, "Cannot use invDXDirect() without UMFPACK. " + OPM_THROW(std::runtime_error, "Cannot use aplyUMFPack() without UMFPACK. " "Reconfigure opm-simulator with SuiteSparse/UMFPACK support and recompile."); #endif // HAVE_UMFPACK } - - // obtain y = D^-1 * x with a BICSSTAB iterative solver template VectorType From ae8e2fb8c2f6a9ba4df0b842b478ac1cd7b58805 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 18 Aug 2020 17:53:37 +0200 Subject: [PATCH 5/9] Removed commented out code from MSWell. --- opm/simulators/wells/MultisegmentWell_impl.hpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 0189b2797..7b9622cb5 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -624,7 +624,6 @@ namespace Opm duneB_.mv(x, Bx); // invDBx = duneD^-1 * Bx_ - // const BVectorWell invDBx = mswellhelpers::invDXDirect(duneD_, Bx); const BVectorWell invDBx = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, Bx); // Ax = Ax - duneC_^T * invDBx @@ -1016,7 +1015,6 @@ namespace Opm { // We assemble the well equations, then we check the convergence, // which is why we do not put the assembleWellEq here. - //const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); const BVectorWell dx_well = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_); updateWellState(dx_well, well_state, deferred_logger); @@ -1114,14 +1112,12 @@ namespace Opm for (int seg = 0; seg < numberOfSegments(); ++seg) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { const int sign = dwells[seg][WFrac] > 0. ? 1 : -1; - // const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit); const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]) * relaxation_factor, dFLimit); primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited; } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { const int sign = dwells[seg][GFrac] > 0. ? 1 : -1; - // const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit); const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]) * relaxation_factor, dFLimit); primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited; } @@ -1132,7 +1128,6 @@ namespace Opm // update the segment pressure { const int sign = dwells[seg][SPres] > 0.? 1 : -1; - //const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change); const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change); primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5); } From 0c0be3a9411f63243d5a80f9605f5325072d3da4 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 18 Aug 2020 17:54:15 +0200 Subject: [PATCH 6/9] Reset duneDSolver only once. --- opm/simulators/wells/MultisegmentWell_impl.hpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 7b9622cb5..ab1ea3c26 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1922,7 +1922,6 @@ namespace Opm for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { duneD_[0][0][SPres][pv_idx] = control_eq.derivative(pv_idx + numEq); } - duneDSolver_.reset(); } @@ -2087,7 +2086,6 @@ namespace Opm if (accelerationalPressureLossConsidered()) { handleAccelerationPressureLoss(seg, well_state); } - duneDSolver_.reset(); } @@ -2182,7 +2180,6 @@ namespace Opm duneD_[seg][seg][SPres][GTotal] -= accelerationPressureLoss.derivative(GTotal + numEq); duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + numEq); duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + numEq); - duneDSolver_.reset(); } @@ -2516,6 +2513,8 @@ namespace Opm duneD_ = 0.0; resWell_ = 0.0; + duneDSolver_.reset(); + well_state.wellVaporizedOilRates()[index_of_well_] = 0.; well_state.wellDissolvedGasRates()[index_of_well_] = 0.; @@ -2659,7 +2658,6 @@ namespace Opm well_state.segPressDropFriction()[seg] + well_state.segPressDropAcceleration()[seg]; } - duneDSolver_.reset(); } @@ -3211,7 +3209,6 @@ namespace Opm duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq); } - duneDSolver_.reset(); } @@ -3253,7 +3250,6 @@ namespace Opm for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq); } - duneDSolver_.reset(); } From 94eb0d5ffe2e4043d090e47f6372fa58ec131146 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 19 Aug 2020 08:11:50 +0200 Subject: [PATCH 7/9] Fixes method name in throw message. --- opm/simulators/wells/MSWellHelpers.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index 8295a481a..c53f7887b 100644 --- a/opm/simulators/wells/MSWellHelpers.hpp +++ b/opm/simulators/wells/MSWellHelpers.hpp @@ -64,7 +64,7 @@ namespace mswellhelpers for (size_t i_block = 0; i_block < y.size(); ++i_block) { for (size_t i_elem = 0; i_elem < y[i_block].size(); ++i_elem) { if (std::isinf(y[i_block][i_elem]) || std::isnan(y[i_block][i_elem]) ) { - OPM_THROW(Opm::NumericalIssue, "nan or inf value found in invDXDirect due to singular matrix"); + OPM_THROW(Opm::NumericalIssue, "nan or inf value found after UMFPack solve due to singular matrix"); } } } From cca38ddc36102520ece3860cf1f7487607243963 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 20 Aug 2020 10:07:55 +0200 Subject: [PATCH 8/9] Fixed spelling in error message. --- opm/simulators/wells/MSWellHelpers.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index c53f7887b..226ca51e0 100644 --- a/opm/simulators/wells/MSWellHelpers.hpp +++ b/opm/simulators/wells/MSWellHelpers.hpp @@ -71,7 +71,7 @@ namespace mswellhelpers return y; #else // this is not thread safe - OPM_THROW(std::runtime_error, "Cannot use aplyUMFPack() without UMFPACK. " + OPM_THROW(std::runtime_error, "Cannot use applyUMFPack() without UMFPACK. " "Reconfigure opm-simulator with SuiteSparse/UMFPACK support and recompile."); #endif // HAVE_UMFPACK } From 49ed3d6468def4f1a7bdac68fc146750e156dee3 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 20 Aug 2020 10:13:53 +0200 Subject: [PATCH 9/9] Removes rest of the commented out code introduced with applyUMFPack. --- opm/simulators/wells/MultisegmentWell_impl.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index ab1ea3c26..f12d3dfac 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -640,7 +640,6 @@ namespace Opm apply(BVector& r) const { // invDrw_ = duneD^-1 * resWell_ - // const BVectorWell invDrw = mswellhelpers::invDXDirect(duneD_, resWell_); const BVectorWell invDrw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_); // r = r - duneC_^T * invDrw duneC_.mmtv(invDrw, r); @@ -1000,7 +999,6 @@ namespace Opm // resWell = resWell - B * x duneB_.mmv(x, resWell); // xw = D^-1 * resWell - //xw = mswellhelpers::invDXDirect(duneD_, resWell); xw = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell); } @@ -2398,7 +2396,6 @@ namespace Opm assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, deferred_logger); - // const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); const BVectorWell dx_well = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, resWell_); if (it > param_.strict_inner_iter_ms_wells_)