OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

add and use class wellModelMatrixAdapter

Browse Source 
The well contribution is substracted in the MatrixAdapter
i.e. Ax - BinvDCx
This commit is contained in:
Tor Harald Sandve
2016-09-07 13:23:52 +02:00
parent a4dcc4b13d
commit 746f9a56cc
2 changed files with 64 additions and 77 deletions
Download Patch File Download Diff File Expand all files Collapse all files

54
opm/autodiff/BlackoilModelEbos.hpp
View File

@@ -98,6 +98,7 @@ namespace Opm {
/// upwind weighting of mobilities.
class BlackoilModelEbos
{
typedef BlackoilModelEbos ThisType;
public:
// --------- Types and enums ---------
typedef AutoDiffBlock<double> ADB;
@@ -362,6 +363,12 @@ namespace Opm {
return linsolver_.iterations();
}
template <class X, class Y>
void applyWellModel(const X& x, Y& y )
{
wellModel().apply(x, y);
}
/// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual.
@@ -374,17 +381,19 @@ namespace Opm {
typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> BVector;
typedef Dune::MatrixAdapter<Mat,BVector,BVector> Operator;
auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
auto& ebosResid = ebosSimulator_.model().linearizer().residual();
Operator opA(ebosJac);
typedef WellModelMatrixAdapter<Mat,BVector,BVector, ThisType> Operator;
Operator opA(ebosJac, const_cast< ThisType& > (*this));
const double relax = 0.9;
typedef Dune::SeqILU0<Mat, BVector, BVector> SeqPreconditioner;
SeqPreconditioner precond(opA.getmat(), relax);
Dune::SeqScalarProduct<BVector> sp;
wellModel().apply(ebosJac, ebosResid);
// apply well residual to the residual.
wellModel().apply(ebosResid);
Dune::BiCGSTABSolver<BVector> linsolve(opA, sp, precond,
0.01,
100,
@@ -393,10 +402,49 @@ namespace Opm {
x = 0.0;
linsolve.apply(x, ebosResid, result);
// recover wells.
xw = 0.0;
wellModel().recoverVariable(x, xw);
}
//=====================================================================
// Implementation for ISTL-matrix based operator
//=====================================================================
/*!
\brief Adapter to turn a matrix into a linear operator.
Adapts a matrix to the assembled linear operator interface
*/
template<class M, class X, class Y, class WellModel>
class WellModelMatrixAdapter : public Dune::MatrixAdapter<M,X,Y>
{
typedef Dune::MatrixAdapter<M,X,Y> BaseType;
public:
//! export types
typedef M matrix_type;
typedef X domain_type;
typedef Y range_type;
typedef typename X::field_type field_type;
//! constructor: just store a reference to a matrix
explicit WellModelMatrixAdapter (const M& A, WellModel& wellMod ) : BaseType( A ), wellMod_( wellMod ) {}
//! apply operator to x: \f$ y = A(x) \f$
virtual void apply (const X& x, Y& y) const
{
BaseType::apply( x, y );
wellMod_.applyWellModel(x, y );
}
private:
WellModel& wellMod_;
};
/** @} end documentation */
/// Apply an update to the primary variables, chopped if appropriate.
/// \param[in] dx updates to apply to primary variables
/// \param[in, out] reservoir_state reservoir state variables

87
opm/autodiff/StandardWellsDense.hpp
View File

@@ -238,87 +238,26 @@ namespace Opm {
}
}
Mat matAdd( const Mat& A, const Mat& B ) const
{
return matBinaryOp( A, B, true );
}
Mat matSubstract( const Mat& A, const Mat& B ) const
{
return matBinaryOp( A, B, false );
}
Mat matBinaryOp( const Mat& A, const Mat& B, const bool add ) const
{
const int avg_cols_per_row = 20;
const double overflow_fraction = 0.4;
Mat res( A.N(), A.M(), avg_cols_per_row,overflow_fraction, Mat::implicit );
assert( A.N() == B.N() && A.M() == B.M() );
// res = A
for( auto rowit = A.begin(), rowEnd = A.end(); rowit != rowEnd; ++rowit )
{
const int rowIdx = rowit.index();
const auto colEnd = rowit->end();
for( auto colit = rowit->begin(); colit != colEnd; ++colit )
{
const int colIdx = colit.index();
res.entry( rowIdx, colIdx ) = (*colit);
}
}
// res += B
for( auto rowit = B.begin(), rowEnd = B.end(); rowit != rowEnd; ++rowit )
{
const int rowIdx = rowit.index();
const auto colEnd = rowit->end();
for( auto colit = rowit->begin(); colit != colEnd; ++colit )
{
const int colIdx = colit.index();
if( add )
{
res.entry( rowIdx, colIdx ) += (*colit);
}
else
{
res.entry( rowIdx, colIdx ) -= (*colit);
}
}
}
res.compress();
return res;
}
void addRhs(BVector& x, Mat& jac) const {
assert(x.size() == rhs.size());
x += rhs_;
// jac = A + duneA
//jac = matAdd( jac, duneA_ );
jac += duneA_;
}
void apply( Mat& A,
BVector& res) const {
// substract Binv(D)rw from r;
void apply( BVector& r) const {
BVector invDrw(invDuneD_.N());
invDuneD_.mv(resWell_,invDrw);
duneB_.mmv(invDrw, r);
}
Mat BmultinvD;
Mat duneA;
Dune::matMultMat(BmultinvD, duneB_ , invDuneD_);
Dune::matMultMat(duneA, BmultinvD, duneC_);
//std::cout << "before" << std::endl;
//std::cout << "A" << std::endl;
//print(A);
//std::cout << "duneA" << std::endl;
//print(duneA);
// A = E - duneA
A = matSubstract( A, duneA );
//A -= duneA;
//std::cout << "after" << std::endl;
//print(A);
BmultinvD.mmv(resWell_, res);
// subtract B*inv(D)*C * x from A*x
void apply(const BVector& x, BVector& Ax) {
BVector Cx(duneC_.N());
duneC_.mv(x, Cx);
BVector invDCx(invDuneD_.N());
invDuneD_.mv(Cx, invDCx);
duneB_.mmv(invDCx,Ax);
}
void recoverVariable(const BVector& x, BVector& xw) const {