From f913ed10fd4ed63d3731f3705fa0ed293819d6a2 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Mon, 15 May 2017 14:46:43 +0200 Subject: [PATCH] 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 | 34 ++++++++++++++++------- opm/autodiff/StandardWellsDense.hpp | 27 ++++++++++++++++++ opm/autodiff/StandardWellsDense_impl.hpp | 35 ++++++++++++++++++++++++ 3 files changed, 86 insertions(+), 10 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index c3435e589..78ebb2cde 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/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() ) - : A_( A ), wellMod_( wellMod ), comm_() + WellModelMatrixAdapter (const M& A, WellModel& wellMod, + 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_ ) - comm_->project( y ); + if( comm_ ) + 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. diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 5c2906a49..837e2d52c 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -58,6 +58,27 @@ #include 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; diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index d1af3d5dd..e29964b05 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -116,6 +116,36 @@ namespace Opm { } + template + void + StandardWellsDense:: + 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 SimulatorReport @@ -243,6 +273,11 @@ namespace Opm { // do the local inversion of D. localInvert( invDuneD_ ); + + if ( param_.matrix_add_well_contributions_ ) + { + addWellContributions( ebosJac ); + } }