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, 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,
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>
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>
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>
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>
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
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 }
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 }
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 }
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 }
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{
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{
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{
981 (this->matLawMgr_, static_cast<int>(phasePos),
982 this->evalPt_.position->cell, pc, isincr);
983}
984
985template<class FluidSystem, class Region>
987PressureTable(const Scalar gravity,
988 const int samplePoints)
989 : gravity_(gravity)
990 , nsample_(samplePoints)
991{
992}
993
994template <class FluidSystem, class Region>
997 : gravity_(rhs.gravity_)
998 , nsample_(rhs.nsample_)
999{
1000 this->copyInPointers(rhs);
1001}
1002
1003template <class FluidSystem, class Region>
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>
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>
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>
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>
1054oilActive() const
1055{
1056 return FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1057}
1058
1059template <class FluidSystem, class Region>
1061gasActive() const
1062{
1063 return FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
1064}
1065
1066template <class FluidSystem, class Region>
1068waterActive() const
1069{
1070 return FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1071}
1072
1073template <class FluidSystem, class Region>
1074typename FluidSystem::Scalar
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
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
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>
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,
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,
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,
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,
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,
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) {
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,
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,
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,
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.) {
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
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:386
const PhaseQuantityValue< Scalar > & 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:596
PhaseSaturations(MaterialLawManager &matLawMgr, const std::vector< Scalar > &swatInit)
Constructor.
Definition InitStateEquil_impl.hpp:572
Definition InitStateEquil.hpp:167
PressureTable & operator=(const PressureTable &rhs)
Assignment operator.
Definition InitStateEquil_impl.hpp:1017
Scalar water(const Scalar depth) const
Evaluate water phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1097
Scalar gas(const Scalar depth) const
Evaluate gas phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1086
bool waterActive() const
Predicate for whether or not water is an active phase.
Definition InitStateEquil_impl.hpp:1068
bool gasActive() const
Predicate for whether or not gas is an active phase.
Definition InitStateEquil_impl.hpp:1061
Scalar oil(const Scalar depth) const
Evaluate oil phase pressure at specified depth.
Definition InitStateEquil_impl.hpp:1076
bool oilActive() const
Predicate for whether or not oil is an active phase.
Definition InitStateEquil_impl.hpp:1054
PressureTable(const Scalar gravity, const int samplePoints=2000)
Constructor.
Definition InitStateEquil_impl.hpp:987
bool water(const PhaseUsage &pu)
Active water predicate.
Definition RegionAttributeHelpers.hpp:309
bool oil(const PhaseUsage &pu)
Active oil predicate.
Definition RegionAttributeHelpers.hpp:322
bool gas(const PhaseUsage &pu)
Active gas predicate.
Definition RegionAttributeHelpers.hpp:335
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Simple set of per-phase (named by primary component) quantities.
Definition InitStateEquil.hpp:338
Evaluation point within a model geometry.
Definition InitStateEquil.hpp:392