Merge pull request #2882 from atgeirr/msw-add-well-contrib

Implement addWellContributions() for multi-segment wells

opm/simulators/linalg/MatrixBlock.hpp

@@ -458,6 +458,37 @@ namespace Opm
 {
 namespace Detail
 {
+    //! calculates ret = A * B
+    template< class TA, class TB, class TC, class PositiveSign >
+    static inline void multMatrixImpl( const TA &A, // n x m
+                                       const TB &B, // n x p
+                                       TC &ret,     // m x p
+                                       const PositiveSign )
+    {
+        typedef typename TA :: size_type size_type;
+        typedef typename TA :: field_type K;
+        assert( A.N() == B.N() );
+        assert( A.M() == ret.N() );
+        assert( B.M() == ret.M() );
+
+        const size_type n = A.N();
+        const size_type m = ret.N();
+        const size_type p = B.M();
+        for( size_type i = 0; i < m; ++i )
+        {
+            for( size_type j = 0; j < p; ++j )
+            {
+                K sum = 0;
+                for( size_type k = 0; k < n; ++k )
+                {
+                    sum += A[ i ][ k ] * B[ k ][ j ];
+                }
+                // set value depending on given sign
+                ret[ i ][ j ] = PositiveSign::value ? sum : -sum;
+            }
+        }
+    }
+
     //! calculates ret = sign * (A^T * B)
     //! TA, TB, and TC are not necessarily FieldMatrix, but those should
     //! follow the Dune::DenseMatrix interface.

opm/simulators/wells/MSWellHelpers.hpp

@@ -72,7 +72,45 @@ namespace mswellhelpers
 #else
         // this is not thread safe
         OPM_THROW(std::runtime_error, "Cannot use applyUMFPack() without UMFPACK. "
-                  "Reconfigure opm-simulator with SuiteSparse/UMFPACK support and recompile.");
+                  "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile.");
+#endif // HAVE_UMFPACK
+    }
+
+
+
+    /// Applies umfpack and checks for singularity
+    template <typename MatrixType, typename VectorType>
+    Dune::Matrix<typename MatrixType::block_type>
+    invertWithUMFPack(const MatrixType& D, std::shared_ptr<Dune::UMFPack<MatrixType> >& linsolver)
+    {
+#if HAVE_UMFPACK
+        const int sz = D.M();
+        const int bsz = D[0][0].M();
+        VectorType e(sz);
+        e = 0.0;
+
+        // Make a full block matrix.
+        Dune::Matrix<typename MatrixType::block_type> inv(sz, sz);
+
+        // Create inverse by passing basis vectors to the solver.
+        for (int ii = 0; ii < sz; ++ii) {
+            for (int jj = 0; jj < bsz; ++jj) {
+                e[ii][jj] = 1.0;
+                auto col = applyUMFPack(D, linsolver, e);
+                for (int cc = 0; cc < sz; ++cc) {
+                    for (int dd = 0; dd < bsz; ++dd) {
+                        inv[cc][ii][dd][jj] = col[cc][dd];
+                    }
+                }
+                e[ii][jj] = 0.0;
+            }
+        }
+
+        return inv;
+#else
+        // this is not thread safe
+        OPM_THROW(std::runtime_error, "Cannot use invertWithUMFPack() without UMFPACK. "
+                  "Reconfigure opm-simulators with SuiteSparse/UMFPACK support and recompile.");
 #endif // HAVE_UMFPACK
     }
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -619,6 +619,11 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     apply(const BVector& x, BVector& Ax) const
     {
+        if ( param_.matrix_add_well_contributions_ )
+        {
+            // Contributions are already in the matrix itself
+            return;
+        }
         BVectorWell Bx(duneB_.N());
 
         duneB_.mv(x, Bx);
@@ -1174,9 +1179,33 @@ namespace Opm
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    addWellContributions(SparseMatrixAdapter& /* jacobian */) const
+    addWellContributions(SparseMatrixAdapter& jacobian) const
     {
-        OPM_THROW(std::runtime_error, "addWellContributions is not supported by multisegment well yet");
+        const auto invDuneD = mswellhelpers::invertWithUMFPack<DiagMatWell, BVectorWell>(duneD_, duneDSolver_);
+
+        // 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.
+        // B and C are (nseg x ncells) block matrices with (4 x 4 blocks).
+        // They have nonzeros at (i, j) only if this well has a
+        // perforation at cell j connected to segment i.  The code
+        // assumes that no cell is connected to more than one segment,
+        // i.e. the columns of B/C have no more than one nonzero.
+        for (int rowC = 0; rowC < duneC_.N(); ++rowC) {
+            for (auto colC = duneC_[rowC].begin(), endC = duneC_[rowC].end(); colC != endC; ++colC) {
+                const auto row_index = colC.index();
+                for (int rowB = 0; rowB < duneB_.N(); ++rowB) {
+                    for (auto colB = duneB_[rowB].begin(), endB = duneB_[rowB].end(); colB != endB; ++colB) {
+                        const auto col_index = colB.index();
+                        OffDiagMatrixBlockWellType tmp1;
+                        Detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type());
+                        typename SparseMatrixAdapter::MatrixBlock tmp2;
+                        Detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type());
+                        jacobian.addToBlock(row_index, col_index, tmp2);
+                    }
+                }
+            }
+        }
     }