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()
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));
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));
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));
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>
252 if (
tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
261 else if (
tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
276 void computeFreeFlux_(TracerEvaluation &
freeFlux,
298 Scalar
A =
scvf.area();
309 void computeSolFlux_(TracerEvaluation &
solFlux,
326 if (
tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
337 else if (
tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
352 Scalar
A =
scvf.area();
364 void assembleTracerEquationVolume(
TrRe&
tr,
372 if (
tr.numTracer() == 0)
379 if (
elemCtx.enableStorageCache()) {
386 Scalar
fVolume1 = computeFreeVolume_(
tr.phaseIdx_,
I1, 1);
387 Scalar
sVolume1 = computeSolutionVolume_(
tr.phaseIdx_,
I1, 1);
417 void assembleTracerEquationFlux(
TrRe&
tr,
424 if (
tr.numTracer() == 0)
427 TracerEvaluation
fFlux;
428 TracerEvaluation
sFlux;
433 dsVol_[
tr.phaseIdx_][
I] +=
sFlux.value() *
dt;
434 dfVol_[
tr.phaseIdx_][
I] +=
fFlux.value() *
dt;
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();
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(),
472 eclWell.getConnections().get(i).segment())] = 0.0;
477 std::vector<Scalar>
wtracer(
tr.numTracer());
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];
508 if (
eclWell.isMultiSegment()) {
509 this->mSwTracerRate_[std::make_tuple(
eclWell.name(),
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];
581 void assembleTracerEquations_()
589 for (
auto&
tr : tbatch) {
590 if (
tr.numTracer() != 0) {
599 const auto&
wellPtrs = simulator_.problem().wellModel().localNonshutWells();
601 for (
auto&
tr : tbatch) {
602 this->assembleTracerEquationWell(
tr, *
wellPtr);
606 ElementContext
elemCtx(simulator_);
607 for (
const auto&
elem :
elements(simulator_.gridView())) {
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);
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) {
642 std::size_t numInteriorFaces =
elemCtx.numInteriorFaces(0);
645 unsigned j =
face.exteriorIndex();
646 unsigned J =
elemCtx.globalSpaceIndex( j, 0);
647 for (
auto&
tr : tbatch) {
653 for (
auto&
tr : tbatch) {
654 this->assembleTracerEquationSource(
tr,
dt,
I);
660 for (
auto&
tr : tbatch) {
661 if (
tr.numTracer() == 0)
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())) {
681 elemCtx.updatePrimaryIntensiveQuantities(0);
682 Scalar extrusionFactor =
elemCtx.intensiveQuantities( 0, 0).extrusionFactor();
683 Scalar
scvVolume =
elemCtx.stencil(0).subControlVolume( 0).volume() * extrusionFactor;
685 for (
auto&
tr : tbatch) {
686 if (
tr.numTracer() == 0)
703 void advanceTracerFields()
705 assembleTracerEquations_();
707 for (
auto&
tr : tbatch) {
708 if (
tr.numTracer() == 0)
713 std::vector<TracerVector>
dx(
tr.concentration_);
718 bool converged = this->linearSolveBatchwise_(*
tr.mat,
dx,
tr.residual_);
720 OpmLog::warning(
"### Tracer model: Linear solver did not converge. ###");
731 if (
tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas())
733 else if (
tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil())
755 const auto&
wellPtrs = simulator_.problem().wellModel().localNonshutWells();
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];
786 this->wellTracerRate_.at(std::make_pair(
eclWell.name(),
this->name(
tr.idx_[
tIdx]))) +=
788 this->wellFreeTracerRate_.at(std::make_pair(
eclWell.name(),
this->wellfname(
tr.idx_[
tIdx]))) +=
790 if (
eclWell.isMultiSegment()) {
791 this->mSwTracerRate_[std::make_tuple(
eclWell.name(),
793 eclWell.getConnections().get(i).segment())] +=
801 this->wellTracerRate_.at(std::make_pair(
eclWell.name(),
this->name(
tr.idx_[
tIdx]))) +=
803 this->wellSolTracerRate_.at(std::make_pair(
eclWell.name(),
this->wellsname(
tr.idx_[
tIdx]))) +=
805 if (
eclWell.isMultiSegment()) {
806 this->mSwTracerRate_[std::make_tuple(
eclWell.name(),
808 eclWell.getConnections().get(i).segment())] +=
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_;
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>
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_;