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

[bugfix] Use correct matrix for scaling.

Browse Source 
For some of the scaling approaches the wrong matrix (dereferenced
nullptr) would have been used which should have resulted segmentation
faults. With this commit we add a method getMatrix() that returns the
correct one and use that for scaling.
This commit is contained in:
Markus Blatt 2020-05-13 12:47:30 +02:00
parent 9f278db5a6
commit 664e582b5f
1 changed files with 14 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -502,6 +502,7 @@ 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;
assert(matrix_);
Operator opA(*matrix_, wellModel, interiorCellNum_); Operator opA(*matrix_, wellModel, interiorCellNum_);
solve( opA, x, *rhs_, *comm_ ); solve( opA, x, *rhs_, *comm_ );
@ -615,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);
@ -870,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);
}; };
@ -1018,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) {
@ -1035,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);
@ -1060,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)
@ -1075,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();
@ -1131,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_;