Add influences by well perforations to matrix.

This is only done upon request and uses the auxiliary module approach
provided by ewoms.

In the case of adding the influences we do not execute applyWellModelScaleAdd
or applyWellModel in the operator

opm/autodiff/BlackoilModelEbos.hpp

@@ -479,14 +479,17 @@ namespace Opm {
             if( isParallel() )
             {
                 typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, true > Operator;
-                Operator opA(ebosJac, const_cast< ThisType& > (*this), istlSolver().parallelInformation() );
+                Operator opA(ebosJac, const_cast< ThisType& > (*this),
+                             param_.matrix_add_well_contributions_,
+                             istlSolver().parallelInformation() );
                 assert( opA.comm() );
                 istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
             }
             else
             {
                 typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, false > Operator;
-                Operator opA(ebosJac, const_cast< ThisType& > (*this) );
+                Operator opA(ebosJac, const_cast< ThisType& > (*this),
+                             param_.matrix_add_well_contributions_ );
                 istlSolver().solve( opA, x, ebosResid );
             }
 
@@ -532,8 +535,11 @@ namespace Opm {
           };
 
           //! constructor: just store a reference to a matrix
-          WellModelMatrixAdapter (const M& A, WellModel& wellMod, const boost::any& parallelInformation = boost::any() )
+          WellModelMatrixAdapter (const M& A, WellModel& wellMod,
-              : A_( A ), wellMod_( wellMod ), comm_()
+                                  bool matrix_add_well_contributions,
+                                  const boost::any& parallelInformation = boost::any() )
+              : A_( A ), wellMod_( wellMod ), comm_(),
+                matrix_add_well_contributions_(matrix_add_well_contributions)
           {
 #if HAVE_MPI
             if( parallelInformation.type() == typeid(ParallelISTLInformation) )
@@ -548,22 +554,29 @@ namespace Opm {
           virtual void apply( const X& x, Y& y ) const
           {
             A_.mv( x, y );
-            // add well model modification to y
+
-            wellMod_.applyWellModelAdd(x, y );
+            if ( ! matrix_add_well_contributions_ )
+            {
+              // add well model modification to y
+              wellMod_.applyWellModelAdd(x, y );
 
 #if HAVE_MPI
-            if( comm_ )
+              if( comm_ )
-              comm_->project( y );
+                comm_->project( y );
 #endif
+            }
           }
 
           // y += \alpha * A * x
           virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
           {
             A_.usmv(alpha,x,y);
-            // add scaled well model modification to y
-            wellMod_.applyWellModelScaleAdd( alpha, x, y );
 
+            if ( ! matrix_add_well_contributions_ )
+            {
+                // add scaled well model modification to y
+                wellMod_.applyWellModelScaleAdd( alpha, x, y );
+            }
 #if HAVE_MPI
             if( comm_ )
               comm_->project( y );
@@ -581,6 +594,7 @@ namespace Opm {
           const matrix_type& A_ ;
           WellModel& wellMod_;
           std::unique_ptr< communication_type > comm_;
+          bool matrix_add_well_contributions_;
         };
 
         /// Apply an update to the primary variables, chopped if appropriate.

opm/autodiff/StandardWellsDense.hpp

@@ -58,6 +58,27 @@
 #include <opm/simulators/WellSwitchingLogger.hpp>
 
 namespace Opm {
+namespace Detail
+{
+    //! calculates ret = A^T * B
+    template< class K, int m, int n, int p >
+    static inline void multMatrixTransposed ( const Dune::FieldMatrix< K, m, n > &A,
+                                              const Dune::FieldMatrix< K, n, p > &B,
+                                              Dune::FieldMatrix< K, m, p > &ret )
+    {
+        typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type;
+
+        for( size_type i = 0; i < m; ++i )
+        {
+            for( size_type j = 0; j < p; ++j )
+            {
+                ret[ i ][ j ] = K( 0 );
+                for( size_type k = 0; k < n; ++k )
+                    ret[ i ][ j ] += A[ k ][ i ] * B[ k ][ j ];
+            }
+        }
+    }
+}
 
 enum WellVariablePositions {
     XvarWell = 0,
@@ -338,6 +359,12 @@ enum WellVariablePositions {
                                      const WellState& well_state,
                                      const int well_number) const;
 
+            /// \! brief Modifies matrix to include influences of the well perforations.
+            ///
+            /// \param mat The linear system with the assembled mass balance
+            ///            equations
+            void addWellContributions(Mat& mat) const;
+
             using WellMapType = typename WellState::WellMapType;
             using WellMapEntryType = typename WellState::mapentry_t;
 

opm/autodiff/StandardWellsDense_impl.hpp

@@ -116,6 +116,36 @@ namespace Opm {
     }
 
 
+    template<class TypeTag>
+    void
+    StandardWellsDense<TypeTag>::
+    addWellContributions(Mat& mat) const
+    {
+        const int nw = wells().number_of_wells;
+
+        for(int w = 0; w < nw; ++w)
+        {
+            for(int perf1 = wells().well_connpos[w] ; perf1 < wells().well_connpos[w+1]; ++perf1) {
+                const auto pi = wells().well_cells[perf1];
+                auto& row = mat[pi];
+
+                for(int perf2 = wells().well_connpos[w] ; perf2 < wells().well_connpos[w+1]; ++perf2) {
+                    const auto pj = wells().well_cells[perf2];
+                    auto col = row.find(pj);
+                    assert(col != row.end());
+                    auto colC = duneC_[w].find(pj);
+                    assert(colC !=  duneC_[w].end());
+                    auto colB = duneB_[w].find(pi);
+                    assert(colB !=  duneB_[w].end());
+
+                    typename Mat::block_type tmp, tmp1;
+                    Dune::FMatrixHelp::multMatrix(invDuneD_[w][w],  (*colC), tmp);
+                    Detail::multMatrixTransposed((*colB), tmp, tmp1);
+                    (*col) -= tmp1;
+                }
+            }
+        }
+    }
 
     template<typename TypeTag>
     SimulatorReport
@@ -243,6 +273,11 @@ namespace Opm {
 
         // do the local inversion of D.
         localInvert( invDuneD_ );
+
+        if ( param_.matrix_add_well_contributions_ )
+        {
+            addWellContributions( ebosJac );
+        }
     }