OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-02 05:49:09 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2613 from blattms/use-correct-matrix

Browse Source 
Use correct matrix for scaling and get rid of unused preconditioner matrix member
This commit is contained in:
Markus Blatt 2020-05-13 14:54:07 +02:00 committed by GitHub
commit 90a0e2b392
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 21 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@ -126,10 +126,9 @@ public:
//! constructor: just store a reference to a matrix //! constructor: just store a reference to a matrix
WellModelMatrixAdapter (const M& A, WellModelMatrixAdapter (const M& A,
const M& A_for_precond,
const WellModel& wellMod, const WellModel& wellMod,
const std::shared_ptr< communication_type >& comm = std::shared_ptr< communication_type >()) const std::shared_ptr< communication_type >& comm = std::shared_ptr< communication_type >())
: A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), comm_(comm) : A_( A ), wellMod_( wellMod ), comm_(comm)
{} {}
@ -160,11 +159,10 @@ public:
#endif #endif
} }
virtual const matrix_type& getmat() const override { return A_for_precond_; } virtual const matrix_type& getmat() const override { return A_; }
protected: protected:
const matrix_type& A_ ; const matrix_type& A_ ;
const matrix_type& A_for_precond_ ;
const WellModel& wellMod_; const WellModel& wellMod_;
std::shared_ptr< communication_type > comm_; std::shared_ptr< communication_type > comm_;
}; };
@ -201,10 +199,9 @@ public:
//! constructor: just store a reference to a matrix //! constructor: just store a reference to a matrix
WellModelGhostLastMatrixAdapter (const M& A, WellModelGhostLastMatrixAdapter (const M& A,
const M& A_for_precond,
const WellModel& wellMod, const WellModel& wellMod,
const size_t interiorSize ) const size_t interiorSize )
: A_( A ), A_for_precond_(A_for_precond), wellMod_( wellMod ), interiorSize_(interiorSize) : A_( A ), wellMod_( wellMod ), interiorSize_(interiorSize)
{} {}
virtual void apply( const X& x, Y& y ) const override virtual void apply( const X& x, Y& y ) const override
@ -238,7 +235,7 @@ public:
ghostLastProject( y ); ghostLastProject( y );
} }
virtual const matrix_type& getmat() const override { return A_for_precond_; } virtual const matrix_type& getmat() const override { return A_; }
protected: protected:
void ghostLastProject(Y& y) const void ghostLastProject(Y& y) const
@ -249,7 +246,6 @@ protected:
} }
const matrix_type& A_ ; const matrix_type& A_ ;
const matrix_type& A_for_precond_ ;
const WellModel& wellMod_; const WellModel& wellMod_;
size_t interiorSize_; size_t interiorSize_;
}; };
@ -380,7 +376,6 @@ protected:
// nothing to clean here // nothing to clean here
void eraseMatrix() { void eraseMatrix() {
matrix_for_preconditioner_.reset();
} }
void prepare(const SparseMatrixAdapter& M, Vector& b) void prepare(const SparseMatrixAdapter& M, Vector& b)
@ -507,7 +502,8 @@ protected:
{ {
if ( ownersFirst_ && !parameters_.linear_solver_use_amg_ && !useFlexible_) { if ( ownersFirst_ && !parameters_.linear_solver_use_amg_ && !useFlexible_) {
typedef WellModelGhostLastMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator; typedef WellModelGhostLastMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
Operator opA(*matrix_, *matrix_, wellModel, interiorCellNum_); assert(matrix_);
Operator opA(*matrix_, wellModel, interiorCellNum_);
solve( opA, x, *rhs_, *comm_ ); solve( opA, x, *rhs_, *comm_ );
} }
@ -516,7 +512,7 @@ protected:
typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator; typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator;
assert (noGhostMat_); assert (noGhostMat_);
copyJacToNoGhost(*matrix_, *noGhostMat_); copyJacToNoGhost(*matrix_, *noGhostMat_);
Operator opA(*noGhostMat_, *noGhostMat_, wellModel, Operator opA(*noGhostMat_, wellModel,
comm_ ); comm_ );
solve( opA, x, *rhs_, *comm_ ); solve( opA, x, *rhs_, *comm_ );
@ -526,7 +522,7 @@ protected:
else else
{ {
typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false > Operator; typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false > Operator;
Operator opA(*matrix_, *matrix_, wellModel); Operator opA(*matrix_, wellModel);
solve( opA, x, *rhs_ ); solve( opA, x, *rhs_ );
} }
@ -620,7 +616,7 @@ protected:
if (!useWellConn_) { if (!useWellConn_) {
simulator_.problem().wellModel().getWellContributions(wellContribs); simulator_.problem().wellModel().getWellContributions(wellContribs);
} }
bdaBridge->solve_system(matrix_, istlb, wellContribs, result); bdaBridge->solve_system(&getMatrix(), istlb, wellContribs, result);
if (result.converged) { if (result.converged) {
// get result vector x from non-Dune backend, iff solve was successful // get result vector x from non-Dune backend, iff solve was successful
bdaBridge->get_result(x); bdaBridge->get_result(x);
@ -875,7 +871,7 @@ protected:
weightsCalculator = weightsCalculator =
[this, transpose, pressureIndex](){ [this, transpose, pressureIndex](){
return Opm::Amg::getQuasiImpesWeights<Matrix, return Opm::Amg::getQuasiImpesWeights<Matrix,
Vector>(*this->matrix_, Vector>(getMatrix(),
pressureIndex, pressureIndex,
transpose); transpose);
}; };
@ -1023,8 +1019,8 @@ protected:
void scaleEquationsAndVariables(Vector& weights) void scaleEquationsAndVariables(Vector& weights)
{ {
// loop over primary variables // loop over primary variables
const auto endi = matrix_->end(); const auto endi = getMatrix().end();
for (auto i = matrix_->begin(); i != endi; ++i) { for (auto i = getMatrix().begin(); i != endi; ++i) {
const auto endj = (*i).end(); const auto endj = (*i).end();
BlockVector& brhs = (*rhs_)[i.index()]; BlockVector& brhs = (*rhs_)[i.index()];
for (auto j = (*i).begin(); j != endj; ++j) { for (auto j = (*i).begin(); j != endj; ++j) {
@ -1040,7 +1036,7 @@ protected:
for (std::size_t ii = 0; ii < brhs.size(); ii++) { for (std::size_t ii = 0; ii < brhs.size(); ii++) {
brhs[ii] *= simulator_.model().eqWeight(i.index(), ii); brhs[ii] *= simulator_.model().eqWeight(i.index(), ii);
} }
if (weights.size() == matrix_->N()) { if (weights.size() == getMatrix().N()) {
BlockVector& bw = weights[i.index()]; BlockVector& bw = weights[i.index()];
for (std::size_t ii = 0; ii < brhs.size(); ii++) { for (std::size_t ii = 0; ii < brhs.size(); ii++) {
bw[ii] /= simulator_.model().eqWeight(i.index(), ii); bw[ii] /= simulator_.model().eqWeight(i.index(), ii);
@ -1065,7 +1061,7 @@ protected:
Vector getQuasiImpesWeights() Vector getQuasiImpesWeights()
{ {
return Amg::getQuasiImpesWeights<Matrix,Vector>(*matrix_, pressureVarIndex, /* transpose=*/ true); return Amg::getQuasiImpesWeights<Matrix,Vector>(getMatrix(), pressureVarIndex, /* transpose=*/ true);
} }
Vector getSimpleWeights(const BlockVector& rhs) Vector getSimpleWeights(const BlockVector& rhs)
@ -1080,8 +1076,8 @@ protected:
void scaleMatrixAndRhs(const Vector& weights) void scaleMatrixAndRhs(const Vector& weights)
{ {
using Block = typename Matrix::block_type; using Block = typename Matrix::block_type;
const auto endi = matrix_->end(); const auto endi = getMatrix().end();
for (auto i = matrix_->begin(); i !=endi; ++i) { for (auto i = getMatrix().begin(); i !=endi; ++i) {
const BlockVector& bweights = weights[i.index()]; const BlockVector& bweights = weights[i.index()];
BlockVector& brhs = (*rhs_)[i.index()]; BlockVector& brhs = (*rhs_)[i.index()];
const auto endj = (*i).end(); const auto endj = (*i).end();
@ -1136,6 +1132,11 @@ protected:
multBlocksVector(b_cp, leftTrans); multBlocksVector(b_cp, leftTrans);
} }
Matrix& getMatrix()
{
return noGhostMat_ ? *noGhostMat_ : *matrix_;
}
const Simulator& simulator_; const Simulator& simulator_;
mutable int iterations_; mutable int iterations_;
mutable bool converged_; mutable bool converged_;
@ -1145,7 +1146,6 @@ protected:
Matrix* matrix_; Matrix* matrix_;
std::unique_ptr<Matrix> noGhostMat_; std::unique_ptr<Matrix> noGhostMat_;
Vector *rhs_; Vector *rhs_;
std::unique_ptr<Matrix> matrix_for_preconditioner_;
std::unique_ptr<FlexibleSolverType> flexibleSolver_; std::unique_ptr<FlexibleSolverType> flexibleSolver_;
std::vector<int> overlapRows_; std::vector<int> overlapRows_;