63 GetPropType<TypeTag, Properties::GridView>,
64 GetPropType<TypeTag, Properties::DofMapper>,
65 GetPropType<TypeTag, Properties::Stencil>,
66 GetPropType<TypeTag, Properties::FluidSystem>,
67 GetPropType<TypeTag, Properties::Scalar>>
85 using TracerEvaluation = DenseAd::Evaluation<Scalar,1>;
87 using TracerMatrix =
typename BaseType::TracerMatrix;
88 using TracerVector =
typename BaseType::TracerVector;
91 enum { numPhases = FluidSystem::numPhases };
92 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
93 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
94 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
98 :
BaseType(simulator.vanguard().gridView(),
99 simulator.vanguard().eclState(),
100 simulator.vanguard().cartesianIndexMapper(),
101 simulator.model().dofMapper(),
102 simulator.vanguard().cellCentroids())
103 , simulator_(simulator)
104 , tbatch({waterPhaseIdx, oilPhaseIdx, gasPhaseIdx})
131 this->
doInit(rst, simulator_.model().numGridDof(),
132 gasPhaseIdx, oilPhaseIdx, waterPhaseIdx);
135 void prepareTracerBatches()
137 for (std::size_t tracerIdx = 0; tracerIdx < this->tracerPhaseIdx_.size(); ++tracerIdx) {
138 if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::waterPhaseIdx) {
139 if (! FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)){
140 throw std::runtime_error(
"Water tracer specified for non-water fluid system:" + this->
name(tracerIdx));
143 wat_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
145 else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::oilPhaseIdx) {
146 if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
147 throw std::runtime_error(
"Oil tracer specified for non-oil fluid system:" + this->
name(tracerIdx));
150 oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
152 else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::gasPhaseIdx) {
153 if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
154 throw std::runtime_error(
"Gas tracer specified for non-gas fluid system:" + this->
name(tracerIdx));
157 gas_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
161 fVol1_[this->tracerPhaseIdx_[tracerIdx]].resize(this->freeTracerConcentration_[tracerIdx].size());
162 sVol1_[this->tracerPhaseIdx_[tracerIdx]].resize(this->solTracerConcentration_[tracerIdx].size());
163 dsVol_[this->tracerPhaseIdx_[tracerIdx]].resize(this->solTracerConcentration_[tracerIdx].size());
164 dfVol_[this->tracerPhaseIdx_[tracerIdx]].resize(this->solTracerConcentration_[tracerIdx].size());
168 TracerMatrix* base = this->tracerMatrix_.get();
169 for (
auto& tr : this->tbatch) {
170 if (tr.numTracer() != 0) {
171 if (this->tracerMatrix_)
172 tr.mat = std::move(this->tracerMatrix_);
174 tr.mat = std::make_unique<TracerMatrix>(*base);
184 updateStorageCache();
195 advanceTracerFields();
202 template <
class Restarter>
212 template <
class Restarter>
216 template<
class Serializer>
217 void serializeOp(Serializer& serializer)
219 serializer(
static_cast<BaseType&
>(*
this));
226 Scalar computeFreeVolume_(
const int tracerPhaseIdx,
227 unsigned globalDofIdx,
230 const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx);
231 const auto& fs = intQuants.fluidState();
234 decay<Scalar>(fs.saturation(tracerPhaseIdx))
235 *decay<Scalar>(fs.invB(tracerPhaseIdx))
236 *decay<Scalar>(intQuants.porosity());
238 return max(phaseVolume, 1e-10);
242 Scalar computeSolutionVolume_(
const int tracerPhaseIdx,
243 unsigned globalDofIdx,
246 const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx);
247 const auto& fs = intQuants.fluidState();
252 if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
254 decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx))
255 * decay<Scalar>(fs.invB(FluidSystem::gasPhaseIdx))
256 * decay<Scalar>(fs.Rv())
257 * decay<Scalar>(intQuants.porosity());
261 else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
263 decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx))
264 * decay<Scalar>(fs.invB(FluidSystem::oilPhaseIdx))
265 * decay<Scalar>(fs.Rs())
266 * decay<Scalar>(intQuants.porosity());
272 return max(phaseVolume, 1e-10);
276 void computeFreeFlux_(TracerEvaluation & freeFlux,
278 const int tracerPhaseIdx,
279 const ElementContext& elemCtx,
283 const auto& stencil = elemCtx.stencil(timeIdx);
284 const auto& scvf = stencil.interiorFace(scvfIdx);
286 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
287 unsigned inIdx = extQuants.interiorIndex();
289 unsigned upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
291 const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
292 const auto& fs = intQuants.fluidState();
295 decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx))
296 * decay<Scalar>(fs.invB(tracerPhaseIdx));
298 Scalar A = scvf.area();
299 if (inIdx == upIdx) {
300 freeFlux = A*v*variable<TracerEvaluation>(1.0, 0);
309 void computeSolFlux_(TracerEvaluation & solFlux,
311 const int tracerPhaseIdx,
312 const ElementContext& elemCtx,
316 const auto& stencil = elemCtx.stencil(timeIdx);
317 const auto& scvf = stencil.interiorFace(scvfIdx);
319 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
320 unsigned inIdx = extQuants.interiorIndex();
326 if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
327 upIdx = extQuants.upstreamIndex(FluidSystem::gasPhaseIdx);
329 const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
330 const auto& fs = intQuants.fluidState();
332 decay<Scalar>(fs.invB(FluidSystem::gasPhaseIdx))
333 * decay<Scalar>(extQuants.volumeFlux(FluidSystem::gasPhaseIdx))
334 * decay<Scalar>(fs.Rv());
337 else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
338 upIdx = extQuants.upstreamIndex(FluidSystem::oilPhaseIdx);
340 const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
341 const auto& fs = intQuants.fluidState();
343 decay<Scalar>(fs.invB(FluidSystem::oilPhaseIdx))
344 * decay<Scalar>(extQuants.volumeFlux(FluidSystem::oilPhaseIdx))
345 * decay<Scalar>(fs.Rs());
352 Scalar A = scvf.area();
353 if (inIdx == upIdx) {
354 solFlux = A*v*variable<TracerEvaluation>(1.0, 0);
364 void assembleTracerEquationVolume(TrRe& tr,
365 const ElementContext& elemCtx,
366 const Scalar scvVolume,
372 if (tr.numTracer() == 0)
376 std::vector<Scalar> storageOfTimeIndex1(tr.numTracer());
377 std::vector<Scalar> fStorageOfTimeIndex1(tr.numTracer());
378 std::vector<Scalar> sStorageOfTimeIndex1(tr.numTracer());
379 if (elemCtx.enableStorageCache()) {
380 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
381 fStorageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I][0];
382 sStorageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I][1];
386 Scalar fVolume1 = computeFreeVolume_(tr.phaseIdx_, I1, 1);
387 Scalar sVolume1 = computeSolutionVolume_(tr.phaseIdx_, I1, 1);
388 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
389 fStorageOfTimeIndex1[tIdx] = fVolume1 * tr.concentration_[tIdx][I1][0];
390 sStorageOfTimeIndex1[tIdx] = sVolume1 * tr.concentration_[tIdx][I1][1];
394 TracerEvaluation fVol = computeFreeVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
395 TracerEvaluation sVol = computeSolutionVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
396 dsVol_[tr.phaseIdx_][I] += sVol.value() * scvVolume - sVol1_[tr.phaseIdx_][I];
397 dfVol_[tr.phaseIdx_][I] += fVol.value() * scvVolume - fVol1_[tr.phaseIdx_][I];
398 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
400 Scalar fStorageOfTimeIndex0 = fVol.value() * tr.concentration_[tIdx][I][0];
401 Scalar fLocalStorage = (fStorageOfTimeIndex0 - fStorageOfTimeIndex1[tIdx]) * scvVolume/dt;
402 tr.residual_[tIdx][I][0] += fLocalStorage;
405 Scalar sStorageOfTimeIndex0 = sVol.value() * tr.concentration_[tIdx][I][1];
406 Scalar sLocalStorage = (sStorageOfTimeIndex0 - sStorageOfTimeIndex1[tIdx]) * scvVolume/dt;
407 tr.residual_[tIdx][I][1] += sLocalStorage;
411 (*tr.mat)[I][I][0][0] += fVol.derivative(0) * scvVolume/dt;
412 (*tr.mat)[I][I][1][1] += sVol.derivative(0) * scvVolume/dt;
417 void assembleTracerEquationFlux(TrRe& tr,
418 const ElementContext& elemCtx,
424 if (tr.numTracer() == 0)
427 TracerEvaluation fFlux;
428 TracerEvaluation sFlux;
431 computeFreeFlux_(fFlux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0);
432 computeSolFlux_(sFlux, isUpS, tr.phaseIdx_, elemCtx, scvfIdx, 0);
433 dsVol_[tr.phaseIdx_][I] += sFlux.value() * dt;
434 dfVol_[tr.phaseIdx_][I] += fFlux.value() * dt;
435 int fGlobalUpIdx = isUpF ? I : J;
436 int sGlobalUpIdx = isUpS ? I : J;
437 for (
int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
439 tr.residual_[tIdx][I][0] += fFlux.value()*tr.concentration_[tIdx][fGlobalUpIdx][0];
440 tr.residual_[tIdx][I][1] += sFlux.value()*tr.concentration_[tIdx][sGlobalUpIdx][1];
445 (*tr.mat)[J][I][0][0] = -fFlux.derivative(0);
446 (*tr.mat)[I][I][0][0] += fFlux.derivative(0);
449 (*tr.mat)[J][I][1][1] = -sFlux.derivative(0);
450 (*tr.mat)[I][I][1][1] += sFlux.derivative(0);
454 template<
class TrRe,
class Well>
455 void assembleTracerEquationWell(TrRe& tr,
458 if (tr.numTracer() == 0)
461 const auto& eclWell = well.wellEcl();
464 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
465 this->wellTracerRate_[std::make_pair(eclWell.name(), this->name(tr.idx_[tIdx]))] = 0.0;
466 this->wellFreeTracerRate_[std::make_pair(eclWell.name(), this->wellfname(tr.idx_[tIdx]))] = 0.0;
467 this->wellSolTracerRate_[std::make_pair(eclWell.name(), this->wellsname(tr.idx_[tIdx]))] = 0.0;
468 if (eclWell.isMultiSegment()) {
469 for (std::size_t i = 0; i < eclWell.getConnections().size(); ++i) {
470 this->mSwTracerRate_[std::make_tuple(eclWell.name(),
471 this->name(tr.idx_[tIdx]),
472 eclWell.getConnections().get(i).segment())] = 0.0;
477 std::vector<Scalar> wtracer(tr.numTracer());
478 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
479 wtracer[tIdx] = this->currentConcentration_(eclWell, this->
name(tr.idx_[tIdx]));
482 Scalar dt = simulator_.timeStepSize();
483 std::size_t well_index = simulator_.problem().wellModel().wellState().index(well.name()).value();
484 const auto& ws = simulator_.problem().wellModel().wellState().well(well_index);
485 for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
486 const auto I = ws.perf_data.cell_index[i];
487 Scalar rate = well.volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
489 if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
490 rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
492 else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
493 rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
499 Scalar rate_f = rate - rate_s;
501 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
503 tr.residual_[tIdx][I][0] -= rate_f*wtracer[tIdx];
506 this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) += rate_f*wtracer[tIdx];
507 this->wellFreeTracerRate_.at(std::make_pair(eclWell.name(),this->wellfname(tr.idx_[tIdx]))) += rate_f*wtracer[tIdx];
508 if (eclWell.isMultiSegment()) {
509 this->mSwTracerRate_[std::make_tuple(eclWell.name(),
510 this->name(tr.idx_[tIdx]),
511 eclWell.getConnections().get(i).segment())] += rate_f*wtracer[tIdx];
514 dfVol_[tr.phaseIdx_][I] -= rate_f * dt;
516 else if (rate_f < 0) {
517 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
519 tr.residual_[tIdx][I][0] -= rate_f * tr.concentration_[tIdx][I][0];
522 dfVol_[tr.phaseIdx_][I] -= rate_f * dt;
525 (*tr.mat)[I][I][0][0] -= rate_f * variable<TracerEvaluation>(1.0, 0).derivative(0);
528 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
530 tr.residual_[tIdx][I][1] -= rate_s * tr.concentration_[tIdx][I][1];
532 dsVol_[tr.phaseIdx_][I] -= rate_s * dt;
535 (*tr.mat)[I][I][1][1] -= rate_s * variable<TracerEvaluation>(1.0, 0).derivative(0);
541 void assembleTracerEquationSource(TrRe& tr,
545 if (tr.numTracer() == 0)
549 if (tr.phaseIdx_ == FluidSystem::waterPhaseIdx ||
550 (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && !FluidSystem::enableDissolvedGas()) ||
551 (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && !FluidSystem::enableVaporizedOil())) {
555 const Scalar& dsVol = dsVol_[tr.phaseIdx_][I];
556 const Scalar& dfVol = dfVol_[tr.phaseIdx_][I];
559 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
561 tr.residual_[tIdx][I][0] -= (dfVol / dt) * tr.concentration_[tIdx][I][0];
562 tr.residual_[tIdx][I][1] += (dfVol / dt) * tr.concentration_[tIdx][I][0];
565 tr.residual_[tIdx][I][0] += (dsVol / dt) * tr.concentration_[tIdx][I][1];
566 tr.residual_[tIdx][I][1] -= (dsVol / dt) * tr.concentration_[tIdx][I][1];
572 (*tr.mat)[I][I][0][0] -= (dfVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
573 (*tr.mat)[I][I][1][0] += (dfVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
576 (*tr.mat)[I][I][0][1] += (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
577 (*tr.mat)[I][I][1][1] -= (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
581 void assembleTracerEquations_()
589 for (
auto& tr : tbatch) {
590 if (tr.numTracer() != 0) {
592 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
593 tr.residual_[tIdx] = 0.0;
599 const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
600 for (
const auto& wellPtr : wellPtrs) {
601 for (
auto& tr : tbatch) {
602 this->assembleTracerEquationWell(tr, *wellPtr);
606 ElementContext elemCtx(simulator_);
607 for (
const auto& elem : elements(simulator_.gridView())) {
608 elemCtx.updateStencil(elem);
610 std::size_t I = elemCtx.globalSpaceIndex( 0, 0);
612 if (elem.partitionType() != Dune::InteriorEntity)
615 for (
auto& tr : tbatch) {
616 if (tr.numTracer() != 0) {
617 (*tr.mat)[I][I][0][0] = 1.;
618 (*tr.mat)[I][I][1][1] = 1.;
623 elemCtx.updateAllIntensiveQuantities();
624 elemCtx.updateAllExtensiveQuantities();
626 Scalar extrusionFactor =
627 elemCtx.intensiveQuantities( 0, 0).extrusionFactor();
628 Valgrind::CheckDefined(extrusionFactor);
629 assert(isfinite(extrusionFactor));
630 assert(extrusionFactor > 0.0);
632 elemCtx.stencil(0).subControlVolume( 0).volume()
634 Scalar dt = elemCtx.simulator().timeStepSize();
636 std::size_t I1 = elemCtx.globalSpaceIndex( 0, 1);
638 for (
auto& tr : tbatch) {
639 this->assembleTracerEquationVolume(tr, elemCtx, scvVolume, dt, I, I1);
642 std::size_t numInteriorFaces = elemCtx.numInteriorFaces(0);
643 for (
unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
644 const auto& face = elemCtx.stencil(0).interiorFace(scvfIdx);
645 unsigned j = face.exteriorIndex();
646 unsigned J = elemCtx.globalSpaceIndex( j, 0);
647 for (
auto& tr : tbatch) {
648 this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J, dt);
653 for (
auto& tr : tbatch) {
654 this->assembleTracerEquationSource(tr, dt, I);
660 for (
auto& tr : tbatch) {
661 if (tr.numTracer() == 0)
663 auto handle = VectorVectorDataHandle<GridView, std::vector<TracerVector>>(tr.residual_,
664 simulator_.gridView());
665 simulator_.gridView().communicate(handle, Dune::InteriorBorder_All_Interface,
666 Dune::ForwardCommunication);
670 void updateStorageCache()
672 for (
auto& tr : tbatch) {
673 if (tr.numTracer() != 0) {
674 tr.concentrationInitial_ = tr.concentration_;
678 ElementContext elemCtx(simulator_);
679 for (
const auto& elem : elements(simulator_.gridView())) {
680 elemCtx.updatePrimaryStencil(elem);
681 elemCtx.updatePrimaryIntensiveQuantities(0);
682 Scalar extrusionFactor = elemCtx.intensiveQuantities( 0, 0).extrusionFactor();
683 Scalar scvVolume = elemCtx.stencil(0).subControlVolume( 0).volume() * extrusionFactor;
684 int globalDofIdx = elemCtx.globalSpaceIndex(0, 0);
685 for (
auto& tr : tbatch) {
686 if (tr.numTracer() == 0)
689 Scalar fVol1 = computeFreeVolume_(tr.phaseIdx_, globalDofIdx, 0);
690 Scalar sVol1 = computeSolutionVolume_(tr.phaseIdx_, globalDofIdx, 0);
691 fVol1_[tr.phaseIdx_][globalDofIdx] = fVol1 * scvVolume;
692 sVol1_[tr.phaseIdx_][globalDofIdx] = sVol1 * scvVolume;
693 dsVol_[tr.phaseIdx_][globalDofIdx] = 0.0;
694 dfVol_[tr.phaseIdx_][globalDofIdx] = 0.0;
695 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
696 tr.storageOfTimeIndex1_[tIdx][globalDofIdx][0] = fVol1 * tr.concentrationInitial_[tIdx][globalDofIdx][0];
697 tr.storageOfTimeIndex1_[tIdx][globalDofIdx][1] = sVol1 * tr.concentrationInitial_[tIdx][globalDofIdx][1];
703 void advanceTracerFields()
705 assembleTracerEquations_();
707 for (
auto& tr : tbatch) {
708 if (tr.numTracer() == 0)
713 std::vector<TracerVector> dx(tr.concentration_);
714 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
718 bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_);
720 OpmLog::warning(
"### Tracer model: Linear solver did not converge. ###");
723 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
724 for (std::size_t globalDofIdx = 0; globalDofIdx < tr.concentration_[tIdx].size(); ++globalDofIdx) {
727 const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, 0);
728 const auto& fs = intQuants.fluidState();
729 Scalar Sf = decay<Scalar>(fs.saturation(tr.phaseIdx_));
731 if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas())
732 Ss = decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx));
733 else if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil())
734 Ss = decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx));
736 const Scalar tol_gas_sat = 1e-6;
737 if ((tr.concentration_[tIdx][globalDofIdx][0] - dx[tIdx][globalDofIdx][0] < 0.0)
739 tr.concentration_[tIdx][globalDofIdx][0] = 0.0;
741 tr.concentration_[tIdx][globalDofIdx][0] -= dx[tIdx][globalDofIdx][0];
742 if (tr.concentration_[tIdx][globalDofIdx][1] - dx[tIdx][globalDofIdx][1] < 0.0
744 tr.concentration_[tIdx][globalDofIdx][1] = 0.0;
746 tr.concentration_[tIdx][globalDofIdx][1] -= dx[tIdx][globalDofIdx][1];
749 this->freeTracerConcentration_[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][0];
750 this->solTracerConcentration_[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][1];
755 const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
756 for (
const auto& wellPtr : wellPtrs) {
757 const auto& eclWell = wellPtr->wellEcl();
760 if (!eclWell.isProducer())
763 Scalar rateWellPos = 0.0;
764 Scalar rateWellNeg = 0.0;
765 std::size_t well_index = simulator_.problem().wellModel().wellState().index(eclWell.name()).value();
766 const auto& ws = simulator_.problem().wellModel().wellState().well(well_index);
767 for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
768 const auto I = ws.perf_data.cell_index[i];
769 Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
772 if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
773 rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
775 else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
776 rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
782 Scalar rate_f = rate - rate_s;
784 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
786 this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) +=
787 rate_f * tr.concentration_[tIdx][I][0];
788 this->wellFreeTracerRate_.at(std::make_pair(eclWell.name(),this->wellfname(tr.idx_[tIdx]))) +=
789 rate_f * tr.concentration_[tIdx][I][0];
790 if (eclWell.isMultiSegment()) {
791 this->mSwTracerRate_[std::make_tuple(eclWell.name(),
792 this->name(tr.idx_[tIdx]),
793 eclWell.getConnections().get(i).segment())] +=
794 rate_f * tr.concentration_[tIdx][I][0];
799 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
801 this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) +=
802 rate_s * tr.concentration_[tIdx][I][1];
803 this->wellSolTracerRate_.at(std::make_pair(eclWell.name(),this->wellsname(tr.idx_[tIdx]))) +=
804 rate_s * tr.concentration_[tIdx][I][1];
805 if (eclWell.isMultiSegment()) {
806 this->mSwTracerRate_[std::make_tuple(eclWell.name(),
807 this->name(tr.idx_[tIdx]),
808 eclWell.getConnections().get(i).segment())] +=
809 rate_s * tr.concentration_[tIdx][I][1];
822 Scalar rateWellTotal = rateWellNeg + rateWellPos;
827 Scalar official_well_rate_total = simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_];
829 rateWellTotal = official_well_rate_total;
831 if (rateWellTotal > rateWellNeg) {
832 const Scalar bucketPrDay = 10.0/(1000.*3600.*24.);
833 const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal/rateWellNeg : 0.0;
834 for (
int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
835 this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) *= factor;
844 Simulator& simulator_;
854 template <
typename TV>
856 std::vector<int> idx_;
858 std::vector<TV> concentrationInitial_;
859 std::vector<TV> concentration_;
860 std::vector<TV> storageOfTimeIndex1_;
861 std::vector<TV> residual_;
862 std::unique_ptr<TracerMatrix> mat;
866 return this->concentrationInitial_ == rhs.concentrationInitial_ &&
867 this->concentration_ == rhs.concentration_;
873 result.idx_ = {1,2,3};
874 result.concentrationInitial_ = {5.0, 6.0};
875 result.concentration_ = {7.0, 8.0};
876 result.storageOfTimeIndex1_ = {9.0, 10.0, 11.0};
877 result.residual_ = {12.0, 13.0};
882 template<
class Serializer>
883 void serializeOp(Serializer& serializer)
885 serializer(concentrationInitial_);
886 serializer(concentration_);
889 TracerBatch(
int phaseIdx = 0) : phaseIdx_(phaseIdx) {}
891 int numTracer()
const {
return idx_.size(); }
893 void addTracer(
const int idx,
const TV & concentration)
895 int numGridDof = concentration.size();
896 idx_.emplace_back(idx);
897 concentrationInitial_.emplace_back(concentration);
898 concentration_.emplace_back(concentration);
899 residual_.emplace_back(numGridDof);
900 storageOfTimeIndex1_.emplace_back(numGridDof);
904 std::array<TracerBatch<TracerVector>,3> tbatch;
908 std::array<std::vector<Scalar>, 3> fVol1_;
909 std::array<std::vector<Scalar>, 3> sVol1_;
910 std::array<std::vector<Scalar>, 3> dsVol_;
911 std::array<std::vector<Scalar>, 3> dfVol_;