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
36#include <opm/grid/utility/ElementChunks.hpp>
37
40
44
45#include <array>
46#include <cstddef>
47#include <memory>
48#include <stdexcept>
49#include <string>
50#include <vector>
51
52namespace Opm::Properties {
53
54template<class TypeTag, class MyTypeTag>
57};
58
59} // namespace Opm::Properties
60
61namespace Opm {
62
68template <class TypeTag>
69class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::Grid>,
70 GetPropType<TypeTag, Properties::GridView>,
71 GetPropType<TypeTag, Properties::DofMapper>,
72 GetPropType<TypeTag, Properties::Stencil>,
73 GetPropType<TypeTag, Properties::FluidSystem>,
74 GetPropType<TypeTag, Properties::Scalar>>
75{
91
92 using TracerEvaluation = DenseAd::Evaluation<Scalar,1>;
93
94 using TracerMatrix = typename BaseType::TracerMatrix;
95 using TracerVector = typename BaseType::TracerVector;
96 using TracerVectorSingle = typename BaseType::TracerVectorSingle;
97
98 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
99 enum { numPhases = FluidSystem::numPhases };
100 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
101 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
102 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
103
104public:
105 explicit TracerModel(Simulator& simulator)
106 : BaseType(simulator.vanguard().gridView(),
107 simulator.vanguard().eclState(),
108 simulator.vanguard().cartesianIndexMapper(),
109 simulator.model().dofMapper(),
110 simulator.vanguard().cellCentroids())
111 , simulator_(simulator)
112 , tbatch({waterPhaseIdx, oilPhaseIdx, gasPhaseIdx})
113 , wat_(tbatch[0])
114 , oil_(tbatch[1])
115 , gas_(tbatch[2])
116 , element_chunks_(simulator.gridView(), Dune::Partitions::all, ThreadManager::maxThreads())
117 { }
118
119
120 /*
121 The initialization of the tracer model is a three step process:
122
123 1. The init() method is called. This will allocate buffers and initialize
124 some phase index stuff. If this is a normal run the initial tracer
125 concentrations will be assigned from the TBLK or TVDPF keywords.
126
127 2. [Restart only:] The tracer concentration are read from the restart
128 file and the concentrations are applied with repeated calls to the
129 setTracerConcentration() method. This is currently done in the
130 eclwriter::beginRestart() method.
131
132 3. Internally the tracer model manages the concentrations in "batches" for
133 the oil, water and gas tracers respectively. The batches should be
134 initialized with the initial concentration, that must be performed
135 after the concentration values have been assigned. This is done in
136 method prepareTracerBatches() called from eclproblem::finishInit().
137 */
138 void init(bool rst)
139 {
140 this->doInit(rst, simulator_.model().numGridDof(),
141 gasPhaseIdx, oilPhaseIdx, waterPhaseIdx);
142 }
143
145 {
146 for (std::size_t tracerIdx = 0; tracerIdx < this->tracerPhaseIdx_.size(); ++tracerIdx) {
147 if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::waterPhaseIdx) {
148 if (! FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)){
149 throw std::runtime_error("Water tracer specified for non-water fluid system: " +
150 this->name(tracerIdx));
151 }
152
153 wat_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
154 }
155 else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::oilPhaseIdx) {
156 if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
157 throw std::runtime_error("Oil tracer specified for non-oil fluid system: " +
158 this->name(tracerIdx));
159 }
160
161 oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
162 }
163 else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::gasPhaseIdx) {
164 if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
165 throw std::runtime_error("Gas tracer specified for non-gas fluid system: " +
166 this->name(tracerIdx));
167 }
168
169 gas_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
170 }
171
172 // resize free and solution volume storages
173 vol1_[0][this->tracerPhaseIdx_[tracerIdx]].
174 resize(this->freeTracerConcentration_[tracerIdx].size());
175 vol1_[1][this->tracerPhaseIdx_[tracerIdx]].
176 resize(this->freeTracerConcentration_[tracerIdx].size());
177 dVol_[0][this->tracerPhaseIdx_[tracerIdx]].
178 resize(this->solTracerConcentration_[tracerIdx].size());
179 dVol_[1][this->tracerPhaseIdx_[tracerIdx]].
180 resize(this->solTracerConcentration_[tracerIdx].size());
181 }
182
183 // will be valid after we move out of tracerMatrix_
184 TracerMatrix* base = this->tracerMatrix_.get();
185 for (auto& tr : this->tbatch) {
186 if (tr.numTracer() != 0) {
187 if (this->tracerMatrix_) {
188 tr.mat = std::move(this->tracerMatrix_);
189 }
190 else {
191 tr.mat = std::make_unique<TracerMatrix>(*base);
192 }
193 }
194 }
195 }
196
198 {
199 if (this->numTracers() == 0) {
200 return;
201 }
202
203 OPM_TIMEBLOCK(tracerUpdateCache);
205 }
206
211 {
212 if (this->numTracers() == 0) {
213 return;
214 }
215
216 OPM_TIMEBLOCK(tracerAdvance);
218 }
219
224 template <class Restarter>
225 void serialize(Restarter&)
226 { /* not implemented */ }
227
234 template <class Restarter>
235 void deserialize(Restarter&)
236 { /* not implemented */ }
237
238 template<class Serializer>
239 void serializeOp(Serializer& serializer)
240 {
241 serializer(static_cast<BaseType&>(*this));
242 serializer(tbatch);
243 }
244
245protected:
247 using BaseType::Free;
248 using BaseType::Solution;
249
250 // compute volume associated with free/solution concentration
251 template<TracerTypeIdx Index>
252 Scalar computeVolume_(const int tracerPhaseIdx,
253 const unsigned globalDofIdx,
254 const unsigned timeIdx) const
255 {
256 const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx);
257 const auto& fs = intQuants.fluidState();
258 constexpr Scalar min_volume = 1e-10;
259
260 if constexpr (Index == Free) {
261 return std::max(decay<Scalar>(fs.saturation(tracerPhaseIdx)) *
262 decay<Scalar>(fs.invB(tracerPhaseIdx)) *
263 decay<Scalar>(intQuants.porosity()),
264 min_volume);
265 } else {
266 // vaporized oil
267 if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
268 return std::max(decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx)) *
269 decay<Scalar>(fs.invB(FluidSystem::gasPhaseIdx)) *
270 decay<Scalar>(fs.Rv()) *
271 decay<Scalar>(intQuants.porosity()),
272 min_volume);
273 }
274
275 // dissolved gas
276 else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
277 return std::max(decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx)) *
278 decay<Scalar>(fs.invB(FluidSystem::oilPhaseIdx)) *
279 decay<Scalar>(fs.Rs()) *
280 decay<Scalar>(intQuants.porosity()),
281 min_volume);
282 }
283
284 return min_volume;
285 }
286 }
287
288 template<TracerTypeIdx Index>
289 std::pair<TracerEvaluation, bool>
290 computeFlux_(const int tracerPhaseIdx,
291 const ElementContext& elemCtx,
292 const unsigned scvfIdx,
293 const unsigned timeIdx) const
294 {
295 const auto& stencil = elemCtx.stencil(timeIdx);
296 const auto& scvf = stencil.interiorFace(scvfIdx);
297
298 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
299 const unsigned inIdx = extQuants.interiorIndex();
300
301 Scalar v;
302 unsigned upIdx;
303
304 if constexpr (Index == Free) {
305 upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
306 const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
307 const auto& fs = intQuants.fluidState();
308 v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx)) *
309 decay<Scalar>(fs.invB(tracerPhaseIdx));
310 } else {
311 if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
312 upIdx = extQuants.upstreamIndex(FluidSystem::gasPhaseIdx);
313
314 const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
315 const auto& fs = intQuants.fluidState();
316 v = decay<Scalar>(fs.invB(FluidSystem::gasPhaseIdx)) *
317 decay<Scalar>(extQuants.volumeFlux(FluidSystem::gasPhaseIdx)) *
318 decay<Scalar>(fs.Rv());
319 }
320 // dissolved gas
321 else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
322 upIdx = extQuants.upstreamIndex(FluidSystem::oilPhaseIdx);
323
324 const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
325 const auto& fs = intQuants.fluidState();
326 v = decay<Scalar>(fs.invB(FluidSystem::oilPhaseIdx)) *
327 decay<Scalar>(extQuants.volumeFlux(FluidSystem::oilPhaseIdx)) *
328 decay<Scalar>(fs.Rs());
329 }
330 else {
331 upIdx = 0;
332 v = 0.0;
333 }
334 }
335
336 const Scalar A = scvf.area();
337 return inIdx == upIdx
338 ? std::pair{A * v * variable<TracerEvaluation>(1.0, 0), true}
339 : std::pair{A * v, false};
340 }
341
342 template<TracerTypeIdx Index, class TrRe>
343 Scalar storage1_(const TrRe& tr,
344 const unsigned tIdx,
345 const unsigned I,
346 const unsigned I1,
347 const bool cache)
348 {
349 if (cache) {
350 return tr.storageOfTimeIndex1_[tIdx][I][Index];
351 } else {
352 return computeVolume_<Index>(tr.phaseIdx_, I1, 1) *
353 tr.concentration_[tIdx][I1][Index];
354 }
355 }
356
357 template<class TrRe>
359 const ElementContext& elemCtx,
360 const Scalar scvVolume,
361 const Scalar dt,
362 unsigned I,
363 unsigned I1)
364
365 {
366 if (tr.numTracer() == 0) {
367 return;
368 }
369
370 const TracerEvaluation fVol = computeVolume_<Free>(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
371 const TracerEvaluation sVol = computeVolume_<Solution>(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
372 dVol_[Solution][tr.phaseIdx_][I] += sVol.value() * scvVolume - vol1_[1][tr.phaseIdx_][I];
373 dVol_[Free][tr.phaseIdx_][I] += fVol.value() * scvVolume - vol1_[0][tr.phaseIdx_][I];
374 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
375 // Free part
376 const Scalar fStorageOfTimeIndex0 = fVol.value() * tr.concentration_[tIdx][I][Free];
377 const Scalar fLocalStorage = (fStorageOfTimeIndex0 - storage1_<Free>(tr, tIdx, I, I1,
378 elemCtx.enableStorageCache())) * scvVolume / dt;
379 tr.residual_[tIdx][I][Free] += fLocalStorage; // residual + flux
380
381 // Solution part
382 const Scalar sStorageOfTimeIndex0 = sVol.value() * tr.concentration_[tIdx][I][Solution];
383 const Scalar sLocalStorage = (sStorageOfTimeIndex0 - storage1_<Solution>(tr, tIdx, I, I1,
384 elemCtx.enableStorageCache())) * scvVolume / dt;
385 tr.residual_[tIdx][I][Solution] += sLocalStorage; // residual + flux
386 }
387
388 // Derivative matrix
389 (*tr.mat)[I][I][Free][Free] += fVol.derivative(0) * scvVolume/dt;
390 (*tr.mat)[I][I][Solution][Solution] += sVol.derivative(0) * scvVolume/dt;
391 }
392
393 template<class TrRe>
395 const ElementContext& elemCtx,
396 unsigned scvfIdx,
397 unsigned I,
398 unsigned J,
399 const Scalar dt)
400 {
401 if (tr.numTracer() == 0) {
402 return;
403 }
404
405 const auto& [fFlux, isUpF] = computeFlux_<Free>(tr.phaseIdx_, elemCtx, scvfIdx, 0);
406 const auto& [sFlux, isUpS] = computeFlux_<Solution>(tr.phaseIdx_, elemCtx, scvfIdx, 0);
407 dVol_[Solution][tr.phaseIdx_][I] += sFlux.value() * dt;
408 dVol_[Free][tr.phaseIdx_][I] += fFlux.value() * dt;
409 const int fGlobalUpIdx = isUpF ? I : J;
410 const int sGlobalUpIdx = isUpS ? I : J;
411 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
412 // Free and solution fluxes
413 tr.residual_[tIdx][I][Free] += fFlux.value()*tr.concentration_[tIdx][fGlobalUpIdx][Free]; // residual + flux
414 tr.residual_[tIdx][I][Solution] += sFlux.value()*tr.concentration_[tIdx][sGlobalUpIdx][Solution]; // residual + flux
415 }
416
417 // Derivative matrix
418 if (isUpF){
419 (*tr.mat)[J][I][Free][Free] = -fFlux.derivative(0);
420 (*tr.mat)[I][I][Free][Free] += fFlux.derivative(0);
421 }
422 if (isUpS) {
423 (*tr.mat)[J][I][Solution][Solution] = -sFlux.derivative(0);
424 (*tr.mat)[I][I][Solution][Solution] += sFlux.derivative(0);
425 }
426 }
427
428 template<class TrRe, class Well>
430 const Well& well)
431 {
432 if (tr.numTracer() == 0) {
433 return;
434 }
435
436 const auto& eclWell = well.wellEcl();
437
438 // Init. well output to zero
439 auto& tracerRate = this->wellTracerRate_[eclWell.seqIndex()];
440 tracerRate.reserve(tr.numTracer());
441 auto& solTracerRate = this->wellSolTracerRate_[eclWell.seqIndex()];
442 solTracerRate.reserve(tr.numTracer());
443 auto& freeTracerRate = this->wellFreeTracerRate_[eclWell.seqIndex()];
444 freeTracerRate.reserve(tr.numTracer());
445 auto* mswTracerRate = eclWell.isMultiSegment()
446 ? &this->mSwTracerRate_[eclWell.seqIndex()]
447 : nullptr;
448 if (mswTracerRate) {
449 mswTracerRate->reserve(tr.numTracer());
450 }
451 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
452 tracerRate.emplace_back(this->name(tr.idx_[tIdx]), 0.0);
453 freeTracerRate.emplace_back(this->wellfname(tr.idx_[tIdx]), 0.0);
454 solTracerRate.emplace_back(this->wellsname(tr.idx_[tIdx]), 0.0);
455 if (eclWell.isMultiSegment()) {
456 auto& wtr = mswTracerRate->emplace_back(this->name(tr.idx_[tIdx]));
457 wtr.rate.reserve(eclWell.getConnections().size());
458 for (std::size_t i = 0; i < eclWell.getConnections().size(); ++i) {
459 wtr.rate.emplace(eclWell.getConnections().get(i).segment(), 0.0);
460 }
461 }
462 }
463
464 std::vector<Scalar> wtracer(tr.numTracer());
465 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
466 wtracer[tIdx] = this->currentConcentration_(eclWell, this->name(tr.idx_[tIdx]));
467 }
468
469 const Scalar dt = simulator_.timeStepSize();
470 const auto& ws = simulator_.problem().wellModel().wellState().well(well.name());
471 for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
472 const auto I = ws.perf_data.cell_index[i];
473 const Scalar rate = well.volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
474 Scalar rate_s;
475 if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
476 rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
477 }
478 else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
479 rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
480 }
481 else {
482 rate_s = 0.0;
483 }
484
485 const Scalar rate_f = rate - rate_s;
486 if (rate_f > 0) {
487 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
488 const Scalar delta = rate_f * wtracer[tIdx];
489 // Injection of free tracer only
490 tr.residual_[tIdx][I][Free] -= delta;
491
492 // Store _injector_ tracer rate for reporting
493 // (can be done here since WTRACER is constant)
494 tracerRate[tIdx].rate += delta;
495 freeTracerRate[tIdx].rate += delta;
496 if (eclWell.isMultiSegment()) {
497 (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += delta;
498 }
499 }
500 dVol_[Free][tr.phaseIdx_][I] -= rate_f * dt;
501 }
502 else if (rate_f < 0) {
503 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
504 const Scalar delta = rate_f * wtracer[tIdx];
505 // Store _injector_ tracer rate for cross-flowing well connections
506 // (can be done here since WTRACER is constant)
507 tracerRate[tIdx].rate += delta;
508 freeTracerRate[tIdx].rate += delta;
509
510 // Production of free tracer
511 tr.residual_[tIdx][I][Free] -= rate_f * tr.concentration_[tIdx][I][Free];
512 }
513 dVol_[Free][tr.phaseIdx_][I] -= rate_f * dt;
514
515 // Derivative matrix for free tracer producer
516 (*tr.mat)[I][I][Free][Free] -= rate_f * variable<TracerEvaluation>(1.0, 0).derivative(0);
517 }
518 if (rate_s < 0) {
519 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
520 // Production of solution tracer
521 tr.residual_[tIdx][I][Solution] -= rate_s * tr.concentration_[tIdx][I][Solution];
522 }
523 dVol_[Solution][tr.phaseIdx_][I] -= rate_s * dt;
524
525 // Derivative matrix for solution tracer producer
526 (*tr.mat)[I][I][Solution][Solution] -= rate_s * variable<TracerEvaluation>(1.0, 0).derivative(0);
527 }
528 }
529 }
530
531 template<class TrRe>
533 const Scalar dt,
534 unsigned I)
535 {
536 if (tr.numTracer() == 0) {
537 return;
538 }
539
540 // Skip if solution tracers do not exist
541 if (tr.phaseIdx_ == FluidSystem::waterPhaseIdx ||
542 (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && !FluidSystem::enableDissolvedGas()) ||
543 (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && !FluidSystem::enableVaporizedOil()))
544 {
545 return;
546 }
547
548 const Scalar& dsVol = dVol_[Solution][tr.phaseIdx_][I];
549 const Scalar& dfVol = dVol_[Free][tr.phaseIdx_][I];
550
551 // Source term determined by sign of dsVol: if dsVol > 0 then ms -> mf, else mf -> ms
552 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
553 if (dsVol >= 0) {
554 const auto delta = (dfVol / dt) * tr.concentration_[tIdx][I][Free];
555 tr.residual_[tIdx][I][Free] -= delta;
556 tr.residual_[tIdx][I][Solution] += delta;
557 }
558 else {
559 const auto delta = (dsVol / dt) * tr.concentration_[tIdx][I][Solution];
560 tr.residual_[tIdx][I][Free] += delta;
561 tr.residual_[tIdx][I][Solution] -= delta;
562 }
563 }
564
565 // Derivative matrix
566 if (dsVol >= 0) {
567 const auto delta = (dfVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
568 (*tr.mat)[I][I][Free][Free] -= delta;
569 (*tr.mat)[I][I][Solution][Free] += delta;
570 }
571 else {
572 const auto delta = (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
573 (*tr.mat)[I][I][Free][Solution] += delta;
574 (*tr.mat)[I][I][Solution][Solution] -= delta;
575 }
576 }
577
579 {
580 // Note that we formulate the equations in terms of a concentration update
581 // (compared to previous time step) and not absolute concentration.
582 // This implies that current concentration (tr.concentration_[][]) contributes
583 // to the rhs both through storage and flux terms.
584 // Compare also advanceTracerFields(...) below.
585
586 DeferredLogger local_deferredLogger{};
588 {
589 OPM_TIMEBLOCK(tracerAssemble);
590 for (auto& tr : tbatch) {
591 if (tr.numTracer() != 0) {
592 (*tr.mat) = 0.0;
593 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
594 tr.residual_[tIdx] = 0.0;
595 }
596 }
597 }
598
599 this->wellTracerRate_.clear();
600 this->wellFreeTracerRate_.clear();
601 this->wellSolTracerRate_.clear();
602
603 // educated guess for new container size
604 const auto num_msw = this->mSwTracerRate_.size();
605 this->mSwTracerRate_.clear();
606
607 // Well terms
608 const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
609 this->wellTracerRate_.reserve(wellPtrs.size());
610 this->wellFreeTracerRate_.reserve(wellPtrs.size());
611 this->wellSolTracerRate_.reserve(wellPtrs.size());
612 this->mSwTracerRate_.reserve(num_msw);
613 for (const auto& wellPtr : wellPtrs) {
614 for (auto& tr : tbatch) {
615 this->assembleTracerEquationWell(tr, *wellPtr);
616 }
617 }
618
619 // Parallel loop over element chunks
620 #ifdef _OPENMP
621 #pragma omp parallel for
622 #endif
623 for (const auto& chunk : element_chunks_) {
624 ElementContext elemCtx(simulator_);
625 const Scalar dt = elemCtx.simulator().timeStepSize();
626
627 for (const auto& elem : chunk) {
628 elemCtx.updateStencil(elem);
629 const std::size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timeIdx=*/0);
630
631 if (elem.partitionType() != Dune::InteriorEntity) {
632 // Dirichlet boundary conditions for parallel matrix
633 // This is safe as each element has a unique I. So each thread
634 // always writes to different memory locations in the shared arrays.
635 for (const auto& tr : tbatch) {
636 if (tr.numTracer() != 0) {
637 (*tr.mat)[I][I][0][0] = 1.;
638 (*tr.mat)[I][I][1][1] = 1.;
639 }
640 }
641 continue;
642 }
643 elemCtx.updateAllIntensiveQuantities();
644 elemCtx.updateAllExtensiveQuantities();
645
646 const Scalar extrusionFactor =
647 elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
648 Valgrind::CheckDefined(extrusionFactor);
649 assert(isfinite(extrusionFactor));
650 assert(extrusionFactor > 0.0);
651 const Scalar scvVolume =
652 elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume() * extrusionFactor;
653 const std::size_t I1 = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timeIdx=*/1);
654
655 // This is safe as each element has a unique I. So each thread
656 // always writes to different memory locations in the shared arrays.
657 for (auto& tr : tbatch) {
658 if (tr.numTracer() == 0) {
659 continue;
660 }
661 this->assembleTracerEquationVolume(tr, elemCtx, scvVolume, dt, I, I1);
662 }
663
664 const std::size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
665 for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
666 const auto& face = elemCtx.stencil(0).interiorFace(scvfIdx);
667 const unsigned j = face.exteriorIndex();
668 const unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
669 for (auto& tr : tbatch) {
670 if (tr.numTracer() == 0) {
671 continue;
672 }
673 this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J, dt);
674 }
675 }
676
677 // Source terms (mass transfer between free and solution tracer)
678 for (auto& tr : tbatch) {
679 if (tr.numTracer() == 0) {
680 continue;
681 }
682 this->assembleTracerEquationSource(tr, dt, I);
683 }
684 }
685 }
686 }
687 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
688 "assembleTracerEquations() failed: ",
689 true, simulator_.gridView().comm())
690
691 // Communicate overlap using grid Communication
692 for (auto& tr : tbatch) {
693 if (tr.numTracer() == 0) {
694 continue;
695 }
697 simulator_.gridView());
698 simulator_.gridView().communicate(handle, Dune::InteriorBorder_All_Interface,
699 Dune::ForwardCommunication);
700 }
701 }
702
703 template<TracerTypeIdx Index, class TrRe>
704 void updateElem(TrRe& tr,
705 const Scalar scvVolume,
706 const unsigned globalDofIdx)
707 {
708 const Scalar vol1 = computeVolume_<Index>(tr.phaseIdx_, globalDofIdx, 0);
709 vol1_[Index][tr.phaseIdx_][globalDofIdx] = vol1 * scvVolume;
710 dVol_[Index][tr.phaseIdx_][globalDofIdx] = 0.0;
711 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
712 tr.storageOfTimeIndex1_[tIdx][globalDofIdx][Index] =
713 vol1 * tr.concentrationInitial_[tIdx][globalDofIdx][Index];
714 }
715 }
716
718 {
719 for (auto& tr : tbatch) {
720 if (tr.numTracer() != 0) {
721 tr.concentrationInitial_ = tr.concentration_;
722 }
723 }
724
725 // Parallel loop over element chunks
726 #ifdef _OPENMP
727 #pragma omp parallel for
728 #endif
729 for (const auto& chunk : element_chunks_) {
730 ElementContext elemCtx(simulator_);
731
732 for (const auto& elem : chunk) {
733 elemCtx.updatePrimaryStencil(elem);
734 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
735 const Scalar extrusionFactor = elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
736 const Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume() * extrusionFactor;
737 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(0, /*timeIdx=*/0);
738
739 for (auto& tr : tbatch) {
740 if (tr.numTracer() == 0) {
741 continue;
742 }
743 // This is safe as each element has a unique globalDofIdx. So each thread
744 // always writes to different memory locations in the shared arrays.
745 updateElem<Free>(tr, scvVolume, globalDofIdx);
746 updateElem<Solution>(tr, scvVolume, globalDofIdx);
747 }
748 }
749 }
750 }
751
752 template<TracerTypeIdx Index, class TrRe>
753 void copyForOutput(TrRe& tr,
754 const std::vector<TracerVector>& dx,
755 const Scalar S,
756 const unsigned tIdx,
757 const unsigned globalDofIdx,
758 std::vector<TracerVectorSingle>& sc)
759 {
760 constexpr Scalar tol_gas_sat = 1e-6;
761 tr.concentration_[tIdx][globalDofIdx][Index] -= dx[tIdx][globalDofIdx][Index];
762 if (tr.concentration_[tIdx][globalDofIdx][Index] < 0.0 || S < tol_gas_sat) {
763 tr.concentration_[tIdx][globalDofIdx][Index] = 0.0;
764 }
765 sc[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][Index];
766 }
767
768 template<TracerTypeIdx Index, class TrRe>
769 void assignRates(const TrRe& tr,
770 const Well& eclWell,
771 const std::size_t i,
772 const std::size_t I,
773 const Scalar rate,
774 std::vector<WellTracerRate<Scalar>>& tracerRate,
775 std::vector<MSWellTracerRate<Scalar>>* mswTracerRate,
776 std::vector<WellTracerRate<Scalar>>& splitRate)
777 {
778 if (rate < 0) {
779 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
780 // Store _producer_ free tracer rate for reporting
781 const Scalar delta = rate * tr.concentration_[tIdx][I][Index];
782 tracerRate[tIdx].rate += delta;
783 splitRate[tIdx].rate += delta;
784 if (eclWell.isMultiSegment()) {
785 (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += delta;
786 }
787 }
788 }
789 }
790
792 {
794
795 for (auto& tr : tbatch) {
796 if (tr.numTracer() == 0) {
797 continue;
798 }
799
800 // Note that we solve for a concentration update (compared to previous time step)
801 // Confer also assembleTracerEquations_(...) above.
802 std::vector<TracerVector> dx(tr.concentration_);
803 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
804 dx[tIdx] = 0.0;
805 }
806
807 const bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_);
808 if (!converged) {
809 OpmLog::warning("### Tracer model: Linear solver did not converge. ###");
810 }
811
812 OPM_TIMEBLOCK(tracerPost);
813
814 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
815 for (std::size_t globalDofIdx = 0; globalDofIdx < tr.concentration_[tIdx].size(); ++globalDofIdx) {
816 // New concetration. Concentrations that are negative or where free/solution phase is not
817 // present are set to zero
818 const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, 0);
819 const auto& fs = intQuants.fluidState();
820 const Scalar Sf = decay<Scalar>(fs.saturation(tr.phaseIdx_));
821 Scalar Ss = 0.0;
822
823 if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
824 Ss = decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx));
825 }
826 else if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
827 Ss = decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx));
828 }
829
830 copyForOutput<Free>(tr, dx, Sf, tIdx, globalDofIdx, this->freeTracerConcentration_);
831 copyForOutput<Solution>(tr, dx, Ss, tIdx, globalDofIdx, this->solTracerConcentration_);
832 }
833 }
834
835 // Store _producer_ tracer rate for reporting
836 const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
837 for (const auto& wellPtr : wellPtrs) {
838 const auto& eclWell = wellPtr->wellEcl();
839
840 // Injection rates already reported during assembly
841 if (!eclWell.isProducer()) {
842 continue;
843 }
844
845 Scalar rateWellPos = 0.0;
846 Scalar rateWellNeg = 0.0;
847 const std::size_t well_index = simulator_.problem().wellModel().wellState().index(eclWell.name()).value();
848 const auto& ws = simulator_.problem().wellModel().wellState().well(well_index);
849 auto& tracerRate = this->wellTracerRate_[eclWell.seqIndex()];
850 auto& freeTracerRate = this->wellFreeTracerRate_[eclWell.seqIndex()];
851 auto& solTracerRate = this->wellSolTracerRate_[eclWell.seqIndex()];
852 auto* mswTracerRate = eclWell.isMultiSegment() ? &this->mSwTracerRate_[eclWell.seqIndex()] : nullptr;
853
854 for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
855 const auto I = ws.perf_data.cell_index[i];
856 const Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
857
858 Scalar rate_s;
859 if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
860 rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
861 }
862 else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
863 rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
864 }
865 else {
866 rate_s = 0.0;
867 }
868
869 const Scalar rate_f = rate - rate_s;
870 assignRates<Free>(tr, eclWell, i, I, rate_f,
871 tracerRate, mswTracerRate, freeTracerRate);
872 assignRates<Solution>(tr, eclWell, i, I, rate_s,
873 tracerRate, mswTracerRate, solTracerRate);
874
875 if (rate < 0) {
876 rateWellNeg += rate;
877 }
878 else {
879 rateWellPos += rate;
880 }
881 }
882
883 // TODO: Some inconsistencies here that perhaps should be clarified.
884 // The "offical" rate as reported below is occasionally significant
885 // different from the sum over connections (as calculated above). Only observed
886 // for small values, neglible for the rate itself, but matters when used to
887 // calculate tracer concentrations.
888 const Scalar official_well_rate_total =
889 simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_];
890
891 const Scalar rateWellTotal = official_well_rate_total;
892
893 if (rateWellTotal > rateWellNeg) { // Cross flow
894 constexpr Scalar bucketPrDay = 10.0 / (1000. * 3600. * 24.); // ... keeps (some) trouble away
895 const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal / rateWellNeg : 0.0;
896 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
897 tracerRate[tIdx].rate *= factor;
898 }
899 }
900 }
901 }
902 }
903
904 Simulator& simulator_;
905
906 // This struct collects tracers of the same type (i.e, transported in same phase).
907 // The idea being that, under the assumption of linearity, tracers of same type can
908 // be solved in concert, having a common system matrix but separate right-hand-sides.
909
910 // Since oil or gas tracers appears in dual compositions when VAPOIL respectively DISGAS
911 // is active, the template argument is intended to support future extension to these
912 // scenarios by supplying an extended vector type.
913
914 template <typename TV>
916 {
917 std::vector<int> idx_;
918 const int phaseIdx_;
919 std::vector<TV> concentrationInitial_;
920 std::vector<TV> concentration_;
921 std::vector<TV> storageOfTimeIndex1_;
922 std::vector<TV> residual_;
923 std::unique_ptr<TracerMatrix> mat;
924
925 bool operator==(const TracerBatch& rhs) const
926 {
927 return this->concentrationInitial_ == rhs.concentrationInitial_ &&
928 this->concentration_ == rhs.concentration_;
929 }
930
932 {
933 TracerBatch<TV> result(4);
934 result.idx_ = {1,2,3};
935 result.concentrationInitial_ = {5.0, 6.0};
936 result.concentration_ = {7.0, 8.0};
937 result.storageOfTimeIndex1_ = {9.0, 10.0, 11.0};
938 result.residual_ = {12.0, 13.0};
939
940 return result;
941 }
942
943 template<class Serializer>
944 void serializeOp(Serializer& serializer)
945 {
946 serializer(concentrationInitial_);
947 serializer(concentration_);
948 }
949
950 TracerBatch(int phaseIdx = 0) : phaseIdx_(phaseIdx) {}
951
952 int numTracer() const
953 { return idx_.size(); }
954
955 void addTracer(const int idx, const TV& concentration)
956 {
957 const int numGridDof = concentration.size();
958 idx_.emplace_back(idx);
959 concentrationInitial_.emplace_back(concentration);
960 concentration_.emplace_back(concentration);
961 residual_.emplace_back(numGridDof);
962 storageOfTimeIndex1_.emplace_back(numGridDof);
963 }
964 };
965
966 std::array<TracerBatch<TracerVector>,numPhases> tbatch;
970 std::array<std::array<std::vector<Scalar>,numPhases>,2> vol1_;
971 std::array<std::array<std::vector<Scalar>,numPhases>,2> dVol_;
972 ElementChunks<GridView, Dune::Partitions::All> element_chunks_;
973};
974
975} // namespace Opm
976
977#endif // OPM_TRACER_MODEL_HPP
#define OPM_END_PARALLEL_TRY_CATCH_LOG(obptc_logger, obptc_prefix, obptc_output, comm)
Catch exception, log, and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:202
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:158
A datahandle sending data located in multiple vectors.
Definition: DeferredLogger.hpp:57
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
static unsigned maxThreads()
Return the maximum number of threads of the current process.
Definition: threadmanager.hpp:66
A class which handles tracers as specified in by ECL.
Definition: TracerModel.hpp:75
void updateStorageCache()
Definition: TracerModel.hpp:717
std::array< std::array< std::vector< Scalar >, numPhases >, 2 > vol1_
Definition: TracerModel.hpp:970
TracerBatch< TracerVector > & oil_
Definition: TracerModel.hpp:968
void advanceTracerFields()
Definition: TracerModel.hpp:791
void init(bool rst)
Definition: TracerModel.hpp:138
void assembleTracerEquationSource(TrRe &tr, const Scalar dt, unsigned I)
Definition: TracerModel.hpp:532
TracerModel(Simulator &simulator)
Definition: TracerModel.hpp:105
void assignRates(const TrRe &tr, const Well &eclWell, const std::size_t i, const std::size_t I, const Scalar rate, std::vector< WellTracerRate< Scalar > > &tracerRate, std::vector< MSWellTracerRate< Scalar > > *mswTracerRate, std::vector< WellTracerRate< Scalar > > &splitRate)
Definition: TracerModel.hpp:769
void copyForOutput(TrRe &tr, const std::vector< TracerVector > &dx, const Scalar S, const unsigned tIdx, const unsigned globalDofIdx, std::vector< TracerVectorSingle > &sc)
Definition: TracerModel.hpp:753
void assembleTracerEquationFlux(TrRe &tr, const ElementContext &elemCtx, unsigned scvfIdx, unsigned I, unsigned J, const Scalar dt)
Definition: TracerModel.hpp:394
void assembleTracerEquationWell(TrRe &tr, const Well &well)
Definition: TracerModel.hpp:429
void updateElem(TrRe &tr, const Scalar scvVolume, const unsigned globalDofIdx)
Definition: TracerModel.hpp:704
Scalar computeVolume_(const int tracerPhaseIdx, const unsigned globalDofIdx, const unsigned timeIdx) const
Definition: TracerModel.hpp:252
typename BaseType::TracerTypeIdx TracerTypeIdx
Definition: TracerModel.hpp:246
void assembleTracerEquationVolume(TrRe &tr, const ElementContext &elemCtx, const Scalar scvVolume, const Scalar dt, unsigned I, unsigned I1)
Definition: TracerModel.hpp:358
Scalar storage1_(const TrRe &tr, const unsigned tIdx, const unsigned I, const unsigned I1, const bool cache)
Definition: TracerModel.hpp:343
void serialize(Restarter &)
This method writes the complete state of all tracer to the hard disk.
Definition: TracerModel.hpp:225
TracerBatch< TracerVector > & gas_
Definition: TracerModel.hpp:969
void assembleTracerEquations_()
Definition: TracerModel.hpp:578
std::array< TracerBatch< TracerVector >, numPhases > tbatch
Definition: TracerModel.hpp:966
std::pair< TracerEvaluation, bool > computeFlux_(const int tracerPhaseIdx, const ElementContext &elemCtx, const unsigned scvfIdx, const unsigned timeIdx) const
Definition: TracerModel.hpp:290
void beginTimeStep()
Definition: TracerModel.hpp:197
void serializeOp(Serializer &serializer)
Definition: TracerModel.hpp:239
void prepareTracerBatches()
Definition: TracerModel.hpp:144
std::array< std::array< std::vector< Scalar >, numPhases >, 2 > dVol_
Definition: TracerModel.hpp:971
ElementChunks< GridView, Dune::Partitions::All > element_chunks_
Definition: TracerModel.hpp:972
TracerBatch< TracerVector > & wat_
Definition: TracerModel.hpp:967
void endTimeStep()
Informs the tracer model that a time step has just been finished.
Definition: TracerModel.hpp:210
Simulator & simulator_
Definition: TracerModel.hpp:904
void deserialize(Restarter &)
This method restores the complete state of the tracer from disk.
Definition: TracerModel.hpp:235
A data handle sending multiple data store in vectors attached to cells.
Definition: VectorVectorDataHandle.hpp:50
int Index
The type of an index of a degree of freedom.
Definition: overlaptypes.hh:44
Definition: blackoilmodel.hh:79
Definition: blackoilboundaryratevector.hh:39
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:233
The Opm property system, traits with inheritance.
Definition: WellTracerRate.hpp:60
Definition: TracerModel.hpp:55
a tag to mark properties as undefined
Definition: propertysystem.hh:38
Definition: TracerModel.hpp:916
const int phaseIdx_
Definition: TracerModel.hpp:918
void addTracer(const int idx, const TV &concentration)
Definition: TracerModel.hpp:955
bool operator==(const TracerBatch &rhs) const
Definition: TracerModel.hpp:925
static TracerBatch serializationTestObject()
Definition: TracerModel.hpp:931
std::vector< TV > concentrationInitial_
Definition: TracerModel.hpp:919
void serializeOp(Serializer &serializer)
Definition: TracerModel.hpp:944
TracerBatch(int phaseIdx=0)
Definition: TracerModel.hpp:950
int numTracer() const
Definition: TracerModel.hpp:952
std::vector< int > idx_
Definition: TracerModel.hpp:917
std::unique_ptr< TracerMatrix > mat
Definition: TracerModel.hpp:923
std::vector< TV > storageOfTimeIndex1_
Definition: TracerModel.hpp:921
std::vector< TV > residual_
Definition: TracerModel.hpp:922
std::vector< TV > concentration_
Definition: TracerModel.hpp:920
Definition: WellTracerRate.hpp:33