Merge pull request #2718 from blattms/reuse-umfpack-ms

Reuse umfpack decomposition for multisegment wells.

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 <typename MatrixType, typename VectorType>
     VectorType
-    invDXDirect(const MatrixType& D, VectorType x)
+    applyUMFPack(const MatrixType& D, std::shared_ptr<Dune::UMFPack<MatrixType> >& linsolver, VectorType x)
     {
 #if HAVE_UMFPACK
+        if (!linsolver)
+        {
+            linsolver.reset(new Dune::UMFPack<MatrixType>(D, 0));
+        }
+
+        // The copy of x seems mandatory for calling UMFPack!
         VectorType y(x.size());
         y = 0.;
 
-        Dune::UMFPack<MatrixType> 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 <typename MatrixType, typename VectorType>
     VectorType

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<Dune::UMFPack<DiagMatWell> > duneDSolver_;
 
         // residuals of the well equations
         mutable BVectorWell resWell_;

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.;