From 4bef9259742f69bc5747279639b2d03e6271235e Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 12 Apr 2024 23:30:19 +0200 Subject: [PATCH] InitStateEquil: use Scalar type from FluidSystem --- opm/simulators/flow/equil/InitStateEquil.cpp | 53 +-- opm/simulators/flow/equil/InitStateEquil.hpp | 166 +++---- .../flow/equil/InitStateEquil_impl.hpp | 405 ++++++++++-------- tests/test_equil.cpp | 2 +- 4 files changed, 345 insertions(+), 281 deletions(-) diff --git a/opm/simulators/flow/equil/InitStateEquil.cpp b/opm/simulators/flow/equil/InitStateEquil.cpp index 53f51afd7..cae9c56ab 100644 --- a/opm/simulators/flow/equil/InitStateEquil.cpp +++ b/opm/simulators/flow/equil/InitStateEquil.cpp @@ -36,34 +36,36 @@ #endif namespace Opm { + +template +using MatLaw = EclMaterialLawManager>; + namespace EQUIL { namespace DeckDependent { -using MatLaw = EclMaterialLawManager>; - -#define INSTANCE_COMP(GridView, Mapper) \ - template class InitialStateComputer, \ +#define INSTANTIATE_COMP(T, GridView, Mapper) \ + template class InitialStateComputer, \ Dune::CpGrid, \ GridView, \ Mapper, \ Dune::CartesianIndexMapper>; \ - template InitialStateComputer, \ + template InitialStateComputer, \ Dune::CpGrid, \ GridView, \ Mapper, \ Dune::CartesianIndexMapper>::\ - InitialStateComputer(MatLaw&, \ - const EclipseState&, \ - const Dune::CpGrid&, \ - const GridView&, \ - const Dune::CartesianIndexMapper&, \ - const double, \ - const int, \ - const bool); + InitialStateComputer(MatLaw&, \ + const EclipseState&, \ + const Dune::CpGrid&, \ + const GridView&, \ + const Dune::CartesianIndexMapper&, \ + const T, \ + const int, \ + const bool); using GridView = Dune::GridView>; using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper; -INSTANCE_COMP(GridView, Mapper) +INSTANTIATE_COMP(double, GridView, Mapper) #if HAVE_DUNE_FEM #if DUNE_VERSION_GTE(DUNE_FEM, 2, 9) @@ -78,23 +80,26 @@ using GridViewFem = Dune::Fem::GridPart2GridViewImpl< false>>; #endif using MapperFem = Dune::MultipleCodimMultipleGeomTypeMapper; -INSTANCE_COMP(GridViewFem, MapperFem) + +INSTANTIATE_COMP(double, GridViewFem, MapperFem) + #endif // HAVE_DUNE_FEM } // namespace DeckDependent namespace Details { - template class PressureTable,EquilReg>; - template void verticalExtent(const std::vector&, - const std::vector>&, - const Parallel::Communication&, - std::array&); +#define INSTANTIATE_TYPE(T) \ + template class PressureTable,EquilReg>; \ + template void verticalExtent(const std::vector&, \ + const std::vector>&, \ + const Parallel::Communication&, \ + std::array&); \ + template class PhaseSaturations,BlackOilFluidSystem, \ + EquilReg,std::size_t>; \ + template std::pair cellZMinMax(const Dune::cpgrid::Entity<0>&); - using MatLaw = EclMaterialLawManager>; - template class PhaseSaturations, - EquilReg,std::size_t>; +INSTANTIATE_TYPE(double) - template std::pair cellZMinMax(const Dune::cpgrid::Entity<0>& element); } } // namespace EQUIL diff --git a/opm/simulators/flow/equil/InitStateEquil.hpp b/opm/simulators/flow/equil/InitStateEquil.hpp index ece0eedbf..f7ea44f6b 100644 --- a/opm/simulators/flow/equil/InitStateEquil.hpp +++ b/opm/simulators/flow/equil/InitStateEquil.hpp @@ -62,96 +62,102 @@ template class EquilReg; namespace Miscibility { template class RsFunction; } namespace Details { -template -class RK4IVP { +template +class RK4IVP +{ public: RK4IVP(const RHS& f, - const std::array& span, - const double y0, + const std::array& span, + const Scalar y0, const int N); - double - operator()(const double x) const; + Scalar operator()(const Scalar x) const; private: int N_; - std::array span_; - std::vector y_; - std::vector f_; + std::array span_; + std::vector y_; + std::vector f_; - double stepsize() const; + Scalar stepsize() const; }; namespace PhasePressODE { template class Water { -using TabulatedFunction = Tabulated1DFunction; + using Scalar = typename FluidSystem::Scalar; + using TabulatedFunction = Tabulated1DFunction; + public: Water(const TabulatedFunction& tempVdTable, const TabulatedFunction& saltVdTable, const int pvtRegionIdx, - const double normGrav); + const Scalar normGrav); - double operator()(const double depth, - const double press) const; + Scalar operator()(const Scalar depth, + const Scalar press) const; private: const TabulatedFunction& tempVdTable_; const TabulatedFunction& saltVdTable_; const int pvtRegionIdx_; - const double g_; + const Scalar g_; - double density(const double depth, - const double press) const; + Scalar density(const Scalar depth, + const Scalar press) const; }; template class Oil { -using TabulatedFunction = Tabulated1DFunction; + using Scalar = typename FluidSystem::Scalar; + using TabulatedFunction = Tabulated1DFunction; + public: Oil(const TabulatedFunction& tempVdTable, const RS& rs, const int pvtRegionIdx, - const double normGrav); + const Scalar normGrav); - double operator()(const double depth, - const double press) const; + Scalar operator()(const Scalar depth, + const Scalar press) const; private: const TabulatedFunction& tempVdTable_; const RS& rs_; const int pvtRegionIdx_; - const double g_; + const Scalar g_; - double density(const double depth, - const double press) const; + Scalar density(const Scalar depth, + const Scalar press) const; }; template class Gas { -using TabulatedFunction = Tabulated1DFunction; + using Scalar = typename FluidSystem::Scalar; + using TabulatedFunction = Tabulated1DFunction; + public: Gas(const TabulatedFunction& tempVdTable, const RV& rv, const RVW& rvw, const int pvtRegionIdx, - const double normGrav); + const Scalar normGrav); - double operator()(const double depth, - const double press) const; + Scalar operator()(const Scalar depth, + const Scalar press) const; private: const TabulatedFunction& tempVdTable_; const RV& rv_; const RVW& rvw_; const int pvtRegionIdx_; - const double g_; + const Scalar g_; - double density(const double depth, - const double press) const; + Scalar density(const Scalar depth, + const Scalar press) const; }; } // namespace PhasePressODE @@ -160,7 +166,8 @@ template class PressureTable { public: - using VSpan = std::array; + using Scalar = typename FluidSystem::Scalar; + using VSpan = std::array; /// Constructor /// @@ -170,7 +177,7 @@ public: /// /// \param[in] samplePoints Number of equally spaced depth sample points /// in each internal phase pressure table. - explicit PressureTable(const double gravity, + explicit PressureTable(const Scalar gravity, const int samplePoints = 2000); /// Copy constructor @@ -220,7 +227,7 @@ public: /// \endcode. /// /// \return Oil phase pressure at specified depth. - double oil(const double depth) const; + Scalar oil(const Scalar depth) const; /// Evaluate gas phase pressure at specified depth. /// @@ -229,7 +236,7 @@ public: /// \endcode. /// /// \return Gas phase pressure at specified depth. - double gas(const double depth) const; + Scalar gas(const Scalar depth) const; /// Evaluate water phase pressure at specified depth. /// @@ -238,7 +245,7 @@ public: /// \endcode. /// /// \return Water phase pressure at specified depth. - double water(const double depth) const; + Scalar water(const Scalar depth) const; private: template @@ -246,8 +253,8 @@ private: { public: struct InitCond { - double depth; - double pressure; + Scalar depth; + Scalar pressure; }; explicit PressureFunction(const ODE& ode, @@ -263,12 +270,12 @@ private: PressureFunction& operator=(PressureFunction&& rhs); - double value(const double depth) const; + Scalar value(const Scalar depth) const; private: enum Direction : std::size_t { Up, Down, NumDir }; - using Distribution = Details::RK4IVP; + using Distribution = Details::RK4IVP; using DistrPtr = std::unique_ptr; InitCond initial_; @@ -292,7 +299,7 @@ private: using Strategy = void (PressureTable::*) (const Region&, const VSpan&); - double gravity_; + Scalar gravity_; int nsample_; std::unique_ptr oil_{}; @@ -327,12 +334,13 @@ private: // =========================================================================== /// Simple set of per-phase (named by primary component) quantities. +template struct PhaseQuantityValue { - double oil{0.0}; - double gas{0.0}; - double water{0.0}; + Scalar oil{0.0}; + Scalar gas{0.0}; + Scalar water{0.0}; - PhaseQuantityValue& axpy(const PhaseQuantityValue& rhs, const double a) + PhaseQuantityValue& axpy(const PhaseQuantityValue& rhs, const Scalar a) { this->oil += a * rhs.oil; this->gas += a * rhs.gas; @@ -341,7 +349,7 @@ struct PhaseQuantityValue { return *this; } - PhaseQuantityValue& operator/=(const double x) + PhaseQuantityValue& operator/=(const Scalar x) { this->oil /= x; this->gas /= x; @@ -377,12 +385,13 @@ template & swatInit); + const std::vector& swatInit); /// Copy constructor. /// @@ -421,7 +430,7 @@ public: /// pertaining to the equilibration region \p reg. /// /// \return Set of phase saturation values defined at particular point. - const PhaseQuantityValue& + const PhaseQuantityValue& deriveSaturations(const Position& x, const Region& reg, const PTable& ptable); @@ -430,7 +439,7 @@ public: /// /// Values associated with evaluation point of previous call to \code /// deriveSaturations() \endcode. - const PhaseQuantityValue& correctedPhasePressures() const + const PhaseQuantityValue& correctedPhasePressures() const { return this->press_; } @@ -448,7 +457,7 @@ private: /// information needed to calculate the capillary pressure values from /// the current set of material laws. using FluidState = ::Opm:: - SimpleModularFluidState& swatInit_; + const std::vector& swatInit_; /// Evaluated phase saturations. - PhaseQuantityValue sat_; + PhaseQuantityValue sat_; /// Saturation-corrected phase pressure values. - PhaseQuantityValue press_; + PhaseQuantityValue press_; /// Current evaluation point. EvaluationPoint evalPt_; @@ -486,7 +495,7 @@ private: FluidState fluidState_; /// Evaluated capillary pressures from current set of material laws. - std::array matLawCapPress_; + std::array matLawCapPress_; /// Capture the input evaluation point information in internal state. /// @@ -547,7 +556,7 @@ private: /// \param[in] pcow O/W capillary pressure value (Po - Pw). /// /// \return Water saturation value. - std::pair applySwatInit(const double pcow); + std::pair applySwatInit(const Scalar pcow); /// Derive water saturation from SWATINIT data. /// @@ -561,7 +570,7 @@ private: /// /// \return Water saturation value. Input value, possibly mollified by /// current set of material laws. - std::pair applySwatInit(const double pc, const double sw); + std::pair applySwatInit(const Scalar pc, const Scalar sw); /// Invoke material law container's capillary pressure calculator on /// current fluid state. @@ -569,15 +578,15 @@ private: /// Extract gas/oil capillary pressure value (Pg - Po) from current /// fluid state. - double materialLawCapPressGasOil() const; + Scalar materialLawCapPressGasOil() const; /// Extract oil/water capillary pressure value (Po - Pw) from current /// fluid state. - double materialLawCapPressOilWater() const; + Scalar materialLawCapPressOilWater() const; /// Extract gas/water capillary pressure value (Pg - Pw) from current /// fluid state. - double materialLawCapPressGasWater() const; + Scalar materialLawCapPressGasWater() const; /// Predicate for whether specific phase has constant capillary pressure /// curve in current cell. @@ -614,7 +623,7 @@ private: /// function of phase saturation. /// /// \return Phase saturation. - double fromDepthTable(const double contactdepth, + Scalar fromDepthTable(const Scalar contactdepth, const PhaseIdx phasePos, const bool isincr) const; @@ -635,7 +644,7 @@ private: /// /// \return Phase saturation at which capillary pressure attains target /// value. - double invertCapPress(const double pc, + Scalar invertCapPress(const Scalar pc, const PhaseIdx phasePos, const bool isincr) const; @@ -660,14 +669,14 @@ private: // =========================================================================== -template +template void verticalExtent(const CellRange& cells, - const std::vector>& cellZMinMax, + const std::vector>& cellZMinMax, const Parallel::Communication& comm, - std::array& span); + std::array& span); -template -std::pair cellZMinMax(const Element& element); +template +std::pair cellZMinMax(const Element& element); } // namespace Details @@ -681,6 +690,7 @@ template::Entity; + using Scalar = typename FluidSystem::Scalar; public: template InitialStateComputer(MaterialLawManager& materialLawManager, @@ -688,11 +698,11 @@ public: const Grid& grid, const GridView& gridView, const CartesianIndexMapper& cartMapper, - const double grav, + const Scalar grav, const int num_pressure_points = 2000, const bool applySwatInit = true); - using Vec = std::vector; + using Vec = std::vector; using PVec = std::vector; // One per phase. const Vec& temperature() const { return temperature_; } @@ -729,7 +739,7 @@ private: const std::vector& rec, MaterialLawManager& materialLawManager, const Comm& comm, - const double grav); + const Scalar grav); template void cellLoop(const CellRange& cells, @@ -737,21 +747,21 @@ private: template void equilibrateCellCentres(const CellRange& cells, - const EquilReg& eqreg, + const EquilReg& eqreg, const PressTable& ptable, PhaseSat& psat); template void equilibrateHorizontal(const CellRange& cells, - const EquilReg& eqreg, + const EquilReg& eqreg, const int acc, const PressTable& ptable, PhaseSat& psat); - std::vector< std::shared_ptr> > rsFunc_; - std::vector< std::shared_ptr> > rvFunc_; - std::vector< std::shared_ptr> > rvwFunc_; - using TabulatedFunction = Tabulated1DFunction; + std::vector< std::shared_ptr> > rsFunc_; + std::vector< std::shared_ptr> > rvFunc_; + std::vector< std::shared_ptr> > rvwFunc_; + using TabulatedFunction = Tabulated1DFunction; std::vector tempVdTable_; std::vector saltVdTable_; std::vector saltpVdTable_; @@ -767,8 +777,8 @@ private: const CartesianIndexMapper& cartesianIndexMapper_; Vec swatInit_; Vec cellCenterDepth_; - std::vector> cellZSpan_; - std::vector> cellZMinMax_; + std::vector> cellZSpan_; + std::vector> cellZMinMax_; int num_pressure_points_; }; diff --git a/opm/simulators/flow/equil/InitStateEquil_impl.hpp b/opm/simulators/flow/equil/InitStateEquil_impl.hpp index 00ed87e47..e384bccf8 100644 --- a/opm/simulators/flow/equil/InitStateEquil_impl.hpp +++ b/opm/simulators/flow/equil/InitStateEquil_impl.hpp @@ -60,14 +60,14 @@ namespace EQUIL { namespace Details { -template +template void verticalExtent(const CellRange& cells, - const std::vector>& cellZMinMax, + const std::vector>& cellZMinMax, const Parallel::Communication& comm, - std::array& span) + std::array& span) { - span[0] = std::numeric_limits::max(); - span[1] = std::numeric_limits::lowest(); + span[0] = std::numeric_limits::max(); + span[1] = std::numeric_limits::lowest(); // Define vertical span as // @@ -85,10 +85,11 @@ void verticalExtent(const CellRange& cells, span[1] = comm.max(span[1]); } -void subdivisionCentrePoints(const double left, - const double right, +template +void subdivisionCentrePoints(const Scalar left, + const Scalar right, const int numIntervals, - std::vector>& subdiv) + std::vector>& subdiv) { const auto h = (right - left) / numIntervals; @@ -101,13 +102,13 @@ void subdivisionCentrePoints(const double left, } } -template -std::vector> +template +std::vector> horizontalSubdivision(const CellID cell, - const std::pair topbot, + const std::pair topbot, const int numIntervals) { - auto subdiv = std::vector>{}; + auto subdiv = std::vector>{}; subdiv.reserve(2 * numIntervals); if (topbot.first > topbot.second) { @@ -123,12 +124,12 @@ horizontalSubdivision(const CellID cell, return subdiv; } -template -double cellCenterDepth(const Element& element) +template +Scalar cellCenterDepth(const Element& element) { typedef typename Element::Geometry Geometry; static constexpr int zCoord = Element::dimension - 1; - double zz = 0.0; + Scalar zz = 0.0; const Geometry& geometry = element.geometry(); const int corners = geometry.corners(); @@ -138,13 +139,13 @@ double cellCenterDepth(const Element& element) return zz/corners; } -template -std::pair cellZSpan(const Element& element) +template +std::pair cellZSpan(const Element& element) { typedef typename Element::Geometry Geometry; static constexpr int zCoord = Element::dimension - 1; - double bot = 0.0; - double top = 0.0; + Scalar bot = 0.0; + Scalar top = 0.0; const Geometry& geometry = element.geometry(); const int corners = geometry.corners(); @@ -157,36 +158,36 @@ std::pair cellZSpan(const Element& element) return std::make_pair(bot/4, top/4); } -template -std::pair cellZMinMax(const Element& element) +template +std::pair cellZMinMax(const Element& element) { typedef typename Element::Geometry Geometry; static constexpr int zCoord = Element::dimension - 1; const Geometry& geometry = element.geometry(); const int corners = geometry.corners(); assert(corners == 8); - auto min = std::numeric_limits::max(); - auto max = std::numeric_limits::lowest(); + auto min = std::numeric_limits::max(); + auto max = std::numeric_limits::lowest(); for (int i=0; i < corners; ++i) { - min = std::min(min, geometry.corner(i)[zCoord]); - max = std::max(max, geometry.corner(i)[zCoord]); + min = std::min(min, static_cast(geometry.corner(i)[zCoord])); + max = std::max(max, static_cast(geometry.corner(i)[zCoord])); } return std::make_pair(min, max); } -template -RK4IVP::RK4IVP(const RHS& f, - const std::array& span, - const double y0, - const int N) +template +RK4IVP::RK4IVP(const RHS& f, + const std::array& span, + const Scalar y0, + const int N) : N_(N) , span_(span) { - const double h = stepsize(); - const double h2 = h / 2; - const double h6 = h / 6; + const Scalar h = stepsize(); + const Scalar h2 = h / 2; + const Scalar h6 = h / 6; y_.reserve(N + 1); f_.reserve(N + 1); @@ -195,39 +196,39 @@ RK4IVP::RK4IVP(const RHS& f, 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 Scalar x = span_[0] + i*h; + const Scalar 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 Scalar k1 = f_[i]; + const Scalar k2 = f(x + h2, y + h2*k1); + const Scalar k3 = f(x + h2, y + h2*k2); + const Scalar 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())); } - assert (y_.size() == std::vector::size_type(N + 1)); + assert (y_.size() == typename std::vector::size_type(N + 1)); } -template -double RK4IVP:: -operator()(const double x) const +template +Scalar RK4IVP:: +operator()(const Scalar x) const { // Dense output (O(h**3)) according to Shampine // (Hermite interpolation) - const double h = stepsize(); + const Scalar h = stepsize(); int i = (x - span_[0]) / h; - const double t = (x - (span_[0] + i*h)) / h; + const Scalar t = (x - (span_[0] + i*h)) / h; // Crude handling of evaluation point outside "span_"; if (i < 0) { i = 0; } if (N_ <= i) { i = N_ - 1; } - const double y0 = y_[i], y1 = y_[i + 1]; - const double f0 = f_[i], f1 = f_[i + 1]; + const Scalar y0 = y_[i], y1 = y_[i + 1]; + const Scalar f0 = f_[i], f1 = f_[i + 1]; - double u = (1 - 2*t) * (y1 - y0); + Scalar u = (1 - 2*t) * (y1 - y0); u += h * ((t - 1)*f0 + t*f1); u *= t * (t - 1); u += (1 - t)*y0 + t*y1; @@ -235,8 +236,8 @@ operator()(const double x) const return u; } -template -double RK4IVP:: +template +Scalar RK4IVP:: stepsize() const { return (span_[1] - span_[0]) / N_; @@ -249,7 +250,7 @@ Water:: Water(const TabulatedFunction& tempVdTable, const TabulatedFunction& saltVdTable, const int pvtRegionIdx, - const double normGrav) + const Scalar normGrav) : tempVdTable_(tempVdTable) , saltVdTable_(saltVdTable) , pvtRegionIdx_(pvtRegionIdx) @@ -258,22 +259,28 @@ Water(const TabulatedFunction& tempVdTable, } template -double Water:: -operator()(const double depth, - const double press) const +typename Water::Scalar +Water:: +operator()(const Scalar depth, + const Scalar press) const { return this->density(depth, press) * g_; } template -double Water:: -density(const double depth, - const double press) const +typename Water::Scalar +Water:: +density(const Scalar depth, + const Scalar press) const { // The initializing algorithm can give depths outside the range due to numerical noise i.e. we extrapolate - double saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true); - double temp = tempVdTable_.eval(depth, /*extrapolate=*/true); - double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0 /*=Rsw*/, saltConcentration); + Scalar saltConcentration = saltVdTable_.eval(depth, /*extrapolate=*/true); + Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true); + Scalar rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, + temp, + press, + Scalar{0.0} /*=Rsw*/, + saltConcentration); rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); return rho; } @@ -283,7 +290,7 @@ Oil:: Oil(const TabulatedFunction& tempVdTable, const RS& rs, const int pvtRegionIdx, - const double normGrav) + const Scalar normGrav) : tempVdTable_(tempVdTable) , rs_(rs) , pvtRegionIdx_(pvtRegionIdx) @@ -292,31 +299,33 @@ Oil(const TabulatedFunction& tempVdTable, } template -double Oil:: -operator()(const double depth, - const double press) const +typename Oil::Scalar +Oil:: +operator()(const Scalar depth, + const Scalar press) const { return this->density(depth, press) * g_; } template -double Oil:: -density(const double depth, - const double press) const +typename Oil::Scalar +Oil:: +density(const Scalar depth, + const Scalar press) const { - const double temp = tempVdTable_.eval(depth, /*extrapolate=*/true); - double rs = 0.0; - if(FluidSystem::enableDissolvedGas()) + const Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true); + Scalar rs = 0.0; + if (FluidSystem::enableDissolvedGas()) rs = rs_(depth, press, temp); - double bOil = 0.0; + Scalar bOil = 0.0; if (rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press)) { bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); } else { bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rs); } - double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); + Scalar rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); if (FluidSystem::enableDissolvedGas()) { rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); } @@ -330,7 +339,7 @@ Gas(const TabulatedFunction& tempVdTable, const RV& rv, const RVW& rvw, const int pvtRegionIdx, - const double normGrav) + const Scalar normGrav) : tempVdTable_(tempVdTable) , rv_(rv) , rvw_(rvw) @@ -340,28 +349,30 @@ Gas(const TabulatedFunction& tempVdTable, } template -double Gas:: -operator()(const double depth, - const double press) const +typename Gas::Scalar +Gas:: +operator()(const Scalar depth, + const Scalar press) const { return this->density(depth, press) * g_; } template -double Gas:: -density(const double depth, - const double press) const +typename Gas::Scalar +Gas:: +density(const Scalar depth, + const Scalar press) const { - const double temp = tempVdTable_.eval(depth, /*extrapolate=*/true); - double rv = 0.0; + const Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true); + Scalar rv = 0.0; if (FluidSystem::enableVaporizedOil()) rv = rv_(depth, press, temp); - double rvw = 0.0; + Scalar rvw = 0.0; if (FluidSystem::enableVaporizedWater()) rvw = rvw_(depth, press, temp); - double bGas = 0.0; + Scalar bGas = 0.0; if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) { if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press) @@ -371,7 +382,7 @@ density(const double depth, } else { bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, rvw); } - double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_) + rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); return rho; @@ -381,9 +392,13 @@ density(const double depth, if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)) { bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); } else { - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, 0.0/*=rvw*/); + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, + temp, + press, + rv, + Scalar{0.0}/*=rvw*/); } - double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); return rho; } @@ -393,16 +408,23 @@ density(const double depth, bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); } else { - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0/*=rv*/, rvw); + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, + temp, + press, + Scalar{0.0} /*=rv*/, + rvw); } - double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); rho += rvw * bGas * FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); return rho; } // immiscible gas - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, 0.0/*=rv*/, 0.0/*=rvw*/); - double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, + press, + Scalar{0.0} /*=rv*/, + Scalar{0.0} /*=rvw*/); + Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); return rho; } @@ -471,10 +493,10 @@ operator=(PressureFunction&& rhs) template template -double +typename PressureTable::Scalar PressureTable:: PressureFunction:: -value(const double depth) const +value(const Scalar depth) const { if (depth < this->initial_.depth) { // Value above initial condition depth. @@ -548,7 +570,7 @@ copyInPointers(const PressureTable& rhs) template PhaseSaturations:: PhaseSaturations(MaterialLawManager& matLawMgr, - const std::vector& swatInit) + const std::vector& swatInit) : matLawMgr_(matLawMgr) , swatInit_ (swatInit) { @@ -569,7 +591,7 @@ PhaseSaturations(const PhaseSaturations& rhs) } template -const PhaseQuantityValue& +const PhaseQuantityValue& PhaseSaturations:: deriveSaturations(const Position& x, const Region& reg, @@ -864,17 +886,17 @@ accountForScaledSaturations() } template -std::pair +std::pair PhaseSaturations:: -applySwatInit(const double pcow) +applySwatInit(const Scalar pcow) { return this->applySwatInit(pcow, this->swatInit_[this->evalPt_.position->cell]); } template -std::pair +std::pair PhaseSaturations:: -applySwatInit(const double pcow, const double sw) +applySwatInit(const Scalar pcow, const Scalar sw) { return this->matLawMgr_.applySwatinit(this->evalPt_.position->cell, pcow, sw); } @@ -892,7 +914,8 @@ computeMaterialLawCapPress() } template -double PhaseSaturations:: +typename FluidSystem::Scalar +PhaseSaturations:: materialLawCapPressGasOil() const { return this->matLawCapPress_[this->oilPos()] @@ -900,7 +923,8 @@ materialLawCapPressGasOil() const } template -double PhaseSaturations:: +typename FluidSystem::Scalar +PhaseSaturations:: materialLawCapPressOilWater() const { return this->matLawCapPress_[this->oilPos()] @@ -908,7 +932,8 @@ materialLawCapPressOilWater() const } template -double PhaseSaturations:: +typename FluidSystem::Scalar +PhaseSaturations:: materialLawCapPressGasWater() const { return this->matLawCapPress_[this->gasPos()] @@ -933,8 +958,9 @@ isOverlappingTransition() const } template -double PhaseSaturations:: -fromDepthTable(const double contactdepth, +typename FluidSystem::Scalar +PhaseSaturations:: +fromDepthTable(const Scalar contactdepth, const PhaseIdx phasePos, const bool isincr) const { @@ -945,8 +971,9 @@ fromDepthTable(const double contactdepth, } template -double PhaseSaturations:: -invertCapPress(const double pc, +typename FluidSystem::Scalar +PhaseSaturations:: +invertCapPress(const Scalar pc, const PhaseIdx phasePos, const bool isincr) const { @@ -957,7 +984,7 @@ invertCapPress(const double pc, template PressureTable:: -PressureTable(const double gravity, +PressureTable(const Scalar gravity, const int samplePoints) : gravity_(gravity) , nsample_(samplePoints) @@ -1044,8 +1071,9 @@ waterActive() const } template -double PressureTable:: -oil(const double depth) const +typename FluidSystem::Scalar +PressureTable:: +oil(const Scalar depth) const { this->checkPtr(this->oil_.get(), "OIL"); @@ -1053,8 +1081,9 @@ oil(const double depth) const } template -double PressureTable:: -gas(const double depth) const +typename FluidSystem::Scalar +PressureTable:: +gas(const Scalar depth) const { this->checkPtr(this->gas_.get(), "GAS"); @@ -1063,8 +1092,9 @@ gas(const double depth) const template -double PressureTable:: -water(const double depth) const +typename FluidSystem::Scalar +PressureTable:: +water(const Scalar depth) const { this->checkPtr(this->wat_.get(), "WATER"); @@ -1310,16 +1340,16 @@ InitialStateComputer(MaterialLawManager& materialLawManager, const Grid& grid, const GridView& gridView, const CartesianIndexMapper& cartMapper, - const double grav, + const Scalar grav, const int num_pressure_points, const bool applySwatInit) : temperature_(grid.size(/*codim=*/0), eclipseState.getTableManager().rtemp()), saltConcentration_(grid.size(/*codim=*/0)), saltSaturation_(grid.size(/*codim=*/0)), pp_(FluidSystem::numPhases, - std::vector(grid.size(/*codim=*/0))), + std::vector(grid.size(/*codim=*/0))), sat_(FluidSystem::numPhases, - std::vector(grid.size(/*codim=*/0))), + std::vector(grid.size(/*codim=*/0))), rs_(grid.size(/*codim=*/0)), rv_(grid.size(/*codim=*/0)), rvw_(grid.size(/*codim=*/0)), @@ -1329,7 +1359,13 @@ InitialStateComputer(MaterialLawManager& materialLawManager, //Check for presence of kw SWATINIT if (applySwatInit) { if (eclipseState.fieldProps().has_double("SWATINIT")) { - swatInit_ = eclipseState.fieldProps().get_double("SWATINIT"); + if constexpr (std::is_same_v) { + swatInit_ = eclipseState.fieldProps().get_double("SWATINIT"); + } else { + const auto& input = eclipseState.fieldProps().get_double("SWATINIT"); + swatInit_.resize(input.size()); + std::copy(input.begin(), input.end(), swatInit_.begin()); + } } } @@ -1349,6 +1385,19 @@ InitialStateComputer(MaterialLawManager& materialLawManager, // Create Rs functions. rsFunc_.reserve(rec.size()); + + auto getArray = [](const std::vector& input) + { + if constexpr (std::is_same_v) { + return input; + } else { + std::vector output; + output.resize(input.size()); + std::copy(input.begin(), input.end(), output.begin()); + return output; + } + }; + if (FluidSystem::enableDissolvedGas()) { for (std::size_t i = 0; i < rec.size(); ++i) { if (eqlmap.cells(i).empty()) { @@ -1360,16 +1409,15 @@ InitialStateComputer(MaterialLawManager& materialLawManager, const TableContainer& rsvdTables = tables.getRsvdTables(); const TableContainer& pbvdTables = tables.getPbvdTables(); if (rsvdTables.size() > 0) { - const RsvdTable& rsvdTable = rsvdTables.getTable(i); - std::vector depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy(); - std::vector rsColumn = rsvdTable.getColumn("RS").vectorCopy(); + auto depthColumn = getArray(rsvdTable.getColumn("DEPTH").vectorCopy()); + auto rsColumn = getArray(rsvdTable.getColumn("RS").vectorCopy()); rsFunc_.push_back(std::make_shared>(pvtIdx, depthColumn, rsColumn)); } else if (pbvdTables.size() > 0) { const PbvdTable& pbvdTable = pbvdTables.getTable(i); - std::vector depthColumn = pbvdTable.getColumn("DEPTH").vectorCopy(); - std::vector pbubColumn = pbvdTable.getColumn("PBUB").vectorCopy(); + auto depthColumn = getArray(pbvdTable.getColumn("DEPTH").vectorCopy()); + auto pbubColumn = getArray(pbvdTable.getColumn("PBUB").vectorCopy()); rsFunc_.push_back(std::make_shared>(pvtIdx, depthColumn, pbubColumn)); @@ -1384,15 +1432,15 @@ InitialStateComputer(MaterialLawManager& materialLawManager, "datum depth must be at the gas-oil-contact. " "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold."); } - const double pContact = rec[i].datumDepthPressure(); - const double TContact = 273.15 + 20; // standard temperature for now + const Scalar pContact = rec[i].datumDepthPressure(); + const Scalar TContact = 273.15 + 20; // standard temperature for now rsFunc_.push_back(std::make_shared>(pvtIdx, pContact, TContact)); } } } else { for (std::size_t i = 0; i < rec.size(); ++i) { - rsFunc_.push_back(std::make_shared>()); + rsFunc_.push_back(std::make_shared>()); } } @@ -1410,14 +1458,14 @@ InitialStateComputer(MaterialLawManager& materialLawManager, if (rvvdTables.size() > 0) { const RvvdTable& rvvdTable = rvvdTables.getTable(i); - std::vector depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy(); - std::vector rvColumn = rvvdTable.getColumn("RV").vectorCopy(); + auto depthColumn = getArray(rvvdTable.getColumn("DEPTH").vectorCopy()); + auto rvColumn = getArray(rvvdTable.getColumn("RV").vectorCopy()); rvFunc_.push_back(std::make_shared>(pvtIdx, depthColumn, rvColumn)); } else if (pdvdTables.size() > 0) { const PdvdTable& pdvdTable = pdvdTables.getTable(i); - std::vector depthColumn = pdvdTable.getColumn("DEPTH").vectorCopy(); - std::vector pdewColumn = pdvdTable.getColumn("PDEW").vectorCopy(); + auto depthColumn = getArray(pdvdTable.getColumn("DEPTH").vectorCopy()); + auto pdewColumn = getArray(pdvdTable.getColumn("PDEW").vectorCopy()); rvFunc_.push_back(std::make_shared>(pvtIdx, depthColumn, pdewColumn)); } else { @@ -1431,19 +1479,19 @@ InitialStateComputer(MaterialLawManager& materialLawManager, "datum depth must be at the gas-oil-contact. " "In EQUIL region "+std::to_string(i + 1)+" (counting from 1), this does not hold."); } - const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); - const double TContact = 273.15 + 20; // standard temperature for now + const Scalar pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); + const Scalar TContact = 273.15 + 20; // standard temperature for now rvFunc_.push_back(std::make_shared>(pvtIdx,pContact, TContact)); } } } else { for (std::size_t i = 0; i < rec.size(); ++i) { - rvFunc_.push_back(std::make_shared>()); + rvFunc_.push_back(std::make_shared>()); } } - rvwFunc_.reserve(rec.size()); + rvwFunc_.reserve(rec.size()); if (FluidSystem::enableVaporizedWater()) { for (std::size_t i = 0; i < rec.size(); ++i) { if (eqlmap.cells(i).empty()) { @@ -1456,8 +1504,8 @@ InitialStateComputer(MaterialLawManager& materialLawManager, if (rvwvdTables.size() > 0) { const RvwvdTable& rvwvdTable = rvwvdTables.getTable(i); - std::vector depthColumn = rvwvdTable.getColumn("DEPTH").vectorCopy(); - std::vector rvwvdColumn = rvwvdTable.getColumn("RVWVD").vectorCopy(); + auto depthColumn = getArray(rvwvdTable.getColumn("DEPTH").vectorCopy()); + auto rvwvdColumn = getArray(rvwvdTable.getColumn("RVWVD").vectorCopy()); rvwFunc_.push_back(std::make_shared>(pvtIdx, depthColumn, rvwvdColumn)); } else { @@ -1466,9 +1514,9 @@ InitialStateComputer(MaterialLawManager& materialLawManager, } else { const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); - if (oilActive){ + if (oilActive) { if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { - rvwFunc_.push_back(std::make_shared>()); + rvwFunc_.push_back(std::make_shared>()); const auto msg = "No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +". \n" "and datum depth is not at the gas-oil-contact. \n" "Rvw is set to 0.0 in all cells. \n"; @@ -1476,8 +1524,8 @@ InitialStateComputer(MaterialLawManager& materialLawManager, } else { // pg = po + Pcgo = po + (pg - po) // for gas-condensate with initial no oil zone: water-oil contact depth (OWC) equal gas-oil contact depth (GOC) - const double pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); - const double TContact = 273.15 + 20; // standard temperature for now + const Scalar pContact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); + const Scalar TContact = 273.15 + 20; // standard temperature for now rvwFunc_.push_back(std::make_shared>(pvtIdx,pContact, TContact)); } } @@ -1485,15 +1533,15 @@ InitialStateComputer(MaterialLawManager& materialLawManager, // two-phase gas-water sytem: water-oil contact depth is taken equal to gas-water contact depth (GWC) // and water-oil capillary pressure (Pcwo) is taken equal to gas-water capillary pressure (Pcgw) at GWC if (rec[i].waterOilContactDepth() != rec[i].datumDepth()) { - rvwFunc_.push_back(std::make_shared>()); + rvwFunc_.push_back(std::make_shared>()); const auto msg = "No explicit RVWVD table is given for EQUIL region " + std::to_string(i + 1) +". \n" - "and datum depth is not at the gas-water-contact. \n" - "Rvw is set to 0.0 in all cells. \n"; + "and datum depth is not at the gas-water-contact. \n" + "Rvw is set to 0.0 in all cells. \n"; OpmLog::warning(msg); } else { // pg = pw + Pcgw = pw + (pg - pw) - const double pContact = rec[i].datumDepthPressure() + rec[i].waterOilContactCapillaryPressure(); - const double TContact = 273.15 + 20; // standard temperature for now + const Scalar pContact = rec[i].datumDepthPressure() + rec[i].waterOilContactCapillaryPressure(); + const Scalar TContact = 273.15 + 20; // standard temperature for now rvwFunc_.push_back(std::make_shared>(pvtIdx,pContact, TContact)); } } @@ -1502,7 +1550,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager, } else { for (std::size_t i = 0; i < rec.size(); ++i) { - rvwFunc_.push_back(std::make_shared>()); + rvwFunc_.push_back(std::make_shared>()); } } @@ -1544,8 +1592,9 @@ updateInitialTemperature_(const EclipseState& eclState, const RMap& reg) tempVdTable_.resize(numEquilReg); const auto& tables = eclState.getTableManager(); if (!tables.hasTables("RTEMPVD")) { - std::vector x = {0.0,1.0}; - std::vector y = {tables.rtemp(),tables.rtemp()}; + std::vector x = {0.0,1.0}; + std::vector y = {static_cast(tables.rtemp()), + static_cast(tables.rtemp())}; for (auto& table : this->tempVdTable_) { table.setXYContainers(x, y); } @@ -1556,7 +1605,7 @@ updateInitialTemperature_(const EclipseState& eclState, const RMap& reg) tempVdTable_[i].setXYContainers(tempvdTable.getDepthColumn(), tempvdTable.getTemperatureColumn()); const auto& cells = reg.cells(i); for (const auto& cell : cells) { - const double depth = cellCenterDepth_[cell]; + const Scalar depth = cellCenterDepth_[cell]; this->temperature_[cell] = tempVdTable_[i].eval(depth, /*extrapolate=*/true); } } @@ -1583,8 +1632,8 @@ updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg) // If no saltvd table is given, we create a trivial table for the density calculations if (saltvdTables.empty()) { - std::vector x = {0.0,1.0}; - std::vector y = {0.0,0.0}; + std::vector x = {0.0,1.0}; + std::vector y = {0.0,0.0}; for (auto& table : this->saltVdTable_) { table.setXYContainers(x, y); } @@ -1595,7 +1644,7 @@ updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg) const auto& cells = reg.cells(i); for (const auto& cell : cells) { - const double depth = cellCenterDepth_[cell]; + const Scalar depth = cellCenterDepth_[cell]; this->saltConcentration_[cell] = saltVdTable_[i].eval(depth, /*extrapolate=*/true); } } @@ -1626,7 +1675,7 @@ updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg) const auto& cells = reg.cells(i); for (const auto& cell : cells) { - const double depth = cellCenterDepth_[cell]; + const Scalar depth = cellCenterDepth_[cell]; this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth, /*extrapolate=*/true); } } @@ -1658,15 +1707,15 @@ updateCellProps_(const GridView& gridView, for (; elemIt != elemEndIt; ++elemIt) { const Element& element = *elemIt; const unsigned int elemIdx = elemMapper.index(element); - cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element); + cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element); const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); - cellZSpan_[elemIdx] = Details::cellZSpan(element); - cellZMinMax_[elemIdx] = Details::cellZMinMax(element); + cellZSpan_[elemIdx] = Details::cellZSpan(element); + cellZMinMax_[elemIdx] = Details::cellZMinMax(element); if (!num_aqu_cells.empty()) { const auto search = num_aqu_cells.find(cartIx); if (search != num_aqu_cells.end()) { const auto* aqu_cell = num_aqu_cells.at(cartIx); - const double depth_change_num_aqu = aqu_cell->depth - cellCenterDepth_[elemIdx]; + const Scalar depth_change_num_aqu = aqu_cell->depth - cellCenterDepth_[elemIdx]; cellCenterDepth_[elemIdx] += depth_change_num_aqu; cellZSpan_[elemIdx].first += depth_change_num_aqu; cellZSpan_[elemIdx].second += depth_change_num_aqu; @@ -1731,7 +1780,7 @@ applyNumericalAquifers_(const GridView& gridView, // if pressure is specified for numerical aquifers, we use these pressure values // for numerical aquifer cells if (aqu_cell->init_pressure) { - const double pres = *(aqu_cell->init_pressure); + const Scalar pres = *(aqu_cell->init_pressure); this->pp_[watPos][elemIdx] = pres; if (FluidSystem::phaseIsActive(gasPos)) { this->pp_[gasPos][elemIdx] = pres; @@ -1780,15 +1829,15 @@ calcPressSatRsRv(const RMap& reg, const std::vector& rec, MaterialLawManager& materialLawManager, const Comm& comm, - const double grav) + const Scalar grav) { using PhaseSat = Details::PhaseSaturations< - MaterialLawManager, FluidSystem, EquilReg, typename RMap::CellId + MaterialLawManager, FluidSystem, EquilReg, typename RMap::CellId >; - auto ptable = Details::PressureTable>{ grav, this->num_pressure_points_ }; + auto ptable = Details::PressureTable>{ grav, this->num_pressure_points_ }; auto psat = PhaseSat { materialLawManager, this->swatInit_ }; - auto vspan = std::array{}; + auto vspan = std::array{}; std::vector regionIsEmpty(rec.size(), 0); for (std::size_t r = 0; r < rec.size(); ++r) { @@ -1867,11 +1916,11 @@ cellLoop(const CellRange& cells, const auto gasActive = FluidSystem::phaseIsActive(gasPos); const auto watActive = FluidSystem::phaseIsActive(watPos); - auto pressures = Details::PhaseQuantityValue{}; - auto saturations = Details::PhaseQuantityValue{}; - auto Rs = 0.0; - auto Rv = 0.0; - auto Rvw = 0.0; + auto pressures = Details::PhaseQuantityValue{}; + auto saturations = Details::PhaseQuantityValue{}; + Scalar Rs = 0.0; + Scalar Rv = 0.0; + Scalar Rvw = 0.0; for (const auto& cell : cells) { eqmethod(cell, pressures, saturations, Rs, Rv, Rvw); @@ -1914,7 +1963,7 @@ void InitialStateComputer:: equilibrateCellCentres(const CellRange& cells, - const EquilReg& eqreg, + const EquilReg& eqreg, const PressTable& ptable, PhaseSat& psat) { @@ -1923,11 +1972,11 @@ equilibrateCellCentres(const CellRange& cells, decltype(std::declval().cell)>>; this->cellLoop(cells, [this, &eqreg, &ptable, &psat] (const CellID cell, - Details::PhaseQuantityValue& pressures, - Details::PhaseQuantityValue& saturations, - double& Rs, - double& Rv, - double& Rvw) -> void + Details::PhaseQuantityValue& pressures, + Details::PhaseQuantityValue& saturations, + Scalar& Rs, + Scalar& Rv, + Scalar& Rvw) -> void { const auto pos = CellPos { cell, cellCenterDepth_[cell] @@ -1961,7 +2010,7 @@ void InitialStateComputer:: equilibrateHorizontal(const CellRange& cells, - const EquilReg& eqreg, + const EquilReg& eqreg, const int acc, const PressTable& ptable, PhaseSat& psat) @@ -1972,16 +2021,16 @@ equilibrateHorizontal(const CellRange& cells, this->cellLoop(cells, [this, acc, &eqreg, &ptable, &psat] (const CellID cell, - Details::PhaseQuantityValue& pressures, - Details::PhaseQuantityValue& saturations, - double& Rs, - double& Rv, - double& Rvw) -> void + Details::PhaseQuantityValue& pressures, + Details::PhaseQuantityValue& saturations, + Scalar& Rs, + Scalar& Rv, + Scalar& Rvw) -> void { pressures .reset(); saturations.reset(); - auto totfrac = 0.0; + Scalar totfrac = 0.0; for (const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) { const auto pos = CellPos { cell, depth }; diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 2d560e188..f688f8a8e 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -137,7 +137,7 @@ static std::vector> cellVerticalExtent(const GridView& for (; elemIt != elemEndIt; ++elemIt) { const auto& element = *elemIt; const unsigned int elemIdx = elemMapper.index(element); - cellZMinMax[elemIdx] = Opm::EQUIL::Details::cellZMinMax(element); + cellZMinMax[elemIdx] = Opm::EQUIL::Details::cellZMinMax(element); } return cellZMinMax; }