28#ifndef OPM_TRACER_MODEL_HPP
29#define OPM_TRACER_MODEL_HPP
31#include <opm/common/OpmLog/OpmLog.hpp>
33#include <opm/models/utils/propertysystem.hh>
47template<
class TypeTag,
class MyTypeTag>
49 using type = UndefinedProperty;
61template <
class TypeTag>
63 GetPropType<TypeTag, Properties::GridView>,
64 GetPropType<TypeTag, Properties::DofMapper>,
65 GetPropType<TypeTag, Properties::Stencil>,
66 GetPropType<TypeTag, Properties::Scalar>>
69 GetPropType<TypeTag, Properties::GridView>,
70 GetPropType<TypeTag, Properties::DofMapper>,
71 GetPropType<TypeTag, Properties::Stencil>,
72 GetPropType<TypeTag, Properties::Scalar>>;
73 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
74 using GridView = GetPropType<TypeTag, Properties::GridView>;
75 using Grid = GetPropType<TypeTag, Properties::Grid>;
76 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
77 using Stencil = GetPropType<TypeTag, Properties::Stencil>;
78 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
79 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
80 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
81 using Indices = GetPropType<TypeTag, Properties::Indices>;
83 using TracerEvaluation = DenseAd::Evaluation<Scalar,1>;
88 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
89 enum { numPhases = FluidSystem::numPhases };
90 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
91 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
92 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
96 :
BaseType(simulator.vanguard().gridView(),
97 simulator.vanguard().eclState(),
98 simulator.vanguard().cartesianIndexMapper(),
99 simulator.model().dofMapper(),
100 simulator.vanguard().cellCentroids())
102 ,
tbatch({waterPhaseIdx, oilPhaseIdx, gasPhaseIdx})
130 gasPhaseIdx, oilPhaseIdx, waterPhaseIdx);
135 for (std::size_t tracerIdx = 0; tracerIdx < this->
tracerPhaseIdx_.size(); ++tracerIdx) {
137 if (! FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)){
138 throw std::runtime_error(
"Water tracer specified for non-water fluid system:" + this->
name(tracerIdx));
143 else if (this->
tracerPhaseIdx_[tracerIdx] == FluidSystem::oilPhaseIdx) {
144 if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
145 throw std::runtime_error(
"Oil tracer specified for non-oil fluid system:" + this->
name(tracerIdx));
147 if (FluidSystem::enableVaporizedOil()) {
148 throw std::runtime_error(
"Oil tracer in combination with kw VAPOIL is not supported: " + this->
name(tracerIdx));
153 else if (this->
tracerPhaseIdx_[tracerIdx] == FluidSystem::gasPhaseIdx) {
154 if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
155 throw std::runtime_error(
"Gas tracer specified for non-gas fluid system:" + this->
name(tracerIdx));
157 if (FluidSystem::enableDissolvedGas()) {
158 throw std::runtime_error(
"Gas tracer in combination with kw DISGAS is not supported: " + this->
name(tracerIdx));
167 for (
auto& tr : this->
tbatch) {
168 if (tr.numTracer() != 0) {
172 tr.mat = std::make_unique<TracerMatrix>(*base);
200 template <
class Restarter>
210 template <
class Restarter>
214 template<
class Serializer>
217 serializer(
static_cast<BaseType&
>(*
this));
224 template <
class LhsEval>
226 const int tracerPhaseIdx,
227 const ElementContext& elemCtx,
231 const auto& intQuants = elemCtx.intensiveQuantities(scvIdx, timeIdx);
232 const auto& fs = intQuants.fluidState();
234 decay<Scalar>(fs.saturation(tracerPhaseIdx))
235 *decay<Scalar>(fs.invB(tracerPhaseIdx))
236 *decay<Scalar>(intQuants.porosity());
239 phaseVolume = max(phaseVolume, 1e-10);
241 if (std::is_same<LhsEval, Scalar>::value)
242 freeVolume = phaseVolume;
244 freeVolume = phaseVolume * variable<LhsEval>(1.0, 0);
251 const int tracerPhaseIdx,
252 const ElementContext& elemCtx,
257 const auto& stencil = elemCtx.stencil(timeIdx);
258 const auto& scvf = stencil.interiorFace(scvfIdx);
260 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
261 unsigned inIdx = extQuants.interiorIndex();
263 unsigned upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
265 const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
266 const auto& fs = intQuants.fluidState();
268 Scalar A = scvf.area();
269 Scalar v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
270 Scalar b = decay<Scalar>(fs.invB(tracerPhaseIdx));
272 if (inIdx == upIdx) {
273 freeFlux = A*v*b*variable<TracerEvaluation>(1.0, 0);
277 freeFlux = A*v*b*1.0;
284 const ElementContext& elemCtx,
285 const Scalar scvVolume,
291 if (tr.numTracer() == 0)
294 std::vector<Scalar> storageOfTimeIndex1(tr.numTracer());
295 if (elemCtx.enableStorageCache()) {
296 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
297 storageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I];
303 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
304 storageOfTimeIndex1[tIdx] = fVolume1*tr.concentrationInitial_[tIdx][I1];
308 TracerEvaluation fVolume;
310 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
311 Scalar storageOfTimeIndex0 = fVolume.value()*tr.concentration_[tIdx][I];
312 Scalar localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1[tIdx]) * scvVolume/dt;
313 tr.residual_[tIdx][I][0] += localStorage;
315 (*tr.mat)[I][I][0][0] += fVolume.derivative(0) * scvVolume/dt;
320 const ElementContext& elemCtx,
325 if (tr.numTracer() == 0)
328 TracerEvaluation flux;
330 computeFlux_(flux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0);
331 int globalUpIdx = isUpF ? I : J;
332 for (
int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
333 tr.residual_[tIdx][I][0] += flux.value()*tr.concentration_[tIdx][globalUpIdx];
336 (*tr.mat)[J][I][0][0] = -flux.derivative(0);
337 (*tr.mat)[I][I][0][0] += flux.derivative(0);
341 template<
class TrRe,
class Well>
345 if (tr.numTracer() == 0)
348 const auto& eclWell = well.wellEcl();
349 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
350 this->
wellTracerRate_[std::make_pair(eclWell.name(), this->name(tr.idx_[tIdx]))] = 0.0;
353 std::vector<double> wtracer(tr.numTracer());
354 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
358 for (
auto& perfData : well.perforationData()) {
359 const auto I = perfData.cell_index;
360 Scalar rate = well.volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
362 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
363 tr.residual_[tIdx][I][0] -= rate*wtracer[tIdx];
365 this->
wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) += rate*wtracer[tIdx];
369 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
370 tr.residual_[tIdx][I][0] -= rate*tr.concentration_[tIdx][I];
372 (*tr.mat)[I][I][0][0] -= rate*variable<TracerEvaluation>(1.0, 0).derivative(0);
386 if (tr.numTracer() != 0) {
388 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx)
389 tr.residual_[tIdx] = 0.0;
394 for (
const auto& elem : elements(
simulator_.gridView())) {
395 elemCtx.updateStencil(elem);
397 std::size_t I = elemCtx.globalSpaceIndex( 0, 0);
399 if (elem.partitionType() != Dune::InteriorEntity)
403 if (tr.numTracer() != 0)
404 (*tr.mat)[I][I][0][0] = 1.;
408 elemCtx.updateAllIntensiveQuantities();
409 elemCtx.updateAllExtensiveQuantities();
411 Scalar extrusionFactor =
412 elemCtx.intensiveQuantities( 0, 0).extrusionFactor();
413 Valgrind::CheckDefined(extrusionFactor);
414 assert(isfinite(extrusionFactor));
415 assert(extrusionFactor > 0.0);
417 elemCtx.stencil(0).subControlVolume( 0).volume()
419 Scalar dt = elemCtx.simulator().timeStepSize();
421 std::size_t I1 = elemCtx.globalSpaceIndex( 0, 1);
427 std::size_t numInteriorFaces = elemCtx.numInteriorFaces(0);
428 for (
unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
429 const auto& face = elemCtx.stencil(0).interiorFace(scvfIdx);
430 unsigned j = face.exteriorIndex();
431 unsigned J = elemCtx.globalSpaceIndex( j, 0);
439 const auto& wellPtrs =
simulator_.problem().wellModel().localNonshutWells();
440 for (
const auto& wellPtr : wellPtrs) {
448 if (tr.numTracer() == 0)
452 simulator_.gridView().communicate(handle, Dune::InteriorBorder_All_Interface,
453 Dune::ForwardCommunication);
460 if (tr.numTracer() != 0)
461 tr.concentrationInitial_ = tr.concentration_;
465 for (
const auto& elem : elements(
simulator_.gridView())) {
466 elemCtx.updatePrimaryStencil(elem);
467 elemCtx.updatePrimaryIntensiveQuantities(0);
468 int globalDofIdx = elemCtx.globalSpaceIndex(0, 0);
470 if (tr.numTracer() == 0)
474 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
475 tr.storageOfTimeIndex1_[tIdx][globalDofIdx] = fVolume*tr.concentrationInitial_[tIdx][globalDofIdx];
486 if (tr.numTracer() == 0)
491 std::vector<TracerVector> dx(tr.concentration_);
492 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx)
497 OpmLog::warning(
"### Tracer model: Linear solver did not converge. ###");
500 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
501 tr.concentration_[tIdx] -= dx[tIdx];
507 const auto& wellPtrs =
simulator_.problem().wellModel().localNonshutWells();
508 for (
const auto& wellPtr : wellPtrs) {
509 const auto& well = wellPtr->wellEcl();
511 if (!well.isProducer())
514 Scalar rateWellPos = 0.0;
515 Scalar rateWellNeg = 0.0;
516 for (
auto& perfData : wellPtr->perforationData()) {
517 const int I = perfData.cell_index;
518 Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
521 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
522 this->
wellTracerRate_.at(std::make_pair(well.name(),this->name(tr.idx_[tIdx]))) += rate*tr.concentration_[tIdx][I];
530 Scalar rateWellTotal = rateWellNeg + rateWellPos;
535 std::size_t well_index =
simulator_.problem().wellModel().wellState().index(well.name()).value();
536 Scalar official_well_rate_total =
simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_];
538 rateWellTotal = official_well_rate_total;
540 if (rateWellTotal > rateWellNeg) {
541 const Scalar bucketPrDay = 10.0/(1000.*3600.*24.);
542 const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal/rateWellNeg : 0.0;
543 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
562 template <
typename TV>
570 std::unique_ptr<TracerMatrix>
mat;
581 result.
idx_ = {1,2,3};
590 template<
class Serializer>
603 int numGridDof = concentration.size();
604 idx_.emplace_back(idx);
612 std::array<TracerBatch<TracerVector>,3>
tbatch;
A datahandle sending data located in multiple vectors.
Definition: GenericTracerModel.hpp:53
std::vector< Dune::BlockVector< Dune::FieldVector< GetPropType< TypeTag, Properties::Scalar >, 1 > > > tracerConcentration_
Definition: GenericTracerModel.hpp:118
std::unique_ptr< TracerMatrix > tracerMatrix_
Definition: GenericTracerModel.hpp:119
std::vector< int > tracerPhaseIdx_
Definition: GenericTracerModel.hpp:117
void doInit(bool rst, std::size_t numGridDof, std::size_t gasPhaseIdx, std::size_t oilPhaseIdx, std::size_t waterPhaseIdx)
Initialize all internal data structures needed by the tracer module.
Definition: GenericTracerModel_impl.hpp:162
const std::string & name(int tracerIdx) const
Return the tracer name.
Definition: GenericTracerModel_impl.hpp:155
Dune::BlockVector< Dune::FieldVector< GetPropType< TypeTag, Properties::Scalar >, 1 > > TracerVector
Definition: GenericTracerModel.hpp:56
std::map< std::pair< std::string, std::string >, double > wellTracerRate_
Definition: GenericTracerModel.hpp:123
bool linearSolveBatchwise_(const TracerMatrix &M, std::vector< TracerVector > &x, std::vector< TracerVector > &b)
Definition: GenericTracerModel_impl.hpp:313
int numTracers() const
Return the number of tracers considered by the tracerModel.
Definition: GenericTracerModel_impl.hpp:134
Dune::BCRSMatrix< Opm::MatrixBlock< GetPropType< TypeTag, Properties::Scalar >, 1, 1 > > TracerMatrix
Definition: GenericTracerModel.hpp:55
double currentConcentration_(const Well &eclWell, const std::string &name) const
Definition: GenericTracerModel_impl.hpp:148
A class which handles tracers as specified in by ECL.
Definition: TracerModel.hpp:67
std::array< TracerBatch< TracerVector >, 3 > tbatch
Definition: TracerModel.hpp:612
void updateStorageCache()
Definition: TracerModel.hpp:457
TracerBatch< TracerVector > & oil_
Definition: TracerModel.hpp:614
void advanceTracerFields()
Definition: TracerModel.hpp:481
void init(bool rst)
Definition: TracerModel.hpp:127
TracerModel(Simulator &simulator)
Definition: TracerModel.hpp:95
void computeFlux_(TracerEvaluation &freeFlux, bool &isUpFree, const int tracerPhaseIdx, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: TracerModel.hpp:249
void assembleTracerEquationWell(TrRe &tr, const Well &well)
Definition: TracerModel.hpp:342
void assembleTracerEquationVolume(TrRe &tr, const ElementContext &elemCtx, const Scalar scvVolume, const Scalar dt, unsigned I, unsigned I1)
Definition: TracerModel.hpp:283
void serialize(Restarter &)
This method writes the complete state of all tracer to the hard disk.
Definition: TracerModel.hpp:201
TracerBatch< TracerVector > & gas_
Definition: TracerModel.hpp:615
void assembleTracerEquations_()
Definition: TracerModel.hpp:377
void computeVolume_(LhsEval &freeVolume, const int tracerPhaseIdx, const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: TracerModel.hpp:225
void beginTimeStep()
Definition: TracerModel.hpp:177
void serializeOp(Serializer &serializer)
Definition: TracerModel.hpp:215
void prepareTracerBatches()
Definition: TracerModel.hpp:133
void assembleTracerEquationFlux(TrRe &tr, const ElementContext &elemCtx, unsigned scvfIdx, unsigned I, unsigned J)
Definition: TracerModel.hpp:319
TracerBatch< TracerVector > & wat_
Definition: TracerModel.hpp:613
void endTimeStep()
Informs the tracer model that a time step has just been finished.
Definition: TracerModel.hpp:188
Simulator & simulator_
Definition: TracerModel.hpp:552
void deserialize(Restarter &)
This method restores the complete state of the tracer from disk.
Definition: TracerModel.hpp:211
A data handle sending multiple data store in vectors attached to cells.
Definition: VectorVectorDataHandle.hpp:49
Definition: AluGridVanguard.hpp:57
Definition: BlackoilPhases.hpp:27
Definition: TracerModel.hpp:48
UndefinedProperty type
Definition: TracerModel.hpp:49
Definition: TracerModel.hpp:563
const int phaseIdx_
Definition: TracerModel.hpp:565
void addTracer(const int idx, const TV &concentration)
Definition: TracerModel.hpp:601
bool operator==(const TracerBatch &rhs) const
Definition: TracerModel.hpp:572
static TracerBatch serializationTestObject()
Definition: TracerModel.hpp:578
std::vector< TV > concentrationInitial_
Definition: TracerModel.hpp:566
void serializeOp(Serializer &serializer)
Definition: TracerModel.hpp:591
TracerBatch(int phaseIdx=0)
Definition: TracerModel.hpp:597
int numTracer() const
Definition: TracerModel.hpp:599
std::vector< int > idx_
Definition: TracerModel.hpp:564
std::unique_ptr< TracerMatrix > mat
Definition: TracerModel.hpp:570
std::vector< TV > storageOfTimeIndex1_
Definition: TracerModel.hpp:568
std::vector< TV > residual_
Definition: TracerModel.hpp:569
std::vector< TV > concentration_
Definition: TracerModel.hpp:567