My Project
Loading...
Searching...
No Matches
InitStateEquil_impl.hpp
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 3 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*/
23#ifndef OPM_INIT_STATE_EQUIL_IMPL_HPP
24#define OPM_INIT_STATE_EQUIL_IMPL_HPP
25
26#include <dune/grid/common/mcmgmapper.hh>
27
28#include <opm/common/OpmLog/OpmLog.hpp>
29
30#include <opm/grid/utility/RegionMapping.hpp>
31
32#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
33#include <opm/input/eclipse/EclipseState/Tables/RsvdTable.hpp>
34#include <opm/input/eclipse/EclipseState/Tables/RvvdTable.hpp>
35#include <opm/input/eclipse/EclipseState/Tables/RvwvdTable.hpp>
36#include <opm/input/eclipse/EclipseState/Tables/PbvdTable.hpp>
37#include <opm/input/eclipse/EclipseState/Tables/PdvdTable.hpp>
38#include <opm/input/eclipse/EclipseState/Tables/SaltvdTable.hpp>
39#include <opm/input/eclipse/EclipseState/Tables/RtempvdTable.hpp>
40
41#include <opm/input/eclipse/EclipseState/Tables/SaltpvdTable.hpp>
42
43#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
44#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
45
48#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
49
50#include <fmt/format.h>
51
52#include <algorithm>
53#include <cassert>
54#include <cstddef>
55#include <limits>
56#include <stdexcept>
57
58namespace Opm {
59namespace EQUIL {
60
61namespace Details {
62
63template <typename CellRange, typename Comm>
64void verticalExtent(const CellRange& cells,
65 const std::vector<std::pair<double, double>>& cellZMinMax,
66 const Comm& comm,
67 std::array<double,2>& span)
68{
69 span[0] = std::numeric_limits<double>::max();
70 span[1] = std::numeric_limits<double>::lowest();
71
72 // Define vertical span as
73 //
74 // [minimum(node depth(cells)), maximum(node depth(cells))]
75 //
76 // Note: The implementation of 'RK4IVP<>' implicitly
77 // imposes the requirement that cell centroids are all
78 // within this vertical span. That requirement is not
79 // checked.
80 for (const auto& cell : cells) {
81 if (cellZMinMax[cell].first < span[0]) { span[0] = cellZMinMax[cell].first; }
82 if (cellZMinMax[cell].second > span[1]) { span[1] = cellZMinMax[cell].second; }
83 }
84 span[0] = comm.min(span[0]);
85 span[1] = comm.max(span[1]);
86}
87
88void subdivisionCentrePoints(const double left,
89 const double right,
90 const int numIntervals,
91 std::vector<std::pair<double, double>>& subdiv)
92{
93 const auto h = (right - left) / numIntervals;
94
95 auto end = left;
96 for (auto i = 0*numIntervals; i < numIntervals; ++i) {
97 const auto start = end;
98 end = left + (i + 1)*h;
99
100 subdiv.emplace_back((start + end) / 2, h);
101 }
102}
103
104template <typename CellID>
105std::vector<std::pair<double, double>>
106horizontalSubdivision(const CellID cell,
107 const std::pair<double, double> topbot,
108 const int numIntervals)
109{
110 auto subdiv = std::vector<std::pair<double, double>>{};
111 subdiv.reserve(2 * numIntervals);
112
113 if (topbot.first > topbot.second) {
114 throw std::out_of_range {
115 "Negative thickness (inverted top/bottom faces) in cell "
116 + std::to_string(cell)
117 };
118 }
119
120 subdivisionCentrePoints(topbot.first, topbot.second,
121 2*numIntervals, subdiv);
122
123 return subdiv;
124}
125
126template <class Element>
127double cellCenterDepth(const Element& element)
128{
129 typedef typename Element::Geometry Geometry;
130 static constexpr int zCoord = Element::dimension - 1;
131 double zz = 0.0;
132
133 const Geometry& geometry = element.geometry();
134 const int corners = geometry.corners();
135 for (int i=0; i < corners; ++i)
136 zz += geometry.corner(i)[zCoord];
137
138 return zz/corners;
139}
140
141template <class Element>
142std::pair<double,double> cellZSpan(const Element& element)
143{
144 typedef typename Element::Geometry Geometry;
145 static constexpr int zCoord = Element::dimension - 1;
146 double bot = 0.0;
147 double top = 0.0;
148
149 const Geometry& geometry = element.geometry();
150 const int corners = geometry.corners();
151 assert(corners == 8);
152 for (int i=0; i < 4; ++i)
153 bot += geometry.corner(i)[zCoord];
154 for (int i=4; i < corners; ++i)
155 top += geometry.corner(i)[zCoord];
156
157 return std::make_pair(bot/4, top/4);
158}
159
160template <class Element>
161std::pair<double,double> cellZMinMax(const Element& element)
162{
163 typedef typename Element::Geometry Geometry;
164 static constexpr int zCoord = Element::dimension - 1;
165 const Geometry& geometry = element.geometry();
166 const int corners = geometry.corners();
167 assert(corners == 8);
168 auto min = std::numeric_limits<double>::max();
169 auto max = std::numeric_limits<double>::lowest();
170
171
172 for (int i=0; i < corners; ++i) {
173 min = std::min(min, geometry.corner(i)[zCoord]);
174 max = std::max(max, geometry.corner(i)[zCoord]);
175 }
176 return std::make_pair(min, max);
177}
178
179template<class RHS>
180RK4IVP<RHS>::RK4IVP(const RHS& f,
181 const std::array<double,2>& span,
182 const double y0,
183 const int N)
184 : N_(N)
185 , span_(span)
186{
187 const double h = stepsize();
188 const double h2 = h / 2;
189 const double h6 = h / 6;
190
191 y_.reserve(N + 1);
192 f_.reserve(N + 1);
193
194 y_.push_back(y0);
195 f_.push_back(f(span_[0], y0));
196
197 for (int i = 0; i < N; ++i) {
198 const double x = span_[0] + i*h;
199 const double y = y_.back();
200
201 const double k1 = f_[i];
202 const double k2 = f(x + h2, y + h2*k1);
203 const double k3 = f(x + h2, y + h2*k2);
204 const double k4 = f(x + h, y + h*k3);
205
206 y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4));
207 f_.push_back(f(x + h, y_.back()));
208 }
209
210 assert (y_.size() == std::vector<double>::size_type(N + 1));
211}
212
213template<class RHS>
214double RK4IVP<RHS>::
215operator()(const double x) const
216{
217 // Dense output (O(h**3)) according to Shampine
218 // (Hermite interpolation)
219 const double h = stepsize();
220 int i = (x - span_[0]) / h;
221 const double t = (x - (span_[0] + i*h)) / h;
222
223 // Crude handling of evaluation point outside "span_";
224 if (i < 0) { i = 0; }
225 if (N_ <= i) { i = N_ - 1; }
226
227 const double y0 = y_[i], y1 = y_[i + 1];
228 const double f0 = f_[i], f1 = f_[i + 1];
229
230 double u = (1 - 2*t) * (y1 - y0);
231 u += h * ((t - 1)*f0 + t*f1);
232 u *= t * (t - 1);
233 u += (1 - t)*y0 + t*y1;
234
235 return u;
236}
237
238template<class RHS>
239double RK4IVP<RHS>::
240stepsize() const
241{
242 return (span_[1] - span_[0]) / N_;
243}
244
245namespace PhasePressODE {
246
247template<class FluidSystem>
248Water<FluidSystem>::
249Water(const TabulatedFunction& tempVdTable,
250 const TabulatedFunction& saltVdTable,
251 const int pvtRegionIdx,
252 const double normGrav)
253 : tempVdTable_(tempVdTable)
254 , saltVdTable_(saltVdTable)
255 , pvtRegionIdx_(pvtRegionIdx)
256 , g_(normGrav)
257{
258}
259
260template<class FluidSystem>
261double Water<FluidSystem>::
262operator()(const double depth,
263 const double press) const
264{
265 return this->density(depth, press) * g_;
266}
267
268template<class FluidSystem>
269double Water<FluidSystem>::
270density(const double depth,
271 const double press) const
272{
273 // The initializing algorithm can give depths outside the range due to numerical noise i.e. we extrapolate
274 double saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true);
275 double temp = tempVdTable_.eval(depth, /*extrapolate=*/true);
276 double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0 /*=Rsw*/, saltConcentration);
277 rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
278 return rho;
279}
280
281template<class FluidSystem, class RS>
282Oil<FluidSystem,RS>::
283Oil(const TabulatedFunction& tempVdTable,
284 const RS& rs,
285 const int pvtRegionIdx,
286 const double normGrav)
287 : tempVdTable_(tempVdTable)
288 , rs_(rs)
289 , pvtRegionIdx_(pvtRegionIdx)
290 , g_(normGrav)
291{
292}
293
294template<class FluidSystem, class RS>
295double Oil<FluidSystem,RS>::
296operator()(const double depth,
297 const double press) const
298{
299 return this->density(depth, press) * g_;
300}
301
302template<class FluidSystem, class RS>
303double Oil<FluidSystem,RS>::
304density(const double depth,
305 const double press) const
306{
307 const double temp = tempVdTable_.eval(depth, /*extrapolate=*/true);
308 double rs = 0.0;
309 if(FluidSystem::enableDissolvedGas())
310 rs = rs_(depth, press, temp);
311
312 double bOil = 0.0;
313 if (rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press)) {
314 bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
315 }
316 else {
317 bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rs);
318 }
319 double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
320 if (FluidSystem::enableDissolvedGas()) {
321 rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
322 }
323
324 return rho;
325}
326
327template<class FluidSystem, class RV, class RVW>
328Gas<FluidSystem,RV,RVW>::
329Gas(const TabulatedFunction& tempVdTable,
330 const RV& rv,
331 const RVW& rvw,
332 const int pvtRegionIdx,
333 const double normGrav)
334 : tempVdTable_(tempVdTable)
335 , rv_(rv)
336 , rvw_(rvw)
337 , pvtRegionIdx_(pvtRegionIdx)
338 , g_(normGrav)
339{
340}
341
342template<class FluidSystem, class RV, class RVW>
343double Gas<FluidSystem,RV,RVW>::
344operator()(const double depth,
345 const double press) const
346{
347 return this->density(depth, press) * g_;
348}
349
350template<class FluidSystem, class RV, class RVW>
351double Gas<FluidSystem,RV,RVW>::
352density(const double depth,
353 const double press) const
354{
355 const double temp = tempVdTable_.eval(depth, /*extrapolate=*/true);
356 double rv = 0.0;
357 if (FluidSystem::enableVaporizedOil())
358 rv = rv_(depth, press, temp);
359
360 double rvw = 0.0;
361 if (FluidSystem::enableVaporizedWater())
362 rvw = rvw_(depth, press, temp);
363
364 double bGas = 0.0;
365
366 if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) {
367 if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)
368 && rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press))
369 {
370 bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
371 } else {
372 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, rvw);
373 }
374 double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
375 rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_)
376 + rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
377 return rho;
378 }
379
380 if (FluidSystem::enableVaporizedOil()){
381 if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)) {
382 bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
383 } else {
384 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, 0.0/*=rvw*/);
385 }
386 double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
387 rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
388 return rho;
389 }
390
391 if (FluidSystem::enableVaporizedWater()){
392 if (rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press)) {
393 bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press);
394 }
395 else {
396 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0/*=rv*/, rvw);
397 }
398 double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
399 rho += rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
400 return rho;
401 }
402
403 // immiscible gas
404 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0/*=rv*/, 0.0/*=rvw*/);
405 double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
406
407 return rho;
408}
409
410}
411
412template<class FluidSystem, class Region>
413template<class ODE>
414PressureTable<FluidSystem,Region>::
415PressureFunction<ODE>::PressureFunction(const ODE& ode,
416 const InitCond& ic,
417 const int nsample,
418 const VSpan& span)
419 : initial_(ic)
420{
421 this->value_[Direction::Up] = std::make_unique<Distribution>
422 (ode, VSpan {{ ic.depth, span[0] }}, ic.pressure, nsample);
423
424 this->value_[Direction::Down] = std::make_unique<Distribution>
425 (ode, VSpan {{ ic.depth, span[1] }}, ic.pressure, nsample);
426}
427
428template<class FluidSystem, class Region>
429template<class ODE>
430PressureTable<FluidSystem,Region>::
431PressureFunction<ODE>::PressureFunction(const PressureFunction& rhs)
432 : initial_(rhs.initial_)
433{
434 this->value_[Direction::Up] =
435 std::make_unique<Distribution>(*rhs.value_[Direction::Up]);
436
437 this->value_[Direction::Down] =
438 std::make_unique<Distribution>(*rhs.value_[Direction::Down]);
439}
440
441template<class FluidSystem, class Region>
442template<class ODE>
443typename PressureTable<FluidSystem,Region>::template PressureFunction<ODE>&
444PressureTable<FluidSystem,Region>::
445PressureFunction<ODE>::
446operator=(const PressureFunction& rhs)
447{
448 this->initial_ = rhs.initial_;
449
450 this->value_[Direction::Up] =
451 std::make_unique<Distribution>(*rhs.value_[Direction::Up]);
452
453 this->value_[Direction::Down] =
454 std::make_unique<Distribution>(*rhs.value_[Direction::Down]);
455
456 return *this;
457}
458
459template<class FluidSystem, class Region>
460template<class ODE>
461typename PressureTable<FluidSystem,Region>::template PressureFunction<ODE>&
462PressureTable<FluidSystem,Region>::
463PressureFunction<ODE>::
464operator=(PressureFunction&& rhs)
465{
466 this->initial_ = rhs.initial_;
467 this->value_ = std::move(rhs.value_);
468
469 return *this;
470}
471
472template<class FluidSystem, class Region>
473template<class ODE>
474double
475PressureTable<FluidSystem,Region>::
476PressureFunction<ODE>::
477value(const double depth) const
478{
479 if (depth < this->initial_.depth) {
480 // Value above initial condition depth.
481 return (*this->value_[Direction::Up])(depth);
482 }
483 else if (depth > this->initial_.depth) {
484 // Value below initial condition depth.
485 return (*this->value_[Direction::Down])(depth);
486 }
487 else {
488 // Value *at* initial condition depth.
489 return this->initial_.pressure;
490 }
491}
492
493
494template<class FluidSystem, class Region>
495template<typename PressFunc>
496void PressureTable<FluidSystem,Region>::
497checkPtr(const PressFunc* phasePress,
498 const std::string& phaseName) const
499{
500 if (phasePress != nullptr) { return; }
501
502 throw std::invalid_argument {
503 "Phase pressure function for \"" + phaseName
504 + "\" most not be null"
505 };
506}
507
508template<class FluidSystem, class Region>
509typename PressureTable<FluidSystem,Region>::Strategy
510PressureTable<FluidSystem,Region>::
511selectEquilibrationStrategy(const Region& reg) const
512{
513 if (!this->oilActive()) {
514 if (reg.datum() > reg.zwoc()) { // Datum in water zone
515 return &PressureTable::equil_WOG;
516 }
517 return &PressureTable::equil_GOW;
518 }
519
520 if (reg.datum() > reg.zwoc()) { // Datum in water zone
521 return &PressureTable::equil_WOG;
522 }
523 else if (reg.datum() < reg.zgoc()) { // Datum in gas zone
524 return &PressureTable::equil_GOW;
525 }
526 else { // Datum in oil zone
527 return &PressureTable::equil_OWG;
528 }
529}
530
531template<class FluidSystem, class Region>
532void PressureTable<FluidSystem,Region>::
533copyInPointers(const PressureTable& rhs)
534{
535 if (rhs.oil_ != nullptr) {
536 this->oil_ = std::make_unique<OPress>(*rhs.oil_);
537 }
538
539 if (rhs.gas_ != nullptr) {
540 this->gas_ = std::make_unique<GPress>(*rhs.gas_);
541 }
542
543 if (rhs.wat_ != nullptr) {
544 this->wat_ = std::make_unique<WPress>(*rhs.wat_);
545 }
546}
547
548template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
550PhaseSaturations(MaterialLawManager& matLawMgr,
551 const std::vector<double>& swatInit)
552 : matLawMgr_(matLawMgr)
553 , swatInit_ (swatInit)
554{
555}
556
557template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
560 : matLawMgr_(rhs.matLawMgr_)
561 , swatInit_ (rhs.swatInit_)
562 , sat_ (rhs.sat_)
563 , press_ (rhs.press_)
564{
565 // Note: We don't need to do anything to the 'fluidState_' here.
566 this->setEvaluationPoint(*rhs.evalPt_.position,
567 *rhs.evalPt_.region,
568 *rhs.evalPt_.ptable);
569}
570
571template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
575 const Region& reg,
576 const PTable& ptable)
577{
578 this->setEvaluationPoint(x, reg, ptable);
579 this->initializePhaseQuantities();
580
581 if (ptable.gasActive()) { this->deriveGasSat(); }
582
583 if (ptable.waterActive()) { this->deriveWaterSat(); }
584
585
586 if (this->isOverlappingTransition()) {
587 this->fixUnphysicalTransition();
588 }
589
590 if (ptable.oilActive()) { this->deriveOilSat(); }
591
592 this->accountForScaledSaturations();
593
594 return this->sat_;
595}
596
597template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
599setEvaluationPoint(const Position& x,
600 const Region& reg,
601 const PTable& ptable)
602{
603 this->evalPt_.position = &x;
604 this->evalPt_.region = &reg;
605 this->evalPt_.ptable = &ptable;
606}
607
608template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
609void PhaseSaturations<MaterialLawManager,FluidSystem,Region,CellID>::
610initializePhaseQuantities()
611{
612 this->sat_.reset();
613 this->press_.reset();
614
615 const auto depth = this->evalPt_.position->depth;
616 const auto& ptable = *this->evalPt_.ptable;
617
618 if (ptable.oilActive()) {
619 this->press_.oil = ptable.oil(depth);
620 }
621
622 if (ptable.gasActive()) {
623 this->press_.gas = ptable.gas(depth);
624 }
625
626 if (ptable.waterActive()) {
627 this->press_.water = ptable.water(depth);
628 }
629}
630
631template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
632void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveOilSat()
633{
634 this->sat_.oil = 1.0 - this->sat_.water - this->sat_.gas;
635}
636
637template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
638void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveGasSat()
639{
640 auto& sg = this->sat_.gas;
641
642 const auto isIncr = true; // dPcgo/dSg >= 0 for all Sg.
643 const auto oilActive = this->evalPt_.ptable->oilActive();
644
645 if (this->isConstCapPress(this->gasPos())) {
646 // Sharp interface between phases. Can derive phase saturation
647 // directly from knowing where 'depth' of evaluation point is
648 // relative to depth of O/G contact.
649 const auto gas_contact = oilActive? this->evalPt_.region->zgoc() : this->evalPt_.region->zwoc();
650 sg = this->fromDepthTable(gas_contact,
651 this->gasPos(), isIncr);
652 }
653 else {
654 // Capillary pressure curve is non-constant, meaning there is a
655 // transition zone between the gas and oil phases. Invert capillary
656 // pressure relation
657 //
658 // Pcgo(Sg) = Pg - Po
659 //
660 // Note that Pcgo is defined to be (Pg - Po), not (Po - Pg).
661 const auto pw = oilActive? this->press_.oil : this->press_.water;
662 const auto pcgo = this->press_.gas - pw;
663 sg = this->invertCapPress(pcgo, this->gasPos(), isIncr);
664 }
665}
666
667template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
668void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::deriveWaterSat()
669{
670 auto& sw = this->sat_.water;
671
672 const auto oilActive = this->evalPt_.ptable->oilActive();
673 if (!oilActive) {
674 // for 2p gas+water we set the water saturation to 1.0 - sg
675 sw = 1.0 - this->sat_.gas;
676 }
677 else {
678 const auto isIncr = false; // dPcow/dSw <= 0 for all Sw.
679
680 if (this->isConstCapPress(this->waterPos())) {
681 // Sharp interface between phases. Can derive phase saturation
682 // directly from knowing where 'depth' of evaluation point is
683 // relative to depth of O/W contact.
684 sw = this->fromDepthTable(this->evalPt_.region->zwoc(),
685 this->waterPos(), isIncr);
686 }
687 else {
688 // Capillary pressure curve is non-constant, meaning there is a
689 // transition zone between the oil and water phases. Invert
690 // capillary pressure relation
691 //
692 // Pcow(Sw) = Po - Pw
693 //
694 // unless the model uses "SWATINIT". In the latter case, pick the
695 // saturation directly from the SWATINIT array of the pertinent
696 // cell.
697 const auto pcow = this->press_.oil - this->press_.water;
698
699 if (this->swatInit_.empty()) {
700 sw = this->invertCapPress(pcow, this->waterPos(), isIncr);
701 }
702 else {
703 auto [swout, newSwatInit] = this->applySwatInit(pcow);
704 if (newSwatInit)
705 sw = this->invertCapPress(pcow, this->waterPos(), isIncr);
706 else {
707 sw = swout;
708 }
709 }
710 }
711 }
712}
713
714template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
715void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
716fixUnphysicalTransition()
717{
718 auto& sg = this->sat_.gas;
719 auto& sw = this->sat_.water;
720
721 // Overlapping gas/oil and oil/water transition zones can lead to
722 // unphysical phase saturations when individual saturations are derived
723 // directly from inverting O/G and O/W capillary pressure curves.
724 //
725 // Recalculate phase saturations using the implied gas/water capillary
726 // pressure: Pg - Pw.
727 const auto pcgw = this->press_.gas - this->press_.water;
728 if (! this->swatInit_.empty()) {
729 // Re-scale Pc to reflect imposed sw for vanishing oil phase. This
730 // seems consistent with ECLIPSE, but fails to honour SWATINIT in
731 // case of non-trivial gas/oil capillary pressure.
732 auto [swout, newSwatInit] = this->applySwatInit(pcgw, sw);
733 if (newSwatInit){
734 const auto isIncr = false; // dPcow/dSw <= 0 for all Sw.
735 sw = this->invertCapPress(pcgw, this->waterPos(), isIncr);
736 }
737 else {
738 sw = swout;
739 }
740 }
741
742 sw = satFromSumOfPcs<FluidSystem>
743 (this->matLawMgr_, this->waterPos(), this->gasPos(),
744 this->evalPt_.position->cell, pcgw);
745 sg = 1.0 - sw;
746
747 this->fluidState_.setSaturation(this->oilPos(), 1.0 - sw - sg);
748 this->fluidState_.setSaturation(this->gasPos(), sg);
749 this->fluidState_.setSaturation(this->waterPos(), this->evalPt_
750 .ptable->waterActive() ? sw : 0.0);
751
752 // Pcgo = Pg - Po => Po = Pg - Pcgo
753 this->computeMaterialLawCapPress();
754 this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil();
755}
756
757template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
758void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
759accountForScaledSaturations()
760{
761 const auto gasActive = this->evalPt_.ptable->gasActive();
762 const auto watActive = this->evalPt_.ptable->waterActive();
763 const auto oilActive = this->evalPt_.ptable->oilActive();
764
765 auto sg = gasActive? this->sat_.gas : 0.0;
766 auto sw = watActive? this->sat_.water : 0.0;
767 auto so = oilActive? this->sat_.oil : 0.0;
768
769 this->fluidState_.setSaturation(this->waterPos(), sw);
770 this->fluidState_.setSaturation(this->oilPos(), so);
771 this->fluidState_.setSaturation(this->gasPos(), sg);
772
773 const auto& scaledDrainageInfo = this->matLawMgr_
774 .oilWaterScaledEpsInfoDrainage(this->evalPt_.position->cell);
775
776 const auto thresholdSat = 1.0e-6;
777 if (watActive && ((sw + thresholdSat) > scaledDrainageInfo.Swu)) {
778 // Water saturation exceeds maximum possible value. Reset oil phase
779 // pressure to that which corresponds to maximum possible water
780 // saturation value.
781 this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swu);
782 if (oilActive) {
783 this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swu);
784 } else if (gasActive) {
785 this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swu);
786 }
787 sw = scaledDrainageInfo.Swu;
788 this->computeMaterialLawCapPress();
789
790 if (oilActive) {
791 // Pcow = Po - Pw => Po = Pw + Pcow
792 this->press_.oil = this->press_.water + this->materialLawCapPressOilWater();
793 } else {
794 // Pcgw = Pg - Pw => Pg = Pw + Pcgw
795 this->press_.gas = this->press_.water + this->materialLawCapPressGasWater();
796 }
797
798 }
799 if (gasActive && ((sg + thresholdSat) > scaledDrainageInfo.Sgu)) {
800 // Gas saturation exceeds maximum possible value. Reset oil phase
801 // pressure to that which corresponds to maximum possible gas
802 // saturation value.
803 this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgu);
804 if (oilActive) {
805 this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgu);
806 } else if (watActive) {
807 this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgu);
808 }
809 sg = scaledDrainageInfo.Sgu;
810 this->computeMaterialLawCapPress();
811
812 if (oilActive) {
813 // Pcgo = Pg - Po => Po = Pg - Pcgo
814 this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil();
815 } else {
816 // Pcgw = Pg - Pw => Pw = Pg - Pcgw
817 this->press_.water = this->press_.gas - this->materialLawCapPressGasWater();
818 }
819 }
820
821 if (watActive && ((sw - thresholdSat) < scaledDrainageInfo.Swl)) {
822 // Water saturation less than minimum possible value in cell. Reset
823 // water phase pressure to that which corresponds to minimum
824 // possible water saturation value.
825 this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swl);
826 if (oilActive) {
827 this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swl);
828 } else if (gasActive) {
829 this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swl);
830 }
831 sw = scaledDrainageInfo.Swl;
832 this->computeMaterialLawCapPress();
833
834 if (oilActive) {
835 // Pcwo = Po - Pw => Pw = Po - Pcow
836 this->press_.water = this->press_.oil - this->materialLawCapPressOilWater();
837 } else {
838 // Pcgw = Pg - Pw => Pw = Pg - Pcgw
839 this->press_.water = this->press_.gas - this->materialLawCapPressGasWater();
840 }
841 }
842
843 if (gasActive && ((sg - thresholdSat) < scaledDrainageInfo.Sgl)) {
844 // Gas saturation less than minimum possible value in cell. Reset
845 // gas phase pressure to that which corresponds to minimum possible
846 // gas saturation.
847 this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgl);
848 if (oilActive) {
849 this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgl);
850 } else if (watActive) {
851 this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgl);
852 }
853 sg = scaledDrainageInfo.Sgl;
854 this->computeMaterialLawCapPress();
855
856 if (oilActive) {
857 // Pcgo = Pg - Po => Pg = Po + Pcgo
858 this->press_.gas = this->press_.oil + this->materialLawCapPressGasOil();
859 } else {
860 // Pcgw = Pg - Pw => Pg = Pw + Pcgw
861 this->press_.gas = this->press_.water + this->materialLawCapPressGasWater();
862 }
863 }
864}
865
866template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
867std::pair<double, bool>
868PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
869applySwatInit(const double pcow)
870{
871 return this->applySwatInit(pcow, this->swatInit_[this->evalPt_.position->cell]);
872}
873
874template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
875std::pair<double, bool>
876PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
877applySwatInit(const double pcow, const double sw)
878{
879 return this->matLawMgr_.applySwatinit(this->evalPt_.position->cell, pcow, sw);
880}
881
882template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
883void PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
884computeMaterialLawCapPress()
885{
886 const auto& matParams = this->matLawMgr_
887 .materialLawParams(this->evalPt_.position->cell);
888
889 this->matLawCapPress_.fill(0.0);
890 MaterialLaw::capillaryPressures(this->matLawCapPress_,
891 matParams, this->fluidState_);
892}
893
894template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
895double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
896materialLawCapPressGasOil() const
897{
898 return this->matLawCapPress_[this->oilPos()]
899 + this->matLawCapPress_[this->gasPos()];
900}
901
902template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
903double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
904materialLawCapPressOilWater() const
905{
906 return this->matLawCapPress_[this->oilPos()]
907 - this->matLawCapPress_[this->waterPos()];
908}
909
910template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
911double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
912materialLawCapPressGasWater() const
913{
914 return this->matLawCapPress_[this->gasPos()]
915 - this->matLawCapPress_[this->waterPos()];
916}
917
918template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
919bool PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
920isConstCapPress(const PhaseIdx phaseIdx) const
921{
922 return isConstPc<FluidSystem>
923 (this->matLawMgr_, phaseIdx, this->evalPt_.position->cell);
924}
925
926template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
927bool PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
928isOverlappingTransition() const
929{
930 return this->evalPt_.ptable->gasActive()
931 && this->evalPt_.ptable->waterActive()
932 && ((this->sat_.gas + this->sat_.water) > 1.0);
933}
934
935template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
936double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
937fromDepthTable(const double contactdepth,
938 const PhaseIdx phasePos,
939 const bool isincr) const
940{
941 return satFromDepth<FluidSystem>
942 (this->matLawMgr_, this->evalPt_.position->depth,
943 contactdepth, static_cast<int>(phasePos),
944 this->evalPt_.position->cell, isincr);
945}
946
947template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
948double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
949invertCapPress(const double pc,
950 const PhaseIdx phasePos,
951 const bool isincr) const
952{
953 return satFromPc<FluidSystem>
954 (this->matLawMgr_, static_cast<int>(phasePos),
955 this->evalPt_.position->cell, pc, isincr);
956}
957
958template<class FluidSystem, class Region>
960PressureTable(const double gravity,
961 const int samplePoints)
962 : gravity_(gravity)
963 , nsample_(samplePoints)
964{
965}
966
967template <class FluidSystem, class Region>
970 : gravity_(rhs.gravity_)
971 , nsample_(rhs.nsample_)
972{
973 this->copyInPointers(rhs);
974}
975
976template <class FluidSystem, class Region>
979 : gravity_(rhs.gravity_)
980 , nsample_(rhs.nsample_)
981 , oil_ (std::move(rhs.oil_))
982 , gas_ (std::move(rhs.gas_))
983 , wat_ (std::move(rhs.wat_))
984{
985}
986
987template <class FluidSystem, class Region>
991{
992 this->gravity_ = rhs.gravity_;
993 this->nsample_ = rhs.nsample_;
994 this->copyInPointers(rhs);
995
996 return *this;
997}
998
999template <class FluidSystem, class Region>
1003{
1004 this->gravity_ = rhs.gravity_;
1005 this->nsample_ = rhs.nsample_;
1006
1007 this->oil_ = std::move(rhs.oil_);
1008 this->gas_ = std::move(rhs.gas_);
1009 this->wat_ = std::move(rhs.wat_);
1010
1011 return *this;
1012}
1013
1014template <class FluidSystem, class Region>
1016equilibrate(const Region& reg,
1017 const VSpan& span)
1018{
1019 // One of the PressureTable::equil_*() member functions.
1020 auto equil = this->selectEquilibrationStrategy(reg);
1021
1022 (this->*equil)(reg, span);
1023}
1024
1025template <class FluidSystem, class Region>
1027oilActive() const
1028{
1029 return FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1030}
1031
1032template <class FluidSystem, class Region>
1034gasActive() const
1035{
1036 return FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
1037}
1038
1039template <class FluidSystem, class Region>
1041waterActive() const
1042{
1043 return FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1044}
1045
1046template <class FluidSystem, class Region>
1048oil(const double depth) const
1049{
1050 this->checkPtr(this->oil_.get(), "OIL");
1051
1052 return this->oil_->value(depth);
1053}
1054
1055template <class FluidSystem, class Region>
1057gas(const double depth) const
1058{
1059 this->checkPtr(this->gas_.get(), "GAS");
1060
1061 return this->gas_->value(depth);
1062}
1063
1064
1065template <class FluidSystem, class Region>
1067water(const double depth) const
1068{
1069 this->checkPtr(this->wat_.get(), "WATER");
1070
1071 return this->wat_->value(depth);
1072}
1073
1074template <class FluidSystem, class Region>
1076equil_WOG(const Region& reg, const VSpan& span)
1077{
1078 // Datum depth in water zone. Calculate phase pressure for water first,
1079 // followed by oil and gas if applicable.
1080
1081 if (! this->waterActive()) {
1082 throw std::invalid_argument {
1083 "Don't know how to interpret EQUIL datum depth in "
1084 "WATER zone in model without active water phase"
1085 };
1086 }
1087
1088 {
1089 const auto ic = typename WPress::InitCond {
1090 reg.datum(), reg.pressure()
1091 };
1092
1093 this->makeWatPressure(ic, reg, span);
1094 }
1095
1096 if (this->oilActive()) {
1097 // Pcow = Po - Pw => Po = Pw + Pcow
1098 const auto ic = typename OPress::InitCond {
1099 reg.zwoc(),
1100 this->water(reg.zwoc()) + reg.pcowWoc()
1101 };
1102
1103 this->makeOilPressure(ic, reg, span);
1104 }
1105
1106 if (this->gasActive() && this->oilActive()) {
1107 // Pcgo = Pg - Po => Pg = Po + Pcgo
1108 const auto ic = typename GPress::InitCond {
1109 reg.zgoc(),
1110 this->oil(reg.zgoc()) + reg.pcgoGoc()
1111 };
1112
1113 this->makeGasPressure(ic, reg, span);
1114 } else if (this->gasActive() && !this->oilActive()) {
1115 // No oil phase set Pg = Pw + Pcgw
1116 const auto ic = typename GPress::InitCond {
1117 reg.zwoc(), // The WOC is really the GWC for gas/water cases
1118 this->water(reg.zwoc()) + reg.pcowWoc() // Pcow(WOC) is really Pcgw(GWC) for gas/water cases
1119 };
1120 this->makeGasPressure(ic, reg, span);
1121 }
1122}
1123
1124template <class FluidSystem, class Region>
1125void PressureTable<FluidSystem, Region>::
1126equil_GOW(const Region& reg, const VSpan& span)
1127{
1128 // Datum depth in gas zone. Calculate phase pressure for gas first,
1129 // followed by oil and water if applicable.
1130
1131 if (! this->gasActive()) {
1132 throw std::invalid_argument {
1133 "Don't know how to interpret EQUIL datum depth in "
1134 "GAS zone in model without active gas phase"
1135 };
1136 }
1137
1138 {
1139 const auto ic = typename GPress::InitCond {
1140 reg.datum(), reg.pressure()
1141 };
1142
1143 this->makeGasPressure(ic, reg, span);
1144 }
1145
1146 if (this->oilActive()) {
1147 // Pcgo = Pg - Po => Po = Pg - Pcgo
1148 const auto ic = typename OPress::InitCond {
1149 reg.zgoc(),
1150 this->gas(reg.zgoc()) - reg.pcgoGoc()
1151 };
1152 this->makeOilPressure(ic, reg, span);
1153 }
1154
1155 if (this->waterActive() && this->oilActive()) {
1156 // Pcow = Po - Pw => Pw = Po - Pcow
1157 const auto ic = typename WPress::InitCond {
1158 reg.zwoc(),
1159 this->oil(reg.zwoc()) - reg.pcowWoc()
1160 };
1161
1162 this->makeWatPressure(ic, reg, span);
1163 } else if (this->waterActive() && !this->oilActive()) {
1164 // No oil phase set Pw = Pg - Pcgw
1165 const auto ic = typename WPress::InitCond {
1166 reg.zwoc(), // The WOC is really the GWC for gas/water cases
1167 this->gas(reg.zwoc()) - reg.pcowWoc() // Pcow(WOC) is really Pcgw(GWC) for gas/water cases
1168 };
1169 this->makeWatPressure(ic, reg, span);
1170 }
1171}
1172
1173template <class FluidSystem, class Region>
1174void PressureTable<FluidSystem, Region>::
1175equil_OWG(const Region& reg, const VSpan& span)
1176{
1177 // Datum depth in oil zone. Calculate phase pressure for oil first,
1178 // followed by gas and water if applicable.
1179
1180 if (! this->oilActive()) {
1181 throw std::invalid_argument {
1182 "Don't know how to interpret EQUIL datum depth in "
1183 "OIL zone in model without active oil phase"
1184 };
1185 }
1186
1187 {
1188 const auto ic = typename OPress::InitCond {
1189 reg.datum(), reg.pressure()
1190 };
1191
1192 this->makeOilPressure(ic, reg, span);
1193 }
1194
1195 if (this->waterActive()) {
1196 // Pcow = Po - Pw => Pw = Po - Pcow
1197 const auto ic = typename WPress::InitCond {
1198 reg.zwoc(),
1199 this->oil(reg.zwoc()) - reg.pcowWoc()
1200 };
1201
1202 this->makeWatPressure(ic, reg, span);
1203 }
1204
1205 if (this->gasActive()) {
1206 // Pcgo = Pg - Po => Pg = Po + Pcgo
1207 const auto ic = typename GPress::InitCond {
1208 reg.zgoc(),
1209 this->oil(reg.zgoc()) + reg.pcgoGoc()
1210 };
1211 this->makeGasPressure(ic, reg, span);
1212 }
1213}
1214
1215template <class FluidSystem, class Region>
1216void PressureTable<FluidSystem, Region>::
1217makeOilPressure(const typename OPress::InitCond& ic,
1218 const Region& reg,
1219 const VSpan& span)
1220{
1221 const auto drho = OilPressODE {
1222 reg.tempVdTable(), reg.dissolutionCalculator(),
1223 reg.pvtIdx(), this->gravity_
1224 };
1225
1226 this->oil_ = std::make_unique<OPress>(drho, ic, this->nsample_, span);
1227}
1228
1229template <class FluidSystem, class Region>
1230void PressureTable<FluidSystem, Region>::
1231makeGasPressure(const typename GPress::InitCond& ic,
1232 const Region& reg,
1233 const VSpan& span)
1234{
1235 const auto drho = GasPressODE {
1236 reg.tempVdTable(), reg.evaporationCalculator(), reg.waterEvaporationCalculator(),
1237 reg.pvtIdx(), this->gravity_
1238 };
1239
1240 this->gas_ = std::make_unique<GPress>(drho, ic, this->nsample_, span);
1241}
1242
1243template <class FluidSystem, class Region>
1244void PressureTable<FluidSystem, Region>::
1245makeWatPressure(const typename WPress::InitCond& ic,
1246 const Region& reg,
1247 const VSpan& span)
1248{
1249 const auto drho = WatPressODE {
1250 reg.tempVdTable(), reg.saltVdTable(), reg.pvtIdx(), this->gravity_
1251 };
1252
1253 this->wat_ = std::make_unique<WPress>(drho, ic, this->nsample_, span);
1254}
1255
1256}
1257
1258namespace DeckDependent {
1259
1260std::vector<EquilRecord>
1261getEquil(const EclipseState& state)
1262{
1263 const auto& init = state.getInitConfig();
1264
1265 if(!init.hasEquil()) {
1266 throw std::domain_error("Deck does not provide equilibration data.");
1267 }
1268
1269 const auto& equil = init.getEquil();
1270 return { equil.begin(), equil.end() };
1271}
1272
1273template<class GridView>
1274std::vector<int>
1275equilnum(const EclipseState& eclipseState,
1276 const GridView& gridview)
1277{
1278 std::vector<int> eqlnum(gridview.size(0), 0);
1279
1280 if (eclipseState.fieldProps().has_int("EQLNUM")) {
1281 const auto& e = eclipseState.fieldProps().get_int("EQLNUM");
1282 std::transform(e.begin(), e.end(), eqlnum.begin(), [](int n){ return n - 1;});
1283 }
1284 OPM_BEGIN_PARALLEL_TRY_CATCH();
1285 const int num_regions = eclipseState.getTableManager().getEqldims().getNumEquilRegions();
1286 if ( std::any_of(eqlnum.begin(), eqlnum.end(), [num_regions](int n){return n >= num_regions;}) ) {
1287 throw std::runtime_error("Values larger than maximum Equil regions " + std::to_string(num_regions) + " provided in EQLNUM");
1288 }
1289 if ( std::any_of(eqlnum.begin(), eqlnum.end(), [](int n){return n < 0;}) ) {
1290 throw std::runtime_error("zero or negative values provided in EQLNUM");
1291 }
1292 OPM_END_PARALLEL_TRY_CATCH("Invalied EQLNUM numbers: ", gridview.comm());
1293
1294 return eqlnum;
1295}
1296
1297template<class FluidSystem,
1298 class Grid,
1299 class GridView,
1300 class ElementMapper,
1301 class CartesianIndexMapper>
1302template<class MaterialLawManager>
1303InitialStateComputer<FluidSystem,
1304 Grid,
1305 GridView,
1306 ElementMapper,
1307 CartesianIndexMapper>::
1308InitialStateComputer(MaterialLawManager& materialLawManager,
1309 const EclipseState& eclipseState,
1310 const Grid& grid,
1311 const GridView& gridView,
1312 const CartesianIndexMapper& cartMapper,
1313 const double grav,
1314 const int num_pressure_points,
1315 const bool applySwatInit)
1316 : temperature_(grid.size(/*codim=*/0), eclipseState.getTableManager().rtemp()),
1317 saltConcentration_(grid.size(/*codim=*/0)),
1318 saltSaturation_(grid.size(/*codim=*/0)),
1319 pp_(FluidSystem::numPhases,
1320 std::vector<double>(grid.size(/*codim=*/0))),
1321 sat_(FluidSystem::numPhases,
1322 std::vector<double>(grid.size(/*codim=*/0))),
1323 rs_(grid.size(/*codim=*/0)),
1324 rv_(grid.size(/*codim=*/0)),
1325 rvw_(grid.size(/*codim=*/0)),
1326 cartesianIndexMapper_(cartMapper),
1327 num_pressure_points_(num_pressure_points)
1328{
1329 //Check for presence of kw SWATINIT
1330 if (applySwatInit) {
1331 if (eclipseState.fieldProps().has_double("SWATINIT")) {
1332 swatInit_ = eclipseState.fieldProps().get_double("SWATINIT");
1333 }
1334 }
1335
1336 // Querry cell depth, cell top-bottom.
1337 // numerical aquifer cells might be specified with different depths.
1338 const auto& num_aquifers = eclipseState.aquifer().numericalAquifers();
1339 updateCellProps_(gridView, num_aquifers);
1340
1341 // Get the equilibration records.
1342 const std::vector<EquilRecord> rec = getEquil(eclipseState);
1343 const auto& tables = eclipseState.getTableManager();
1344 // Create (inverse) region mapping.
1345 const RegionMapping<> eqlmap(equilnum(eclipseState, grid));
1346 const int invalidRegion = -1;
1347 regionPvtIdx_.resize(rec.size(), invalidRegion);
1348 setRegionPvtIdx(eclipseState, eqlmap);
1349
1350 // Create Rs functions.
1351 rsFunc_.reserve(rec.size());
1352 if (FluidSystem::enableDissolvedGas()) {
1353 for (std::size_t i = 0; i < rec.size(); ++i) {
1354 if (eqlmap.cells(i).empty()) {
1355 rsFunc_.push_back(std::shared_ptr<Miscibility::RsVD<FluidSystem>>());
1356 continue;
1357 }
1358 const int pvtIdx = regionPvtIdx_[i];
1359 if (!rec[i].liveOilInitConstantRs()) {
1360 const TableContainer& rsvdTables = tables.getRsvdTables();
1361 const TableContainer& pbvdTables = tables.getPbvdTables();
1362 if (rsvdTables.size() > 0) {
1363
1364 const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i);
1365 std::vector<double> depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy();
1366 std::vector<double> rsColumn = rsvdTable.getColumn("RS").vectorCopy();
1367 rsFunc_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
1368 depthColumn, rsColumn));
1369 } else if (pbvdTables.size() > 0) {
1370 const PbvdTable& pbvdTable = pbvdTables.getTable<PbvdTable>(i);
1371 std::vector<double> depthColumn = pbvdTable.getColumn("DEPTH").vectorCopy();
1372 std::vector<double> pbubColumn = pbvdTable.getColumn("PBUB").vectorCopy();
1373 rsFunc_.push_back(std::make_shared<Miscibility::PBVD<FluidSystem>>(pvtIdx,
1374 depthColumn, pbubColumn));
1375
1376 } else {
1377 throw std::runtime_error("Cannot initialise: RSVD or PBVD table not available.");
1378 }
1379
1380 }
1381 else {
1382 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
1383 throw std::runtime_error("Cannot initialise: when no explicit RSVD table is given, \n"
1384 "datum depth must be at the gas-oil-contact. "
1385 "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold.");
1386 }
1387 const double pContact = rec[i].datumDepthPressure();
1388 const double TContact = 273.15 + 20; // standard temperature for now
1389 rsFunc_.push_back(std::make_shared<Miscibility::RsSatAtContact<FluidSystem>>(pvtIdx, pContact, TContact));
1390 }
1391 }
1392 }
1393 else {
1394 for (std::size_t i = 0; i < rec.size(); ++i) {
1395 rsFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1396 }
1397 }
1398
1399 rvFunc_.reserve(rec.size());
1400 if (FluidSystem::enableVaporizedOil()) {
1401 for (std::size_t i = 0; i < rec.size(); ++i) {
1402 if (eqlmap.cells(i).empty()) {
1403 rvFunc_.push_back(std::shared_ptr<Miscibility::RvVD<FluidSystem>>());
1404 continue;
1405 }
1406 const int pvtIdx = regionPvtIdx_[i];
1407 if (!rec[i].wetGasInitConstantRv()) {
1408 const TableContainer& rvvdTables = tables.getRvvdTables();
1409 const TableContainer& pdvdTables = tables.getPdvdTables();
1410
1411 if (rvvdTables.size() > 0) {
1412 const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i);
1413 std::vector<double> depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy();
1414 std::vector<double> rvColumn = rvvdTable.getColumn("RV").vectorCopy();
1415 rvFunc_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
1416 depthColumn, rvColumn));
1417 } else if (pdvdTables.size() > 0) {
1418 const PdvdTable& pdvdTable = pdvdTables.getTable<PdvdTable>(i);
1419 std::vector<double> depthColumn = pdvdTable.getColumn("DEPTH").vectorCopy();
1420 std::vector<double> pdewColumn = pdvdTable.getColumn("PDEW").vectorCopy();
1421 rvFunc_.push_back(std::make_shared<Miscibility::PDVD<FluidSystem>>(pvtIdx,
1422 depthColumn, pdewColumn));
1423 } else {
1424 throw std::runtime_error("Cannot initialise: RVVD or PDCD table not available.");
1425 }
1426 }
1427 else {
1428 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
1429 throw std::runtime_error(
1430 "Cannot initialise: when no explicit RVVD table is given, \n"
1431 "datum depth must be at the gas-oil-contact. "
1432 "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold.");
1433 }
1434 const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
1435 const double TContact = 273.15 + 20; // standard temperature for now
1436 rvFunc_.push_back(std::make_shared<Miscibility::RvSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
1437 }
1438 }
1439 }
1440 else {
1441 for (std::size_t i = 0; i < rec.size(); ++i) {
1442 rvFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1443 }
1444 }
1445
1446 rvwFunc_.reserve(rec.size());
1447 if (FluidSystem::enableVaporizedWater()) {
1448 for (std::size_t i = 0; i < rec.size(); ++i) {
1449 if (eqlmap.cells(i).empty()) {
1450 rvwFunc_.push_back(std::shared_ptr<Miscibility::RvwVD<FluidSystem>>());
1451 continue;
1452 }
1453 const int pvtIdx = regionPvtIdx_[i];
1454 if (!rec[i].humidGasInitConstantRvw()) {
1455 const TableContainer& rvwvdTables = tables.getRvwvdTables();
1456
1457 if (rvwvdTables.size() > 0) {
1458 const RvwvdTable& rvwvdTable = rvwvdTables.getTable<RvwvdTable>(i);
1459 std::vector<double> depthColumn = rvwvdTable.getColumn("DEPTH").vectorCopy();
1460 std::vector<double> rvwvdColumn = rvwvdTable.getColumn("RVWVD").vectorCopy();
1461 rvwFunc_.push_back(std::make_shared<Miscibility::RvwVD<FluidSystem>>(pvtIdx,
1462 depthColumn, rvwvdColumn));
1463 } else {
1464 throw std::runtime_error("Cannot initialise: RVWVD table not available.");
1465 }
1466 }
1467 else {
1468 const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1469 if (oilActive){
1470 if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
1471 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1472 const auto msg = "No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +". \n"
1473 "and datum depth is not at the gas-oil-contact. \n"
1474 "Rvw is set to 0.0 in all cells. \n";
1475 OpmLog::warning(msg);
1476 } else {
1477 // pg = po + Pcgo = po + (pg - po)
1478 // for gas-condensate with initial no oil zone: water-oil contact depth (OWC) equal gas-oil contact depth (GOC)
1479 const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
1480 const double TContact = 273.15 + 20; // standard temperature for now
1481 rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
1482 }
1483 }
1484 else {
1485 // two-phase gas-water sytem: water-oil contact depth is taken equal to gas-water contact depth (GWC)
1486 // and water-oil capillary pressure (Pcwo) is taken equal to gas-water capillary pressure (Pcgw) at GWC
1487 if (rec[i].waterOilContactDepth() != rec[i].datumDepth()) {
1488 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1489 const auto msg = "No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +". \n"
1490 "and datum depth is not at the gas-water-contact. \n"
1491 "Rvw is set to 0.0 in all cells. \n";
1492 OpmLog::warning(msg);
1493 } else {
1494 // pg = pw + Pcgw = pw + (pg - pw)
1495 const double pContact = rec[i].datumDepthPressure() + rec[i].waterOilContactCapillaryPressure();
1496 const double TContact = 273.15 + 20; // standard temperature for now
1497 rvwFunc_.push_back(std::make_shared<Miscibility::RvwSatAtContact<FluidSystem>>(pvtIdx,pContact, TContact));
1498 }
1499 }
1500 }
1501 }
1502 }
1503 else {
1504 for (std::size_t i = 0; i < rec.size(); ++i) {
1505 rvwFunc_.push_back(std::make_shared<Miscibility::NoMixing>());
1506 }
1507 }
1508
1509
1510 // EXTRACT the initial temperature
1511 updateInitialTemperature_(eclipseState, eqlmap);
1512
1513 // EXTRACT the initial salt concentration
1514 updateInitialSaltConcentration_(eclipseState, eqlmap);
1515
1516 // EXTRACT the initial salt saturation
1517 updateInitialSaltSaturation_(eclipseState, eqlmap);
1518
1519 // Compute pressures, saturations, rs and rv factors.
1520 const auto& comm = grid.comm();
1521 calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav);
1522
1523 // modify the pressure and saturation for numerical aquifer cells
1524 applyNumericalAquifers_(gridView, num_aquifers, eclipseState.runspec().co2Storage() || eclipseState.runspec().h2Storage());
1525
1526 // Modify oil pressure in no-oil regions so that the pressures of present phases can
1527 // be recovered from the oil pressure and capillary relations.
1528}
1529
1530template<class FluidSystem,
1531 class Grid,
1532 class GridView,
1533 class ElementMapper,
1534 class CartesianIndexMapper>
1535template<class RMap>
1536void InitialStateComputer<FluidSystem,
1537 Grid,
1538 GridView,
1539 ElementMapper,
1540 CartesianIndexMapper>::
1541updateInitialTemperature_(const EclipseState& eclState, const RMap& reg)
1542{
1543 const int numEquilReg = rsFunc_.size();
1544 tempVdTable_.resize(numEquilReg);
1545 const auto& tables = eclState.getTableManager();
1546 if (!tables.hasTables("RTEMPVD")) {
1547 std::vector<double> x = {0.0,1.0};
1548 std::vector<double> y = {tables.rtemp(),tables.rtemp()};
1549 for (auto& table : this->tempVdTable_) {
1550 table.setXYContainers(x, y);
1551 }
1552 } else {
1553 const TableContainer& tempvdTables = tables.getRtempvdTables();
1554 for (std::size_t i = 0; i < tempvdTables.size(); ++i) {
1555 const RtempvdTable& tempvdTable = tempvdTables.getTable<RtempvdTable>(i);
1556 tempVdTable_[i].setXYContainers(tempvdTable.getDepthColumn(), tempvdTable.getTemperatureColumn());
1557 const auto& cells = reg.cells(i);
1558 for (const auto& cell : cells) {
1559 const double depth = cellCenterDepth_[cell];
1560 this->temperature_[cell] = tempVdTable_[i].eval(depth, /*extrapolate=*/true);
1561 }
1562 }
1563 }
1564}
1565
1566template<class FluidSystem,
1567 class Grid,
1568 class GridView,
1569 class ElementMapper,
1570 class CartesianIndexMapper>
1571template<class RMap>
1572void InitialStateComputer<FluidSystem,
1573 Grid,
1574 GridView,
1575 ElementMapper,
1576 CartesianIndexMapper>::
1577updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg)
1578{
1579 const int numEquilReg = rsFunc_.size();
1580 saltVdTable_.resize(numEquilReg);
1581 const auto& tables = eclState.getTableManager();
1582 const TableContainer& saltvdTables = tables.getSaltvdTables();
1583
1584 // If no saltvd table is given, we create a trivial table for the density calculations
1585 if (saltvdTables.empty()) {
1586 std::vector<double> x = {0.0,1.0};
1587 std::vector<double> y = {0.0,0.0};
1588 for (auto& table : this->saltVdTable_) {
1589 table.setXYContainers(x, y);
1590 }
1591 } else {
1592 for (std::size_t i = 0; i < saltvdTables.size(); ++i) {
1593 const SaltvdTable& saltvdTable = saltvdTables.getTable<SaltvdTable>(i);
1594 saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn());
1595
1596 const auto& cells = reg.cells(i);
1597 for (const auto& cell : cells) {
1598 const double depth = cellCenterDepth_[cell];
1599 this->saltConcentration_[cell] = saltVdTable_[i].eval(depth, /*extrapolate=*/true);
1600 }
1601 }
1602 }
1603}
1604
1605template<class FluidSystem,
1606 class Grid,
1607 class GridView,
1608 class ElementMapper,
1609 class CartesianIndexMapper>
1610template<class RMap>
1611void InitialStateComputer<FluidSystem,
1612 Grid,
1613 GridView,
1614 ElementMapper,
1615 CartesianIndexMapper>::
1616updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg)
1617{
1618 const int numEquilReg = rsFunc_.size();
1619 saltpVdTable_.resize(numEquilReg);
1620 const auto& tables = eclState.getTableManager();
1621 const TableContainer& saltpvdTables = tables.getSaltpvdTables();
1622
1623 for (std::size_t i = 0; i < saltpvdTables.size(); ++i) {
1624 const SaltpvdTable& saltpvdTable = saltpvdTables.getTable<SaltpvdTable>(i);
1625 saltpVdTable_[i].setXYContainers(saltpvdTable.getDepthColumn(), saltpvdTable.getSaltpColumn());
1626
1627 const auto& cells = reg.cells(i);
1628 for (const auto& cell : cells) {
1629 const double depth = cellCenterDepth_[cell];
1630 this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth, /*extrapolate=*/true);
1631 }
1632 }
1633}
1634
1635
1636template<class FluidSystem,
1637 class Grid,
1638 class GridView,
1639 class ElementMapper,
1640 class CartesianIndexMapper>
1641void InitialStateComputer<FluidSystem,
1642 Grid,
1643 GridView,
1644 ElementMapper,
1645 CartesianIndexMapper>::
1646updateCellProps_(const GridView& gridView,
1647 const NumericalAquifers& aquifer)
1648{
1649 ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
1650 int numElements = gridView.size(/*codim=*/0);
1651 cellCenterDepth_.resize(numElements);
1652 cellZSpan_.resize(numElements);
1653 cellZMinMax_.resize(numElements);
1654
1655 auto elemIt = gridView.template begin</*codim=*/0>();
1656 const auto& elemEndIt = gridView.template end</*codim=*/0>();
1657 const auto num_aqu_cells = aquifer.allAquiferCells();
1658 for (; elemIt != elemEndIt; ++elemIt) {
1659 const Element& element = *elemIt;
1660 const unsigned int elemIdx = elemMapper.index(element);
1661 cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element);
1662 const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
1663 cellZSpan_[elemIdx] = Details::cellZSpan(element);
1664 cellZMinMax_[elemIdx] = Details::cellZMinMax(element);
1665 if (!num_aqu_cells.empty()) {
1666 const auto search = num_aqu_cells.find(cartIx);
1667 if (search != num_aqu_cells.end()) {
1668 const auto* aqu_cell = num_aqu_cells.at(cartIx);
1669 const double depth_change_num_aqu = aqu_cell->depth - cellCenterDepth_[elemIdx];
1670 cellCenterDepth_[elemIdx] += depth_change_num_aqu;
1671 cellZSpan_[elemIdx].first += depth_change_num_aqu;
1672 cellZSpan_[elemIdx].second += depth_change_num_aqu;
1673 cellZMinMax_[elemIdx].first += depth_change_num_aqu;
1674 cellZMinMax_[elemIdx].second += depth_change_num_aqu;
1675 }
1676 }
1677 }
1678}
1679
1680template<class FluidSystem,
1681 class Grid,
1682 class GridView,
1683 class ElementMapper,
1684 class CartesianIndexMapper>
1685void InitialStateComputer<FluidSystem,
1686 Grid,
1687 GridView,
1688 ElementMapper,
1689 CartesianIndexMapper>::
1690applyNumericalAquifers_(const GridView& gridView,
1691 const NumericalAquifers& aquifer,
1692 const bool co2store_or_h2store)
1693{
1694 const auto num_aqu_cells = aquifer.allAquiferCells();
1695 if (num_aqu_cells.empty()) return;
1696
1697 // Check if water phase is active, or in the case of CO2STORE and H2STORE, water is modelled as oil phase
1698 bool oil_as_brine = co2store_or_h2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1699 const auto watPos = oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx;
1700 if (!FluidSystem::phaseIsActive(watPos)){
1701 throw std::logic_error { "Water phase has to be active for numerical aquifer case" };
1702 }
1703
1704 ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
1705 auto elemIt = gridView.template begin</*codim=*/0>();
1706 const auto& elemEndIt = gridView.template end</*codim=*/0>();
1707 const auto oilPos = FluidSystem::oilPhaseIdx;
1708 const auto gasPos = FluidSystem::gasPhaseIdx;
1709 for (; elemIt != elemEndIt; ++elemIt) {
1710 const Element& element = *elemIt;
1711 const unsigned int elemIdx = elemMapper.index(element);
1712 const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
1713 const auto search = num_aqu_cells.find(cartIx);
1714 if (search != num_aqu_cells.end()) {
1715 // numerical aquifer cells are filled with water initially
1716 this->sat_[watPos][elemIdx] = 1.;
1717
1718 if (!co2store_or_h2store && FluidSystem::phaseIsActive(oilPos)) {
1719 this->sat_[oilPos][elemIdx] = 0.;
1720 }
1721
1722 if (FluidSystem::phaseIsActive(gasPos)) {
1723 this->sat_[gasPos][elemIdx] = 0.;
1724 }
1725 const auto* aqu_cell = num_aqu_cells.at(cartIx);
1726 const auto msg = fmt::format("FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL "
1727 "AQUIFER {}, WATER SATURATION IS SET TO BE UNITY",
1728 aqu_cell->I+1, aqu_cell->J+1, aqu_cell->K+1, aqu_cell->aquifer_id);
1729 OpmLog::info(msg);
1730
1731 // if pressure is specified for numerical aquifers, we use these pressure values
1732 // for numerical aquifer cells
1733 if (aqu_cell->init_pressure) {
1734 const double pres = *(aqu_cell->init_pressure);
1735 this->pp_[watPos][elemIdx] = pres;
1736 if (FluidSystem::phaseIsActive(gasPos)) {
1737 this->pp_[gasPos][elemIdx] = pres;
1738 }
1739 if (FluidSystem::phaseIsActive(oilPos)) {
1740 this->pp_[oilPos][elemIdx] = pres;
1741 }
1742 }
1743 }
1744 }
1745}
1746
1747template<class FluidSystem,
1748 class Grid,
1749 class GridView,
1750 class ElementMapper,
1751 class CartesianIndexMapper>
1752template<class RMap>
1753void InitialStateComputer<FluidSystem,
1754 Grid,
1755 GridView,
1756 ElementMapper,
1757 CartesianIndexMapper>::
1758setRegionPvtIdx(const EclipseState& eclState, const RMap& reg)
1759{
1760 const auto& pvtnumData = eclState.fieldProps().get_int("PVTNUM");
1761
1762 for (const auto& r : reg.activeRegions()) {
1763 const auto& cells = reg.cells(r);
1764 regionPvtIdx_[r] = pvtnumData[*cells.begin()] - 1;
1765 }
1766}
1767
1768template<class FluidSystem,
1769 class Grid,
1770 class GridView,
1771 class ElementMapper,
1772 class CartesianIndexMapper>
1773template<class RMap, class MaterialLawManager, class Comm>
1774void InitialStateComputer<FluidSystem,
1775 Grid,
1776 GridView,
1777 ElementMapper,
1778 CartesianIndexMapper>::
1779calcPressSatRsRv(const RMap& reg,
1780 const std::vector<EquilRecord>& rec,
1781 MaterialLawManager& materialLawManager,
1782 const Comm& comm,
1783 const double grav)
1784{
1785 using PhaseSat = Details::PhaseSaturations<
1786 MaterialLawManager, FluidSystem, EquilReg, typename RMap::CellId
1787 >;
1788
1789 auto ptable = Details::PressureTable<FluidSystem, EquilReg>{ grav, this->num_pressure_points_ };
1790 auto psat = PhaseSat { materialLawManager, this->swatInit_ };
1791 auto vspan = std::array<double, 2>{};
1792
1793 std::vector<int> regionIsEmpty(rec.size(), 0);
1794 for (std::size_t r = 0; r < rec.size(); ++r) {
1795 const auto& cells = reg.cells(r);
1796
1797 Details::verticalExtent(cells, cellZMinMax_, comm, vspan);
1798
1799 const auto acc = rec[r].initializationTargetAccuracy();
1800 if (acc > 0) {
1801 throw std::runtime_error {
1802 "Cannot initialise model: Positive item 9 is not supported "
1803 "in EQUIL keyword, record " + std::to_string(r + 1)
1804 };
1805 }
1806
1807 if (cells.empty()) {
1808 regionIsEmpty[r] = 1;
1809 continue;
1810 }
1811
1812 const auto eqreg = EquilReg {
1813 rec[r], this->rsFunc_[r], this->rvFunc_[r], this->rvwFunc_[r], this->tempVdTable_[r], this->saltVdTable_[r], this->regionPvtIdx_[r]
1814 };
1815
1816 // Ensure gas/oil and oil/water contacts are within the span for the
1817 // phase pressure calculation.
1818 vspan[0] = std::min(vspan[0], std::min(eqreg.zgoc(), eqreg.zwoc()));
1819 vspan[1] = std::max(vspan[1], std::max(eqreg.zgoc(), eqreg.zwoc()));
1820
1821 ptable.equilibrate(eqreg, vspan);
1822
1823 if (acc == 0) {
1824 // Centre-point method
1825 this->equilibrateCellCentres(cells, eqreg, ptable, psat);
1826 }
1827 else if (acc < 0) {
1828 // Horizontal subdivision
1829 this->equilibrateHorizontal(cells, eqreg, -acc,
1830 ptable, psat);
1831 } else {
1832 // Horizontal subdivision with titled fault blocks
1833 // the simulator throw a few line above for the acc > 0 case
1834 // i.e. we should not reach here.
1835 assert(false);
1836 }
1837 }
1838 comm.min(regionIsEmpty.data(),regionIsEmpty.size());
1839 if (comm.rank() == 0) {
1840 for (std::size_t r = 0; r < rec.size(); ++r) {
1841 if (regionIsEmpty[r]) //region is empty on all partitions
1842 OpmLog::warning("Equilibration region " + std::to_string(r + 1)
1843 + " has no active cells");
1844 }
1845 }
1846}
1847
1848template<class FluidSystem,
1849 class Grid,
1850 class GridView,
1851 class ElementMapper,
1852 class CartesianIndexMapper>
1853template<class CellRange, class EquilibrationMethod>
1854void InitialStateComputer<FluidSystem,
1855 Grid,
1856 GridView,
1857 ElementMapper,
1858 CartesianIndexMapper>::
1859cellLoop(const CellRange& cells,
1860 EquilibrationMethod&& eqmethod)
1861{
1862 const auto oilPos = FluidSystem::oilPhaseIdx;
1863 const auto gasPos = FluidSystem::gasPhaseIdx;
1864 const auto watPos = FluidSystem::waterPhaseIdx;
1865
1866 const auto oilActive = FluidSystem::phaseIsActive(oilPos);
1867 const auto gasActive = FluidSystem::phaseIsActive(gasPos);
1868 const auto watActive = FluidSystem::phaseIsActive(watPos);
1869
1870 auto pressures = Details::PhaseQuantityValue{};
1871 auto saturations = Details::PhaseQuantityValue{};
1872 auto Rs = 0.0;
1873 auto Rv = 0.0;
1874 auto Rvw = 0.0;
1875
1876 for (const auto& cell : cells) {
1877 eqmethod(cell, pressures, saturations, Rs, Rv, Rvw);
1878
1879 if (oilActive) {
1880 this->pp_ [oilPos][cell] = pressures.oil;
1881 this->sat_[oilPos][cell] = saturations.oil;
1882 }
1883
1884 if (gasActive) {
1885 this->pp_ [gasPos][cell] = pressures.gas;
1886 this->sat_[gasPos][cell] = saturations.gas;
1887 }
1888
1889 if (watActive) {
1890 this->pp_ [watPos][cell] = pressures.water;
1891 this->sat_[watPos][cell] = saturations.water;
1892 }
1893
1894 if (oilActive && gasActive) {
1895 this->rs_[cell] = Rs;
1896 this->rv_[cell] = Rv;
1897 }
1898
1899 if (watActive && gasActive) {
1900 this->rvw_[cell] = Rvw;
1901 }
1902 }
1903}
1904
1905template<class FluidSystem,
1906 class Grid,
1907 class GridView,
1908 class ElementMapper,
1909 class CartesianIndexMapper>
1910template<class CellRange, class PressTable, class PhaseSat>
1911void InitialStateComputer<FluidSystem,
1912 Grid,
1913 GridView,
1914 ElementMapper,
1915 CartesianIndexMapper>::
1916equilibrateCellCentres(const CellRange& cells,
1917 const EquilReg& eqreg,
1918 const PressTable& ptable,
1919 PhaseSat& psat)
1920{
1921 using CellPos = typename PhaseSat::Position;
1922 using CellID = std::remove_cv_t<std::remove_reference_t<
1923 decltype(std::declval<CellPos>().cell)>>;
1924 this->cellLoop(cells, [this, &eqreg, &ptable, &psat]
1925 (const CellID cell,
1926 Details::PhaseQuantityValue& pressures,
1927 Details::PhaseQuantityValue& saturations,
1928 double& Rs,
1929 double& Rv,
1930 double& Rvw) -> void
1931 {
1932 const auto pos = CellPos {
1933 cell, cellCenterDepth_[cell]
1934 };
1935
1936 saturations = psat.deriveSaturations(pos, eqreg, ptable);
1937 pressures = psat.correctedPhasePressures();
1938
1939 const auto temp = this->temperature_[cell];
1940
1941 Rs = eqreg.dissolutionCalculator()
1942 (pos.depth, pressures.oil, temp, saturations.gas);
1943
1944 Rv = eqreg.evaporationCalculator()
1945 (pos.depth, pressures.gas, temp, saturations.oil);
1946
1947 Rvw = eqreg.waterEvaporationCalculator()
1948 (pos.depth, pressures.gas, temp, saturations.water);
1949 });
1950}
1951
1952template<class FluidSystem,
1953 class Grid,
1954 class GridView,
1955 class ElementMapper,
1956 class CartesianIndexMapper>
1957template<class CellRange, class PressTable, class PhaseSat>
1958void InitialStateComputer<FluidSystem,
1959 Grid,
1960 GridView,
1961 ElementMapper,
1962 CartesianIndexMapper>::
1963equilibrateHorizontal(const CellRange& cells,
1964 const EquilReg& eqreg,
1965 const int acc,
1966 const PressTable& ptable,
1967 PhaseSat& psat)
1968{
1969 using CellPos = typename PhaseSat::Position;
1970 using CellID = std::remove_cv_t<std::remove_reference_t<
1971 decltype(std::declval<CellPos>().cell)>>;
1972
1973 this->cellLoop(cells, [this, acc, &eqreg, &ptable, &psat]
1974 (const CellID cell,
1975 Details::PhaseQuantityValue& pressures,
1976 Details::PhaseQuantityValue& saturations,
1977 double& Rs,
1978 double& Rv,
1979 double& Rvw) -> void
1980 {
1981 pressures .reset();
1982 saturations.reset();
1983
1984 auto totfrac = 0.0;
1985 for (const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) {
1986 const auto pos = CellPos { cell, depth };
1987
1988 saturations.axpy(psat.deriveSaturations(pos, eqreg, ptable), frac);
1989 pressures .axpy(psat.correctedPhasePressures(), frac);
1990
1991 totfrac += frac;
1992 }
1993
1994 if (totfrac > 0.) {
1995 saturations /= totfrac;
1996 pressures /= totfrac;
1997 } else {
1998 // Fall back to centre point method for zero-thickness cells.
1999 const auto pos = CellPos {
2000 cell, cellCenterDepth_[cell]
2001 };
2002
2003 saturations = psat.deriveSaturations(pos, eqreg, ptable);
2004 pressures = psat.correctedPhasePressures();
2005 }
2006
2007 const auto temp = this->temperature_[cell];
2008 const auto cz = cellCenterDepth_[cell];
2009
2010 Rs = eqreg.dissolutionCalculator()
2011 (cz, pressures.oil, temp, saturations.gas);
2012
2013 Rv = eqreg.evaporationCalculator()
2014 (cz, pressures.gas, temp, saturations.oil);
2015
2016 Rvw = eqreg.waterEvaporationCalculator()
2017 (cz, pressures.gas, temp, saturations.water);
2018 });
2019}
2020
2021}
2022} // namespace EQUIL
2023} // namespace Opm
2024
2025#endif // OPM_INIT_STATE_EQUIL_IMPL_HPP
Auxiliary routines that to solve the ODEs that emerge from the hydrostatic equilibrium problem.
Routines that actually solve the ODEs that emerge from the hydrostatic equilibrium problem.
Calculator for phase saturations.
Definition InitStateEquil.hpp:376
const PhaseQuantityValue & deriveSaturations(const Position &x, const Region &reg, const PTable &ptable)
Calculate phase saturations at particular point of the simulation model geometry.
Definition InitStateEquil_impl.hpp:574
PhaseSaturations(MaterialLawManager &matLawMgr, const std::vector< double > &swatInit)
Constructor.
Definition InitStateEquil_impl.hpp:550
Definition InitStateEquil.hpp:159
PressureTable & operator=(const PressureTable &rhs)
Assignment operator.
Definition InitStateEquil_impl.hpp:990
double gas(const double depth) const
Evaluate gas phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1057
bool waterActive() const
Predicate for whether or not water is an active phase.
Definition InitStateEquil_impl.hpp:1041
bool gasActive() const
Predicate for whether or not gas is an active phase.
Definition InitStateEquil_impl.hpp:1034
bool oilActive() const
Predicate for whether or not oil is an active phase.
Definition InitStateEquil_impl.hpp:1027
PressureTable(const double gravity, const int samplePoints=2000)
Constructor.
Definition InitStateEquil_impl.hpp:960
double water(const double depth) const
Evaluate water phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1067
double oil(const double depth) const
Evaluate oil phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1048
bool water(const PhaseUsage &pu)
Active water predicate.
Definition RegionAttributeHelpers.hpp:308
bool oil(const PhaseUsage &pu)
Active oil predicate.
Definition RegionAttributeHelpers.hpp:321
bool gas(const PhaseUsage &pu)
Active gas predicate.
Definition RegionAttributeHelpers.hpp:334
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27
Simple set of per-phase (named by primary component) quantities.
Definition InitStateEquil.hpp:328
Evaluation point within a model geometry.
Definition InitStateEquil.hpp:381