From 69b1a5ca5a9453ea6055b30251f8d5aeb03db81f Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 2 Jan 2018 12:43:56 +0100 Subject: [PATCH] equil init: formating fixes to make it more consistent with the rest of ewoms in particular, this removes excessive whitespace usage. --- ebos/equil/EquilibrationHelpers.hpp | 186 +++++++++-------- ebos/equil/RegionMapping.hpp | 25 +-- ebos/equil/initStateEquil.hpp | 306 ++++++++++++++-------------- 3 files changed, 258 insertions(+), 259 deletions(-) diff --git a/ebos/equil/EquilibrationHelpers.hpp b/ebos/equil/EquilibrationHelpers.hpp index 0db905bac..ebc0ce21b 100644 --- a/ebos/equil/EquilibrationHelpers.hpp +++ b/ebos/equil/EquilibrationHelpers.hpp @@ -61,14 +61,14 @@ template struct PcEq; - template - inline double satFromPc(const MaterialLawManager& materialLawManager, + template + double satFromPc(const MaterialLawManager& materialLawManager, const int phase, const int cell, const double targetPc, const bool increasing = false) template - inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager, + double satFromSumOfPcs(const MaterialLawManager& materialLawManager, const int phase1, const int phase2, const int cell, @@ -78,8 +78,7 @@ ---- end of synopsis of EquilibrationHelpers.hpp ---- */ -namespace Ewoms -{ +namespace Ewoms { /** * Types and routines that collectively implement a basic * ECLIPSE-style equilibration-based initialisation scheme. @@ -143,7 +142,8 @@ public: /** * Type that implements "no phase mixing" policy. */ -class NoMixing : public RsFunction { +class NoMixing : public RsFunction +{ public: /** * Function call. @@ -178,7 +178,8 @@ public: * typically taken from keyword 'RSVD'. */ template -class RsVD : public RsFunction { +class RsVD : public RsFunction +{ public: /** * Constructor. @@ -192,8 +193,7 @@ public: const std::vector& rs) : pvtRegionIdx_(pvtRegionIdx) , rsVsDepth_(depth, rs) - { - } + {} /** * Function call. @@ -210,15 +210,15 @@ public: * \return Dissolved gas-oil ratio (RS) at depth @c * depth and pressure @c press. */ - double - operator()(const double depth, - const double press, - const double temp, - const double satGas = 0.0) const + double operator()(const double depth, + const double press, + const double temp, + const double satGas = 0.0) const { if (satGas > 0.0) { return satRs(press, temp); - } else { + } + else { if (rsVsDepth_.xMin() > depth) return rsVsDepth_.valueAt(0); else if (rsVsDepth_.xMax() < depth) @@ -246,7 +246,8 @@ private: * typically taken from keyword 'RVVD'. */ template -class RvVD : public RsFunction { +class RvVD : public RsFunction +{ public: /** * Constructor. @@ -260,8 +261,7 @@ public: const std::vector& rv) : pvtRegionIdx_(pvtRegionIdx) , rvVsDepth_(depth, rv) - { - } + {} /** * Function call. @@ -278,15 +278,15 @@ public: * \return Vaporized oil-gas ratio (RV) at depth @c * depth and pressure @c press. */ - double - operator()(const double depth, - const double press, - const double temp, - const double satOil = 0.0 ) const + double operator()(const double depth, + const double press, + const double temp, + const double satOil = 0.0) const { if (std::abs(satOil) > 1e-16) { return satRv(press, temp); - } else { + } + else { if (rvVsDepth_.xMin() > depth) return rvVsDepth_.valueAt(0); else if (rvVsDepth_.xMax() < depth) @@ -323,7 +323,8 @@ private: * contact, and decreasing above the contact. */ template -class RsSatAtContact : public RsFunction { +class RsSatAtContact : public RsFunction +{ public: /** * Constructor. @@ -353,15 +354,15 @@ public: * \return Dissolved gas-oil ratio (RS) at depth @c * depth and pressure @c press. */ - double - operator()(const double /* depth */, - const double press, - const double temp, - const double satGas = 0.0) const + double operator()(const double /* depth */, + const double press, + const double temp, + const double satGas = 0.0) const { if (satGas > 0.0) { return satRs(press, temp); - } else { + } + else { return std::min(satRs(press, temp), rsSatContact_); } } @@ -392,7 +393,8 @@ private: * contact, and decreasing above the contact. */ template -class RvSatAtContact : public RsFunction { +class RvSatAtContact : public RsFunction +{ public: /** * Constructor. @@ -422,15 +424,15 @@ public: * \return Dissolved oil-gas ratio (RV) at depth @c * depth and pressure @c press. */ - double - operator()(const double /*depth*/, - const double press, - const double temp, - const double satOil = 0.0) const + double operator()(const double /*depth*/, + const double press, + const double temp, + const double satOil = 0.0) const { if (satOil > 0.0) { return satRv(press, temp); - } else { + } + else { return std::min(satRv(press, temp), rvSatContact_); } } @@ -460,13 +462,14 @@ private: * declared as * * std::vector - * operator()(const double press, - * const std::vector& svol ) + * operator()(const double press, + * const std::vector& svol) * * that calculates the phase densities of all phases in @c * svol at fluid pressure @c press. */ -class EquilReg { +class EquilReg +{ public: /** * Constructor. @@ -484,8 +487,7 @@ public: , rs_ (rs) , rv_ (rv) , pvtIdx_ (pvtIdx) - { - } + {} /** * Type of dissolved gas-oil ratio calculator. @@ -575,9 +577,8 @@ struct PcEq phase_(phase), cell_(cell), targetPc_(targetPc) - { + {} - } double operator()(double s) const { const auto& matParams = materialLawManager_.materialLawParams(cell_); @@ -602,53 +603,47 @@ private: }; template -double minSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) { +double minSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) +{ const auto& scaledDrainageInfo = materialLawManager.oilWaterScaledEpsInfoDrainage(cell); // Find minimum and maximum saturations. switch(phase) { - case FluidSystem::waterPhaseIdx : - { + case FluidSystem::waterPhaseIdx: return scaledDrainageInfo.Swl; - } - case FluidSystem::gasPhaseIdx : - { + + case FluidSystem::gasPhaseIdx: return scaledDrainageInfo.Sgl; - } - case FluidSystem::oilPhaseIdx : - { + + case FluidSystem::oilPhaseIdx: OPM_THROW(std::runtime_error, "Min saturation not implemented for oil phase."); - break; - } - default: OPM_THROW(std::runtime_error, "Unknown phaseIdx ."); + + default: + OPM_THROW(std::runtime_error, "Unknown phaseIdx ."); } return -1.0; } template -double maxSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) { +double maxSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) +{ const auto& scaledDrainageInfo = materialLawManager.oilWaterScaledEpsInfoDrainage(cell); // Find minimum and maximum saturations. switch(phase) { - case FluidSystem::waterPhaseIdx : - { + case FluidSystem::waterPhaseIdx: return scaledDrainageInfo.Swu; - break; - } - case FluidSystem::gasPhaseIdx : - { + + case FluidSystem::gasPhaseIdx: return scaledDrainageInfo.Sgu; - break; - } - case FluidSystem::oilPhaseIdx : - { + + case FluidSystem::oilPhaseIdx: OPM_THROW(std::runtime_error, "Max saturation not implemented for oil phase."); - break; - } - default: OPM_THROW(std::runtime_error, "Unknown phaseIdx ."); + + default: + OPM_THROW(std::runtime_error, "Unknown phaseIdx ."); } return -1.0; } @@ -656,12 +651,12 @@ double maxSaturations(const MaterialLawManager& materialLawManager, const int ph /// Compute saturation of some phase corresponding to a given /// capillary pressure. -template -inline double satFromPc(const MaterialLawManager& materialLawManager, - const int phase, - const int cell, - const double targetPc, - const bool increasing = false) +template +double satFromPc(const MaterialLawManager& materialLawManager, + const int phase, + const int cell, + const double targetPc, + const bool increasing = false) { // Find minimum and maximum saturations. double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); @@ -730,8 +725,8 @@ struct PcEqSum phase2_(phase2), cell_(cell), targetPc_(targetPc) - { - } + {} + double operator()(double s) const { const auto& matParams = materialLawManager_.materialLawParams(cell_); @@ -767,11 +762,11 @@ private: /// capillary pressure, where the capillary pressure function /// is given as a sum of two other functions. template -inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager, - const int phase1, - const int phase2, - const int cell, - const double targetPc) +double satFromSumOfPcs(const MaterialLawManager& materialLawManager, + const int phase1, + const int phase2, + const int cell, + const double targetPc) { // Find minimum and maximum saturations. double s0 = minSaturations(materialLawManager, phase1, cell); @@ -824,19 +819,20 @@ inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager, /// Compute saturation from depth. Used for constant capillary pressure function template -inline double satFromDepth(const MaterialLawManager& materialLawManager, - const double cellDepth, - const double contactDepth, - const int phase, - const int cell, - const bool increasing = false) +double satFromDepth(const MaterialLawManager& materialLawManager, + const double cellDepth, + const double contactDepth, + const int phase, + const int cell, + const bool increasing = false) { const double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); const double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); - if (cellDepth < contactDepth){ + if (cellDepth < contactDepth) { return s0; - } else { + } + else { return s1; } @@ -844,9 +840,9 @@ inline double satFromDepth(const MaterialLawManager& materialLawManager, /// Return true if capillary pressure function is constant template -inline bool isConstPc(const MaterialLawManager& materialLawManager, - const int phase, - const int cell) +bool isConstPc(const MaterialLawManager& materialLawManager, + const int phase, + const int cell) { // Create the equation f(s) = pc(s); const PcEq f(materialLawManager, phase, cell, 0); diff --git a/ebos/equil/RegionMapping.hpp b/ebos/equil/RegionMapping.hpp index 857eade49..77cca4681 100644 --- a/ebos/equil/RegionMapping.hpp +++ b/ebos/equil/RegionMapping.hpp @@ -26,8 +26,7 @@ #include #include -namespace Ewoms -{ +namespace Ewoms { /** * Forward and reverse mappings between cells and @@ -40,8 +39,9 @@ namespace Ewoms * 'valueType', 'size_type', and * 'const_iterator'. */ -template < class Region = std::vector > -class RegionMapping { +template > +class RegionMapping +{ public: /** * Constructor. @@ -80,7 +80,8 @@ public: typedef CellIter iterator; typedef CellIter const_iterator; - Range() {}; + Range() + {}; Range(const CellIter& beg, const CellIter& en) : begin_(beg) @@ -116,8 +117,8 @@ public: * \param[in] c Active cell * \return Region to which @c c belongs. */ - RegionId - region(const CellId c) const { return reg_[c]; } + RegionId region(const CellId c) const + { return reg_[c]; } const std::vector& activeRegions() const @@ -133,8 +134,8 @@ public: * \return Range of active cells in region @c r. Empty if @c r is * not an active region. */ - Range - cells(const RegionId r) const { + Range cells(const RegionId r) const + { const auto id = rev_.binid.find(r); if (id == rev_.binid.end()) { @@ -192,7 +193,7 @@ private: for (decltype(p.size()) i = 1, sz = p.size(); i < sz; ++i) { p[0] += p[i]; - p[i] = p[0] - p[i]; + p[i] = p[0] - p[i]; } assert (p[0] == static_cast(reg.size())); @@ -201,8 +202,8 @@ private: { CellId i = 0; for (const auto& r : reg) { - auto& pos = p[ binid[r] + 1 ]; - c[ pos++ ] = i++; + auto& pos = p[binid[r] + 1]; + c[pos++] = i++; } } diff --git a/ebos/equil/initStateEquil.hpp b/ebos/equil/initStateEquil.hpp index 25600eb18..e883ea5fe 100644 --- a/ebos/equil/initStateEquil.hpp +++ b/ebos/equil/initStateEquil.hpp @@ -70,14 +70,14 @@ namespace Details { template class RK4IVP { public: - RK4IVP(const RHS& f , + RK4IVP(const RHS& f, const std::array& span, - const double y0 , - const int N ) + const double y0, + const int N) : N_(N) , span_(span) { - const double h = stepsize(); + const double h = stepsize(); const double h2 = h / 2; const double h6 = h / 6; @@ -88,13 +88,13 @@ public: f_.push_back(f(span_[0], y0)); for (int i = 0; i < N; ++i) { - const double x = span_[0] + i*h; - const double y = y_.back(); + const double x = span_[0] + i*h; + const double y = y_.back(); const double k1 = f_[i]; const double k2 = f(x + h2, y + h2*k1); const double k3 = f(x + h2, y + h2*k2); - const double k4 = f(x + h , y + h*k3); + const double k4 = f(x + h, y + h*k3); y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4)); f_.push_back(f(x + h, y_.back())); @@ -109,7 +109,7 @@ public: // Dense output (O(h**3)) according to Shampine // (Hermite interpolation) const double h = stepsize(); - int i = (x - span_[0]) / h; + int i = (x - span_[0]) / h; const double t = (x - (span_[0] + i*h)) / h; // Crude handling of evaluation point outside "span_"; @@ -128,7 +128,7 @@ public: } private: - int N_; + int N_; std::array span_; std::vector y_; std::vector f_; @@ -139,16 +139,16 @@ private: namespace PhasePressODE { template -class Water { +class Water +{ public: - Water(const double temp, - const int pvtRegionIdx, - const double normGrav) + Water(const double temp, + const int pvtRegionIdx, + const double normGrav) : temp_(temp) , pvtRegionIdx_(pvtRegionIdx) , g_(normGrav) - { - } + {} double operator()(const double /* depth */, @@ -158,9 +158,9 @@ public: } private: - const double temp_; - const int pvtRegionIdx_; - const double g_; + const double temp_; + const int pvtRegionIdx_; + const double g_; double density(const double press) const @@ -172,18 +172,18 @@ private: }; template -class Oil { +class Oil +{ public: - Oil(const double temp, - const RS& rs, - const int pvtRegionIdx, - const double normGrav) + Oil(const double temp, + const RS& rs, + const int pvtRegionIdx, + const double normGrav) : temp_(temp) , rs_(rs) , pvtRegionIdx_(pvtRegionIdx) , g_(normGrav) - { - } + {} double operator()(const double depth, @@ -193,10 +193,10 @@ public: } private: - const double temp_; - const RS& rs_; - const int pvtRegionIdx_; - const double g_; + const double temp_; + const RS& rs_; + const int pvtRegionIdx_; + const double g_; double density(const double depth, @@ -204,9 +204,10 @@ private: { double rs = rs_(depth, press, temp_); double bOil = 0.0; - if ( !FluidSystem::enableDissolvedGas() || rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press) ) { + if (!FluidSystem::enableDissolvedGas() || rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press)) { bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - } else { + } + else { bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rs); } double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); @@ -219,18 +220,18 @@ private: }; template -class Gas { +class Gas +{ public: - Gas(const double temp, - const RV& rv, - const int pvtRegionIdx, - const double normGrav) + Gas(const double temp, + const RV& rv, + const int pvtRegionIdx, + const double normGrav) : temp_(temp) , rv_(rv) , pvtRegionIdx_(pvtRegionIdx) , g_(normGrav) - { - } + {} double operator()(const double depth, @@ -240,10 +241,10 @@ public: } private: - const double temp_; - const RV& rv_; - const int pvtRegionIdx_; - const double g_; + const double temp_; + const RV& rv_; + const int pvtRegionIdx_; + const double g_; double density(const double depth, @@ -251,9 +252,10 @@ private: { double rv = rv_(depth, press, temp_); double bGas = 0.0; - if ( !FluidSystem::enableVaporizedOil() || rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press) ) { + if (!FluidSystem::enableVaporizedOil() || rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press)) { bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - } else { + } + else { bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv); } double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); @@ -271,12 +273,11 @@ namespace PhasePressure { template -void -assign(const Grid& grid , - const std::array& f , - const double split, - const CellRange& cells, - std::vector& p ) +void assign(const Grid& grid, + const std::array& f , + const double split, + const CellRange& cells, + std::vector& p) { enum { up = 0, down = 1 }; @@ -297,14 +298,13 @@ template -void -water(const Grid& grid , - const Region& reg , - const std::array& span , - const double grav , - double& poWoc, - const CellRange& cells , - std::vector& press ) +void water(const Grid& grid, + const Region& reg, + const std::array& span , + const double grav, + double& poWoc, + const CellRange& cells, + std::vector& press) { using PhasePressODE::Water; typedef Water ODE; @@ -317,12 +317,13 @@ water(const Grid& grid , if (reg.datum() > reg.zwoc()) {//Datum in water zone z0 = reg.datum(); p0 = reg.pressure(); - } else { + } + else { z0 = reg.zwoc(); p0 = poWoc - reg.pcowWoc(); // Water pressure at contact } - std::array up = {{ z0, span[0] }}; + std::array up = {{ z0, span[0] }}; std::array down = {{ z0, span[1] }}; typedef Details::RK4IVP WPress; @@ -346,15 +347,14 @@ template -void -oil(const Grid& grid , - const Region& reg , - const std::array& span , - const double grav , - const CellRange& cells , - std::vector& press , - double& poWoc, - double& poGoc) +void oil(const Grid& grid, + const Region& reg, + const std::array& span , + const double grav, + const CellRange& cells, + std::vector& press, + double& poWoc, + double& poGoc) { using PhasePressODE::Oil; typedef Oil ODE; @@ -368,15 +368,17 @@ oil(const Grid& grid , if (reg.datum() > reg.zwoc()) {//Datum in water zone, poWoc given z0 = reg.zwoc(); p0 = poWoc; - } else if (reg.datum() < reg.zgoc()) {//Datum in gas zone, poGoc given + } + else if (reg.datum() < reg.zgoc()) {//Datum in gas zone, poGoc given z0 = reg.zgoc(); p0 = poGoc; - } else { //Datum in oil zone + } + else { //Datum in oil zone z0 = reg.datum(); p0 = reg.pressure(); } - std::array up = {{ z0, span[0] }}; + std::array up = {{ z0, span[0] }}; std::array down = {{ z0, span[1] }}; typedef Details::RK4IVP OPress; @@ -405,14 +407,13 @@ template -void -gas(const Grid& grid , - const Region& reg , - const std::array& span , - const double grav , - double& poGoc, - const CellRange& cells , - std::vector& press ) +void gas(const Grid& grid, + const Region& reg, + const std::array& span , + const double grav, + double& poGoc, + const CellRange& cells, + std::vector& press) { using PhasePressODE::Gas; typedef Gas ODE; @@ -426,12 +427,13 @@ gas(const Grid& grid , if (reg.datum() < reg.zgoc()) {//Datum in gas zone z0 = reg.datum(); p0 = reg.pressure(); - } else { + } + else { z0 = reg.zgoc(); p0 = poGoc + reg.pcgoGoc(); // Gas pressure at contact } - std::array up = {{ z0, span[0] }}; + std::array up = {{ z0, span[0] }}; std::array down = {{ z0, span[1] }}; typedef Details::RK4IVP GPress; @@ -456,13 +458,12 @@ template -void -equilibrateOWG(const Grid& grid, - const Region& reg, - const double grav, - const std::array& span, - const CellRange& cells, - std::vector< std::vector >& press) +void equilibrateOWG(const Grid& grid, + const Region& reg, + const double grav, + const std::array& span, + const CellRange& cells, + std::vector< std::vector >& press) { const bool water = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); @@ -477,53 +478,55 @@ equilibrateOWG(const Grid& grid, if (water) { PhasePressure::water(grid, reg, span, grav, poWoc, - cells, press[ waterpos ]); + cells, press[waterpos]); } if (oil) { PhasePressure::oil(grid, reg, span, grav, cells, - press[ oilpos ], poWoc, poGoc); + press[oilpos], poWoc, poGoc); } if (gas) { PhasePressure::gas(grid, reg, span, grav, poGoc, - cells, press[ gaspos ]); + cells, press[gaspos]); } - } else if (reg.datum() < reg.zgoc()) { // Datum in gas zone + } + else if (reg.datum() < reg.zgoc()) { // Datum in gas zone double poWoc = -1; double poGoc = -1; if (gas) { PhasePressure::gas(grid, reg, span, grav, poGoc, - cells, press[ gaspos ]); + cells, press[gaspos]); } if (oil) { PhasePressure::oil(grid, reg, span, grav, cells, - press[ oilpos ], poWoc, poGoc); + press[oilpos], poWoc, poGoc); } if (water) { PhasePressure::water(grid, reg, span, grav, poWoc, - cells, press[ waterpos ]); + cells, press[waterpos]); } - } else { // Datum in oil zone + } + else { // Datum in oil zone double poWoc = -1; double poGoc = -1; if (oil) { PhasePressure::oil(grid, reg, span, grav, cells, - press[ oilpos ], poWoc, poGoc); + press[oilpos], poWoc, poGoc); } if (water) { PhasePressure::water(grid, reg, span, grav, poWoc, - cells, press[ waterpos ]); + cells, press[waterpos]); } if (gas) { PhasePressure::gas(grid, reg, span, grav, poGoc, - cells, press[ gaspos ]); + cells, press[gaspos]); } } } @@ -568,14 +571,14 @@ equilibrateOWG(const Grid& grid, * equilibration region. */ template -std::vector< std::vector > -phasePressures(const Grid& grid, - const Region& reg, - const CellRange& cells, - const double grav = Opm::unit::gravity) +std::vector< std::vector> +phasePressures(const Grid& grid, + const Region& reg, + const CellRange& cells, + const double grav = Opm::unit::gravity) { std::array span = - {{ std::numeric_limits::max() , + {{ std::numeric_limits::max(), -std::numeric_limits::max() }}; // Symm. about 0. int ncell = 0; @@ -642,9 +645,9 @@ template std::vector -temperature(const Grid& /* G */, - const Region& /* reg */, - const CellRange& cells) +temperature(const Grid& /* G */, + const Region& /* reg */, + const CellRange& cells) { // use the standard temperature for everything for now return std::vector(cells.size(), 273.15 + 20.0); @@ -685,10 +688,10 @@ temperature(const Grid& /* G */, * one saturation value per cell in the region. */ template -std::vector< std::vector > -phaseSaturations(const Grid& grid, - const Region& reg, - const CellRange& cells, +std::vector< std::vector> +phaseSaturations(const Grid& grid, + const Region& reg, + const CellRange& cells, MaterialLawManager& materialLawManager, const std::vector swatInit, std::vector< std::vector >& phasePressures) @@ -733,17 +736,18 @@ phaseSaturations(const Grid& grid, double sw = 0.0; if (water) { if (isConstPc(materialLawManager,FluidSystem::waterPhaseIdx, cell)){ - const double cellDepth = Opm::UgGridHelpers::cellCenterDepth(grid, - cell); + const double cellDepth = Opm::UgGridHelpers::cellCenterDepth(grid, + cell); sw = satFromDepth(materialLawManager,cellDepth,reg.zwoc(),waterpos,cell,false); phaseSaturations[waterpos][localIndex] = sw; } - else{ + else { const double pcov = phasePressures[oilpos][localIndex] - phasePressures[waterpos][localIndex]; if (swatInit.empty()) { // Invert Pc to find sw sw = satFromPc(materialLawManager, waterpos, cell, pcov); phaseSaturations[waterpos][localIndex] = sw; - } else { // Scale Pc to reflect imposed sw + } + else { // Scale Pc to reflect imposed sw sw = swatInit[cell]; sw = materialLawManager.applySwatinit(cell, pcov, sw); phaseSaturations[waterpos][localIndex] = sw; @@ -753,12 +757,12 @@ phaseSaturations(const Grid& grid, double sg = 0.0; if (gas) { if (isConstPc(materialLawManager,FluidSystem::gasPhaseIdx,cell)){ - const double cellDepth = Opm::UgGridHelpers::cellCenterDepth(grid, - cell); + const double cellDepth = Opm::UgGridHelpers::cellCenterDepth(grid, + cell); sg = satFromDepth(materialLawManager,cellDepth,reg.zgoc(),gaspos,cell,true); phaseSaturations[gaspos][localIndex] = sg; } - else{ + else { // Note that pcog is defined to be (pg - po), not (po - pg). const double pcog = phasePressures[gaspos][localIndex] - phasePressures[oilpos][localIndex]; const double increasing = true; // pcog(sg) expected to be increasing function @@ -782,7 +786,7 @@ phaseSaturations(const Grid& grid, sg = 1.0 - sw; phaseSaturations[waterpos][localIndex] = sw; phaseSaturations[gaspos][localIndex] = sg; - if ( water ) { + if (water) { fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw); } else { @@ -815,12 +819,13 @@ phaseSaturations(const Grid& grid, } fluidState.setSaturation(FluidSystem::oilPhaseIdx, so); - if (water && sw > scaledDrainageInfo.Swu-thresholdSat ) { + if (water && sw > scaledDrainageInfo.Swu-thresholdSat) { fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swu); MaterialLaw::capillaryPressures(pC, matParams, fluidState); double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx]; phasePressures[oilpos][localIndex] = phasePressures[waterpos][localIndex] + pcWat; - } else if (gas && sg > scaledDrainageInfo.Sgu-thresholdSat) { + } + else if (gas && sg > scaledDrainageInfo.Sgu-thresholdSat) { fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgu); MaterialLaw::capillaryPressures(pC, matParams, fluidState); double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx]; @@ -879,13 +884,12 @@ std::vector computeRs(const Grid& grid, } namespace DeckDependent { -inline -std::vector +inline std::vector getEquil(const Opm::EclipseState& state) { const auto& init = state.getInitConfig(); - if( !init.hasEquil() ) { + if(!init.hasEquil()) { OPM_THROW(std::domain_error, "Deck does not provide equilibration data."); } @@ -894,10 +898,9 @@ getEquil(const Opm::EclipseState& state) } template -inline std::vector equilnum(const Opm::EclipseState& eclipseState, - const Grid& grid ) + const Grid& grid) { std::vector eqlnum; if (eclipseState.get3DProperties().hasDeckIntGridProperty("EQLNUM")) { @@ -931,10 +934,9 @@ public: template InitialStateComputer(MaterialLawManager& materialLawManager, const Opm::EclipseState& eclipseState, - const Grid& grid , + const Grid& grid, const double grav = Opm::unit::gravity, - const bool applySwatInit = true - ) + const bool applySwatInit = true) : pp_(FluidSystem::numPhases, std::vector(grid.size(/*codim=*/0))), sat_(FluidSystem::numPhases, @@ -968,22 +970,22 @@ public: if (FluidSystem::enableDissolvedGas()) { const Opm::TableContainer& rsvdTables = tables.getRsvdTables(); for (size_t i = 0; i < rec.size(); ++i) { - if (eqlmap.cells(i).empty()) - { + if (eqlmap.cells(i).empty()) { rsFunc_.push_back(std::shared_ptr>()); continue; } const int pvtIdx = regionPvtIdx_[i]; if (!rec[i].liveOilInitConstantRs()) { - if (rsvdTables.size() <= 0 ) { + if (rsvdTables.size() <= 0) { OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table not available."); } const Opm::RsvdTable& rsvdTable = rsvdTables.getTable(i); std::vector depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy(); std::vector rsColumn = rsvdTable.getColumn("RS").vectorCopy(); rsFunc_.push_back(std::make_shared>(pvtIdx, - depthColumn , rsColumn)); - } else { + depthColumn, rsColumn)); + } + else { if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { OPM_THROW(std::runtime_error, "Cannot initialise: when no explicit RSVD table is given, \n" @@ -995,7 +997,8 @@ public: rsFunc_.push_back(std::make_shared>(pvtIdx, pContact, TContact)); } } - } else { + } + else { for (size_t i = 0; i < rec.size(); ++i) { rsFunc_.push_back(std::make_shared()); } @@ -1005,8 +1008,7 @@ public: if (FluidSystem::enableVaporizedOil()) { const Opm::TableContainer& rvvdTables = tables.getRvvdTables(); for (size_t i = 0; i < rec.size(); ++i) { - if (eqlmap.cells(i).empty()) - { + if (eqlmap.cells(i).empty()) { rvFunc_.push_back(std::shared_ptr>()); continue; } @@ -1020,9 +1022,10 @@ public: std::vector depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy(); std::vector rvColumn = rvvdTable.getColumn("RV").vectorCopy(); rvFunc_.push_back(std::make_shared>(pvtIdx, - depthColumn , rvColumn)); + depthColumn, rvColumn)); - } else { + } + else { if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { OPM_THROW(std::runtime_error, "Cannot initialise: when no explicit RVVD table is given, \n" @@ -1031,10 +1034,11 @@ public: } const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); const double TContact = 273.15 + 20; // standard temperature for now - rvFunc_.push_back(std::make_shared>(pvtIdx ,pContact, TContact)); + rvFunc_.push_back(std::make_shared>(pvtIdx,pContact, TContact)); } } - } else { + } + else { for (size_t i = 0; i < rec.size(); ++i) { rvFunc_.push_back(std::make_shared()); } @@ -1094,17 +1098,15 @@ private: } template - void - calcPressSatRsRv(const RMap& reg , - const std::vector< Opm::EquilRecord >& rec , - MaterialLawManager& materialLawManager, - const Grid& grid , - const double grav) + void calcPressSatRsRv(const RMap& reg, + const std::vector< Opm::EquilRecord >& rec, + MaterialLawManager& materialLawManager, + const Grid& grid, + const double grav) { for (const auto& r : reg.activeRegions()) { const auto& cells = reg.cells(r); - if (cells.empty()) - { + if (cells.empty()) { Opm::OpmLog::warning("Equilibration region " + std::to_string(r + 1) + " has no active cells"); continue; @@ -1143,7 +1145,7 @@ private: auto c = cells.begin(); const auto e = cells.end(); for (; c != e; ++c, ++s) { - destination[*c] = *s; + destination[*c] =*s; } }