67 GetPropType<TypeTag, Properties::GridView>,
68 GetPropType<TypeTag, Properties::DofMapper>,
69 GetPropType<TypeTag, Properties::Stencil>,
70 GetPropType<TypeTag, Properties::FluidSystem>,
71 GetPropType<TypeTag, Properties::Scalar>>
89 using TracerEvaluation = DenseAd::Evaluation<Scalar,1>;
91 using TracerMatrix =
typename BaseType::TracerMatrix;
92 using TracerVector =
typename BaseType::TracerVector;
93 using TracerVectorSingle =
typename BaseType::TracerVectorSingle;
96 enum { numPhases = FluidSystem::numPhases };
97 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
98 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
99 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
103 :
BaseType(simulator.vanguard().gridView(),
104 simulator.vanguard().eclState(),
105 simulator.vanguard().cartesianIndexMapper(),
106 simulator.model().dofMapper(),
107 simulator.vanguard().cellCentroids())
108 , simulator_(simulator)
109 , tbatch({waterPhaseIdx, oilPhaseIdx, gasPhaseIdx})
136 this->
doInit(rst, simulator_.model().numGridDof(),
137 gasPhaseIdx, oilPhaseIdx, waterPhaseIdx);
140 void prepareTracerBatches()
143 if (this->tracerPhaseIdx_[
tracerIdx] == FluidSystem::waterPhaseIdx) {
144 if (! FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)){
145 throw std::runtime_error(
"Water tracer specified for non-water fluid system: " +
151 else if (this->tracerPhaseIdx_[
tracerIdx] == FluidSystem::oilPhaseIdx) {
152 if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
153 throw std::runtime_error(
"Oil tracer specified for non-oil fluid system: " +
159 else if (this->tracerPhaseIdx_[
tracerIdx] == FluidSystem::gasPhaseIdx) {
160 if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
161 throw std::runtime_error(
"Gas tracer specified for non-gas fluid system: " +
169 vol1_[0][this->tracerPhaseIdx_[
tracerIdx]].
170 resize(this->freeTracerConcentration_[
tracerIdx].size());
171 vol1_[1][this->tracerPhaseIdx_[
tracerIdx]].
172 resize(this->freeTracerConcentration_[
tracerIdx].size());
173 dVol_[0][this->tracerPhaseIdx_[
tracerIdx]].
174 resize(this->solTracerConcentration_[
tracerIdx].size());
175 dVol_[1][this->tracerPhaseIdx_[
tracerIdx]].
176 resize(this->solTracerConcentration_[
tracerIdx].size());
180 TracerMatrix*
base = this->tracerMatrix_.get();
181 for (
auto&
tr : this->tbatch) {
182 if (
tr.numTracer() != 0) {
183 if (this->tracerMatrix_) {
184 tr.mat = std::move(this->tracerMatrix_);
187 tr.mat = std::make_unique<TracerMatrix>(*
base);
200 updateStorageCache();
213 advanceTracerFields();
220 template <
class Restarter>
230 template <
class Restarter>
234 template<
class Serializer>
243 using BaseType::Free;
244 using BaseType::Solution;
247 template<TracerTypeIdx Index>
249 const unsigned globalDofIdx,
252 const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx,
timeIdx);
253 const auto& fs = intQuants.fluidState();
256 if constexpr (Index == Free) {
263 if (
tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
264 return std::max(
decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx)) *
272 else if (
tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
273 return std::max(
decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx)) *
284 template<TracerTypeIdx Index>
285 std::pair<TracerEvaluation, bool>
287 const ElementContext& elemCtx,
291 const auto& stencil = elemCtx.stencil(
timeIdx);
300 if constexpr (
Index == Free) {
302 const auto& intQuants = elemCtx.intensiveQuantities(
upIdx,
timeIdx);
303 const auto& fs = intQuants.fluidState();
307 if (
tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
310 const auto& intQuants = elemCtx.intensiveQuantities(
upIdx,
timeIdx);
311 const auto& fs = intQuants.fluidState();
317 else if (
tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
320 const auto& intQuants = elemCtx.intensiveQuantities(
upIdx,
timeIdx);
321 const auto& fs = intQuants.fluidState();
332 const Scalar
A =
scvf.area();
335 : std::pair{
A *
v,
false};
338 template<TracerTypeIdx Index,
class TrRe>
339 Scalar storage1_(
const TrRe&
tr,
354 void assembleTracerEquationVolume(
TrRe&
tr,
355 const ElementContext& elemCtx,
362 if (
tr.numTracer() == 0) {
368 dVol_[Solution][
tr.phaseIdx_][I] +=
sVol.value() *
scvVolume - vol1_[1][
tr.phaseIdx_][I];
369 dVol_[Free][
tr.phaseIdx_][I] +=
fVol.value() *
scvVolume - vol1_[0][
tr.phaseIdx_][I];
390 void assembleTracerEquationFlux(
TrRe&
tr,
391 const ElementContext& elemCtx,
397 if (
tr.numTracer() == 0) {
403 dVol_[Solution][
tr.phaseIdx_][I] +=
sFlux.value() *
dt;
404 dVol_[Free][
tr.phaseIdx_][I] +=
fFlux.value() *
dt;
415 (*
tr.mat)[J][I][Free][Free] = -
fFlux.derivative(0);
416 (*
tr.mat)[I][I][Free][Free] +=
fFlux.derivative(0);
419 (*
tr.mat)[J][I][Solution][Solution] = -
sFlux.derivative(0);
420 (*
tr.mat)[I][I][Solution][Solution] +=
sFlux.derivative(0);
424 template<
class TrRe,
class Well>
425 void assembleTracerEquationWell(
TrRe&
tr,
428 if (
tr.numTracer() == 0) {
432 const auto&
eclWell = well.wellEcl();
442 ? &this->mSwTracerRate_[
eclWell.seqIndex()]
451 if (
eclWell.isMultiSegment()) {
453 wtr.rate.reserve(
eclWell.getConnections().size());
454 for (std::size_t i = 0; i <
eclWell.getConnections().size(); ++i) {
455 wtr.rate.emplace(
eclWell.getConnections().get(i).segment(), 0.0);
460 std::vector<Scalar>
wtracer(
tr.numTracer());
465 const Scalar
dt = simulator_.timeStepSize();
466 const auto&
ws = simulator_.problem().wellModel().wellState().well(well.name());
467 for (std::size_t i = 0; i <
ws.perf_data.size(); ++i) {
468 const auto I =
ws.perf_data.cell_index[i];
469 const Scalar rate = well.volumetricSurfaceRateForConnection(I,
tr.phaseIdx_);
471 if (
tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
472 rate_s =
ws.perf_data.phase_mixing_rates[i][
ws.vaporized_oil];
474 else if (
tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
475 rate_s =
ws.perf_data.phase_mixing_rates[i][
ws.dissolved_gas];
492 if (
eclWell.isMultiSegment()) {
493 (*mswTracerRate)[
tIdx].rate[
eclWell.getConnections().get(i).segment()] +=
delta;
519 dVol_[Solution][
tr.phaseIdx_][I] -=
rate_s *
dt;
528 void assembleTracerEquationSource(
TrRe&
tr,
532 if (
tr.numTracer() == 0) {
537 if (
tr.phaseIdx_ == FluidSystem::waterPhaseIdx ||
538 (
tr.phaseIdx_ == FluidSystem::gasPhaseIdx && !FluidSystem::enableDissolvedGas()) ||
539 (
tr.phaseIdx_ == FluidSystem::oilPhaseIdx && !FluidSystem::enableVaporizedOil()))
544 const Scalar&
dsVol = dVol_[Solution][
tr.phaseIdx_][I];
545 const Scalar&
dfVol = dVol_[Free][
tr.phaseIdx_][I];
564 (*
tr.mat)[I][I][Free][Free] -=
delta;
565 (*
tr.mat)[I][I][Solution][Free] +=
delta;
569 (*
tr.mat)[I][I][Free][Solution] +=
delta;
570 (*
tr.mat)[I][I][Solution][Solution] -=
delta;
574 void assembleTracerEquations_()
583 OPM_BEGIN_PARALLEL_TRY_CATCH()
586 for (
auto&
tr : tbatch) {
587 if (
tr.numTracer() != 0) {
595 this->wellTracerRate_.clear();
596 this->wellFreeTracerRate_.clear();
597 this->wellSolTracerRate_.clear();
600 const auto num_msw = this->mSwTracerRate_.size();
601 this->mSwTracerRate_.clear();
604 const auto&
wellPtrs = simulator_.problem().wellModel().localNonshutWells();
605 this->wellTracerRate_.reserve(
wellPtrs.size());
606 this->wellFreeTracerRate_.reserve(
wellPtrs.size());
607 this->wellSolTracerRate_.reserve(
wellPtrs.size());
608 this->mSwTracerRate_.reserve(
num_msw);
610 for (
auto&
tr : tbatch) {
611 this->assembleTracerEquationWell(
tr, *
wellPtr);
615 ElementContext elemCtx(simulator_);
616 const Scalar
dt = elemCtx.simulator().timeStepSize();
617 for (
const auto&
elem :
elements(simulator_.gridView())) {
618 elemCtx.updateStencil(
elem);
620 const std::size_t I = elemCtx.globalSpaceIndex( 0, 0);
622 if (
elem.partitionType() != Dune::InteriorEntity) {
624 for (
const auto&
tr : tbatch) {
625 if (
tr.numTracer() != 0) {
626 (*
tr.mat)[I][I][0][0] = 1.;
627 (*
tr.mat)[I][I][1][1] = 1.;
632 elemCtx.updateAllIntensiveQuantities();
633 elemCtx.updateAllExtensiveQuantities();
635 const Scalar extrusionFactor =
636 elemCtx.intensiveQuantities( 0, 0).extrusionFactor();
637 Valgrind::CheckDefined(extrusionFactor);
639 assert(extrusionFactor > 0.0);
641 elemCtx.stencil(0).subControlVolume( 0).volume()
644 const std::size_t
I1 = elemCtx.globalSpaceIndex( 0, 1);
646 for (
auto&
tr : tbatch) {
650 const std::size_t numInteriorFaces = elemCtx.numInteriorFaces(0);
652 const auto&
face = elemCtx.stencil(0).interiorFace(
scvfIdx);
653 const unsigned j =
face.exteriorIndex();
654 const unsigned J = elemCtx.globalSpaceIndex( j, 0);
655 for (
auto&
tr : tbatch) {
656 this->assembleTracerEquationFlux(
tr, elemCtx,
scvfIdx, I, J,
dt);
661 for (
auto&
tr : tbatch) {
662 this->assembleTracerEquationSource(
tr,
dt, I);
667 "assembleTracerEquations() failed: ",
668 true, simulator_.gridView().comm())
672 if (
tr.numTracer() == 0) {
676 simulator_.gridView());
677 simulator_.gridView().communicate(handle, Dune::InteriorBorder_All_Interface,
678 Dune::ForwardCommunication);
682 template<TracerTypeIdx Index,
class TrRe>
685 const unsigned globalDofIdx)
689 dVol_[
Index][
tr.phaseIdx_][globalDofIdx] = 0.0;
691 tr.storageOfTimeIndex1_[
tIdx][globalDofIdx][
Index] =
696 void updateStorageCache()
698 for (
auto&
tr : tbatch) {
699 if (
tr.numTracer() != 0) {
700 tr.concentrationInitial_ =
tr.concentration_;
704 ElementContext elemCtx(simulator_);
705 for (
const auto&
elem :
elements(simulator_.gridView())) {
706 elemCtx.updatePrimaryStencil(
elem);
707 elemCtx.updatePrimaryIntensiveQuantities(0);
708 const Scalar extrusionFactor = elemCtx.intensiveQuantities( 0, 0).extrusionFactor();
709 const Scalar
scvVolume = elemCtx.stencil(0).subControlVolume( 0).volume() * extrusionFactor;
710 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(0, 0);
711 for (
auto&
tr : tbatch) {
712 if (
tr.numTracer() == 0) {
721 template<TracerTypeIdx Index,
class TrRe>
722 void copyForOutput(
TrRe&
tr,
723 const std::vector<TracerVector>&
dx,
726 const unsigned globalDofIdx,
727 std::vector<TracerVectorSingle>&
sc)
732 tr.concentration_[
tIdx][globalDofIdx][
Index] = 0.0;
737 template<TracerTypeIdx Index,
class TrRe>
738 void assignRates(
const TrRe&
tr,
753 if (
eclWell.isMultiSegment()) {
754 (*mswTracerRate)[
tIdx].rate[
eclWell.getConnections().get(i).segment()] +=
delta;
760 void advanceTracerFields()
762 assembleTracerEquations_();
764 for (
auto&
tr : tbatch) {
765 if (
tr.numTracer() == 0) {
771 std::vector<TracerVector>
dx(
tr.concentration_);
776 const bool converged = this->linearSolveBatchwise_(*
tr.mat,
dx,
tr.residual_);
778 OpmLog::warning(
"### Tracer model: Linear solver did not converge. ###");
784 for (std::size_t globalDofIdx = 0; globalDofIdx <
tr.concentration_[
tIdx].size(); ++globalDofIdx) {
787 const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, 0);
788 const auto& fs = intQuants.fluidState();
792 if (
tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
795 else if (
tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
805 const auto&
wellPtrs = simulator_.problem().wellModel().localNonshutWells();
816 const std::size_t well_index = simulator_.problem().wellModel().wellState().index(
eclWell.name()).value();
817 const auto&
ws = simulator_.problem().wellModel().wellState().well(well_index);
823 for (std::size_t i = 0; i <
ws.perf_data.size(); ++i) {
824 const auto I =
ws.perf_data.cell_index[i];
825 const Scalar rate =
wellPtr->volumetricSurfaceRateForConnection(I,
tr.phaseIdx_);
828 if (
tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
829 rate_s =
ws.perf_data.phase_mixing_rates[i][
ws.vaporized_oil];
831 else if (
tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
832 rate_s =
ws.perf_data.phase_mixing_rates[i][
ws.dissolved_gas];
858 simulator_.problem().wellModel().wellState().well(well_index).surface_rates[
tr.phaseIdx_];
863 constexpr Scalar
bucketPrDay = 10.0 / (1000. * 3600. * 24.);
873 Simulator& simulator_;
883 template <
typename TV>
886 std::vector<int> idx_;
888 std::vector<TV> concentrationInitial_;
889 std::vector<TV> concentration_;
890 std::vector<TV> storageOfTimeIndex1_;
891 std::vector<TV> residual_;
892 std::unique_ptr<TracerMatrix> mat;
896 return this->concentrationInitial_ == rhs.concentrationInitial_ &&
897 this->concentration_ == rhs.concentration_;
904 result.concentrationInitial_ = {5.0, 6.0};
905 result.concentration_ = {7.0, 8.0};
906 result.storageOfTimeIndex1_ = {9.0, 10.0, 11.0};
907 result.residual_ = {12.0, 13.0};
912 template<
class Serializer>
921 int numTracer()
const
922 {
return idx_.size(); }
924 void addTracer(
const int idx,
const TV& concentration)
926 const int numGridDof = concentration.size();
927 idx_.emplace_back(
idx);
928 concentrationInitial_.emplace_back(concentration);
929 concentration_.emplace_back(concentration);
930 residual_.emplace_back(numGridDof);
931 storageOfTimeIndex1_.emplace_back(numGridDof);
935 std::array<TracerBatch<TracerVector>,numPhases> tbatch;
939 std::array<std::array<std::vector<Scalar>,numPhases>,2> vol1_;
940 std::array<std::array<std::vector<Scalar>,numPhases>,2> dVol_;