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>;
85 using TracerMatrix =
typename BaseType::TracerMatrix;
86 using TracerVector =
typename BaseType::TracerVector;
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())
101 , simulator_(simulator)
102 , tbatch({waterPhaseIdx, oilPhaseIdx, gasPhaseIdx})
129 this->
doInit(rst, simulator_.model().numGridDof(),
130 gasPhaseIdx, oilPhaseIdx, waterPhaseIdx);
133 void prepareTracerBatches()
135 for (std::size_t tracerIdx = 0; tracerIdx < this->tracerPhaseIdx_.size(); ++tracerIdx) {
136 if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::waterPhaseIdx) {
137 if (! FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)){
138 throw std::runtime_error(
"Water tracer specified for non-water fluid system:" + this->
name(tracerIdx));
141 wat_.addTracer(tracerIdx, this->tracerConcentration_[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));
151 oil_.addTracer(tracerIdx, this->tracerConcentration_[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));
161 gas_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
166 TracerMatrix* base = this->tracerMatrix_.get();
167 for (
auto& tr : this->tbatch) {
168 if (tr.numTracer() != 0) {
169 if (this->tracerMatrix_)
170 tr.mat = std::move(this->tracerMatrix_);
172 tr.mat = std::make_unique<TracerMatrix>(*base);
182 updateStorageCache();
193 advanceTracerFields();
200 template <
class Restarter>
210 template <
class Restarter>
214 template<
class Serializer>
215 void serializeOp(Serializer& serializer)
217 serializer(
static_cast<BaseType&
>(*
this));
224 template <
class LhsEval>
225 void computeVolume_(LhsEval& freeVolume,
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);
249 void computeFlux_(TracerEvaluation & freeFlux,
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;
283 void assembleTracerEquationVolume(TrRe& tr,
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];
302 computeVolume_(fVolume1, tr.phaseIdx_, elemCtx, 0, 1);
303 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
304 storageOfTimeIndex1[tIdx] = fVolume1*tr.concentrationInitial_[tIdx][I1];
308 TracerEvaluation fVolume;
309 computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, 0);
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;
319 void assembleTracerEquationFlux(TrRe& tr,
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>
342 void assembleTracerEquationWell(TrRe& tr,
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) {
355 wtracer[tIdx] = this->currentConcentration_(eclWell, this->
name(tr.idx_[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);
377 void assembleTracerEquations_()
385 for (
auto& tr : tbatch) {
386 if (tr.numTracer() != 0) {
388 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx)
389 tr.residual_[tIdx] = 0.0;
393 ElementContext elemCtx(simulator_);
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)
402 for (
auto& tr : tbatch) {
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);
423 for (
auto& tr : tbatch) {
424 this->assembleTracerEquationVolume(tr, elemCtx, scvVolume, dt, I, I1);
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);
432 for (
auto& tr : tbatch) {
433 this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J);
439 const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
440 for (
const auto& wellPtr : wellPtrs) {
441 for (
auto& tr : tbatch) {
442 this->assembleTracerEquationWell(tr, *wellPtr);
447 for (
auto& tr : tbatch) {
448 if (tr.numTracer() == 0)
450 auto handle = VectorVectorDataHandle<GridView, std::vector<TracerVector>>(tr.residual_,
451 simulator_.gridView());
452 simulator_.gridView().communicate(handle, Dune::InteriorBorder_All_Interface,
453 Dune::ForwardCommunication);
457 void updateStorageCache()
459 for (
auto& tr : tbatch) {
460 if (tr.numTracer() != 0)
461 tr.concentrationInitial_ = tr.concentration_;
464 ElementContext elemCtx(simulator_);
465 for (
const auto& elem : elements(simulator_.gridView())) {
466 elemCtx.updatePrimaryStencil(elem);
467 elemCtx.updatePrimaryIntensiveQuantities(0);
468 int globalDofIdx = elemCtx.globalSpaceIndex(0, 0);
469 for (
auto& tr : tbatch) {
470 if (tr.numTracer() == 0)
473 computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, 0);
474 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
475 tr.storageOfTimeIndex1_[tIdx][globalDofIdx] = fVolume*tr.concentrationInitial_[tIdx][globalDofIdx];
481 void advanceTracerFields()
483 assembleTracerEquations_();
485 for (
auto& tr : tbatch) {
486 if (tr.numTracer() == 0)
491 std::vector<TracerVector> dx(tr.concentration_);
492 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx)
495 bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_);
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];
503 this->tracerConcentration_[tr.idx_[tIdx]] = tr.concentration_[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) {
544 this->wellTracerRate_.at(std::make_pair(well.name(),this->name(tr.idx_[tIdx]))) *= factor;
552 Simulator& simulator_;
562 template <
typename TV>
564 std::vector<int> idx_;
566 std::vector<TV> concentrationInitial_;
567 std::vector<TV> concentration_;
568 std::vector<TV> storageOfTimeIndex1_;
569 std::vector<TV> residual_;
570 std::unique_ptr<TracerMatrix> mat;
574 return this->concentrationInitial_ == rhs.concentrationInitial_ &&
575 this->concentration_ == rhs.concentration_;
581 result.idx_ = {1,2,3};
582 result.concentrationInitial_ = {5.0, 6.0};
583 result.concentration_ = {7.0, 8.0};
584 result.storageOfTimeIndex1_ = {9.0, 10.0, 11.0};
585 result.residual_ = {12.0, 13.0};
590 template<
class Serializer>
591 void serializeOp(Serializer& serializer)
593 serializer(concentrationInitial_);
594 serializer(concentration_);
597 TracerBatch(
int phaseIdx = 0) : phaseIdx_(phaseIdx) {}
599 int numTracer()
const {
return idx_.size(); }
601 void addTracer(
const int idx,
const TV & concentration)
603 int numGridDof = concentration.size();
604 idx_.emplace_back(idx);
605 concentrationInitial_.emplace_back(concentration);
606 concentration_.emplace_back(concentration);
607 storageOfTimeIndex1_.emplace_back(numGridDof);
608 residual_.emplace_back(numGridDof);
612 std::array<TracerBatch<TracerVector>,3> tbatch;