added: StandardWellEquations::extract(WellContributions&)

this adds the well matrices to a WellContributions object.
this is the core of StandardWellEval::addWellContributions.
use the new method in the implementation.
This commit is contained in:
Arne Morten Kvarving 2022-11-11 21:41:24 +01:00
parent 40640e9da6
commit 65efe8e1fc
7 changed files with 67 additions and 58 deletions

View File

@ -1215,7 +1215,7 @@ namespace Opm {
// maybe WellInterface could implement addWellContribution()
auto derived_std = std::dynamic_pointer_cast<StandardWell<TypeTag>>(well);
if (derived_std) {
derived_std->addWellContribution(wellContribs);
derived_std->linSys().extract(derived_std->numStaticWellEq, wellContribs);
} else {
auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
if (derived_ms) {

View File

@ -24,6 +24,7 @@
#include <opm/common/Exceptions.hpp>
#include <opm/simulators/linalg/bda/WellContributions.hpp>
#include <opm/simulators/linalg/matrixblock.hh>
namespace Opm
@ -178,6 +179,60 @@ recoverSolutionWell(const BVector& x, BVectorWell& xw) const
invDuneD_.mv(resWell, xw);
}
template<class Scalar, int numEq>
void StandardWellEquations<Scalar,numEq>::
extract(const int numStaticWellEq,
WellContributions& wellContribs) const
{
std::vector<int> colIndices;
std::vector<double> nnzValues;
colIndices.reserve(duneB_.nonzeroes());
nnzValues.reserve(duneB_.nonzeroes() * numStaticWellEq * numEq);
// duneC
for (auto colC = duneC_[0].begin(),
endC = duneC_[0].end(); colC != endC; ++colC )
{
colIndices.emplace_back(colC.index());
for (int i = 0; i < numStaticWellEq; ++i) {
for (int j = 0; j < numEq; ++j) {
nnzValues.emplace_back((*colC)[i][j]);
}
}
}
wellContribs.addMatrix(WellContributions::MatrixType::C,
colIndices.data(), nnzValues.data(), duneC_.nonzeroes());
// invDuneD
colIndices.clear();
nnzValues.clear();
colIndices.emplace_back(0);
for (int i = 0; i < numStaticWellEq; ++i)
{
for (int j = 0; j < numStaticWellEq; ++j) {
nnzValues.emplace_back(invDuneD_[0][0][i][j]);
}
}
wellContribs.addMatrix(WellContributions::MatrixType::D,
colIndices.data(), nnzValues.data(), 1);
// duneB
colIndices.clear();
nnzValues.clear();
for (auto colB = duneB_[0].begin(),
endB = duneB_[0].end(); colB != endB; ++colB )
{
colIndices.emplace_back(colB.index());
for (int i = 0; i < numStaticWellEq; ++i) {
for (int j = 0; j < numEq; ++j) {
nnzValues.emplace_back((*colB)[i][j]);
}
}
}
wellContribs.addMatrix(WellContributions::MatrixType::B,
colIndices.data(), nnzValues.data(), duneB_.nonzeroes());
}
#define INSTANCE(N) \
template class StandardWellEquations<double,N>;

View File

@ -34,6 +34,7 @@ namespace Opm
{
class ParallelWellInfo;
class WellContributions;
template<class Scalar, int numEq>
class StandardWellEquations
@ -89,6 +90,10 @@ public:
//! \details xw = inv(D)*(rw - C*x)
void recoverSolutionWell(const BVector& x, BVectorWell& xw) const;
//! \brief Add the matrices of this well to the WellContributions object.
void extract(const int numStaticWellEq,
WellContributions& wellContribs) const;
// two off-diagonal matrices
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;

View File

@ -1044,57 +1044,6 @@ init(std::vector<double>& perf_depth,
baseif_.numPerfs(), baseif_.cells());
}
template<class FluidSystem, class Indices, class Scalar>
void
StandardWellEval<FluidSystem,Indices,Scalar>::
addWellContribution(WellContributions& wellContribs) const
{
std::vector<int> colIndices;
std::vector<double> nnzValues;
colIndices.reserve(this->linSys_.duneB_.nonzeroes());
nnzValues.reserve(this->linSys_.duneB_.nonzeroes()*numStaticWellEq * Indices::numEq);
// duneC
for (auto colC = this->linSys_.duneC_[0].begin(),
endC = this->linSys_.duneC_[0].end(); colC != endC; ++colC )
{
colIndices.emplace_back(colC.index());
for (int i = 0; i < numStaticWellEq; ++i) {
for (int j = 0; j < Indices::numEq; ++j) {
nnzValues.emplace_back((*colC)[i][j]);
}
}
}
wellContribs.addMatrix(WellContributions::MatrixType::C, colIndices.data(), nnzValues.data(), this->linSys_.duneC_.nonzeroes());
// invDuneD
colIndices.clear();
nnzValues.clear();
colIndices.emplace_back(0);
for (int i = 0; i < numStaticWellEq; ++i)
{
for (int j = 0; j < numStaticWellEq; ++j) {
nnzValues.emplace_back(this->linSys_.invDuneD_[0][0][i][j]);
}
}
wellContribs.addMatrix(WellContributions::MatrixType::D, colIndices.data(), nnzValues.data(), 1);
// duneB
colIndices.clear();
nnzValues.clear();
for (auto colB = this->linSys_.duneB_[0].begin(),
endB = this->linSys_.duneB_[0].end(); colB != endB; ++colB )
{
colIndices.emplace_back(colB.index());
for (int i = 0; i < numStaticWellEq; ++i) {
for (int j = 0; j < Indices::numEq; ++j) {
nnzValues.emplace_back((*colB)[i][j]);
}
}
}
wellContribs.addMatrix(WellContributions::MatrixType::B, colIndices.data(), nnzValues.data(), this->linSys_.duneB_.nonzeroes());
}
template<class FluidSystem, class Indices, class Scalar>
unsigned int StandardWellEval<FluidSystem,Indices,Scalar>::
getNumBlocks() const

View File

@ -53,7 +53,11 @@ protected:
static constexpr int numWellControlEq = 1;
// number of the well equations that will always be used
// based on the solution strategy, there might be other well equations be introduced
public:
static constexpr int numStaticWellEq = numWellConservationEq + numWellControlEq;
protected:
// the index for Bhp in primary variables and also the index of well control equation
// they both will be the last one in their respective system.
// TODO: we should have indices for the well equations and well primary variables separately
@ -97,9 +101,6 @@ public:
using Eval = DenseAd::Evaluation<Scalar, Indices::numEq>;
using BVectorWell = typename StandardWellGeneric<Scalar>::BVectorWell;
/// add the contribution (C, D^-1, B matrices) of this Well to the WellContributions object
void addWellContribution(WellContributions& wellContribs) const;
//! \brief Returns a const reference to equation system.
const StandardWellEquations<Scalar,Indices::numEq>& linSys() const
{ return linSys_; }

View File

@ -71,7 +71,6 @@ class WellInterface : public WellInterfaceIndices<GetPropType<TypeTag, Propertie
GetPropType<TypeTag, Properties::Scalar>>
{
public:
using ModelParameters = BlackoilModelParametersEbos<TypeTag>;
using Grid = GetPropType<TypeTag, Properties::Grid>;