added: MultisegmentWellEquations::apply(x,Ax)

this applies the equation system to a vector.
use the new method in the well implementation.
This commit is contained in:
Arne Morten Kvarving
2022-11-11 21:41:24 +01:00
parent 8fe6b3968e
commit f2acbccc1a
3 changed files with 24 additions and 13 deletions

View File

@@ -22,6 +22,7 @@
#include <config.h> #include <config.h>
#include <opm/simulators/wells/MultisegmentWellEquations.hpp> #include <opm/simulators/wells/MultisegmentWellEquations.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/MultisegmentWellGeneric.hpp> #include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
namespace Opm { namespace Opm {
@@ -114,6 +115,21 @@ void MultisegmentWellEquations<Scalar,numWellEq,numEq>::clear()
duneDSolver_.reset(); duneDSolver_.reset();
} }
template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
apply(const BVector& x, BVector& Ax) const
{
BVectorWell Bx(duneB_.N());
duneB_.mv(x, Bx);
// invDBx = duneD^-1 * Bx_
const BVectorWell invDBx = mswellhelpers::applyUMFPack(duneD_, duneDSolver_, Bx);
// Ax = Ax - duneC_^T * invDBx
duneC_.mmtv(invDBx,Ax);
}
#define INSTANCE(numWellEq, numEq) \ #define INSTANCE(numWellEq, numEq) \
template class MultisegmentWellEquations<double,numWellEq,numEq>; template class MultisegmentWellEquations<double,numWellEq,numEq>;

View File

@@ -74,7 +74,9 @@ public:
//! \brief Set all coefficients to 0. //! \brief Set all coefficients to 0.
void clear(); void clear();
// TODO, the following should go to a class for computing purpose //! \brief Apply linear operator to vector.
void apply(const BVector& x, BVector& Ax) const;
// two off-diagonal matrices // two off-diagonal matrices
OffDiagMatWell duneB_; OffDiagMatWell duneB_;
OffDiagMatWell duneC_; OffDiagMatWell duneC_;

View File

@@ -196,23 +196,16 @@ namespace Opm
MultisegmentWell<TypeTag>:: MultisegmentWell<TypeTag>::
apply(const BVector& x, BVector& Ax) const apply(const BVector& x, BVector& Ax) const
{ {
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return; if (!this->isOperableAndSolvable() && !this->wellIsStopped()) {
return;
}
if ( this->param_.matrix_add_well_contributions_ ) if (this->param_.matrix_add_well_contributions_) {
{
// Contributions are already in the matrix itself // Contributions are already in the matrix itself
return; return;
} }
BVectorWell Bx(this->linSys_.duneB_.N());
this->linSys_.duneB_.mv(x, Bx); this->linSys_.apply(x, Ax);
// invDBx = duneD^-1 * Bx_
const BVectorWell invDBx = mswellhelpers::applyUMFPack(this->linSys_.duneD_,
this->linSys_.duneDSolver_, Bx);
// Ax = Ax - duneC_^T * invDBx
this->linSys_.duneC_.mmtv(invDBx,Ax);
} }