[cleanup][Wells] Consequently use SparseMatrixAdapter in addWellContributions.

This commit is contained in:
Robert Kloefkorn 2019-10-02 12:48:29 +02:00
parent 28cf1c17be
commit f708286111
8 changed files with 60 additions and 59 deletions

View File

@ -1,5 +1,6 @@
/* /*
Copyright 2016 IRIS AS Copyright 2016 IRIS AS
Copyright 2019 NORCE
This file is part of the Open Porous Media project (OPM). This file is part of the Open Porous Media project (OPM).
@ -457,25 +458,57 @@ namespace Opm
{ {
namespace Detail namespace Detail
{ {
//! calculates ret = A^T * B //! calculates ret = sign * (A^T * B)
template< class K, int m, int n, int p > //! TA, TB, and TC are not necessarily FieldMatrix, but those should
static inline void multMatrixTransposed(const Dune::FieldMatrix< K, n, m >& A, //! follow the Dune::DenseMatrix interface.
const Dune::FieldMatrix< K, n, p >& B, template< class TA, class TB, class TC, class PositiveSign >
Dune::FieldMatrix< K, m, p >& ret) static inline void multMatrixTransposedImpl ( const TA &A, // n x m
const TB &B, // n x p
TC &ret, // m x p
const PositiveSign )
{ {
typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type; 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 i = 0; i < m; ++i )
{ {
for( size_type j = 0; j < p; ++j ) for( size_type j = 0; j < p; ++j )
{ {
ret[ i ][ j ] = K( 0 ); K sum = 0;
for( size_type k = 0; k < n; ++k ) for( size_type k = 0; k < n; ++k )
ret[ i ][ j ] += A[ k ][ i ] * B[ k ][ j ]; {
sum += A[ k ][ i ] * B[ k ][ j ];
}
// set value depending on given sign
ret[ i ][ j ] = PositiveSign::value ? sum : -sum;
} }
} }
} }
//! calculates ret = A^T * B
template <class DenseMatrixA, class DenseMatrixB, class DenseMatrixC>
static inline void multMatrixTransposed(const DenseMatrixA& A,
const DenseMatrixB& B,
DenseMatrixC& ret)
{
multMatrixTransposedImpl( A, B, ret, std::true_type() );
}
//! calculates ret = -A^T * B
template <class DenseMatrixA, class DenseMatrixB, class DenseMatrixC>
static inline void negativeMultMatrixTransposed(const DenseMatrixA& A,
const DenseMatrixB& B,
DenseMatrixC& ret)
{
multMatrixTransposedImpl( A, B, ret, std::false_type() );
}
//! calculates ret = A * B //! calculates ret = A * B
template< class K> template< class K>
static inline void multMatrix(const Dune::DynamicMatrix<K>& A, static inline void multMatrix(const Dune::DynamicMatrix<K>& A,
@ -504,32 +537,6 @@ namespace Detail
} }
} }
//! calculates ret = A^T * B
template< class K, int m, int p >
static inline void multMatrixTransposed(const Dune::DynamicMatrix<K>& A,
const Dune::DynamicMatrix<K>& B,
Dune::FieldMatrix< K, m, p>& ret )
{
typedef typename Dune::DynamicMatrix<K> :: size_type size_type;
// A is a tranpose matrix
const size_type n = A.rows();
assert(m == A.cols() );
assert(n == B.rows() );
assert(p == B.cols() );
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 ];
}
}
}
} // namespace Detail } // namespace Detail
} // namespace Opm } // namespace Opm

View File

@ -100,7 +100,6 @@ namespace Opm {
#else #else
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType; typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
#endif #endif
typedef typename SparseMatrixAdapter::IstlMatrix Mat;
typedef Opm::BlackOilPolymerModule<TypeTag> PolymerModule; typedef Opm::BlackOilPolymerModule<TypeTag> PolymerModule;
@ -124,7 +123,7 @@ namespace Opm {
void applyInitial() void applyInitial()
{} {}
void linearize(SparseMatrixAdapter& mat , GlobalEqVector& res); void linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res);
void postSolve(GlobalEqVector& deltaX) void postSolve(GlobalEqVector& deltaX)
{ {
@ -218,10 +217,10 @@ namespace Opm {
const SimulatorReport& lastReport() const; const SimulatorReport& lastReport() const;
void addWellContributions(Mat& mat) const void addWellContributions(SparseMatrixAdapter& jacobian) const
{ {
for ( const auto& well: well_container_ ) { for ( const auto& well: well_container_ ) {
well->addWellContributions(mat); well->addWellContributions(jacobian);
} }
} }

View File

@ -119,7 +119,7 @@ namespace Opm {
template<typename TypeTag> template<typename TypeTag>
void void
BlackoilWellModel<TypeTag>:: BlackoilWellModel<TypeTag>::
linearize(SparseMatrixAdapter& mat , GlobalEqVector& res) linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res)
{ {
if (!localWellsActive()) if (!localWellsActive())
return; return;
@ -135,7 +135,7 @@ namespace Opm {
} }
for (const auto& well: well_container_) { for (const auto& well: well_container_) {
well->addWellContributions(mat.istlMatrix()); well->addWellContributions(jacobian);
// applying the well residual to reservoir residuals // applying the well residual to reservoir residuals
// r = r - duneC_^T * invDuneD_ * resWell_ // r = r - duneC_^T * invDuneD_ * resWell_

View File

@ -42,6 +42,7 @@ namespace Opm
using typename Base::MaterialLaw; using typename Base::MaterialLaw;
using typename Base::Indices; using typename Base::Indices;
using typename Base::RateConverterType; using typename Base::RateConverterType;
using typename Base::SparseMatrixAdapter;
/// the number of reservior equations /// the number of reservior equations
@ -73,7 +74,6 @@ namespace Opm
using typename Base::Scalar; using typename Base::Scalar;
/// the matrix and vector types for the reservoir /// the matrix and vector types for the reservoir
using typename Base::Mat;
using typename Base::BVector; using typename Base::BVector;
using typename Base::Eval; using typename Base::Eval;
@ -152,7 +152,7 @@ namespace Opm
const WellState& well_state, const WellState& well_state,
Opm::DeferredLogger& deferred_logger) override; // should be const? Opm::DeferredLogger& deferred_logger) override; // should be const?
virtual void addWellContributions(Mat& mat) const override; virtual void addWellContributions(SparseMatrixAdapter& jacobian) const override;
/// number of segments for this well /// number of segments for this well
/// int number_of_segments_; /// int number_of_segments_;

View File

@ -885,7 +885,7 @@ namespace Opm
template<typename TypeTag> template<typename TypeTag>
void void
MultisegmentWell<TypeTag>:: MultisegmentWell<TypeTag>::
addWellContributions(Mat& /* mat */) const addWellContributions(SparseMatrixAdapter& /* jacobian */) const
{ {
OPM_THROW(std::runtime_error, "addWellContributions is not supported by multisegment well yet"); OPM_THROW(std::runtime_error, "addWellContributions is not supported by multisegment well yet");
} }

View File

@ -58,6 +58,7 @@ namespace Opm
using typename Base::ModelParameters; using typename Base::ModelParameters;
using typename Base::Indices; using typename Base::Indices;
using typename Base::RateConverterType; using typename Base::RateConverterType;
using typename Base::SparseMatrixAdapter;
using Base::numEq; using Base::numEq;
@ -110,7 +111,6 @@ namespace Opm
using Base::Oil; using Base::Oil;
using Base::Gas; using Base::Gas;
using typename Base::Mat;
using typename Base::BVector; using typename Base::BVector;
using typename Base::Eval; using typename Base::Eval;
@ -191,7 +191,7 @@ namespace Opm
const WellState& well_state, const WellState& well_state,
Opm::DeferredLogger& deferred_logger) override; // should be const? Opm::DeferredLogger& deferred_logger) override; // should be const?
virtual void addWellContributions(Mat& mat) const override; virtual void addWellContributions(SparseMatrixAdapter& mat) const override;
/// \brief Wether the Jacobian will also have well contributions in it. /// \brief Wether the Jacobian will also have well contributions in it.
virtual bool jacobianContainsWellContributions() const override virtual bool jacobianContainsWellContributions() const override

View File

@ -2751,9 +2751,12 @@ namespace Opm
} }
} }
template<typename TypeTag> template<typename TypeTag>
void void
StandardWell<TypeTag>::addWellContributions(Mat& mat) const StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian) const
{ {
// We need to change matrx A as follows // We need to change matrx A as follows
// A -= C^T D^-1 B // A -= C^T D^-1 B
@ -2761,24 +2764,17 @@ namespace Opm
// B and C have 1 row, nc colums and nonzero // B and C have 1 row, nc colums and nonzero
// at (0,j) only if this well has a perforation at cell j. // at (0,j) only if this well has a perforation at cell j.
typename SparseMatrixAdapter::MatrixBlock tmpMat;
Dune::DynamicMatrix<Scalar> tmp;
for ( auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC ) for ( auto colC = duneC_[0].begin(), endC = duneC_[0].end(); colC != endC; ++colC )
{ {
const auto row_index = colC.index(); const auto row_index = colC.index();
auto& row = mat[row_index];
auto col = row.begin();
for ( auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB ) for ( auto colB = duneB_[0].begin(), endB = duneB_[0].end(); colB != endB; ++colB )
{ {
const auto col_index = colB.index();
// Move col to index col_index
while ( col != row.end() && col.index() < col_index ) ++col;
assert(col != row.end() && col.index() == col_index);
Dune::DynamicMatrix<Scalar> tmp;
Detail::multMatrix(invDuneD_[0][0], (*colB), tmp); Detail::multMatrix(invDuneD_[0][0], (*colB), tmp);
typename Mat::block_type tmp1; Detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat);
Detail::multMatrixTransposed((*colC), tmp, tmp1); jacobian.addToBlock( row_index, colB.index(), tmpMat );
(*col) -= tmp1;
} }
} }
} }

View File

@ -88,7 +88,6 @@ namespace Opm
typedef Dune::FieldVector<Scalar, numEq > VectorBlockType; typedef Dune::FieldVector<Scalar, numEq > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType; typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
typedef typename SparseMatrixAdapter::IstlMatrix Mat;
typedef Dune::BlockVector<VectorBlockType> BVector; typedef Dune::BlockVector<VectorBlockType> BVector;
typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval; typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
@ -206,7 +205,7 @@ namespace Opm
void calculateReservoirRates(WellState& well_state) const; void calculateReservoirRates(WellState& well_state) const;
// Add well contributions to matrix // Add well contributions to matrix
virtual void addWellContributions(Mat&) const = 0; virtual void addWellContributions(SparseMatrixAdapter&) const = 0;
void addCellRates(RateVector& rates, int cellIdx) const; void addCellRates(RateVector& rates, int cellIdx) const;