diff --git a/opm/simulators/wells/StandardWellEquations.cpp b/opm/simulators/wells/StandardWellEquations.cpp index 67e0b1464..91d9876e5 100644 --- a/opm/simulators/wells/StandardWellEquations.cpp +++ b/opm/simulators/wells/StandardWellEquations.cpp @@ -130,6 +130,17 @@ void StandardWellEquations::apply(const BVector& x, BVector& Ax) c duneC_.mmtv(invDBx, Ax); } +template +void StandardWellEquations::apply(BVector& r) const +{ + assert(invDrw_.size() == invDuneD_.N()); + + // invDrw_ = invDuneD_ * resWell_ + invDuneD_.mv(resWell_, invDrw_); + // r = r - duneC_^T * invDrw_ + duneC_.mmtv(invDrw_, r); +} + #define INSTANCE(N) \ template class StandardWellEquations; diff --git a/opm/simulators/wells/StandardWellEquations.hpp b/opm/simulators/wells/StandardWellEquations.hpp index 4bf09e8e2..f63baca3d 100644 --- a/opm/simulators/wells/StandardWellEquations.hpp +++ b/opm/simulators/wells/StandardWellEquations.hpp @@ -76,6 +76,9 @@ public: //! \brief Apply linear operator to vector. void apply(const BVector& x, BVector& Ax) const; + //! \brief Apply linear operator to vector. + void apply(BVector& r) const; + // two off-diagonal matrices OffDiagMatWell duneB_; OffDiagMatWell duneC_; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 47d18e284..16a33b1e9 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -1736,12 +1736,7 @@ namespace Opm { if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return; - assert(this->linSys_.invDrw_.size() == this->linSys_.invDuneD_.N()); - - // invDrw_ = invDuneD_ * resWell_ - this->linSys_.invDuneD_.mv(this->linSys_.resWell_, this->linSys_.invDrw_); - // r = r - duneC_^T * invDrw_ - this->linSys_.duneC_.mmtv(this->linSys_.invDrw_, r); + this->linSys_.apply(r); } template