From 65efe8e1fcc0ece9987bfd9404359491d164ff36 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 11 Nov 2022 21:41:24 +0100 Subject: [PATCH] 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. --- opm/simulators/wells/BlackoilWellModel.hpp | 2 +- .../wells/BlackoilWellModel_impl.hpp | 4 +- .../wells/StandardWellEquations.cpp | 55 +++++++++++++++++++ .../wells/StandardWellEquations.hpp | 5 ++ opm/simulators/wells/StandardWellEval.cpp | 51 ----------------- opm/simulators/wells/StandardWellEval.hpp | 7 ++- opm/simulators/wells/WellInterface.hpp | 1 - 7 files changed, 67 insertions(+), 58 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index b84ce54a7..52920d4b3 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -320,7 +320,7 @@ namespace Opm { Simulator& ebosSimulator_; // a vector of all the wells. - std::vector well_container_{}; + std::vector well_container_{}; std::vector is_cell_perforated_{}; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 2de06057a..95d23a832 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1213,9 +1213,9 @@ namespace Opm { for(unsigned int i = 0; i < well_container_.size(); i++){ auto& well = well_container_[i]; // maybe WellInterface could implement addWellContribution() - auto derived_std = std::dynamic_pointer_cast >(well); + auto derived_std = std::dynamic_pointer_cast>(well); if (derived_std) { - derived_std->addWellContribution(wellContribs); + derived_std->linSys().extract(derived_std->numStaticWellEq, wellContribs); } else { auto derived_ms = std::dynamic_pointer_cast >(well); if (derived_ms) { diff --git a/opm/simulators/wells/StandardWellEquations.cpp b/opm/simulators/wells/StandardWellEquations.cpp index 197686fd3..63360af11 100644 --- a/opm/simulators/wells/StandardWellEquations.cpp +++ b/opm/simulators/wells/StandardWellEquations.cpp @@ -24,6 +24,7 @@ #include +#include #include namespace Opm @@ -178,6 +179,60 @@ recoverSolutionWell(const BVector& x, BVectorWell& xw) const invDuneD_.mv(resWell, xw); } +template +void StandardWellEquations:: +extract(const int numStaticWellEq, + WellContributions& wellContribs) const +{ + std::vector colIndices; + std::vector 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; diff --git a/opm/simulators/wells/StandardWellEquations.hpp b/opm/simulators/wells/StandardWellEquations.hpp index e182362dd..ce5625e08 100644 --- a/opm/simulators/wells/StandardWellEquations.hpp +++ b/opm/simulators/wells/StandardWellEquations.hpp @@ -34,6 +34,7 @@ namespace Opm { class ParallelWellInfo; +class WellContributions; template 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_; diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp index 374d08088..df398afa5 100644 --- a/opm/simulators/wells/StandardWellEval.cpp +++ b/opm/simulators/wells/StandardWellEval.cpp @@ -1044,57 +1044,6 @@ init(std::vector& perf_depth, baseif_.numPerfs(), baseif_.cells()); } -template -void -StandardWellEval:: -addWellContribution(WellContributions& wellContribs) const -{ - std::vector colIndices; - std::vector 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 unsigned int StandardWellEval:: getNumBlocks() const diff --git a/opm/simulators/wells/StandardWellEval.hpp b/opm/simulators/wells/StandardWellEval.hpp index 50a853263..3a18f8a37 100644 --- a/opm/simulators/wells/StandardWellEval.hpp +++ b/opm/simulators/wells/StandardWellEval.hpp @@ -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; using BVectorWell = typename StandardWellGeneric::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& linSys() const { return linSys_; } diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 1a2389c26..eb8dc6d5e 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -71,7 +71,6 @@ class WellInterface : public WellInterfaceIndices> { public: - using ModelParameters = BlackoilModelParametersEbos; using Grid = GetPropType;