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 simulator_.problem().wellModel().summaryState());
468 }
469
470 const Scalar dt = simulator_.timeStepSize();
471 const auto& ws = simulator_.problem().wellModel().wellState().well(well.name());
472 for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
473 const auto I = ws.perf_data.cell_index[i];
474 const Scalar rate = well.volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
475 Scalar rate_s;
476 if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
477 rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
478 }
479 else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
480 rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
481 }
482 else {
483 rate_s = 0.0;
484 }
485
486 const Scalar rate_f = rate - rate_s;
487 if (rate_f > 0) {
488 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
489 const Scalar delta = rate_f * wtracer[tIdx];
490 // Injection of free tracer only
491 tr.residual_[tIdx][I][Free] -= delta;
492
493 // Store _injector_ tracer rate for reporting
494 // (can be done here since WTRACER is constant)
495 tracerRate[tIdx].rate += delta;
496 freeTracerRate[tIdx].rate += delta;
497 if (eclWell.isMultiSegment()) {
498 (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += delta;
499 }
500 }
501 dVol_[Free][tr.phaseIdx_][I] -= rate_f * dt;
502 }
503 else if (rate_f < 0) {
504 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
505 const Scalar delta = rate_f * wtracer[tIdx];
506 // Store _injector_ tracer rate for cross-flowing well connections
507 // (can be done here since WTRACER is constant)
508 tracerRate[tIdx].rate += delta;
509 freeTracerRate[tIdx].rate += delta;
510
511 // Production of free tracer
512 tr.residual_[tIdx][I][Free] -= rate_f * tr.concentration_[tIdx][I][Free];
513 }
514 dVol_[Free][tr.phaseIdx_][I] -= rate_f * dt;
515
516 // Derivative matrix for free tracer producer
517 (*tr.mat)[I][I][Free][Free] -= rate_f * variable<TracerEvaluation>(1.0, 0).derivative(0);
518 }
519 if (rate_s < 0) {
520 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
521 // Production of solution tracer
522 tr.residual_[tIdx][I][Solution] -= rate_s * tr.concentration_[tIdx][I][Solution];
523 }
524 dVol_[Solution][tr.phaseIdx_][I] -= rate_s * dt;
525
526 // Derivative matrix for solution tracer producer
527 (*tr.mat)[I][I][Solution][Solution] -= rate_s * variable<TracerEvaluation>(1.0, 0).derivative(0);
528 }
529 }
530 }
531
532 template<class TrRe>
534 const Scalar dt,
535 unsigned I)
536 {
537 if (tr.numTracer() == 0) {
538 return;
539 }
540
541 // Skip if solution tracers do not exist
542 if (tr.phaseIdx_ == FluidSystem::waterPhaseIdx ||
543 (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && !FluidSystem::enableDissolvedGas()) ||
544 (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && !FluidSystem::enableVaporizedOil()))
545 {
546 return;
547 }
548
549 const Scalar& dsVol = dVol_[Solution][tr.phaseIdx_][I];
550 const Scalar& dfVol = dVol_[Free][tr.phaseIdx_][I];
551
552 // Source term determined by sign of dsVol: if dsVol > 0 then ms -> mf, else mf -> ms
553 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
554 if (dsVol >= 0) {
555 const auto delta = (dfVol / dt) * tr.concentration_[tIdx][I][Free];
556 tr.residual_[tIdx][I][Free] -= delta;
557 tr.residual_[tIdx][I][Solution] += delta;
558 }
559 else {
560 const auto delta = (dsVol / dt) * tr.concentration_[tIdx][I][Solution];
561 tr.residual_[tIdx][I][Free] += delta;
562 tr.residual_[tIdx][I][Solution] -= delta;
563 }
564 }
565
566 // Derivative matrix
567 if (dsVol >= 0) {
568 const auto delta = (dfVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
569 (*tr.mat)[I][I][Free][Free] -= delta;
570 (*tr.mat)[I][I][Solution][Free] += delta;
571 }
572 else {
573 const auto delta = (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
574 (*tr.mat)[I][I][Free][Solution] += delta;
575 (*tr.mat)[I][I][Solution][Solution] -= delta;
576 }
577 }
578
580 {
581 // Note that we formulate the equations in terms of a concentration update
582 // (compared to previous time step) and not absolute concentration.
583 // This implies that current concentration (tr.concentration_[][]) contributes
584 // to the rhs both through storage and flux terms.
585 // Compare also advanceTracerFields(...) below.
586
587 DeferredLogger local_deferredLogger{};
589 {
590 OPM_TIMEBLOCK(tracerAssemble);
591 for (auto& tr : tbatch) {
592 if (tr.numTracer() != 0) {
593 (*tr.mat) = 0.0;
594 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
595 tr.residual_[tIdx] = 0.0;
596 }
597 }
598 }
599
600 this->wellTracerRate_.clear();
601 this->wellFreeTracerRate_.clear();
602 this->wellSolTracerRate_.clear();
603
604 // educated guess for new container size
605 const auto num_msw = this->mSwTracerRate_.size();
606 this->mSwTracerRate_.clear();
607
608 // Well terms
609 const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
610 this->wellTracerRate_.reserve(wellPtrs.size());
611 this->wellFreeTracerRate_.reserve(wellPtrs.size());
612 this->wellSolTracerRate_.reserve(wellPtrs.size());
613 this->mSwTracerRate_.reserve(num_msw);
614 for (const auto& wellPtr : wellPtrs) {
615 for (auto& tr : tbatch) {
616 this->assembleTracerEquationWell(tr, *wellPtr);
617 }
618 }
619
620 // Parallel loop over element chunks
621 #ifdef _OPENMP
622 #pragma omp parallel for
623 #endif
624 for (const auto& chunk : element_chunks_) {
625 ElementContext elemCtx(simulator_);
626 const Scalar dt = elemCtx.simulator().timeStepSize();
627
628 for (const auto& elem : chunk) {
629 elemCtx.updateStencil(elem);
630 const std::size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timeIdx=*/0);
631
632 if (elem.partitionType() != Dune::InteriorEntity) {
633 // Dirichlet boundary conditions for parallel matrix
634 // This is safe as each element has a unique I. So each thread
635 // always writes to different memory locations in the shared arrays.
636 for (const auto& tr : tbatch) {
637 if (tr.numTracer() != 0) {
638 (*tr.mat)[I][I][0][0] = 1.;
639 (*tr.mat)[I][I][1][1] = 1.;
640 }
641 }
642 continue;
643 }
644 elemCtx.updateAllIntensiveQuantities();
645 elemCtx.updateAllExtensiveQuantities();
646
647 const Scalar extrusionFactor =
648 elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
649 Valgrind::CheckDefined(extrusionFactor);
650 assert(isfinite(extrusionFactor));
651 assert(extrusionFactor > 0.0);
652 const Scalar scvVolume =
653 elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume() * extrusionFactor;
654 const std::size_t I1 = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timeIdx=*/1);
655
656 // This is safe as each element has a unique I. So each thread
657 // always writes to different memory locations in the shared arrays.
658 for (auto& tr : tbatch) {
659 if (tr.numTracer() == 0) {
660 continue;
661 }
662 this->assembleTracerEquationVolume(tr, elemCtx, scvVolume, dt, I, I1);
663 }
664
665 const std::size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
666 for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
667 const auto& face = elemCtx.stencil(0).interiorFace(scvfIdx);
668 const unsigned j = face.exteriorIndex();
669 const unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
670 for (auto& tr : tbatch) {
671 if (tr.numTracer() == 0) {
672 continue;
673 }
674 this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J, dt);
675 }
676 }
677
678 // Source terms (mass transfer between free and solution tracer)
679 for (auto& tr : tbatch) {
680 if (tr.numTracer() == 0) {
681 continue;
682 }
683 this->assembleTracerEquationSource(tr, dt, I);
684 }
685 }
686 }
687 }
688 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
689 "assembleTracerEquations() failed: ",
690 true, simulator_.gridView().comm())
691
692 // Communicate overlap using grid Communication
693 for (auto& tr : tbatch) {
694 if (tr.numTracer() == 0) {
695 continue;
696 }
698 simulator_.gridView());
699 simulator_.gridView().communicate(handle, Dune::InteriorBorder_All_Interface,
700 Dune::ForwardCommunication);
701 }
702 }
703
704 template<TracerTypeIdx Index, class TrRe>
705 void updateElem(TrRe& tr,
706 const Scalar scvVolume,
707 const unsigned globalDofIdx)
708 {
709 const Scalar vol1 = computeVolume_<Index>(tr.phaseIdx_, globalDofIdx, 0);
710 vol1_[Index][tr.phaseIdx_][globalDofIdx] = vol1 * scvVolume;
711 dVol_[Index][tr.phaseIdx_][globalDofIdx] = 0.0;
712 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
713 tr.storageOfTimeIndex1_[tIdx][globalDofIdx][Index] =
714 vol1 * tr.concentrationInitial_[tIdx][globalDofIdx][Index];
715 }
716 }
717
719 {
720 for (auto& tr : tbatch) {
721 if (tr.numTracer() != 0) {
722 tr.concentrationInitial_ = tr.concentration_;
723 }
724 }
725
726 // Parallel loop over element chunks
727 #ifdef _OPENMP
728 #pragma omp parallel for
729 #endif
730 for (const auto& chunk : element_chunks_) {
731 ElementContext elemCtx(simulator_);
732
733 for (const auto& elem : chunk) {
734 elemCtx.updatePrimaryStencil(elem);
735 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
736 const Scalar extrusionFactor = elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
737 const Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume() * extrusionFactor;
738 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(0, /*timeIdx=*/0);
739
740 for (auto& tr : tbatch) {
741 if (tr.numTracer() == 0) {
742 continue;
743 }
744 // This is safe as each element has a unique globalDofIdx. So each thread
745 // always writes to different memory locations in the shared arrays.
746 updateElem<Free>(tr, scvVolume, globalDofIdx);
747 updateElem<Solution>(tr, scvVolume, globalDofIdx);
748 }
749 }
750 }
751 }
752
753 template<TracerTypeIdx Index, class TrRe>
754 void copyForOutput(TrRe& tr,
755 const std::vector<TracerVector>& dx,
756 const Scalar S,
757 const unsigned tIdx,
758 const unsigned globalDofIdx,
759 std::vector<TracerVectorSingle>& sc)
760 {
761 constexpr Scalar tol_gas_sat = 1e-6;
762 tr.concentration_[tIdx][globalDofIdx][Index] -= dx[tIdx][globalDofIdx][Index];
763 if (tr.concentration_[tIdx][globalDofIdx][Index] < 0.0 || S < tol_gas_sat) {
764 tr.concentration_[tIdx][globalDofIdx][Index] = 0.0;
765 }
766 sc[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][Index];
767 }
768
769 template<TracerTypeIdx Index, class TrRe>
770 void assignRates(const TrRe& tr,
771 const Well& eclWell,
772 const std::size_t i,
773 const std::size_t I,
774 const Scalar rate,
775 std::vector<WellTracerRate<Scalar>>& tracerRate,
776 std::vector<MSWellTracerRate<Scalar>>* mswTracerRate,
777 std::vector<WellTracerRate<Scalar>>& splitRate)
778 {
779 if (rate < 0) {
780 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
781 // Store _producer_ free tracer rate for reporting
782 const Scalar delta = rate * tr.concentration_[tIdx][I][Index];
783 tracerRate[tIdx].rate += delta;
784 splitRate[tIdx].rate += delta;
785 if (eclWell.isMultiSegment()) {
786 (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += delta;
787 }
788 }
789 }
790 }
791
793 {
795
796 for (auto& tr : tbatch) {
797 if (tr.numTracer() == 0) {
798 continue;
799 }
800
801 // Note that we solve for a concentration update (compared to previous time step)
802 // Confer also assembleTracerEquations_(...) above.
803 std::vector<TracerVector> dx(tr.concentration_);
804 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
805 dx[tIdx] = 0.0;
806 }
807
808 const bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_);
809 if (!converged) {
810 OpmLog::warning("### Tracer model: Linear solver did not converge. ###");
811 }
812
813 OPM_TIMEBLOCK(tracerPost);
814
815 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
816 for (std::size_t globalDofIdx = 0; globalDofIdx < tr.concentration_[tIdx].size(); ++globalDofIdx) {
817 // New concetration. Concentrations that are negative or where free/solution phase is not
818 // present are set to zero
819 const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, 0);
820 const auto& fs = intQuants.fluidState();
821 const Scalar Sf = decay<Scalar>(fs.saturation(tr.phaseIdx_));
822 Scalar Ss = 0.0;
823
824 if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
825 Ss = decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx));
826 }
827 else if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
828 Ss = decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx));
829 }
830
831 copyForOutput<Free>(tr, dx, Sf, tIdx, globalDofIdx, this->freeTracerConcentration_);
832 copyForOutput<Solution>(tr, dx, Ss, tIdx, globalDofIdx, this->solTracerConcentration_);
833 }
834 }
835
836 // Store _producer_ tracer rate for reporting
837 const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
838 for (const auto& wellPtr : wellPtrs) {
839 const auto& eclWell = wellPtr->wellEcl();
840
841 // Injection rates already reported during assembly
842 if (!eclWell.isProducer()) {
843 continue;
844 }
845
846 Scalar rateWellPos = 0.0;
847 Scalar rateWellNeg = 0.0;
848 const std::size_t well_index = simulator_.problem().wellModel().wellState().index(eclWell.name()).value();
849 const auto& ws = simulator_.problem().wellModel().wellState().well(well_index);
850 auto& tracerRate = this->wellTracerRate_[eclWell.seqIndex()];
851 auto& freeTracerRate = this->wellFreeTracerRate_[eclWell.seqIndex()];
852 auto& solTracerRate = this->wellSolTracerRate_[eclWell.seqIndex()];
853 auto* mswTracerRate = eclWell.isMultiSegment() ? &this->mSwTracerRate_[eclWell.seqIndex()] : nullptr;
854
855 for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
856 const auto I = ws.perf_data.cell_index[i];
857 const Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
858
859 Scalar rate_s;
860 if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
861 rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
862 }
863 else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
864 rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
865 }
866 else {
867 rate_s = 0.0;
868 }
869
870 const Scalar rate_f = rate - rate_s;
871 assignRates<Free>(tr, eclWell, i, I, rate_f,
872 tracerRate, mswTracerRate, freeTracerRate);
873 assignRates<Solution>(tr, eclWell, i, I, rate_s,
874 tracerRate, mswTracerRate, solTracerRate);
875
876 if (rate < 0) {
877 rateWellNeg += rate;
878 }
879 else {
880 rateWellPos += rate;
881 }
882 }
883
884 // TODO: Some inconsistencies here that perhaps should be clarified.
885 // The "offical" rate as reported below is occasionally significant
886 // different from the sum over connections (as calculated above). Only observed
887 // for small values, neglible for the rate itself, but matters when used to
888 // calculate tracer concentrations.
889 const Scalar official_well_rate_total =
890 simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_];
891
892 const Scalar rateWellTotal = official_well_rate_total;
893
894 if (rateWellTotal > rateWellNeg) { // Cross flow
895 constexpr Scalar bucketPrDay = 10.0 / (1000. * 3600. * 24.); // ... keeps (some) trouble away
896 const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal / rateWellNeg : 0.0;
897 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
898 tracerRate[tIdx].rate *= factor;
899 }
900 }
901 }
902 }
903 }
904
905 Simulator& simulator_;
906
907 // This struct collects tracers of the same type (i.e, transported in same phase).
908 // The idea being that, under the assumption of linearity, tracers of same type can
909 // be solved in concert, having a common system matrix but separate right-hand-sides.
910
911 // Since oil or gas tracers appears in dual compositions when VAPOIL respectively DISGAS
912 // is active, the template argument is intended to support future extension to these
913 // scenarios by supplying an extended vector type.
914
915 template <typename TV>
917 {
918 std::vector<int> idx_;
919 const int phaseIdx_;
920 std::vector<TV> concentrationInitial_;
921 std::vector<TV> concentration_;
922 std::vector<TV> storageOfTimeIndex1_;
923 std::vector<TV> residual_;
924 std::unique_ptr<TracerMatrix> mat;
925
926 bool operator==(const TracerBatch& rhs) const
927 {
928 return this->concentrationInitial_ == rhs.concentrationInitial_ &&
929 this->concentration_ == rhs.concentration_;
930 }
931
933 {
934 TracerBatch<TV> result(4);
935 result.idx_ = {1,2,3};
936 result.concentrationInitial_ = {5.0, 6.0};
937 result.concentration_ = {7.0, 8.0};
938 result.storageOfTimeIndex1_ = {9.0, 10.0, 11.0};
939 result.residual_ = {12.0, 13.0};
940
941 return result;
942 }
943
944 template<class Serializer>
945 void serializeOp(Serializer& serializer)
946 {
947 serializer(concentrationInitial_);
948 serializer(concentration_);
949 }
950
951 TracerBatch(int phaseIdx = 0) : phaseIdx_(phaseIdx) {}
952
953 int numTracer() const
954 { return idx_.size(); }
955
956 void addTracer(const int idx, const TV& concentration)
957 {
958 const int numGridDof = concentration.size();
959 idx_.emplace_back(idx);
960 concentrationInitial_.emplace_back(concentration);
961 concentration_.emplace_back(concentration);
962 residual_.emplace_back(numGridDof);
963 storageOfTimeIndex1_.emplace_back(numGridDof);
964 }
965 };
966
967 std::array<TracerBatch<TracerVector>,numPhases> tbatch;
971 std::array<std::array<std::vector<Scalar>,numPhases>,2> vol1_;
972 std::array<std::array<std::vector<Scalar>,numPhases>,2> dVol_;
973 ElementChunks<GridView, Dune::Partitions::All> element_chunks_;
974};
975
976} // namespace Opm
977
978#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:227
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:718
std::array< std::array< std::vector< Scalar >, numPhases >, 2 > vol1_
Definition: TracerModel.hpp:971
TracerBatch< TracerVector > & oil_
Definition: TracerModel.hpp:969
void advanceTracerFields()
Definition: TracerModel.hpp:792
void init(bool rst)
Definition: TracerModel.hpp:138
void assembleTracerEquationSource(TrRe &tr, const Scalar dt, unsigned I)
Definition: TracerModel.hpp:533
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:770
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:754
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:705
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:970
void assembleTracerEquations_()
Definition: TracerModel.hpp:579
std::array< TracerBatch< TracerVector >, numPhases > tbatch
Definition: TracerModel.hpp:967
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:972
ElementChunks< GridView, Dune::Partitions::All > element_chunks_
Definition: TracerModel.hpp:973
TracerBatch< TracerVector > & wat_
Definition: TracerModel.hpp:968
void endTimeStep()
Informs the tracer model that a time step has just been finished.
Definition: TracerModel.hpp:210
Simulator & simulator_
Definition: TracerModel.hpp:905
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: blackoilbioeffectsmodules.hh:43
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:917
const int phaseIdx_
Definition: TracerModel.hpp:919
void addTracer(const int idx, const TV &concentration)
Definition: TracerModel.hpp:956
bool operator==(const TracerBatch &rhs) const
Definition: TracerModel.hpp:926
static TracerBatch serializationTestObject()
Definition: TracerModel.hpp:932
std::vector< TV > concentrationInitial_
Definition: TracerModel.hpp:920
void serializeOp(Serializer &serializer)
Definition: TracerModel.hpp:945
TracerBatch(int phaseIdx=0)
Definition: TracerModel.hpp:951
int numTracer() const
Definition: TracerModel.hpp:953
std::vector< int > idx_
Definition: TracerModel.hpp:918
std::unique_ptr< TracerMatrix > mat
Definition: TracerModel.hpp:924
std::vector< TV > storageOfTimeIndex1_
Definition: TracerModel.hpp:922
std::vector< TV > residual_
Definition: TracerModel.hpp:923
std::vector< TV > concentration_
Definition: TracerModel.hpp:921
Definition: WellTracerRate.hpp:33