OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request from osae/tracer

Browse Source 
Tracer summary curves, collect well rates.
This commit is contained in:
Bård Skaflestad 2021-06-25 17:48:40 +02:00 committed by GitHub
commit 71de144399
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 343 additions and 90 deletions

45
ebos/eclgenerictracermodel.cc
View File

@ -133,14 +133,14 @@ doInit(bool enabled, size_t numGridDof,
//TBLK keyword
if (!tracer.concentration.empty()){
int tblkDatasize = tracer.concentration.size();
if (!tracer.free_concentration.empty()){
int tblkDatasize = tracer.free_concentration.size();
if (tblkDatasize < cartMapper_.cartesianSize()){
throw std::runtime_error("Wrong size of TBLK for" + tracer.name);
}
for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
tracerConcentration_[tracerIdx][globalDofIdx] = tracer.concentration[cartDofIdx];
tracerConcentration_[tracerIdx][globalDofIdx] = tracer.free_concentration[cartDofIdx];
}
}
//TVDPF keyword
@ -150,7 +150,7 @@ doInit(bool enabled, size_t numGridDof,
for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
const auto& center = eclGrid.getCellCenter(cartDofIdx);
tracerConcentration_[tracerIdx][globalDofIdx] = tracer.tvdpf.evaluate("TRACER_CONCENTRATION", center[2]);
tracerConcentration_[tracerIdx][globalDofIdx] = tracer.free_tvdp.evaluate("TRACER_CONCENTRATION", center[2]);
}
}
++tracerIdx;
@ -243,6 +243,43 @@ linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
return result.converged;
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
bool EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b)
{
#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
Dune::FMatrixPrecision<Scalar>::set_singular_limit(1.e-30);
Dune::FMatrixPrecision<Scalar>::set_absolute_limit(1.e-30);
#endif
Scalar tolerance = 1e-2;
int maxIter = 100;
int verbosity = 0;
using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
TracerOperator tracerOperator(M);
TracerScalarProduct tracerScalarProduct;
TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
TracerSolver solver (tracerOperator, tracerScalarProduct,
tracerPreconditioner, tolerance, maxIter,
verbosity);
bool converged = true;
for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
x[nrhs] = 0.0;
Dune::InverseOperatorResult result;
solver.apply(x[nrhs], b[nrhs], result);
converged = (converged && result.converged);
}
// return the result of the solver
return converged;
}
#if HAVE_DUNE_FEM
template class EclGenericTracerModel<Dune::CpGrid,
Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,

12
ebos/eclgenerictracermodel.hh
View File

@ -68,6 +68,12 @@ public:
*/
Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const;
/*!
* \brief Return well tracer rates
*/
const std::map<std::pair<std::string, std::string>, double>&
getWellTracerRates() const {return wellTracerRate_;}
protected:
EclGenericTracerModel(const GridView& gridView,
const EclipseState& eclState,
@ -85,6 +91,8 @@ protected:
bool linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b);
bool linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b);
const GridView& gridView_;
const EclipseState& eclState_;
const CartesianIndexMapper& cartMapper_;
@ -98,6 +106,10 @@ protected:
TracerVector tracerResidual_;
std::vector<int> cartToGlobal_;
std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> storageOfTimeIndex1_;
// <wellName, tracerIdx> -> wellRate
std::map<std::pair<std::string, std::string>, double> wellTracerRate_;
};
} // namespace Opm

341
ebos/ecltracermodel.hh
View File

@ -77,6 +77,8 @@ class EclTracerModel : public EclGenericTracerModel<GetPropType<TypeTag, Propert
using TracerEvaluation = DenseAd::Evaluation<Scalar,1>;
using TracerVector = typename BaseType::TracerVector;
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
enum { numPhases = FluidSystem::numPhases };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
@ -90,6 +92,9 @@ public:
simulator.vanguard().cartesianIndexMapper(),
simulator.model().dofMapper())
, simulator_(simulator)
, wat_(TracerBatch<TracerVector>(waterPhaseIdx))
, oil_(TracerBatch<TracerVector>(oilPhaseIdx))
, gas_(TracerBatch<TracerVector>(gasPhaseIdx))
{ }
@ -101,6 +106,8 @@ public:
bool enabled = EWOMS_GET_PARAM(TypeTag, bool, EnableTracerModel);
this->doInit(enabled, simulator_.model().numGridDof(),
gasPhaseIdx, oilPhaseIdx, waterPhaseIdx);
prepareTracerBatches();
}
void beginTimeStep()
@ -108,21 +115,9 @@ public:
if (this->numTracers()==0)
return;
this->tracerConcentrationInitial_ = this->tracerConcentration_;
// compute storageCache
ElementContext elemCtx(simulator_);
auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
elemCtx.updateAll(*elemIt);
int globalDofIdx = elemCtx.globalSpaceIndex(0, 0);
for (int tracerIdx = 0; tracerIdx < this->numTracers(); ++ tracerIdx){
Scalar storageOfTimeIndex1;
computeStorage_(storageOfTimeIndex1, elemCtx, 0, /*timIdx=*/0, tracerIdx);
this->storageOfTimeIndex1_[tracerIdx][globalDofIdx] = storageOfTimeIndex1;
}
}
updateStorageCache(wat_);
updateStorageCache(oil_);
updateStorageCache(gas_);
}
/*!
@ -133,19 +128,9 @@ public:
if (this->numTracers()==0)
return;
for (int tracerIdx = 0; tracerIdx < this->numTracers(); ++ tracerIdx){
typename BaseType::TracerVector dx(this->tracerResidual_.size());
// Newton step (currently the system is linear, converge in one iteration)
for (int iter = 0; iter < 5; ++ iter){
linearize_(tracerIdx);
this->linearSolve_(*this->tracerMatrix_, dx, this->tracerResidual_);
this->tracerConcentration_[tracerIdx] -= dx;
if (dx.two_norm()<1e-2)
break;
}
}
advanceTracerFields(wat_);
advanceTracerFields(oil_);
advanceTracerFields(gas_);
}
/*!
@ -167,40 +152,39 @@ public:
{ /* not implemented */ }
protected:
// evaluate storage term for all tracers in a single cell
template <class LhsEval>
void computeStorage_(LhsEval& tracerStorage,
const ElementContext& elemCtx,
unsigned scvIdx,
unsigned timeIdx,
const int tracerIdx)
{
int globalDofIdx = elemCtx.globalSpaceIndex(scvIdx, timeIdx);
// evaluate water storage volume(s) in a single cell
template <class LhsEval>
void computeVolume_(LhsEval& freeVolume,
const int tracerPhaseIdx,
const ElementContext& elemCtx,
unsigned scvIdx,
unsigned timeIdx)
{
const auto& intQuants = elemCtx.intensiveQuantities(scvIdx, timeIdx);
const auto& fs = intQuants.fluidState();
Scalar phaseVolume =
decay<Scalar>(fs.saturation(this->tracerPhaseIdx_[tracerIdx]))
*decay<Scalar>(fs.invB(this->tracerPhaseIdx_[tracerIdx]))
decay<Scalar>(fs.saturation(tracerPhaseIdx))
*decay<Scalar>(fs.invB(tracerPhaseIdx))
*decay<Scalar>(intQuants.porosity());
// avoid singular matrix if no water is present.
phaseVolume = max(phaseVolume, 1e-10);
if (std::is_same<LhsEval, Scalar>::value)
tracerStorage = phaseVolume * this->tracerConcentrationInitial_[tracerIdx][globalDofIdx];
freeVolume = phaseVolume;
else
tracerStorage =
phaseVolume
* variable<LhsEval>(this->tracerConcentration_[tracerIdx][globalDofIdx][0], 0);
freeVolume = phaseVolume * variable<LhsEval>(1.0, 0);
}
// evaluate the tracerflux over one face
void computeFlux_(TracerEvaluation & tracerFlux,
// evaluate the flux(es) over one face
void computeFlux_(TracerEvaluation & freeFlux,
bool & isUpFree,
const int tracerPhaseIdx,
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx,
const int tracerIdx)
unsigned timeIdx)
{
const auto& stencil = elemCtx.stencil(timeIdx);
@ -209,33 +193,41 @@ protected:
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
unsigned inIdx = extQuants.interiorIndex();
const int tracerPhaseIdx = this->tracerPhaseIdx_[tracerIdx];
unsigned upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
int globalUpIdx = elemCtx.globalSpaceIndex(upIdx, timeIdx);
const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
const auto& fs = intQuants.fluidState();
Scalar A = scvf.area();
Scalar v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
Scalar b = decay<Scalar>(fs.invB(this->tracerPhaseIdx_[tracerIdx]));
Scalar c = this->tracerConcentration_[tracerIdx][globalUpIdx];
if (inIdx == upIdx)
tracerFlux = A*v*b*variable<TracerEvaluation>(c, 0);
else
tracerFlux = A*v*b*c;
Scalar b = decay<Scalar>(fs.invB(tracerPhaseIdx));
if (inIdx == upIdx) {
freeFlux = A*v*b*variable<TracerEvaluation>(1.0, 0);
isUpFree = true;
}
else {
freeFlux = A*v*b*1.0;
isUpFree = false;
}
}
void linearize_(int tracerIdx)
template <class TrRe>
void assembleTracerEquations_(TrRe & tr)
{
(*this->tracerMatrix_) = 0.0;
this->tracerResidual_ = 0.0;
if (tr.numTracer() == 0)
return;
// Note that we formulate the equations in terms of a concentration update
// (compared to previous time step) and not absolute concentration.
// This implies that current concentration (tr.concentration_[][]) contributes
// to the rhs both through storrage and flux terms.
// Compare also advanceTracerFields(...) below.
(*this->tracerMatrix_) = 0.0;
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx)
tr.residual_[tIdx] = 0.0;
size_t numGridDof = simulator_.model().numGridDof();
std::vector<double> volumes(numGridDof, 0.0);
ElementContext elemCtx(simulator_);
auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
@ -253,42 +245,68 @@ protected:
Scalar dt = elemCtx.simulator().timeStepSize();
size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/0);
volumes[I] = scvVolume;
TracerEvaluation localStorage;
TracerEvaluation storageOfTimeIndex0;
Scalar storageOfTimeIndex1;
computeStorage_(storageOfTimeIndex0, elemCtx, 0, /*timIdx=*/0, tracerIdx);
if (elemCtx.enableStorageCache())
storageOfTimeIndex1 = this->storageOfTimeIndex1_[tracerIdx][I];
else
computeStorage_(storageOfTimeIndex1, elemCtx, 0, /*timIdx=*/1, tracerIdx);
size_t I1 = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/1);
Scalar storageOfTimeIndex1[tr.numTracer()];
if (elemCtx.enableStorageCache()) {
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
storageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I];
}
}
else {
Scalar fVolume1;
computeVolume_(fVolume1, tr.phaseIdx_, elemCtx, 0, /*timIdx=*/1);
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
storageOfTimeIndex1[tIdx] = fVolume1*tr.concentrationInitial_[tIdx][I1];
}
}
TracerEvaluation fVolume;
computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, /*timIdx=*/0);
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
Scalar storageOfTimeIndex0 = fVolume.value()*tr.concentration_[tIdx][I];
Scalar localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1[tIdx]) * scvVolume/dt;
tr.residual_[tIdx][I][0] += localStorage; //residual + flux
}
(*this->tracerMatrix_)[I][I][0][0] += fVolume.derivative(0) * scvVolume/dt;
localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1) * scvVolume/dt;
this->tracerResidual_[I][0] += localStorage.value(); //residual + flux
(*this->tracerMatrix_)[I][I][0][0] = localStorage.derivative(0);
size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
TracerEvaluation flux;
const auto& face = elemCtx.stencil(0).interiorFace(scvfIdx);
unsigned j = face.exteriorIndex();
unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
computeFlux_(flux, elemCtx, scvfIdx, 0, tracerIdx);
this->tracerResidual_[I][0] += flux.value(); //residual + flux
(*this->tracerMatrix_)[J][I][0][0] = -flux.derivative(0);
(*this->tracerMatrix_)[I][J][0][0] = flux.derivative(0);
bool isUpF;
computeFlux_(flux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0);
int globalUpIdx = isUpF ? I : J;
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
tr.residual_[tIdx][I][0] += flux.value()*tr.concentration_[tIdx][globalUpIdx]; //residual + flux
}
if (isUpF) {
(*this->tracerMatrix_)[J][I][0][0] = -flux.derivative(0);
(*this->tracerMatrix_)[I][I][0][0] += flux.derivative(0);
}
}
}
// Wells
// Wells terms
const int episodeIdx = simulator_.episodeIndex();
const auto& wells = simulator_.vanguard().schedule().getWells(episodeIdx);
for (const auto& well : wells) {
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
this->wellTracerRate_[std::make_pair(well.name(),this->tracerNames_[tr.idx_[tIdx]])] = 0.0;
}
if (well.getStatus() == Well::Status::SHUT)
continue;
const double wtracer = well.getTracerProperties().getConcentration(this->tracerNames_[tracerIdx]);
double wtracer[tr.numTracer()];
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
wtracer[tIdx] = well.getTracerProperties().getConcentration(this->tracerNames_[tr.idx_[tIdx]]);
}
std::array<int, 3> cartesianCoordinate;
for (auto& connection : well.getConnections()) {
@ -300,16 +318,169 @@ protected:
cartesianCoordinate[2] = connection.getK();
const size_t cartIdx = simulator_.vanguard().cartesianIndex(cartesianCoordinate);
const int I = this->cartToGlobal_[cartIdx];
Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, this->tracerPhaseIdx_[tracerIdx]);
if (rate > 0)
this->tracerResidual_[I][0] -= rate*wtracer;
else if (rate < 0)
this->tracerResidual_[I][0] -= rate*this->tracerConcentration_[tracerIdx][I];
Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
if (rate > 0) {
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
tr.residual_[tIdx][I][0] -= rate*wtracer[tIdx];
// Store _injector_ tracer rate for reporting
this->wellTracerRate_.at(std::make_pair(well.name(),this->tracerNames_[tr.idx_[tIdx]])) += rate*wtracer[tIdx];
}
}
else if (rate < 0) {
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
tr.residual_[tIdx][I][0] -= rate*tr.concentration_[tIdx][I];
}
(*this->tracerMatrix_)[I][I][0][0] -= rate*variable<TracerEvaluation>(1.0, 0).derivative(0);
}
}
}
}
template <class TrRe>
void updateStorageCache(TrRe & tr)
{
if (tr.numTracer() == 0)
return;
tr.concentrationInitial_ = tr.concentration_;
ElementContext elemCtx(simulator_);
auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
elemCtx.updateAll(*elemIt);
int globalDofIdx = elemCtx.globalSpaceIndex(0, /*timIdx=*/0);
Scalar fVolume;
computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, /*timIdx=*/0);
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
tr.storageOfTimeIndex1_[tIdx][globalDofIdx] = fVolume*tr.concentrationInitial_[tIdx][globalDofIdx];
}
}
}
template <class TrRe>
void advanceTracerFields(TrRe & tr)
{
if (tr.numTracer() == 0)
return;
// Note that we solve for a concentration update (compared to previous time step)
// Confer also assembleTracerEquations_(...) above.
std::vector<TracerVector> dx(tr.concentration_);
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx)
dx[tIdx] = 0.0;
assembleTracerEquations_(tr);
bool converged = this->linearSolveBatchwise_(*this->tracerMatrix_, dx, tr.residual_);
if (!converged)
std::cout << "### Tracer model: Warning, linear solver did not converge. ###" << std::endl;
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
tr.concentration_[tIdx] -= dx[tIdx];
// Tracer concentrations for restart report
this->tracerConcentration_[tr.idx_[tIdx]] = tr.concentration_[tIdx];
}
// Store _producer_ tracer rate for reporting
const int episodeIdx = simulator_.episodeIndex();
const auto& wells = simulator_.vanguard().schedule().getWells(episodeIdx);
for (const auto& well : wells) {
if (well.getStatus() == Well::Status::SHUT)
continue;
std::array<int, 3> cartesianCoordinate;
for (auto& connection : well.getConnections()) {
if (connection.state() == Connection::State::SHUT)
continue;
cartesianCoordinate[0] = connection.getI();
cartesianCoordinate[1] = connection.getJ();
cartesianCoordinate[2] = connection.getK();
const size_t cartIdx = simulator_.vanguard().cartesianIndex(cartesianCoordinate);
const int I = this->cartToGlobal_[cartIdx];
Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
if (rate < 0 && well.isProducer()) { //Injection rates already reported during assembly
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
this->wellTracerRate_.at(std::make_pair(well.name(),this->tracerNames_[tr.idx_[tIdx]])) += rate*tr.concentration_[tIdx][I];
}
}
}
}
}
void prepareTracerBatches()
{
for (size_t tracerIdx=0; tracerIdx<this->tracerPhaseIdx_.size(); ++tracerIdx) {
if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::waterPhaseIdx) {
if (! FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)){
throw std::runtime_error("Water tracer specified for non-water fluid system:" + this->tracerNames_[tracerIdx]);
}
wat_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
}
else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::oilPhaseIdx) {
if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
throw std::runtime_error("Oil tracer specified for non-oil fluid system:" + this->tracerNames_[tracerIdx]);
}
if (FluidSystem::enableVaporizedOil()) {
throw std::runtime_error("Oil tracer in combination with kw VAPOIL is not supported: " + this->tracerNames_[tracerIdx]);
}
oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
}
else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::gasPhaseIdx) {
if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
throw std::runtime_error("Gas tracer specified for non-gas fluid system:" + this->tracerNames_[tracerIdx]);
}
if (FluidSystem::enableDissolvedGas()) {
throw std::runtime_error("Gas tracer in combination with kw DISGAS is not supported: " + this->tracerNames_[tracerIdx]);
}
gas_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
}
}
}
Simulator& simulator_;
// This struct collects tracers of the same type (i.e, transported in same phase).
// The idea beeing that, under the assumption of linearity, tracers of same type can
// be solved in concert, having a common system matrix but separate right-hand-sides.
// Since oil or gas tracers appears in dual compositions when VAPOIL respectively DISGAS
// is active, the template argument is intended to support future extension to these
// scenarios by supplying an extended vector type.
template <typename TV>
struct TracerBatch {
std::vector<int> idx_;
const int phaseIdx_;
std::vector<TV> concentrationInitial_;
std::vector<TV> concentration_;
std::vector<TV> storageOfTimeIndex1_;
std::vector<TV> residual_;
TracerBatch(int phaseIdx) : phaseIdx_(phaseIdx) {}
int numTracer() const {return idx_.size(); }
void addTracer(const int idx, const TV & concentration)
{
int numGridDof = concentration.size();
idx_.emplace_back(idx);
concentrationInitial_.emplace_back(concentration);
concentration_.emplace_back(concentration);
storageOfTimeIndex1_.emplace_back(numGridDof);
residual_.emplace_back(numGridDof);
}
};
TracerBatch<TracerVector> wat_;
TracerBatch<TracerVector> oil_;
TracerBatch<TracerVector> gas_;
};
} // namespace Opm

6
opm/simulators/wells/BlackoilWellModel.hpp
View File

@ -229,6 +229,8 @@ namespace Opm {
return this->wasDynamicallyShutThisTimeStep(well_ndex);
});
this->assignWellTracerRates(wsrpt);
this->assignWellGuideRates(wsrpt);
this->assignShutConnections(wsrpt, this->reportStepIndex());
@ -395,7 +397,9 @@ namespace Opm {
const int pvtreg,
std::vector<double>& resv_coeff) override;
void computeWellTemperature();
void computeWellTemperature();
void assignWellTracerRates(data::Wells& wsrpt) const;
private:
BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& pu);

29
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -1736,4 +1736,33 @@ namespace Opm {
this->wellState().update_temperature(wellID, weighted_temperature/total_weight);
}
}
template <typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assignWellTracerRates(data::Wells& wsrpt) const
{
const auto & wellTracerRates = ebosSimulator_.problem().tracerModel().getWellTracerRates();
if (wellTracerRates.empty())
return; // no tracers
for (const auto& wTR : wellTracerRates) {
std::string wellName = wTR.first.first;
auto xwPos = wsrpt.find(wellName);
if (xwPos == wsrpt.end()) { // No well results.
continue;
}
std::string tracerName = wTR.first.second;
double rate = wTR.second;
xwPos->second.rates.set(data::Rates::opt::tracer, rate, tracerName);
}
}
} // namespace Opm