My Project
Loading...
Searching...
No Matches
TracerModel.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef OPM_TRACER_MODEL_HPP
29#define OPM_TRACER_MODEL_HPP
30
31#include <opm/common/OpmLog/OpmLog.hpp>
32#include <opm/common/TimingMacros.hpp>
33
34#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
35
37
39#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
41
42#include <array>
43#include <cstddef>
44#include <memory>
45#include <stdexcept>
46#include <string>
47#include <vector>
48
49namespace Opm::Properties {
50
51template<class TypeTag, class MyTypeTag>
53 using type = UndefinedProperty;
54};
55
56} // namespace Opm::Properties
57
58namespace Opm {
59
65template <class TypeTag>
66class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::Grid>,
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>>
72{
88
89 using TracerEvaluation = DenseAd::Evaluation<Scalar,1>;
90
91 using TracerMatrix = typename BaseType::TracerMatrix;
92 using TracerVector = typename BaseType::TracerVector;
93 using TracerVectorSingle = typename BaseType::TracerVectorSingle;
94
96 enum { numPhases = FluidSystem::numPhases };
97 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
98 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
99 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
100
101public:
102 explicit TracerModel(Simulator& simulator)
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})
110 , wat_(tbatch[0])
111 , oil_(tbatch[1])
112 , gas_(tbatch[2])
113 { }
114
115
116 /*
117 The initialization of the tracer model is a three step process:
118
119 1. The init() method is called. This will allocate buffers and initialize
120 some phase index stuff. If this is a normal run the initial tracer
121 concentrations will be assigned from the TBLK or TVDPF keywords.
122
123 2. [Restart only:] The tracer concentration are read from the restart
124 file and the concentrations are applied with repeated calls to the
125 setTracerConcentration() method. This is currently done in the
126 eclwriter::beginRestart() method.
127
128 3. Internally the tracer model manages the concentrations in "batches" for
129 the oil, water and gas tracers respectively. The batches should be
130 initialized with the initial concentration, that must be performed
131 after the concentration values have been assigned. This is done in
132 method prepareTracerBatches() called from eclproblem::finishInit().
133 */
134 void init(bool rst)
135 {
136 this->doInit(rst, simulator_.model().numGridDof(),
137 gasPhaseIdx, oilPhaseIdx, waterPhaseIdx);
138 }
139
140 void prepareTracerBatches()
141 {
142 for (std::size_t tracerIdx = 0; tracerIdx < this->tracerPhaseIdx_.size(); ++tracerIdx) {
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: " +
146 this->name(tracerIdx));
147 }
148
149 wat_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
150 }
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: " +
154 this->name(tracerIdx));
155 }
156
157 oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
158 }
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: " +
162 this->name(tracerIdx));
163 }
164
165 gas_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
166 }
167
168 // resize free and solution volume storages
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());
177 }
178
179 // will be valid after we move out of tracerMatrix_
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_);
185 }
186 else {
187 tr.mat = std::make_unique<TracerMatrix>(*base);
188 }
189 }
190 }
191 }
192
193 void beginTimeStep()
194 {
195 if (this->numTracers() == 0) {
196 return;
197 }
198
200 updateStorageCache();
201 }
202
207 {
208 if (this->numTracers() == 0) {
209 return;
210 }
211
213 advanceTracerFields();
214 }
215
220 template <class Restarter>
222 { /* not implemented */ }
223
230 template <class Restarter>
232 { /* not implemented */ }
233
234 template<class Serializer>
235 void serializeOp(Serializer& serializer)
236 {
237 serializer(static_cast<BaseType&>(*this));
238 serializer(tbatch);
239 }
240
241protected:
242 using TracerTypeIdx = typename BaseType::TracerTypeIdx;
243 using BaseType::Free;
244 using BaseType::Solution;
245
246 // compute volume associated with free/solution concentration
247 template<TracerTypeIdx Index>
248 Scalar computeVolume_(const int tracerPhaseIdx,
249 const unsigned globalDofIdx,
250 const unsigned timeIdx) const
251 {
252 const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx);
253 const auto& fs = intQuants.fluidState();
254 constexpr Scalar min_volume = 1e-10;
255
256 if constexpr (Index == Free) {
257 return std::max(decay<Scalar>(fs.saturation(tracerPhaseIdx)) *
258 decay<Scalar>(fs.invB(tracerPhaseIdx)) *
259 decay<Scalar>(intQuants.porosity()),
260 min_volume);
261 } else {
262 // vaporized oil
263 if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
264 return std::max(decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx)) *
265 decay<Scalar>(fs.invB(FluidSystem::gasPhaseIdx)) *
266 decay<Scalar>(fs.Rv()) *
267 decay<Scalar>(intQuants.porosity()),
268 min_volume);
269 }
270
271 // dissolved gas
272 else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
273 return std::max(decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx)) *
274 decay<Scalar>(fs.invB(FluidSystem::oilPhaseIdx)) *
275 decay<Scalar>(fs.Rs()) *
276 decay<Scalar>(intQuants.porosity()),
277 min_volume);
278 }
279
280 return min_volume;
281 }
282 }
283
284 template<TracerTypeIdx Index>
285 std::pair<TracerEvaluation, bool>
286 computeFlux_(const int tracerPhaseIdx,
287 const ElementContext& elemCtx,
288 const unsigned scvfIdx,
289 const unsigned timeIdx) const
290 {
291 const auto& stencil = elemCtx.stencil(timeIdx);
292 const auto& scvf = stencil.interiorFace(scvfIdx);
293
294 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
295 const unsigned inIdx = extQuants.interiorIndex();
296
297 Scalar v;
298 unsigned upIdx;
299
300 if constexpr (Index == Free) {
301 upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
302 const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
303 const auto& fs = intQuants.fluidState();
304 v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx)) *
306 } else {
307 if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
308 upIdx = extQuants.upstreamIndex(FluidSystem::gasPhaseIdx);
309
310 const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
311 const auto& fs = intQuants.fluidState();
312 v = decay<Scalar>(fs.invB(FluidSystem::gasPhaseIdx)) *
313 decay<Scalar>(extQuants.volumeFlux(FluidSystem::gasPhaseIdx)) *
314 decay<Scalar>(fs.Rv());
315 }
316 // dissolved gas
317 else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
318 upIdx = extQuants.upstreamIndex(FluidSystem::oilPhaseIdx);
319
320 const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
321 const auto& fs = intQuants.fluidState();
322 v = decay<Scalar>(fs.invB(FluidSystem::oilPhaseIdx)) *
323 decay<Scalar>(extQuants.volumeFlux(FluidSystem::oilPhaseIdx)) *
324 decay<Scalar>(fs.Rs());
325 }
326 else {
327 upIdx = 0;
328 v = 0.0;
329 }
330 }
331
332 const Scalar A = scvf.area();
333 return inIdx == upIdx
334 ? std::pair{A * v * variable<TracerEvaluation>(1.0, 0), true}
335 : std::pair{A * v, false};
336 }
337
338 template<TracerTypeIdx Index, class TrRe>
339 Scalar storage1_(const TrRe& tr,
340 const unsigned tIdx,
341 const unsigned I,
342 const unsigned I1,
343 const bool cache)
344 {
345 if (cache) {
346 return tr.storageOfTimeIndex1_[tIdx][I][Index];
347 } else {
348 return computeVolume_<Index>(tr.phaseIdx_, I1, 1) *
349 tr.concentration_[tIdx][I1][Index];
350 }
351 }
352
353 template<class TrRe>
354 void assembleTracerEquationVolume(TrRe& tr,
355 const ElementContext& elemCtx,
356 const Scalar scvVolume,
357 const Scalar dt,
358 unsigned I,
359 unsigned I1)
360
361 {
362 if (tr.numTracer() == 0) {
363 return;
364 }
365
366 const TracerEvaluation fVol = computeVolume_<Free>(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
367 const TracerEvaluation sVol = computeVolume_<Solution>(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 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];
370 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
371 // Free part
372 const Scalar fStorageOfTimeIndex0 = fVol.value() * tr.concentration_[tIdx][I][Free];
374 elemCtx.enableStorageCache())) * scvVolume / dt;
375 tr.residual_[tIdx][I][Free] += fLocalStorage; // residual + flux
376
377 // Solution part
378 const Scalar sStorageOfTimeIndex0 = sVol.value() * tr.concentration_[tIdx][I][Solution];
380 elemCtx.enableStorageCache())) * scvVolume / dt;
381 tr.residual_[tIdx][I][Solution] += sLocalStorage; // residual + flux
382 }
383
384 // Derivative matrix
385 (*tr.mat)[I][I][Free][Free] += fVol.derivative(0) * scvVolume/dt;
386 (*tr.mat)[I][I][Solution][Solution] += sVol.derivative(0) * scvVolume/dt;
387 }
388
389 template<class TrRe>
390 void assembleTracerEquationFlux(TrRe& tr,
391 const ElementContext& elemCtx,
392 unsigned scvfIdx,
393 unsigned I,
394 unsigned J,
395 const Scalar dt)
396 {
397 if (tr.numTracer() == 0) {
398 return;
399 }
400
401 const auto& [fFlux, isUpF] = computeFlux_<Free>(tr.phaseIdx_, elemCtx, scvfIdx, 0);
402 const auto& [sFlux, isUpS] = computeFlux_<Solution>(tr.phaseIdx_, elemCtx, scvfIdx, 0);
403 dVol_[Solution][tr.phaseIdx_][I] += sFlux.value() * dt;
404 dVol_[Free][tr.phaseIdx_][I] += fFlux.value() * dt;
405 const int fGlobalUpIdx = isUpF ? I : J;
406 const int sGlobalUpIdx = isUpS ? I : J;
407 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
408 // Free and solution fluxes
409 tr.residual_[tIdx][I][Free] += fFlux.value()*tr.concentration_[tIdx][fGlobalUpIdx][Free]; // residual + flux
410 tr.residual_[tIdx][I][Solution] += sFlux.value()*tr.concentration_[tIdx][sGlobalUpIdx][Solution]; // residual + flux
411 }
412
413 // Derivative matrix
414 if (isUpF){
415 (*tr.mat)[J][I][Free][Free] = -fFlux.derivative(0);
416 (*tr.mat)[I][I][Free][Free] += fFlux.derivative(0);
417 }
418 if (isUpS) {
419 (*tr.mat)[J][I][Solution][Solution] = -sFlux.derivative(0);
420 (*tr.mat)[I][I][Solution][Solution] += sFlux.derivative(0);
421 }
422 }
423
424 template<class TrRe, class Well>
425 void assembleTracerEquationWell(TrRe& tr,
426 const Well& well)
427 {
428 if (tr.numTracer() == 0) {
429 return;
430 }
431
432 const auto& eclWell = well.wellEcl();
433
434 // Init. well output to zero
435 auto& tracerRate = this->wellTracerRate_[eclWell.seqIndex()];
436 tracerRate.reserve(tr.numTracer());
437 auto& solTracerRate = this->wellTracerRate_[eclWell.seqIndex()];
438 solTracerRate.reserve(tr.numTracer());
439 auto& freeTracerRate = this->wellFreeTracerRate_[eclWell.seqIndex()];
440 freeTracerRate.reserve(tr.numTracer());
441 auto* mswTracerRate = eclWell.isMultiSegment()
442 ? &this->mSwTracerRate_[eclWell.seqIndex()]
443 : nullptr;
444 if (mswTracerRate) {
445 mswTracerRate->reserve(tr.numTracer());
446 }
447 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
448 tracerRate.emplace_back(this->name(tr.idx_[tIdx]), 0.0);
449 freeTracerRate.emplace_back(this->wellfname(tr.idx_[tIdx]), 0.0);
450 solTracerRate.emplace_back(this->wellsname(tr.idx_[tIdx]), 0.0);
451 if (eclWell.isMultiSegment()) {
452 auto& wtr = mswTracerRate->emplace_back(this->name(tr.idx_[tIdx]));
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);
456 }
457 }
458 }
459
460 std::vector<Scalar> wtracer(tr.numTracer());
461 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
462 wtracer[tIdx] = this->currentConcentration_(eclWell, this->name(tr.idx_[tIdx]));
463 }
464
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_);
470 Scalar rate_s;
471 if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
472 rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
473 }
474 else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
475 rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
476 }
477 else {
478 rate_s = 0.0;
479 }
480
481 const Scalar rate_f = rate - rate_s;
482 if (rate_f > 0) {
483 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
484 const Scalar delta = rate_f * wtracer[tIdx];
485 // Injection of free tracer only
486 tr.residual_[tIdx][I][Free] -= delta;
487
488 // Store _injector_ tracer rate for reporting
489 // (can be done here since WTRACER is constant)
490 tracerRate[tIdx].rate += delta;
491 freeTracerRate[tIdx].rate += delta;
492 if (eclWell.isMultiSegment()) {
493 (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += delta;
494 }
495 }
496 dVol_[Free][tr.phaseIdx_][I] -= rate_f * dt;
497 }
498 else if (rate_f < 0) {
499 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
500 const Scalar delta = rate_f * wtracer[tIdx];
501 // Store _injector_ tracer rate for cross-flowing well connections
502 // (can be done here since WTRACER is constant)
503 tracerRate[tIdx].rate += delta;
504 freeTracerRate[tIdx].rate += delta;
505
506 // Production of free tracer
507 tr.residual_[tIdx][I][Free] -= rate_f * tr.concentration_[tIdx][I][Free];
508 }
509 dVol_[Free][tr.phaseIdx_][I] -= rate_f * dt;
510
511 // Derivative matrix for free tracer producer
512 (*tr.mat)[I][I][Free][Free] -= rate_f * variable<TracerEvaluation>(1.0, 0).derivative(0);
513 }
514 if (rate_s < 0) {
515 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
516 // Production of solution tracer
517 tr.residual_[tIdx][I][Solution] -= rate_s * tr.concentration_[tIdx][I][Solution];
518 }
519 dVol_[Solution][tr.phaseIdx_][I] -= rate_s * dt;
520
521 // Derivative matrix for solution tracer producer
522 (*tr.mat)[I][I][Solution][Solution] -= rate_s * variable<TracerEvaluation>(1.0, 0).derivative(0);
523 }
524 }
525 }
526
527 template<class TrRe>
528 void assembleTracerEquationSource(TrRe& tr,
529 const Scalar dt,
530 unsigned I)
531 {
532 if (tr.numTracer() == 0) {
533 return;
534 }
535
536 // Skip if solution tracers do not exist
537 if (tr.phaseIdx_ == FluidSystem::waterPhaseIdx ||
538 (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && !FluidSystem::enableDissolvedGas()) ||
539 (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && !FluidSystem::enableVaporizedOil()))
540 {
541 return;
542 }
543
544 const Scalar& dsVol = dVol_[Solution][tr.phaseIdx_][I];
545 const Scalar& dfVol = dVol_[Free][tr.phaseIdx_][I];
546
547 // Source term determined by sign of dsVol: if dsVol > 0 then ms -> mf, else mf -> ms
548 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
549 if (dsVol >= 0) {
550 const auto delta = (dfVol / dt) * tr.concentration_[tIdx][I][Free];
551 tr.residual_[tIdx][I][Free] -= delta;
552 tr.residual_[tIdx][I][Solution] += delta;
553 }
554 else {
555 const auto delta = (dsVol / dt) * tr.concentration_[tIdx][I][Solution];
556 tr.residual_[tIdx][I][Free] += delta;
557 tr.residual_[tIdx][I][Solution] -= delta;
558 }
559 }
560
561 // Derivative matrix
562 if (dsVol >= 0) {
563 const auto delta = (dfVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
564 (*tr.mat)[I][I][Free][Free] -= delta;
565 (*tr.mat)[I][I][Solution][Free] += delta;
566 }
567 else {
568 const auto delta = (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
569 (*tr.mat)[I][I][Free][Solution] += delta;
570 (*tr.mat)[I][I][Solution][Solution] -= delta;
571 }
572 }
573
574 void assembleTracerEquations_()
575 {
576 // Note that we formulate the equations in terms of a concentration update
577 // (compared to previous time step) and not absolute concentration.
578 // This implies that current concentration (tr.concentration_[][]) contributes
579 // to the rhs both through storage and flux terms.
580 // Compare also advanceTracerFields(...) below.
581
582 DeferredLogger local_deferredLogger{};
583 OPM_BEGIN_PARALLEL_TRY_CATCH()
584 {
586 for (auto& tr : tbatch) {
587 if (tr.numTracer() != 0) {
588 (*tr.mat) = 0.0;
589 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
590 tr.residual_[tIdx] = 0.0;
591 }
592 }
593 }
594
595 this->wellTracerRate_.clear();
596 this->wellFreeTracerRate_.clear();
597 this->wellSolTracerRate_.clear();
598
599 // educated guess for new container size
600 const auto num_msw = this->mSwTracerRate_.size();
601 this->mSwTracerRate_.clear();
602
603 // Well terms
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);
609 for (const auto& wellPtr : wellPtrs) {
610 for (auto& tr : tbatch) {
611 this->assembleTracerEquationWell(tr, *wellPtr);
612 }
613 }
614
615 ElementContext elemCtx(simulator_);
616 const Scalar dt = elemCtx.simulator().timeStepSize();
617 for (const auto& elem : elements(simulator_.gridView())) {
618 elemCtx.updateStencil(elem);
619
620 const std::size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timeIdx=*/0);
621
622 if (elem.partitionType() != Dune::InteriorEntity) {
623 // Dirichlet boundary conditions needed for the parallel matrix
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.;
628 }
629 }
630 continue;
631 }
632 elemCtx.updateAllIntensiveQuantities();
633 elemCtx.updateAllExtensiveQuantities();
634
635 const Scalar extrusionFactor =
636 elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
637 Valgrind::CheckDefined(extrusionFactor);
638 assert(isfinite(extrusionFactor));
639 assert(extrusionFactor > 0.0);
640 const Scalar scvVolume =
641 elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume()
642 * extrusionFactor;
643
644 const std::size_t I1 = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timeIdx=*/1);
645
646 for (auto& tr : tbatch) {
647 this->assembleTracerEquationVolume(tr, elemCtx, scvVolume, dt, I, I1);
648 }
649
650 const std::size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
651 for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
652 const auto& face = elemCtx.stencil(0).interiorFace(scvfIdx);
653 const unsigned j = face.exteriorIndex();
654 const unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
655 for (auto& tr : tbatch) {
656 this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J, dt);
657 }
658 }
659
660 // Source terms (mass transfer between free and solution tracer)
661 for (auto& tr : tbatch) {
662 this->assembleTracerEquationSource(tr, dt, I);
663 }
664 }
665 }
666 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
667 "assembleTracerEquations() failed: ",
668 true, simulator_.gridView().comm())
669
670 // Communicate overlap using grid Communication
671 for (auto& tr : tbatch) {
672 if (tr.numTracer() == 0) {
673 continue;
674 }
676 simulator_.gridView());
677 simulator_.gridView().communicate(handle, Dune::InteriorBorder_All_Interface,
678 Dune::ForwardCommunication);
679 }
680 }
681
682 template<TracerTypeIdx Index, class TrRe>
683 void updateElem(TrRe& tr,
684 const Scalar scvVolume,
685 const unsigned globalDofIdx)
686 {
687 const Scalar vol1 = computeVolume_<Index>(tr.phaseIdx_, globalDofIdx, 0);
688 vol1_[Index][tr.phaseIdx_][globalDofIdx] = vol1 * scvVolume;
689 dVol_[Index][tr.phaseIdx_][globalDofIdx] = 0.0;
690 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
691 tr.storageOfTimeIndex1_[tIdx][globalDofIdx][Index] =
692 vol1 * tr.concentrationInitial_[tIdx][globalDofIdx][Index];
693 }
694 }
695
696 void updateStorageCache()
697 {
698 for (auto& tr : tbatch) {
699 if (tr.numTracer() != 0) {
700 tr.concentrationInitial_ = tr.concentration_;
701 }
702 }
703
704 ElementContext elemCtx(simulator_);
705 for (const auto& elem : elements(simulator_.gridView())) {
706 elemCtx.updatePrimaryStencil(elem);
707 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
708 const Scalar extrusionFactor = elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
709 const Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume() * extrusionFactor;
710 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(0, /*timeIdx=*/0);
711 for (auto& tr : tbatch) {
712 if (tr.numTracer() == 0) {
713 continue;
714 }
715 updateElem<Free>(tr, scvVolume, globalDofIdx);
716 updateElem<Solution>(tr, scvVolume, globalDofIdx);
717 }
718 }
719 }
720
721 template<TracerTypeIdx Index, class TrRe>
722 void copyForOutput(TrRe& tr,
723 const std::vector<TracerVector>& dx,
724 const Scalar S,
725 const unsigned tIdx,
726 const unsigned globalDofIdx,
727 std::vector<TracerVectorSingle>& sc)
728 {
729 constexpr Scalar tol_gas_sat = 1e-6;
730 tr.concentration_[tIdx][globalDofIdx][Index] -= dx[tIdx][globalDofIdx][Index];
731 if (tr.concentration_[tIdx][globalDofIdx][Index] < 0.0 || S < tol_gas_sat) {
732 tr.concentration_[tIdx][globalDofIdx][Index] = 0.0;
733 }
734 sc[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][Index];
735 }
736
737 template<TracerTypeIdx Index, class TrRe>
738 void assignRates(const TrRe& tr,
739 const Well& eclWell,
740 const std::size_t i,
741 const std::size_t I,
742 const Scalar rate,
745 std::vector<WellTracerRate<Scalar>>& splitRate)
746 {
747 if (rate < 0) {
748 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
749 // Store _producer_ free tracer rate for reporting
750 const Scalar delta = rate * tr.concentration_[tIdx][I][Index];
751 tracerRate[tIdx].rate += delta;
752 splitRate[tIdx].rate += delta;
753 if (eclWell.isMultiSegment()) {
754 (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += delta;
755 }
756 }
757 }
758 }
759
760 void advanceTracerFields()
761 {
762 assembleTracerEquations_();
763
764 for (auto& tr : tbatch) {
765 if (tr.numTracer() == 0) {
766 continue;
767 }
768
769 // Note that we solve for a concentration update (compared to previous time step)
770 // Confer also assembleTracerEquations_(...) above.
771 std::vector<TracerVector> dx(tr.concentration_);
772 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
773 dx[tIdx] = 0.0;
774 }
775
776 const bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_);
777 if (!converged) {
778 OpmLog::warning("### Tracer model: Linear solver did not converge. ###");
779 }
780
782
783 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
784 for (std::size_t globalDofIdx = 0; globalDofIdx < tr.concentration_[tIdx].size(); ++globalDofIdx) {
785 // New concetration. Concentrations that are negative or where free/solution phase is not
786 // present are set to zero
787 const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, 0);
788 const auto& fs = intQuants.fluidState();
789 const Scalar Sf = decay<Scalar>(fs.saturation(tr.phaseIdx_));
790 Scalar Ss = 0.0;
791
792 if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
793 Ss = decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx));
794 }
795 else if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
796 Ss = decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx));
797 }
798
799 copyForOutput<Free>(tr, dx, Sf, tIdx, globalDofIdx, this->freeTracerConcentration_);
800 copyForOutput<Solution>(tr, dx, Ss, tIdx, globalDofIdx, this->solTracerConcentration_);
801 }
802 }
803
804 // Store _producer_ tracer rate for reporting
805 const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
806 for (const auto& wellPtr : wellPtrs) {
807 const auto& eclWell = wellPtr->wellEcl();
808
809 // Injection rates already reported during assembly
810 if (!eclWell.isProducer()) {
811 continue;
812 }
813
814 Scalar rateWellPos = 0.0;
815 Scalar rateWellNeg = 0.0;
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);
818 auto& tracerRate = this->wellTracerRate_[eclWell.seqIndex()];
819 auto& freeTracerRate = this->wellFreeTracerRate_[eclWell.seqIndex()];
820 auto& solTracerRate = this->wellSolTracerRate_[eclWell.seqIndex()];
821 auto* mswTracerRate = eclWell.isMultiSegment() ? &this->mSwTracerRate_[eclWell.seqIndex()] : nullptr;
822
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_);
826
827 Scalar rate_s;
828 if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
829 rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
830 }
831 else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
832 rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
833 }
834 else {
835 rate_s = 0.0;
836 }
837
838 const Scalar rate_f = rate - rate_s;
843
844 if (rate < 0) {
845 rateWellNeg += rate;
846 }
847 else {
848 rateWellPos += rate;
849 }
850 }
851
852 // TODO: Some inconsistencies here that perhaps should be clarified.
853 // The "offical" rate as reported below is occasionally significant
854 // different from the sum over connections (as calculated above). Only observed
855 // for small values, neglible for the rate itself, but matters when used to
856 // calculate tracer concentrations.
857 const Scalar official_well_rate_total =
858 simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_];
859
861
862 if (rateWellTotal > rateWellNeg) { // Cross flow
863 constexpr Scalar bucketPrDay = 10.0 / (1000. * 3600. * 24.); // ... keeps (some) trouble away
864 const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal / rateWellNeg : 0.0;
865 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
866 tracerRate[tIdx].rate *= factor;
867 }
868 }
869 }
870 }
871 }
872
873 Simulator& simulator_;
874
875 // This struct collects tracers of the same type (i.e, transported in same phase).
876 // The idea being that, under the assumption of linearity, tracers of same type can
877 // be solved in concert, having a common system matrix but separate right-hand-sides.
878
879 // Since oil or gas tracers appears in dual compositions when VAPOIL respectively DISGAS
880 // is active, the template argument is intended to support future extension to these
881 // scenarios by supplying an extended vector type.
882
883 template <typename TV>
885 {
886 std::vector<int> idx_;
887 const int phaseIdx_;
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;
893
894 bool operator==(const TracerBatch& rhs) const
895 {
896 return this->concentrationInitial_ == rhs.concentrationInitial_ &&
897 this->concentration_ == rhs.concentration_;
898 }
899
900 static TracerBatch serializationTestObject()
901 {
903 result.idx_ = {1,2,3};
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};
908
909 return result;
910 }
911
912 template<class Serializer>
913 void serializeOp(Serializer& serializer)
914 {
915 serializer(concentrationInitial_);
916 serializer(concentration_);
917 }
918
919 TracerBatch(int phaseIdx = 0) : phaseIdx_(phaseIdx) {}
920
921 int numTracer() const
922 { return idx_.size(); }
923
924 void addTracer(const int idx, const TV& concentration)
925 {
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);
932 }
933 };
934
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_;
941};
942
943} // namespace Opm
944
945#endif // OPM_TRACER_MODEL_HPP
A class which handles tracers as specified in by ECL.
A datahandle sending data located in multiple vectors.
Definition GenericTracerModel.hpp:56
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:225
A class which handles tracers as specified in by ECL.
Definition TracerModel.hpp:72
void serialize(Restarter &)
This method writes the complete state of all tracer to the hard disk.
Definition TracerModel.hpp:221
void endTimeStep()
Informs the tracer model that a time step has just been finished.
Definition TracerModel.hpp:206
void deserialize(Restarter &)
This method restores the complete state of the tracer from disk.
Definition TracerModel.hpp:231
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition propertysystem.hh:235
int Index
The type of an index of a degree of freedom.
Definition overlaptypes.hh:44
The Opm property system, traits with inheritance.
Definition TracerModel.hpp:52
a tag to mark properties as undefined
Definition propertysystem.hh:40
Definition TracerModel.hpp:885