added: MultisegmentWellEquations::extract(WellContributions&)

this adds the well matrices to a WellContributions object.
use the new method in the implementation.
This commit is contained in:
Arne Morten Kvarving 2022-11-11 21:41:24 +01:00
parent d50aaf8ed4
commit 4a2fcd5f09
5 changed files with 80 additions and 74 deletions

View File

@ -1219,7 +1219,7 @@ namespace Opm {
} else {
auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
if (derived_ms) {
derived_ms->addWellContribution(wellContribs);
derived_ms->linSys().extract(wellContribs);
} else {
OpmLog::warning("Warning unknown type of well");
}

View File

@ -26,6 +26,8 @@
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/linalg/bda/WellContributions.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
@ -179,6 +181,73 @@ recoverSolutionWell(const BVector& x, BVectorWell& xw) const
xw = mswellhelpers::applyUMFPack(*duneDSolver_, resWell);
}
template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
extract(WellContributions& wellContribs) const
{
unsigned int Mb = duneB_.N(); // number of blockrows in duneB_, duneC_ and duneD_
unsigned int BnumBlocks = duneB_.nonzeroes();
unsigned int DnumBlocks = duneD_.nonzeroes();
// duneC
std::vector<unsigned int> Ccols;
std::vector<double> Cvals;
Ccols.reserve(BnumBlocks);
Cvals.reserve(BnumBlocks * numEq * numWellEq);
for (auto rowC = duneC_.begin(); rowC != duneC_.end(); ++rowC) {
for (auto colC = rowC->begin(), endC = rowC->end(); colC != endC; ++colC) {
Ccols.emplace_back(colC.index());
for (int i = 0; i < numWellEq; ++i) {
for (int j = 0; j < numEq; ++j) {
Cvals.emplace_back((*colC)[i][j]);
}
}
}
}
// duneD
Dune::UMFPack<DiagMatWell> umfpackMatrix(duneD_, 0);
double* Dvals = umfpackMatrix.getInternalMatrix().getValues();
auto* Dcols = umfpackMatrix.getInternalMatrix().getColStart();
auto* Drows = umfpackMatrix.getInternalMatrix().getRowIndex();
// duneB
std::vector<unsigned int> Bcols;
std::vector<unsigned int> Brows;
std::vector<double> Bvals;
Bcols.reserve(BnumBlocks);
Brows.reserve(Mb+1);
Bvals.reserve(BnumBlocks * numEq * numWellEq);
Brows.emplace_back(0);
unsigned int sumBlocks = 0;
for (auto rowB = duneB_.begin(); rowB != duneB_.end(); ++rowB) {
int sizeRow = 0;
for (auto colB = rowB->begin(), endB = rowB->end(); colB != endB; ++colB) {
Bcols.emplace_back(colB.index());
for (int i = 0; i < numWellEq; ++i) {
for (int j = 0; j < numEq; ++j) {
Bvals.emplace_back((*colB)[i][j]);
}
}
sizeRow++;
}
sumBlocks += sizeRow;
Brows.emplace_back(sumBlocks);
}
wellContribs.addMultisegmentWellContribution(numEq,
numWellEq,
Mb,
Bvals,
Bcols,
Brows,
DnumBlocks,
Dvals,
Dcols,
Drows,
Cvals);
}
#define INSTANCE(numWellEq, numEq) \
template class MultisegmentWellEquations<double,numWellEq,numEq>;

View File

@ -37,6 +37,7 @@ namespace Opm
{
template<class Scalar> class MultisegmentWellGeneric;
class WellContributions;
template<class Scalar, int numWellEq, int numEq>
class MultisegmentWellEquations
@ -90,6 +91,9 @@ 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(WellContributions& wellContribs) const;
// two off-diagonal matrices
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;

View File

@ -30,7 +30,6 @@
#include <opm/models/blackoil/blackoilonephaseindices.hh>
#include <opm/models/blackoil/blackoiltwophaseindices.hh>
#include <opm/simulators/linalg/bda/WellContributions.hpp>
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
@ -1785,74 +1784,6 @@ updateUpwindingSegments()
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
addWellContribution(WellContributions& wellContribs) const
{
unsigned int Mb = linSys_.duneB_.N(); // number of blockrows in duneB_, duneC_ and duneD_
unsigned int BnumBlocks = linSys_.duneB_.nonzeroes();
unsigned int DnumBlocks = linSys_.duneD_.nonzeroes();
// duneC
std::vector<unsigned int> Ccols;
std::vector<double> Cvals;
Ccols.reserve(BnumBlocks);
Cvals.reserve(BnumBlocks * Indices::numEq * numWellEq);
for (auto rowC = linSys_.duneC_.begin(); rowC != linSys_.duneC_.end(); ++rowC) {
for (auto colC = rowC->begin(), endC = rowC->end(); colC != endC; ++colC) {
Ccols.emplace_back(colC.index());
for (int i = 0; i < numWellEq; ++i) {
for (int j = 0; j < Indices::numEq; ++j) {
Cvals.emplace_back((*colC)[i][j]);
}
}
}
}
// duneD
Dune::UMFPack<typename Equations::DiagMatWell> umfpackMatrix(linSys_.duneD_, 0);
double *Dvals = umfpackMatrix.getInternalMatrix().getValues();
auto *Dcols = umfpackMatrix.getInternalMatrix().getColStart();
auto *Drows = umfpackMatrix.getInternalMatrix().getRowIndex();
// duneB
std::vector<unsigned int> Bcols;
std::vector<unsigned int> Brows;
std::vector<double> Bvals;
Bcols.reserve(BnumBlocks);
Brows.reserve(Mb+1);
Bvals.reserve(BnumBlocks * Indices::numEq * numWellEq);
Brows.emplace_back(0);
unsigned int sumBlocks = 0;
for (auto rowB = linSys_.duneB_.begin(); rowB != linSys_.duneB_.end(); ++rowB) {
int sizeRow = 0;
for (auto colB = rowB->begin(), endB = rowB->end(); colB != endB; ++colB) {
Bcols.emplace_back(colB.index());
for (int i = 0; i < numWellEq; ++i) {
for (int j = 0; j < Indices::numEq; ++j) {
Bvals.emplace_back((*colB)[i][j]);
}
}
sizeRow++;
}
sumBlocks += sizeRow;
Brows.emplace_back(sumBlocks);
}
wellContribs.addMultisegmentWellContribution(Indices::numEq,
numWellEq,
Mb,
Bvals,
Bcols,
Brows,
DnumBlocks,
Dvals,
Dcols,
Drows,
Cvals);
}
#define INSTANCE(...) \
template class MultisegmentWellEval<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;

View File

@ -51,10 +51,6 @@ class WellState;
template<typename FluidSystem, typename Indices, typename Scalar>
class MultisegmentWellEval : public MultisegmentWellGeneric<Scalar>
{
public:
/// add the contribution (C, D, B matrices) of this Well to the WellContributions object
void addWellContribution(WellContributions& wellContribs) const;
protected:
// TODO: for now, not considering the polymer, solvent and so on to simplify the development process.
@ -98,6 +94,12 @@ protected:
using EvalWell = DenseAd::Evaluation<double, /*size=*/Indices::numEq + numWellEq>;
using Eval = DenseAd::Evaluation<Scalar, /*size=*/Indices::numEq>;
public:
//! \brief Returns a const reference to equation system.
const Equations& linSys() const
{ return linSys_; }
protected:
MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif);
void initMatrixAndVectors(const int num_cells);