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
This commit is contained in:
Markus Blatt 2017-05-15 14:46:43 +02:00
parent d0fcd1a8e9
commit 799cbb4b62
7 changed files with 124 additions and 12 deletions

View File

@ -358,8 +358,7 @@ namespace Opm {
{
// This might be dangerous?!
ebosSimulator_.model().clearAuxiliaryModules();
const auto* wells = wellModel().wellsPointer();
auto auxMod = std::make_shared<WellConnectionAuxiliaryModule<TypeTag> >(wells);
auto auxMod = std::make_shared<WellConnectionAuxiliaryModule<TypeTag> >( wellModel().wells() );
ebosSimulator_.model().addAuxiliaryModule(auxMod);
}
@ -498,14 +497,17 @@ namespace Opm {
if( isParallel() )
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, true > Operator;
Operator opA(ebosJac, wellModel(), istlSolver().parallelInformation() );
Operator opA(ebosJac, wellModel(),
param_.matrix_add_well_contributions_,
istlSolver().parallelInformation() );
assert( opA.comm() );
istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
}
else
{
typedef WellModelMatrixAdapter< Mat, BVector, BVector, BlackoilWellModel<TypeTag>, false > Operator;
Operator opA(ebosJac, wellModel());
Operator opA(ebosJac, wellModel(),
param_.matrix_add_well_contributions_ );
istlSolver().solve( opA, x, ebosResid );
}
}
@ -552,8 +554,11 @@ namespace Opm {
#endif
//! constructor: just store a reference to a matrix
WellModelMatrixAdapter (const M& A, const WellModel& wellMod, const boost::any& parallelInformation = boost::any() )
: A_( A ), wellMod_( wellMod ), comm_()
WellModelMatrixAdapter (const M& A, const 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) )
@ -568,8 +573,12 @@ namespace Opm {
virtual void apply( const X& x, Y& y ) const
{
A_.mv( x, y );
// add well model modification to y
wellMod_.apply(x, y );
if ( ! matrix_add_well_contributions_ )
{
// add well model modification to y
wellMod_.apply(x, y );
}
#if HAVE_MPI
if( comm_ )
@ -581,8 +590,12 @@ namespace Opm {
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_.applyScaleAdd( alpha, x, y );
if ( ! matrix_add_well_contributions_ )
{
// add scaled well model modification to y
wellMod_.applyScaleAdd( alpha, x, y );
}
#if HAVE_MPI
if( comm_ )
@ -601,6 +614,7 @@ namespace Opm {
const matrix_type& A_ ;
const WellModel& wellMod_;
std::unique_ptr< communication_type > comm_;
bool matrix_add_well_contributions_;
};
/// Apply an update to the primary variables, chopped if appropriate.

View File

@ -61,9 +61,14 @@
namespace Opm {
template<typename TypeTag>
class BlackoilModelEbos;
/// Class for handling the blackoil well model.
template<typename TypeTag>
class BlackoilWellModel {
// Needs acces to wells_ecl_
friend class BlackoilModelEbos<TypeTag>;
public:
// --------- Types ---------
typedef WellStateFullyImplicitBlackoil WellState;
@ -84,6 +89,15 @@ namespace Opm {
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::BlockVector<VectorBlockType> BVector;
#if DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2 , 5, 1)
// 3x3 matrix block inversion was unstable from at least 2.3 until and
// including 2.5.0
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
#else
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
#endif
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
// For the conversion between the surface volume rate and resrevoir voidage rate
@ -143,6 +157,13 @@ namespace Opm {
const SimulatorReport& lastReport() 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;
protected:
Simulator& ebosSimulator_;

View File

@ -295,8 +295,16 @@ namespace Opm {
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
addWellContributions(Mat& mat) const
{
for(const auto& well: well_container_)
{
well->addWellContributions(mat);
}
}
// applying the well residual to reservoir residuals
// r = r - duneC_^T * invDuneD_ * resWell_

View File

@ -157,6 +157,7 @@ namespace Opm
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state); // should be const?
virtual void addWellContributions(Mat& mat) const;
protected:
// protected functions from the Base class

View File

@ -22,6 +22,28 @@
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, n, m > &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 ];
}
}
}
}
template<typename TypeTag>
StandardWell<TypeTag>::
StandardWell(const Well* well, const int time_step, const Wells* wells,
@ -667,6 +689,11 @@ namespace Opm
// do the local inversion of D.
invDuneD_[0][0].invert();
if ( param_.matrix_add_well_contributions_ )
{
addWellContributions( ebosJac );
}
}
@ -2073,4 +2100,30 @@ namespace Opm
mob[waterCompIdx] /= shear_factor;
}
}
template<typename TypeTag>
void
StandardWell<TypeTag>::addWellContributions(Mat& mat) const
{
auto colC = duneC_[0].begin();
for(int perf1 = 0; perf1 < number_of_perforations_; ++perf1, ++colC) {
const auto row_index = well_cells_[perf1];
auto& row = mat[row_index];
auto colB = duneB_[0].begin();
assert(colC.index() == row_index);
for(int perf2 = 0 ; perf2 < number_of_perforations_; ++perf2, ++colB) {
const auto col_index = well_cells_[perf2];
auto col = row.find(col_index);
assert(col != row.end());
assert(colB.index() == col_index);
Dune::FieldMatrix<Scalar, numWellEq, numEq> tmp;
typename Mat::block_type tmp1;
Dune::FMatrixHelp::multMatrix(invDuneD_[0][0], (*colB), tmp);
Detail::multMatrixTransposed((*colC), tmp, tmp1);
(*col) -= tmp1;
}
}
}
}

View File

@ -82,7 +82,13 @@ namespace Opm
typedef double Scalar;
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
#if DUNE_VERSION_NEWER_REV(DUNE_ISTL, 2 , 5, 1)
// 3x3 matrix block inversion was unstable from at least 2.3 until and
// including 2.5.0
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
#else
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
#endif
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
@ -202,6 +208,8 @@ namespace Opm
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state) = 0; // should be const?
virtual void addWellContributions(Mat& mat) const;
// updating the voidage rates in well_state when requested
void calculateReservoirRates(WellState& well_state) const;

View File

@ -844,6 +844,13 @@ namespace Opm
return 1.0;
}
template<typename TypeTag>
void
WellInterface<TypeTag>::addWellContributions(Mat& mat) const
{
OPM_THROW(NotImplemented, "This well class does not support adding well contributions"
<< "to the matrix");
}