diff --git a/opm/simulators/wells/MSWellHelpers.hpp b/opm/simulators/wells/MSWellHelpers.hpp index a3b0cfcc2..226ca51e0 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,45 +37,47 @@ namespace Opm { namespace mswellhelpers { - // obtain y = D^-1 * x with a direct solver + + /// Applies umfpack and checks for singularity template VectorType - invDXDirect(const MatrixType& D, VectorType x) + applyUMFPack(const MatrixType& D, std::shared_ptr >& linsolver, VectorType x) { #if HAVE_UMFPACK + if (!linsolver) + { + linsolver.reset(new Dune::UMFPack(D, 0)); + } + + // 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; // 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 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"); } } } - return y; #else // this is not thread safe - OPM_THROW(std::runtime_error, "Cannot use invDXDirect() without UMFPACK. " + OPM_THROW(std::runtime_error, "Cannot use applyUMFPack() 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 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..f12d3dfac 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -624,7 +624,7 @@ 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 duneC_.mmtv(invDBx,Ax); @@ -640,7 +640,7 @@ 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); } @@ -999,7 +999,7 @@ 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); } @@ -1013,7 +1013,7 @@ 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); } @@ -1110,14 +1110,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; } @@ -1128,7 +1126,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); } @@ -2399,7 +2396,7 @@ 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_) relax_convergence = true; @@ -2513,6 +2510,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.;