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