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