From 072750fedc09550d13bc092ca6343dacf7d94b99 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 27 Mar 2026 15:01:18 +0100 Subject: [PATCH 1/2] refactoring InitStateEquil so that it can get ready for compositional equilibration. --- CMakeLists_files.cmake | 2 + opm/simulators/flow/equil/InitStateEquil.hpp | 722 +------ .../flow/equil/InitStateEquilBase.hpp | 765 ++++++++ .../flow/equil/InitStateEquilBase_impl.hpp | 1727 +++++++++++++++++ .../flow/equil/InitStateEquil_impl.hpp | 366 +--- 5 files changed, 2555 insertions(+), 1027 deletions(-) create mode 100644 opm/simulators/flow/equil/InitStateEquilBase.hpp create mode 100644 opm/simulators/flow/equil/InitStateEquilBase_impl.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index ad147010f22..2aa1db0d7b8 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -1035,6 +1035,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/VtkTracerModule.hpp opm/simulators/flow/equil/EquilibrationHelpers.hpp opm/simulators/flow/equil/EquilibrationHelpers_impl.hpp + opm/simulators/flow/equil/InitStateEquilBase.hpp + opm/simulators/flow/equil/InitStateEquilBase_impl.hpp opm/simulators/flow/equil/InitStateEquil.hpp opm/simulators/flow/equil/InitStateEquil_impl.hpp opm/simulators/flow/rescoup/ReservoirCouplingEnabled.hpp diff --git a/opm/simulators/flow/equil/InitStateEquil.hpp b/opm/simulators/flow/equil/InitStateEquil.hpp index 36edfd4aab8..87e5dacb088 100644 --- a/opm/simulators/flow/equil/InitStateEquil.hpp +++ b/opm/simulators/flow/equil/InitStateEquil.hpp @@ -23,676 +23,22 @@ /** * \file * - * \brief Routines that actually solve the ODEs that emerge from the hydrostatic - * equilibrium problem + * \brief Black-oil specific routines for hydrostatic equilibrium-based + * initialisation, including dissolved gas (Rs), vaporized oil (Rv), + * and vaporized water (Rvw). */ #ifndef OPM_INIT_STATE_EQUIL_HPP #define OPM_INIT_STATE_EQUIL_HPP -#include - -#include -#include - -#include - -#include -#include -#include -#include -#include -#include +#include namespace Opm { -class EclipseState; -class EquilRecord; -class NumericalAquifers; - -/** - * Types and routines that collectively implement a basic - * ECLIPSE-style equilibration-based initialisation scheme. - * - * This namespace is intentionally nested to avoid name clashes - * with other parts of OPM. - */ namespace EQUIL { -template struct CellCornerData { - std::array X; - std::array Y; - std::array Z; - CellCornerData() = default; - - CellCornerData(const std::array& x, - const std::array& y, - const std::array& z) - : X(x), Y(y), Z(z) - {} -}; - template class EquilReg; namespace Miscibility { template class RsFunction; } -namespace Details { -template -class RK4IVP -{ -public: - RK4IVP(const RHS& f, - const std::array& span, - const Scalar y0, - const int N); - - Scalar operator()(const Scalar x) const; - -private: - int N_; - std::array span_; - std::vector y_; - std::vector f_; - - Scalar stepsize() const; -}; - -namespace PhasePressODE { -template -class Water -{ - using Scalar = typename FluidSystem::Scalar; - using TabulatedFunction = Tabulated1DFunction; - -public: - Water(const TabulatedFunction& tempVdTable, - const TabulatedFunction& saltVdTable, - const int pvtRegionIdx, - const Scalar normGrav); - - Scalar operator()(const Scalar depth, - const Scalar press) const; - -private: - const TabulatedFunction& tempVdTable_; - const TabulatedFunction& saltVdTable_; - const int pvtRegionIdx_; - const Scalar g_; - - Scalar density(const Scalar depth, - const Scalar press) const; -}; - -template -class Oil -{ - using Scalar = typename FluidSystem::Scalar; - using TabulatedFunction = Tabulated1DFunction; - -public: - Oil(const TabulatedFunction& tempVdTable, - const RS& rs, - const int pvtRegionIdx, - const Scalar normGrav); - - Scalar operator()(const Scalar depth, - const Scalar press) const; - -private: - const TabulatedFunction& tempVdTable_; - const RS& rs_; - const int pvtRegionIdx_; - const Scalar g_; - - Scalar density(const Scalar depth, - const Scalar press) const; -}; - -template -class Gas -{ - using Scalar = typename FluidSystem::Scalar; - using TabulatedFunction = Tabulated1DFunction; - -public: - Gas(const TabulatedFunction& tempVdTable, - const RV& rv, - const RVW& rvw, - const int pvtRegionIdx, - const Scalar normGrav); - - Scalar operator()(const Scalar depth, - const Scalar press) const; - -private: - const TabulatedFunction& tempVdTable_; - const RV& rv_; - const RVW& rvw_; - const int pvtRegionIdx_; - const Scalar g_; - - Scalar density(const Scalar depth, - const Scalar press) const; -}; - -} // namespace PhasePressODE - -template -class PressureTable -{ -public: - using Scalar = typename FluidSystem::Scalar; - using VSpan = std::array; - - /// Constructor - /// - /// \param[in] gravity Norm of gravity vector (acceleration strength due - /// to gravity). Normally the standardised value at Tellus equator - /// (9.80665 m/s^2). - /// - /// \param[in] samplePoints Number of equally spaced depth sample points - /// in each internal phase pressure table. - explicit PressureTable(const Scalar gravity, - const int samplePoints = 2000); - - /// Copy constructor - /// - /// \param[in] rhs Source object for copy initialization. - PressureTable(const PressureTable& rhs); - - /// Move constructor - /// - /// \param[in,out] rhs Source object for move initialization. On output, - /// left in a moved-from ("valid but unspecified") state. Internal - /// pointers in \p rhs are null (\c unique_ptr guarantee). - PressureTable(PressureTable&& rhs); - - /// Assignment operator - /// - /// \param[in] rhs Source object. - /// - /// \return \code *this \endcode. - PressureTable& operator=(const PressureTable& rhs); - - /// Move-assignment operator - /// - /// \param[in] rhs Source object. On output, left in a moved-from ("valid - /// but unspecified") state. Internal pointers in \p rhs are null (\c - /// unique_ptr guarantee). - /// - /// \return \code *this \endcode. - PressureTable& operator=(PressureTable&& rhs); - - void equilibrate(const Region& reg, - const VSpan& span); - - /// Predicate for whether or not oil is an active phase - bool oilActive() const; - - /// Predicate for whether or not gas is an active phase - bool gasActive() const; - - /// Predicate for whether or not water is an active phase - bool waterActive() const; - - /// Evaluate oil phase pressure at specified depth. - /// - /// \param[in] depth Depth of evaluation point. Should generally be - /// within the \c span from the previous call to \code equilibrate() - /// \endcode. - /// - /// \return Oil phase pressure at specified depth. - Scalar oil(const Scalar depth) const; - - /// Evaluate gas phase pressure at specified depth. - /// - /// \param[in] depth Depth of evaluation point. Should generally be - /// within the \c span from the previous call to \code equilibrate() - /// \endcode. - /// - /// \return Gas phase pressure at specified depth. - Scalar gas(const Scalar depth) const; - - /// Evaluate water phase pressure at specified depth. - /// - /// \param[in] depth Depth of evaluation point. Should generally be - /// within the \c span from the previous call to \code equilibrate() - /// \endcode. - /// - /// \return Water phase pressure at specified depth. - Scalar water(const Scalar depth) const; - -private: - template - class PressureFunction - { - public: - struct InitCond { - Scalar depth; - Scalar pressure; - }; - - explicit PressureFunction(const ODE& ode, - const InitCond& ic, - const int nsample, - const VSpan& span); - - PressureFunction(const PressureFunction& rhs); - - PressureFunction(PressureFunction&& rhs) = default; - - PressureFunction& operator=(const PressureFunction& rhs); - - PressureFunction& operator=(PressureFunction&& rhs); - - Scalar value(const Scalar depth) const; - - private: - enum Direction : std::size_t { Up, Down, NumDir }; - - using Distribution = Details::RK4IVP; - using DistrPtr = std::unique_ptr; - - InitCond initial_; - std::array value_; - }; - - using OilPressODE = PhasePressODE::Oil< - FluidSystem, typename Region::CalcDissolution - >; - - using GasPressODE = PhasePressODE::Gas< - FluidSystem, typename Region::CalcEvaporation, typename Region::CalcWaterEvaporation - >; - - using WatPressODE = PhasePressODE::Water; - - using OPress = PressureFunction; - using GPress = PressureFunction; - using WPress = PressureFunction; - - using Strategy = void (PressureTable::*) - (const Region&, const VSpan&); - - Scalar gravity_; - int nsample_; - - std::unique_ptr oil_{}; - std::unique_ptr gas_{}; - std::unique_ptr wat_{}; - - template - void checkPtr(const PressFunc* phasePress, - const std::string& phaseName) const; - - Strategy selectEquilibrationStrategy(const Region& reg) const; - - void copyInPointers(const PressureTable& rhs); - - void equil_WOG(const Region& reg, const VSpan& span); - void equil_GOW(const Region& reg, const VSpan& span); - void equil_OWG(const Region& reg, const VSpan& span); - - void makeOilPressure(const typename OPress::InitCond& ic, - const Region& reg, - const VSpan& span); - - void makeGasPressure(const typename GPress::InitCond& ic, - const Region& reg, - const VSpan& span); - - void makeWatPressure(const typename WPress::InitCond& ic, - const Region& reg, - const VSpan& span); -}; - -// =========================================================================== - -/// Simple set of per-phase (named by primary component) quantities. -template -struct PhaseQuantityValue { - Scalar oil{0.0}; - Scalar gas{0.0}; - Scalar water{0.0}; - - PhaseQuantityValue& axpy(const PhaseQuantityValue& rhs, const Scalar a) - { - this->oil += a * rhs.oil; - this->gas += a * rhs.gas; - this->water += a * rhs.water; - - return *this; - } - - PhaseQuantityValue& operator/=(const Scalar x) - { - this->oil /= x; - this->gas /= x; - this->water /= x; - - return *this; - } - - void reset() - { - this->oil = this->gas = this->water = 0.0; - } -}; - -/// Calculator for phase saturations -/// -/// Computes saturation values at arbitrary depths. -/// -/// \tparam MaterialLawManager Container for material laws. Typically a -/// specialization of the \code Opm::EclMaterialLawManager<> \endcode -/// template. -/// -/// \tparam FluidSystem An OPM fluid system type. Typically a -/// specialization of the \code Opm::BlackOilFluidSystem<> \endcode -/// template. -/// -/// \tparam Region Representation of an equilibration region. Typically -/// \code Opm::EQUIL::EquilReg \endcode from the equilibrationhelpers. -/// -/// \tparam CellID Representation an equilibration region's cell IDs. -/// Typically \code std::size_t \endcode. -template -class PhaseSaturations -{ -public: - using Scalar = typename FluidSystem::Scalar; - /// Evaluation point within a model geometry. - /// - /// Associates a particular depth to specific cell. - struct Position { - CellID cell; - Scalar depth; - }; - - /// Convenience type alias - using PTable = PressureTable; - - /// Constructor - /// - /// \param[in,out] matLawMgr Read/write reference to a material law - /// container. Mutated by member functions. - /// - /// \param[in] swatInit Initial water saturation array (from SWATINIT - /// data). Empty if SWATINIT is not used in this simulation model. - explicit PhaseSaturations(MaterialLawManager& matLawMgr, - const std::vector& swatInit); - - /// Copy constructor. - /// - /// \param[in] rhs Source object. - PhaseSaturations(const PhaseSaturations& rhs); - - /// Disabled assignment operator. - PhaseSaturations& operator=(const PhaseSaturations&) = delete; - - /// Disabled move-assignment operator. - PhaseSaturations& operator=(PhaseSaturations&&) = delete; - - /// Calculate phase saturations at particular point of the simulation - /// model geometry. - /// - /// \param[in] x Specific geometric point (depth within a specific cell). - /// - /// \param[in] reg Equilibration information for a single equilibration - /// region; notably contact depths. - /// - /// \param[in] ptable Previously equilibrated phase pressure table - /// pertaining to the equilibration region \p reg. - /// - /// \return Set of phase saturation values defined at particular point. - const PhaseQuantityValue& - deriveSaturations(const Position& x, - const Region& reg, - const PTable& ptable); - - /// Retrieve saturation-corrected phase pressures - /// - /// Values associated with evaluation point of previous call to \code - /// deriveSaturations() \endcode. - const PhaseQuantityValue& correctedPhasePressures() const - { - return this->press_; - } - -private: - /// Convenience amalgamation of the deriveSaturations() input state. - /// These values are almost always used in concert. - struct EvaluationPoint { - const Position* position{nullptr}; - const Region* region {nullptr}; - const PTable* ptable {nullptr}; - }; - - /// Simplified fluid state object that contains only the pieces of - /// information needed to calculate the capillary pressure values from - /// the current set of material laws. - using FluidState = ::Opm:: - SimpleModularFluidState; - - /// Convenience type alias. - using MaterialLaw = typename MaterialLawManager::MaterialLaw; - - /// Fluid system's representation of phase indices. - using PhaseIdx = std::remove_cv_t< - std::remove_reference_t - >; - - /// Read/write reference to client's material law container. - MaterialLawManager& matLawMgr_; - - /// Client's SWATINIT data. - const std::vector& swatInit_; - - /// Evaluated phase saturations. - PhaseQuantityValue sat_; - - /// Saturation-corrected phase pressure values. - PhaseQuantityValue press_; - - /// Current evaluation point. - EvaluationPoint evalPt_; - - /// Capillary pressure fluid state. - FluidState fluidState_; - - /// Evaluated capillary pressures from current set of material laws. - std::array matLawCapPress_; - - /// Capture the input evaluation point information in internal state. - /// - /// \param[in] x Specific geometric point (depth within a specific cell). - /// - /// \param[in] reg Equilibration information for a single equilibration - /// region; notably contact depths. - /// - /// \param[in] ptable Previously equilibrated phase pressure table - /// pertaining to the equilibration region \p reg. - void setEvaluationPoint(const Position& x, - const Region& reg, - const PTable& ptable); - - /// Initialize phase saturation and phase pressure values. - /// - /// Looks up phase pressure values from the input pressure table. - void initializePhaseQuantities(); - - /// Derive phase saturation for oil. - /// - /// Calculated as 1 - Sw - Sg. - void deriveOilSat(); - - /// Derive phase saturation for gas. - /// - /// Inverts capillary pressure curve if non-constant or uses a simple - /// depth consideration with respect to G/O contact depth otherwise. - void deriveGasSat(); - - /// Derive phase saturation for water. - /// - /// Uses input data if simulation model is defined in terms of SWATINIT. - /// Otherwise, inverts capillary pressure curve if non-constant or uses - /// a simple depth consideration with respect to the O/W contact depth - /// if capillary pressure curve is constant within the current cell. - void deriveWaterSat(); - - /// Correct phase saturation and pressure values to account for - /// overlapping transition zones between G/O and O/W systems. - void fixUnphysicalTransition(); - - /// Re-adjust phase pressure values to account for phase saturations - /// outside permissible ranges. - void accountForScaledSaturations(); - - // -------------------------------------------------------------------- - // Note: Function 'applySwatInit' is non-const because the overload set - // needs to mutate the 'matLawMgr_'. - // -------------------------------------------------------------------- - - /// Derive water saturation from SWATINIT data. - /// - /// Uses SWATINIT array data from current cell directly. Also updates - /// the material law container's internal notion of the maximum - /// attainable O/W capillary pressure value. - /// - /// \param[in] pcow O/W capillary pressure value (Po - Pw). - /// - /// \return Water saturation value. - std::pair applySwatInit(const Scalar pcow); - - /// Derive water saturation from SWATINIT data. - /// - /// Uses explicitly passed-in saturation value. Also updates the - /// material law container's internal notion of the maximum attainable - /// O/W capillary pressure value. - /// - /// \param[in] pc x/W capillary pressure value (Px - Pw; x in {O, G}). - /// - /// \param[in] sw Water saturation value. - /// - /// \return Water saturation value. Input value, possibly mollified by - /// current set of material laws. - std::pair applySwatInit(const Scalar pc, const Scalar sw); - - /// Invoke material law container's capillary pressure calculator on - /// current fluid state. - void computeMaterialLawCapPress(); - - /// Extract gas/oil capillary pressure value (Pg - Po) from current - /// fluid state. - Scalar materialLawCapPressGasOil() const; - - /// Extract oil/water capillary pressure value (Po - Pw) from current - /// fluid state. - Scalar materialLawCapPressOilWater() const; - - /// Extract gas/water capillary pressure value (Pg - Pw) from current - /// fluid state. - Scalar materialLawCapPressGasWater() const; - - /// Predicate for whether specific phase has constant capillary pressure - /// curve in current cell. - /// - /// \param[in] phaseIdx Phase. Typically gas or water. - /// - /// \return Whether or not \p phaseIdx has constant capillary pressure - /// curve in current cell. - bool isConstCapPress(const PhaseIdx phaseIdx) const; - - /// Predicate for whether or not the G/O and O/W transition zones - /// overlap in the current cell. - /// - /// This is the case when inverting the capillary pressure curves - /// produces a negative oil saturation--i.e., when Sg + Sw > 1. - bool isOverlappingTransition() const; - - /// Derive phase saturation value from simple depth consideration. - /// - /// Assumes that the pertinent capillary pressure curve is constant - /// (typically zero) in the current cell--i.e., that there is a sharp - /// interface between the two phases. - /// - /// \param[in] contactdepth Depth of relevant phase separation contact. - /// - /// \param[in] Position of phase in three-phase enumeration. Typically - /// \code gasPos() \endcode or \code waterPos() \endcode. - /// - /// \param[in] isincr Whether the capillary pressure curve is normally - /// increasing as a function of phase saturation (e.g., Pcgo(Sg) = Pg - /// - Po) or if the curve is normally decreasing as a function of - /// increasing phase saturation (e.g., Pcow(Sw) = Po - Pw). True for - /// capillary pressure functions that are normally increasing as a - /// function of phase saturation. - /// - /// \return Phase saturation. - Scalar fromDepthTable(const Scalar contactdepth, - const PhaseIdx phasePos, - const bool isincr) const; - - /// Derive phase saturation by inverting non-constant capillary pressure - /// curve. - /// - /// \param[in] pc Target capillary pressure value. - /// - /// \param[in] Position of phase in three-phase enumeration. Typically - /// \code gasPos() \endcode or \code waterPos() \endcode. - /// - /// \param[in] isincr Whether the capillary pressure curve is normally - /// increasing as a function of phase saturation (e.g., Pcgo(Sg) = Pg - /// - Po) or if the curve is normally decreasing as a function of - /// increasing phase saturation (e.g., Pcow(Sw) = Po - Pw). True for - /// capillary pressure functions that are normally increasing as a - /// function of phase saturation. - /// - /// \return Phase saturation at which capillary pressure attains target - /// value. - Scalar invertCapPress(const Scalar pc, - const PhaseIdx phasePos, - const bool isincr) const; - - /// Position of oil in fluid system's three-phase enumeration. - PhaseIdx oilPos() const - { - return FluidSystem::oilPhaseIdx; - } - - /// Position of gas in fluid system's three-phase enumeration. - PhaseIdx gasPos() const - { - return FluidSystem::gasPhaseIdx; - } - - /// Position of water in fluid system's three-phase enumeration. - PhaseIdx waterPos() const - { - return FluidSystem::waterPhaseIdx; - } -}; - -// =========================================================================== - -template -void verticalExtent(const CellRange& cells, - const std::vector>& cellZMinMax, - const Parallel::Communication& comm, - std::array& span); - -template -std::pair cellZMinMax(const Element& element); - -} // namespace Details - namespace DeckDependent { template class InitialStateComputer + : public InitialStateComputerBase { - using Element = typename GridView::template Codim<0>::Entity; + using Base = InitialStateComputerBase; using Scalar = typename FluidSystem::Scalar; + public: + using typename Base::Vec; + using typename Base::PVec; + template InitialStateComputer(MaterialLawManager& materialLawManager, const EclipseState& eclipseState, @@ -715,38 +68,17 @@ class InitialStateComputer const int num_pressure_points = 2000, const bool applySwatInit = true); - using Vec = std::vector; - using PVec = std::vector; // One per phase. + using Base::temperature; + using Base::saltConcentration; + using Base::saltSaturation; + using Base::press; + using Base::saturation; - const Vec& temperature() const { return temperature_; } - const Vec& saltConcentration() const { return saltConcentration_; } - const Vec& saltSaturation() const { return saltSaturation_; } - const PVec& press() const { return pp_; } - const PVec& saturation() const { return sat_; } const Vec& rs() const { return rs_; } const Vec& rv() const { return rv_; } const Vec& rvw() const { return rvw_; } private: - template - void updateInitialTemperature_(const EclipseState& eclState, const RMap& reg); - - template - void updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg); - - template - void updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg); - - void updateCellProps_(const GridView& gridView, - const NumericalAquifers& aquifer); - - void applyNumericalAquifers_(const GridView& gridView, - const NumericalAquifers& aquifer, - const bool co2store_or_h2store); - - template - void setRegionPvtIdx(const EclipseState& eclState, const RMap& reg); - template void calcPressSatRsRv(const RMap& reg, const std::vector& rec, @@ -773,7 +105,7 @@ class InitialStateComputer PhaseSat& psat); template - void equilibrateTiltedFaultBlock(const CellRange& cells, + void equilibrateTiltedFaultBlock(const CellRange& cells, const EquilReg& eqreg, const GridView& gridView, const int numLevels, const PressTable& ptable, PhaseSat& psat); @@ -787,27 +119,9 @@ class InitialStateComputer 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_; - std::vector regionPvtIdx_; - Vec temperature_; - Vec saltConcentration_; - Vec saltSaturation_; - PVec pp_; - PVec sat_; Vec rs_; Vec rv_; Vec rvw_; - const CartesianIndexMapper& cartesianIndexMapper_; - Vec swatInit_; - Vec cellCenterDepth_; - std::vector> cellCenterXY_; - std::vector> cellZSpan_; - std::vector> cellZMinMax_; - std::vector> cellCorners_; - int num_pressure_points_; }; } // namespace DeckDependent diff --git a/opm/simulators/flow/equil/InitStateEquilBase.hpp b/opm/simulators/flow/equil/InitStateEquilBase.hpp new file mode 100644 index 00000000000..edcb4d5d7bf --- /dev/null +++ b/opm/simulators/flow/equil/InitStateEquilBase.hpp @@ -0,0 +1,765 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/** + * \file + * + * \brief Base routines for hydrostatic equilibrium-based initialisation, + * common to both black-oil and compositional models. + */ +#ifndef OPM_INIT_STATE_EQUIL_BASE_HPP +#define OPM_INIT_STATE_EQUIL_BASE_HPP + +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { + +class EclipseState; +class EquilRecord; +class NumericalAquifers; + +/** + * Types and routines that collectively implement a basic + * ECLIPSE-style equilibration-based initialisation scheme. + * + * This namespace is intentionally nested to avoid name clashes + * with other parts of OPM. + */ +namespace EQUIL { + +template struct CellCornerData { + std::array X; + std::array Y; + std::array Z; + CellCornerData() = default; + + CellCornerData(const std::array& x, + const std::array& y, + const std::array& z) + : X(x), Y(y), Z(z) + {} +}; + +namespace Details { +template +class RK4IVP +{ +public: + RK4IVP(const RHS& f, + const std::array& span, + const Scalar y0, + const int N); + + Scalar operator()(const Scalar x) const; + +private: + int N_; + std::array span_; + std::vector y_; + std::vector f_; + + Scalar stepsize() const; +}; + +namespace PhasePressODE { +template +class Water +{ + using Scalar = typename FluidSystem::Scalar; + using TabulatedFunction = Tabulated1DFunction; + +public: + Water(const TabulatedFunction& tempVdTable, + const TabulatedFunction& saltVdTable, + const int pvtRegionIdx, + const Scalar normGrav); + + Scalar operator()(const Scalar depth, + const Scalar press) const; + +private: + const TabulatedFunction& tempVdTable_; + const TabulatedFunction& saltVdTable_; + const int pvtRegionIdx_; + const Scalar g_; + + Scalar density(const Scalar depth, + const Scalar press) const; +}; + +template +class Oil +{ + using Scalar = typename FluidSystem::Scalar; + using TabulatedFunction = Tabulated1DFunction; + +public: + Oil(const TabulatedFunction& tempVdTable, + const RS& rs, + const int pvtRegionIdx, + const Scalar normGrav); + + Scalar operator()(const Scalar depth, + const Scalar press) const; + +private: + const TabulatedFunction& tempVdTable_; + const RS& rs_; + const int pvtRegionIdx_; + const Scalar g_; + + Scalar density(const Scalar depth, + const Scalar press) const; +}; + +template +class Gas +{ + using Scalar = typename FluidSystem::Scalar; + using TabulatedFunction = Tabulated1DFunction; + +public: + Gas(const TabulatedFunction& tempVdTable, + const RV& rv, + const RVW& rvw, + const int pvtRegionIdx, + const Scalar normGrav); + + Scalar operator()(const Scalar depth, + const Scalar press) const; + +private: + const TabulatedFunction& tempVdTable_; + const RV& rv_; + const RVW& rvw_; + const int pvtRegionIdx_; + const Scalar g_; + + Scalar density(const Scalar depth, + const Scalar press) const; +}; + +} // namespace PhasePressODE + +template +class PressureTable +{ +public: + using Scalar = typename FluidSystem::Scalar; + using VSpan = std::array; + + /// Constructor + /// + /// \param[in] gravity Norm of gravity vector (acceleration strength due + /// to gravity). Normally the standardised value at Tellus equator + /// (9.80665 m/s^2). + /// + /// \param[in] samplePoints Number of equally spaced depth sample points + /// in each internal phase pressure table. + explicit PressureTable(const Scalar gravity, + const int samplePoints = 2000); + + /// Copy constructor + /// + /// \param[in] rhs Source object for copy initialization. + PressureTable(const PressureTable& rhs); + + /// Move constructor + /// + /// \param[in,out] rhs Source object for move initialization. On output, + /// left in a moved-from ("valid but unspecified") state. Internal + /// pointers in \p rhs are null (\c unique_ptr guarantee). + PressureTable(PressureTable&& rhs); + + /// Assignment operator + /// + /// \param[in] rhs Source object. + /// + /// \return \code *this \endcode. + PressureTable& operator=(const PressureTable& rhs); + + /// Move-assignment operator + /// + /// \param[in] rhs Source object. On output, left in a moved-from ("valid + /// but unspecified") state. Internal pointers in \p rhs are null (\c + /// unique_ptr guarantee). + /// + /// \return \code *this \endcode. + PressureTable& operator=(PressureTable&& rhs); + + void equilibrate(const Region& reg, + const VSpan& span); + + /// Predicate for whether or not oil is an active phase + bool oilActive() const; + + /// Predicate for whether or not gas is an active phase + bool gasActive() const; + + /// Predicate for whether or not water is an active phase + bool waterActive() const; + + /// Evaluate oil phase pressure at specified depth. + /// + /// \param[in] depth Depth of evaluation point. Should generally be + /// within the \c span from the previous call to \code equilibrate() + /// \endcode. + /// + /// \return Oil phase pressure at specified depth. + Scalar oil(const Scalar depth) const; + + /// Evaluate gas phase pressure at specified depth. + /// + /// \param[in] depth Depth of evaluation point. Should generally be + /// within the \c span from the previous call to \code equilibrate() + /// \endcode. + /// + /// \return Gas phase pressure at specified depth. + Scalar gas(const Scalar depth) const; + + /// Evaluate water phase pressure at specified depth. + /// + /// \param[in] depth Depth of evaluation point. Should generally be + /// within the \c span from the previous call to \code equilibrate() + /// \endcode. + /// + /// \return Water phase pressure at specified depth. + Scalar water(const Scalar depth) const; + +private: + template + class PressureFunction + { + public: + struct InitCond { + Scalar depth; + Scalar pressure; + }; + + explicit PressureFunction(const ODE& ode, + const InitCond& ic, + const int nsample, + const VSpan& span); + + PressureFunction(const PressureFunction& rhs); + + PressureFunction(PressureFunction&& rhs) = default; + + PressureFunction& operator=(const PressureFunction& rhs); + + PressureFunction& operator=(PressureFunction&& rhs); + + Scalar value(const Scalar depth) const; + + private: + enum Direction : std::size_t { Up, Down, NumDir }; + + using Distribution = Details::RK4IVP; + using DistrPtr = std::unique_ptr; + + InitCond initial_; + std::array value_; + }; + + using OilPressODE = PhasePressODE::Oil< + FluidSystem, typename Region::CalcDissolution + >; + + using GasPressODE = PhasePressODE::Gas< + FluidSystem, typename Region::CalcEvaporation, typename Region::CalcWaterEvaporation + >; + + using WatPressODE = PhasePressODE::Water; + + using OPress = PressureFunction; + using GPress = PressureFunction; + using WPress = PressureFunction; + + using Strategy = void (PressureTable::*) + (const Region&, const VSpan&); + + Scalar gravity_; + int nsample_; + + std::unique_ptr oil_{}; + std::unique_ptr gas_{}; + std::unique_ptr wat_{}; + + template + void checkPtr(const PressFunc* phasePress, + const std::string& phaseName) const; + + Strategy selectEquilibrationStrategy(const Region& reg) const; + + void copyInPointers(const PressureTable& rhs); + + void equil_WOG(const Region& reg, const VSpan& span); + void equil_GOW(const Region& reg, const VSpan& span); + void equil_OWG(const Region& reg, const VSpan& span); + + void makeOilPressure(const typename OPress::InitCond& ic, + const Region& reg, + const VSpan& span); + + void makeGasPressure(const typename GPress::InitCond& ic, + const Region& reg, + const VSpan& span); + + void makeWatPressure(const typename WPress::InitCond& ic, + const Region& reg, + const VSpan& span); +}; + +// =========================================================================== + +/// Simple set of per-phase (named by primary component) quantities. +template +struct PhaseQuantityValue { + Scalar oil{0.0}; + Scalar gas{0.0}; + Scalar water{0.0}; + + PhaseQuantityValue& axpy(const PhaseQuantityValue& rhs, const Scalar a) + { + this->oil += a * rhs.oil; + this->gas += a * rhs.gas; + this->water += a * rhs.water; + + return *this; + } + + PhaseQuantityValue& operator/=(const Scalar x) + { + this->oil /= x; + this->gas /= x; + this->water /= x; + + return *this; + } + + void reset() + { + this->oil = this->gas = this->water = 0.0; + } +}; + +/// Calculator for phase saturations +/// +/// Computes saturation values at arbitrary depths. +/// +/// \tparam MaterialLawManager Container for material laws. Typically a +/// specialization of the \code Opm::EclMaterialLawManager<> \endcode +/// template. +/// +/// \tparam FluidSystem An OPM fluid system type. Typically a +/// specialization of the \code Opm::BlackOilFluidSystem<> \endcode +/// template. +/// +/// \tparam Region Representation of an equilibration region. Typically +/// \code Opm::EQUIL::EquilReg \endcode from the equilibrationhelpers. +/// +/// \tparam CellID Representation an equilibration region's cell IDs. +/// Typically \code std::size_t \endcode. +template +class PhaseSaturations +{ +public: + using Scalar = typename FluidSystem::Scalar; + /// Evaluation point within a model geometry. + /// + /// Associates a particular depth to specific cell. + struct Position { + CellID cell; + Scalar depth; + }; + + /// Convenience type alias + using PTable = PressureTable; + + /// Constructor + /// + /// \param[in,out] matLawMgr Read/write reference to a material law + /// container. Mutated by member functions. + /// + /// \param[in] swatInit Initial water saturation array (from SWATINIT + /// data). Empty if SWATINIT is not used in this simulation model. + explicit PhaseSaturations(MaterialLawManager& matLawMgr, + const std::vector& swatInit); + + /// Copy constructor. + /// + /// \param[in] rhs Source object. + PhaseSaturations(const PhaseSaturations& rhs); + + /// Disabled assignment operator. + PhaseSaturations& operator=(const PhaseSaturations&) = delete; + + /// Disabled move-assignment operator. + PhaseSaturations& operator=(PhaseSaturations&&) = delete; + + /// Calculate phase saturations at particular point of the simulation + /// model geometry. + /// + /// \param[in] x Specific geometric point (depth within a specific cell). + /// + /// \param[in] reg Equilibration information for a single equilibration + /// region; notably contact depths. + /// + /// \param[in] ptable Previously equilibrated phase pressure table + /// pertaining to the equilibration region \p reg. + /// + /// \return Set of phase saturation values defined at particular point. + const PhaseQuantityValue& + deriveSaturations(const Position& x, + const Region& reg, + const PTable& ptable); + + /// Retrieve saturation-corrected phase pressures + /// + /// Values associated with evaluation point of previous call to \code + /// deriveSaturations() \endcode. + const PhaseQuantityValue& correctedPhasePressures() const + { + return this->press_; + } + +private: + /// Convenience amalgamation of the deriveSaturations() input state. + /// These values are almost always used in concert. + struct EvaluationPoint { + const Position* position{nullptr}; + const Region* region {nullptr}; + const PTable* ptable {nullptr}; + }; + + /// Simplified fluid state object that contains only the pieces of + /// information needed to calculate the capillary pressure values from + /// the current set of material laws. + using FluidState = ::Opm:: + SimpleModularFluidState; + + /// Convenience type alias. + using MaterialLaw = typename MaterialLawManager::MaterialLaw; + + /// Fluid system's representation of phase indices. + using PhaseIdx = std::remove_cv_t< + std::remove_reference_t + >; + + /// Read/write reference to client's material law container. + MaterialLawManager& matLawMgr_; + + /// Client's SWATINIT data. + const std::vector& swatInit_; + + /// Evaluated phase saturations. + PhaseQuantityValue sat_; + + /// Saturation-corrected phase pressure values. + PhaseQuantityValue press_; + + /// Current evaluation point. + EvaluationPoint evalPt_; + + /// Capillary pressure fluid state. + FluidState fluidState_; + + /// Evaluated capillary pressures from current set of material laws. + std::array matLawCapPress_; + + /// Capture the input evaluation point information in internal state. + /// + /// \param[in] x Specific geometric point (depth within a specific cell). + /// + /// \param[in] reg Equilibration information for a single equilibration + /// region; notably contact depths. + /// + /// \param[in] ptable Previously equilibrated phase pressure table + /// pertaining to the equilibration region \p reg. + void setEvaluationPoint(const Position& x, + const Region& reg, + const PTable& ptable); + + /// Initialize phase saturation and phase pressure values. + /// + /// Looks up phase pressure values from the input pressure table. + void initializePhaseQuantities(); + + /// Derive phase saturation for oil. + /// + /// Calculated as 1 - Sw - Sg. + void deriveOilSat(); + + /// Derive phase saturation for gas. + /// + /// Inverts capillary pressure curve if non-constant or uses a simple + /// depth consideration with respect to G/O contact depth otherwise. + void deriveGasSat(); + + /// Derive phase saturation for water. + /// + /// Uses input data if simulation model is defined in terms of SWATINIT. + /// Otherwise, inverts capillary pressure curve if non-constant or uses + /// a simple depth consideration with respect to the O/W contact depth + /// if capillary pressure curve is constant within the current cell. + void deriveWaterSat(); + + /// Correct phase saturation and pressure values to account for + /// overlapping transition zones between G/O and O/W systems. + void fixUnphysicalTransition(); + + /// Re-adjust phase pressure values to account for phase saturations + /// outside permissible ranges. + void accountForScaledSaturations(); + + // -------------------------------------------------------------------- + // Note: Function 'applySwatInit' is non-const because the overload set + // needs to mutate the 'matLawMgr_'. + // -------------------------------------------------------------------- + + /// Derive water saturation from SWATINIT data. + /// + /// Uses SWATINIT array data from current cell directly. Also updates + /// the material law container's internal notion of the maximum + /// attainable O/W capillary pressure value. + /// + /// \param[in] pcow O/W capillary pressure value (Po - Pw). + /// + /// \return Water saturation value. + std::pair applySwatInit(const Scalar pcow); + + /// Derive water saturation from SWATINIT data. + /// + /// Uses explicitly passed-in saturation value. Also updates the + /// material law container's internal notion of the maximum attainable + /// O/W capillary pressure value. + /// + /// \param[in] pc x/W capillary pressure value (Px - Pw; x in {O, G}). + /// + /// \param[in] sw Water saturation value. + /// + /// \return Water saturation value. Input value, possibly mollified by + /// current set of material laws. + std::pair applySwatInit(const Scalar pc, const Scalar sw); + + /// Invoke material law container's capillary pressure calculator on + /// current fluid state. + void computeMaterialLawCapPress(); + + /// Extract gas/oil capillary pressure value (Pg - Po) from current + /// fluid state. + Scalar materialLawCapPressGasOil() const; + + /// Extract oil/water capillary pressure value (Po - Pw) from current + /// fluid state. + Scalar materialLawCapPressOilWater() const; + + /// Extract gas/water capillary pressure value (Pg - Pw) from current + /// fluid state. + Scalar materialLawCapPressGasWater() const; + + /// Predicate for whether specific phase has constant capillary pressure + /// curve in current cell. + /// + /// \param[in] phaseIdx Phase. Typically gas or water. + /// + /// \return Whether or not \p phaseIdx has constant capillary pressure + /// curve in current cell. + bool isConstCapPress(const PhaseIdx phaseIdx) const; + + /// Predicate for whether or not the G/O and O/W transition zones + /// overlap in the current cell. + /// + /// This is the case when inverting the capillary pressure curves + /// produces a negative oil saturation--i.e., when Sg + Sw > 1. + bool isOverlappingTransition() const; + + /// Derive phase saturation value from simple depth consideration. + /// + /// Assumes that the pertinent capillary pressure curve is constant + /// (typically zero) in the current cell--i.e., that there is a sharp + /// interface between the two phases. + /// + /// \param[in] contactdepth Depth of relevant phase separation contact. + /// + /// \param[in] Position of phase in three-phase enumeration. Typically + /// \code gasPos() \endcode or \code waterPos() \endcode. + /// + /// \param[in] isincr Whether the capillary pressure curve is normally + /// increasing as a function of phase saturation (e.g., Pcgo(Sg) = Pg + /// - Po) or if the curve is normally decreasing as a function of + /// increasing phase saturation (e.g., Pcow(Sw) = Po - Pw). True for + /// capillary pressure functions that are normally increasing as a + /// function of phase saturation. + /// + /// \return Phase saturation. + Scalar fromDepthTable(const Scalar contactdepth, + const PhaseIdx phasePos, + const bool isincr) const; + + /// Derive phase saturation by inverting non-constant capillary pressure + /// curve. + /// + /// \param[in] pc Target capillary pressure value. + /// + /// \param[in] Position of phase in three-phase enumeration. Typically + /// \code gasPos() \endcode or \code waterPos() \endcode. + /// + /// \param[in] isincr Whether the capillary pressure curve is normally + /// increasing as a function of phase saturation (e.g., Pcgo(Sg) = Pg + /// - Po) or if the curve is normally decreasing as a function of + /// increasing phase saturation (e.g., Pcow(Sw) = Po - Pw). True for + /// capillary pressure functions that are normally increasing as a + /// function of phase saturation. + /// + /// \return Phase saturation at which capillary pressure attains target + /// value. + Scalar invertCapPress(const Scalar pc, + const PhaseIdx phasePos, + const bool isincr) const; + + /// Position of oil in fluid system's three-phase enumeration. + PhaseIdx oilPos() const + { + return FluidSystem::oilPhaseIdx; + } + + /// Position of gas in fluid system's three-phase enumeration. + PhaseIdx gasPos() const + { + return FluidSystem::gasPhaseIdx; + } + + /// Position of water in fluid system's three-phase enumeration. + PhaseIdx waterPos() const + { + return FluidSystem::waterPhaseIdx; + } +}; + +// =========================================================================== + +template +void verticalExtent(const CellRange& cells, + const std::vector>& cellZMinMax, + const Parallel::Communication& comm, + std::array& span); + +template +std::pair cellZMinMax(const Element& element); + +} // namespace Details + +namespace DeckDependent { + +template +class InitialStateComputerBase +{ +protected: + using Element = typename GridView::template Codim<0>::Entity; + using Scalar = typename FluidSystem::Scalar; + using TabulatedFunction = Tabulated1DFunction; + +public: + using Vec = std::vector; + using PVec = std::vector; // One per phase. + + const Vec& temperature() const { return temperature_; } + const Vec& saltConcentration() const { return saltConcentration_; } + const Vec& saltSaturation() const { return saltSaturation_; } + const PVec& press() const { return pp_; } + const PVec& saturation() const { return sat_; } + +protected: + InitialStateComputerBase(const int numCells, + const Scalar defaultTemperature, + const CartesianIndexMapper& cartMapper, + const int num_pressure_points); + + template + void updateInitialTemperature_(const EclipseState& eclState, const RMap& reg); + + template + void updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg); + + template + void updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg); + + void updateCellProps_(const GridView& gridView, + const NumericalAquifers& aquifer); + + void applyNumericalAquifers_(const GridView& gridView, + const NumericalAquifers& aquifer, + const bool co2store_or_h2store); + + template + void setRegionPvtIdx(const EclipseState& eclState, const RMap& reg); + + std::vector tempVdTable_; + std::vector saltVdTable_; + std::vector saltpVdTable_; + std::vector regionPvtIdx_; + Vec temperature_; + Vec saltConcentration_; + Vec saltSaturation_; + PVec pp_; + PVec sat_; + const CartesianIndexMapper& cartesianIndexMapper_; + Vec swatInit_; + Vec cellCenterDepth_; + std::vector> cellCenterXY_; + std::vector> cellZSpan_; + std::vector> cellZMinMax_; + std::vector> cellCorners_; + int num_pressure_points_; +}; + +} // namespace DeckDependent +} // namespace EQUIL +} // namespace Opm + +#endif // OPM_INIT_STATE_EQUIL_BASE_HPP diff --git a/opm/simulators/flow/equil/InitStateEquilBase_impl.hpp b/opm/simulators/flow/equil/InitStateEquilBase_impl.hpp new file mode 100644 index 00000000000..b8976df605f --- /dev/null +++ b/opm/simulators/flow/equil/InitStateEquilBase_impl.hpp @@ -0,0 +1,1727 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +#ifndef OPM_INIT_STATE_EQUIL_BASE_IMPL_HPP +#define OPM_INIT_STATE_EQUIL_BASE_IMPL_HPP + +#include + +#include + +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace EQUIL { + +namespace Details { + +template +void verticalExtent(const CellRange& cells, + const std::vector>& cellZMinMax, + const Parallel::Communication& comm, + std::array& span) +{ + span[0] = std::numeric_limits::max(); + span[1] = std::numeric_limits::lowest(); + + // Define vertical span as + // + // [minimum(node depth(cells)), maximum(node depth(cells))] + // + // Note: The implementation of 'RK4IVP<>' implicitly + // imposes the requirement that cell centroids are all + // within this vertical span. That requirement is not + // checked. + for (const auto& cell : cells) { + if (cellZMinMax[cell].first < span[0]) { span[0] = cellZMinMax[cell].first; } + if (cellZMinMax[cell].second > span[1]) { span[1] = cellZMinMax[cell].second; } + } + span[0] = comm.min(span[0]); + span[1] = comm.max(span[1]); +} + +template +void subdivisionCentrePoints(const Scalar left, + const Scalar right, + const int numIntervals, + std::vector>& subdiv) +{ + const auto h = (right - left) / numIntervals; + + auto end = left; + for (auto i = 0*numIntervals; i < numIntervals; ++i) { + const auto start = end; + end = left + (i + 1)*h; + + subdiv.emplace_back((start + end) / 2, h); + } +} + +template +std::vector> +horizontalSubdivision(const CellID cell, + const std::pair topbot, + const int numIntervals) +{ + auto subdiv = std::vector>{}; + subdiv.reserve(2 * numIntervals); + + if (topbot.first > topbot.second) { + throw std::out_of_range { + "Negative thickness (inverted top/bottom faces) in cell " + + std::to_string(cell) + }; + } + + subdivisionCentrePoints(topbot.first, topbot.second, + 2*numIntervals, subdiv); + + return subdiv; +} + +template +Scalar cellCenterDepth(const Element& element) +{ + typedef typename Element::Geometry Geometry; + static constexpr int zCoord = Element::dimension - 1; + Scalar zz = 0.0; + + const Geometry& geometry = element.geometry(); + const int corners = geometry.corners(); + for (int i=0; i < corners; ++i) + zz += geometry.corner(i)[zCoord]; + + return zz/corners; +} + +template +std::pair cellCenterXY(const Element& element) +{ + typedef typename Element::Geometry Geometry; + static constexpr int xCoord = Element::dimension - 3; + static constexpr int yCoord = Element::dimension - 2; + Scalar yy = 0.0; + Scalar xx = 0.0; + + + const Geometry& geometry = element.geometry(); + const int corners = geometry.corners(); + for (int i=0; i < corners; ++i) { + xx += geometry.corner(i)[xCoord]; + yy += geometry.corner(i)[yCoord]; + } + return std::make_pair(xx/corners, yy/corners); +} + +template +std::pair cellZSpan(const Element& element) +{ + typedef typename Element::Geometry Geometry; + static constexpr int zCoord = Element::dimension - 1; + Scalar bot = 0.0; + Scalar top = 0.0; + + const Geometry& geometry = element.geometry(); + const int corners = geometry.corners(); + assert(corners == 8); + for (int i=0; i < 4; ++i) + bot += geometry.corner(i)[zCoord]; + for (int i=4; i < corners; ++i) + top += geometry.corner(i)[zCoord]; + + return std::make_pair(bot/4, top/4); +} + +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(); + + + for (int i=0; i < corners; ++i) { + 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 +void computeBlockDip(const CellCornerData& cellCorners, + Scalar& dipAngle, Scalar& dipAzimuth) +{ + const auto& Xc = cellCorners.X; + const auto& Yc = cellCorners.Y; + const auto& Zc = cellCorners.Z; + + Scalar v1x = Xc[1] - Xc[0]; + Scalar v1y = Yc[1] - Yc[0]; + Scalar v1z = Zc[1] - Zc[0]; + + Scalar v2x = Xc[2] - Xc[0]; + Scalar v2y = Yc[2] - Yc[0]; + Scalar v2z = Zc[2] - Zc[0]; + + // Cross product to get normal vector + Scalar nx = v1y * v2z - v1z * v2y; + Scalar ny = v1z * v2x - v1x * v2z; + Scalar nz = v1x * v2y - v1y * v2x; + + // Normalize the normal vector + Scalar norm = std::hypot(nx, ny, nz); + + if (norm > 1e-10) { + nx /= norm; + ny /= norm; + nz /= norm; + + // Dip angle is the angle between normal and vertical (0,0,1) + dipAngle = std::acos(std::abs(nz)); + + // Dip azimuth (direction of dip) + if (std::abs(nx) > 1e-10 || std::abs(ny) > 1e-10) { + dipAzimuth = std::atan2(ny, nx); + // Convert to 0-2π range + dipAzimuth = std::fmod(dipAzimuth + 2*std::numbers::pi_v, 2*std::numbers::pi_v); + } else { + dipAzimuth = 0.0; // Vertical cell + } + + // Clamp dip angle to reasonable values + const Scalar maxDip = std::numbers::pi_v/2 - static_cast(1e-6); + dipAngle = std::min(dipAngle, maxDip); + } else { + // Degenerate cell - assume horizontal + dipAngle = 0.0; + dipAzimuth = 0.0; + } +} + +template +CellCornerData getCellCornerXY(const Element& element) +{ + typedef typename Element::Geometry Geometry; + const Geometry& geometry = element.geometry(); + static constexpr int zCoord = Element::dimension - 1; + static constexpr int yCoord = Element::dimension - 2; + static constexpr int xCoord = Element::dimension - 3; + const int corners = geometry.corners(); + assert(corners == 8); + std::array X {}; + std::array Y {}; + std::array Z {}; + // Get all 8 corners of the hexahedral cell (maybe expensive) + for (int i = 0; i < corners; ++i) { + auto corner = geometry.corner(i); + X[i] = corner[xCoord]; + Y[i] = corner[yCoord]; + Z[i] = corner[zCoord]; + } + + return CellCornerData{X, Y, Z}; +} + +template +Scalar calculateTrueVerticalDepth(Scalar z, Scalar x, Scalar y, + Scalar dipAngle, Scalar dipAzimuth, + const std::array& referencePoint) +{ + // For True Vertical Depth calculation: + // TVD = reference_depth + (z - reference_z) * cos(dipAngle) + // + lateral_distance * sin(dipAngle) * cos(azimuth_difference) + + // Calculate lateral displacement from reference point + Scalar dx = x - referencePoint[0]; + Scalar dy = y - referencePoint[1]; + Scalar dz = z - referencePoint[2]; + + // If no dip, TVD is simply the depth + if (std::abs(dipAngle) < 1e-10) { + return referencePoint[2] + dz; + } + + // Calculate the direction from reference point to current point + Scalar pointAzimuth = std::atan2(dy, dx); + + // Calculate the angle between dip direction and point direction + Scalar azimuthDiff = pointAzimuth - dipAzimuth; + + // Calculate lateral distance + Scalar lateralDist = std::hypot(dx, dy); + + // Project lateral distance onto dip direction + Scalar lateralInDipDir = lateralDist * std::cos(azimuthDiff); + + // True Vertical Depth calculation + // TVD increases with depth (more negative z means deeper) + // For a dipping plane: TVD = vertical_component + dip_component + Scalar tvd = referencePoint[2] + dz * std::cos(dipAngle) + lateralInDipDir * std::sin(dipAngle); + + return tvd; +} + +template +RK4IVP::RK4IVP(const RHS& f, + const std::array& span, + const Scalar y0, + const int N) + : N_(N) + , span_(span) +{ + const Scalar h = stepsize(); + const Scalar h2 = h / 2; + const Scalar h6 = h / 6; + + y_.reserve(N + 1); + f_.reserve(N + 1); + + y_.push_back(y0); + f_.push_back(f(span_[0], y0)); + + for (int i = 0; i < N; ++i) { + const Scalar x = span_[0] + i*h; + const Scalar y = y_.back(); + + 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() == typename std::vector::size_type(N + 1)); +} + +template +Scalar RK4IVP:: +operator()(const Scalar x) const +{ + // Dense output (O(h**3)) according to Shampine + // (Hermite interpolation) + const Scalar h = stepsize(); + int i = (x - span_[0]) / 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 Scalar y0 = y_[i], y1 = y_[i + 1]; + const Scalar f0 = f_[i], f1 = f_[i + 1]; + + Scalar u = (1 - 2*t) * (y1 - y0); + u += h * ((t - 1)*f0 + t*f1); + u *= t * (t - 1); + u += (1 - t)*y0 + t*y1; + + return u; +} + +template +Scalar RK4IVP:: +stepsize() const +{ + return (span_[1] - span_[0]) / N_; +} + +namespace PhasePressODE { + +template +Water:: +Water(const TabulatedFunction& tempVdTable, + const TabulatedFunction& saltVdTable, + const int pvtRegionIdx, + const Scalar normGrav) + : tempVdTable_(tempVdTable) + , saltVdTable_(saltVdTable) + , pvtRegionIdx_(pvtRegionIdx) + , g_(normGrav) +{ +} + +template +typename Water::Scalar +Water:: +operator()(const Scalar depth, + const Scalar press) const +{ + return this->density(depth, press) * g_; +} + +template +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 + 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; +} + +template +Oil:: +Oil(const TabulatedFunction& tempVdTable, + const RS& rs, + const int pvtRegionIdx, + const Scalar normGrav) + : tempVdTable_(tempVdTable) + , rs_(rs) + , pvtRegionIdx_(pvtRegionIdx) + , g_(normGrav) +{ +} + +template +typename Oil::Scalar +Oil:: +operator()(const Scalar depth, + const Scalar press) const +{ + return this->density(depth, press) * g_; +} + +template +typename Oil::Scalar +Oil:: +density(const Scalar depth, + const Scalar press) const +{ + const Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true); + Scalar rs = 0.0; + if (FluidSystem::enableDissolvedGas()) + rs = rs_(depth, press, temp); + + 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); + } + Scalar rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); + if (FluidSystem::enableDissolvedGas()) { + rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + } + + return rho; +} + +template +Gas:: +Gas(const TabulatedFunction& tempVdTable, + const RV& rv, + const RVW& rvw, + const int pvtRegionIdx, + const Scalar normGrav) + : tempVdTable_(tempVdTable) + , rv_(rv) + , rvw_(rvw) + , pvtRegionIdx_(pvtRegionIdx) + , g_(normGrav) +{ +} + +template +typename Gas::Scalar +Gas:: +operator()(const Scalar depth, + const Scalar press) const +{ + return this->density(depth, press) * g_; +} + +template +typename Gas::Scalar +Gas:: +density(const Scalar depth, + const Scalar press) const +{ + const Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true); + Scalar rv = 0.0; + if (FluidSystem::enableVaporizedOil()) + rv = rv_(depth, press, temp); + + Scalar rvw = 0.0; + if (FluidSystem::enableVaporizedWater()) + rvw = rvw_(depth, press, temp); + + Scalar bGas = 0.0; + + if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) { + if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press) + && rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press)) + { + bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); + } else { + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, rvw); + } + 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; + } + + if (FluidSystem::enableVaporizedOil()){ + if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)) { + bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); + } else { + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, + temp, + press, + rv, + Scalar{0.0}/*=rvw*/); + } + Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); + return rho; + } + + if (FluidSystem::enableVaporizedWater()){ + if (rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press)) { + bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); + } + else { + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, + temp, + press, + Scalar{0.0} /*=rv*/, + rvw); + } + 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, + Scalar{0.0} /*=rv*/, + Scalar{0.0} /*=rvw*/); + Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + + return rho; +} + +} + +template +template +PressureTable:: +PressureFunction::PressureFunction(const ODE& ode, + const InitCond& ic, + const int nsample, + const VSpan& span) + : initial_(ic) +{ + this->value_[Direction::Up] = std::make_unique + (ode, VSpan {{ ic.depth, span[0] }}, ic.pressure, nsample); + + this->value_[Direction::Down] = std::make_unique + (ode, VSpan {{ ic.depth, span[1] }}, ic.pressure, nsample); +} + +template +template +PressureTable:: +PressureFunction::PressureFunction(const PressureFunction& rhs) + : initial_(rhs.initial_) +{ + this->value_[Direction::Up] = + std::make_unique(*rhs.value_[Direction::Up]); + + this->value_[Direction::Down] = + std::make_unique(*rhs.value_[Direction::Down]); +} + +template +template +typename PressureTable::template PressureFunction& +PressureTable:: +PressureFunction:: +operator=(const PressureFunction& rhs) +{ + this->initial_ = rhs.initial_; + + this->value_[Direction::Up] = + std::make_unique(*rhs.value_[Direction::Up]); + + this->value_[Direction::Down] = + std::make_unique(*rhs.value_[Direction::Down]); + + return *this; +} + +template +template +typename PressureTable::template PressureFunction& +PressureTable:: +PressureFunction:: +operator=(PressureFunction&& rhs) +{ + this->initial_ = rhs.initial_; + this->value_ = std::move(rhs.value_); + + return *this; +} + +template +template +typename PressureTable::Scalar +PressureTable:: +PressureFunction:: +value(const Scalar depth) const +{ + if (depth < this->initial_.depth) { + // Value above initial condition depth. + return (*this->value_[Direction::Up])(depth); + } + else if (depth > this->initial_.depth) { + // Value below initial condition depth. + return (*this->value_[Direction::Down])(depth); + } + else { + // Value *at* initial condition depth. + return this->initial_.pressure; + } +} + + +template +template +void PressureTable:: +checkPtr(const PressFunc* phasePress, + const std::string& phaseName) const +{ + if (phasePress != nullptr) { return; } + + throw std::invalid_argument { + "Phase pressure function for \"" + phaseName + + "\" most not be null" + }; +} + +template +typename PressureTable::Strategy +PressureTable:: +selectEquilibrationStrategy(const Region& reg) const +{ + if (!this->oilActive()) { + if (reg.datum() > reg.zwoc()) { // Datum in water zone + return &PressureTable::equil_WOG; + } + return &PressureTable::equil_GOW; + } + + if (reg.datum() > reg.zwoc()) { // Datum in water zone + return &PressureTable::equil_WOG; + } + else if (reg.datum() < reg.zgoc()) { // Datum in gas zone + return &PressureTable::equil_GOW; + } + else { // Datum in oil zone + return &PressureTable::equil_OWG; + } +} + +template +void PressureTable:: +copyInPointers(const PressureTable& rhs) +{ + if (rhs.oil_ != nullptr) { + this->oil_ = std::make_unique(*rhs.oil_); + } + + if (rhs.gas_ != nullptr) { + this->gas_ = std::make_unique(*rhs.gas_); + } + + if (rhs.wat_ != nullptr) { + this->wat_ = std::make_unique(*rhs.wat_); + } +} + +template +PhaseSaturations:: +PhaseSaturations(MaterialLawManager& matLawMgr, + const std::vector& swatInit) + : matLawMgr_(matLawMgr) + , swatInit_ (swatInit) +{ +} + +template +PhaseSaturations:: +PhaseSaturations(const PhaseSaturations& rhs) + : matLawMgr_(rhs.matLawMgr_) + , swatInit_ (rhs.swatInit_) + , sat_ (rhs.sat_) + , press_ (rhs.press_) +{ + // Note: We don't need to do anything to the 'fluidState_' here. + this->setEvaluationPoint(*rhs.evalPt_.position, + *rhs.evalPt_.region, + *rhs.evalPt_.ptable); +} + +template +const PhaseQuantityValue& +PhaseSaturations:: +deriveSaturations(const Position& x, + const Region& reg, + const PTable& ptable) +{ + this->setEvaluationPoint(x, reg, ptable); + this->initializePhaseQuantities(); + + if (ptable.gasActive()) { this->deriveGasSat(); } + + if (ptable.waterActive()) { this->deriveWaterSat(); } + + + if (this->isOverlappingTransition()) { + this->fixUnphysicalTransition(); + } + + if (ptable.oilActive()) { this->deriveOilSat(); } + + this->accountForScaledSaturations(); + + return this->sat_; +} + +template +void PhaseSaturations:: +setEvaluationPoint(const Position& x, + const Region& reg, + const PTable& ptable) +{ + this->evalPt_.position = &x; + this->evalPt_.region = ® + this->evalPt_.ptable = &ptable; +} + +template +void PhaseSaturations:: +initializePhaseQuantities() +{ + this->sat_.reset(); + this->press_.reset(); + + const auto depth = this->evalPt_.position->depth; + const auto& ptable = *this->evalPt_.ptable; + + if (ptable.oilActive()) { + this->press_.oil = ptable.oil(depth); + } + + if (ptable.gasActive()) { + this->press_.gas = ptable.gas(depth); + } + + if (ptable.waterActive()) { + this->press_.water = ptable.water(depth); + } +} + +template +void PhaseSaturations::deriveOilSat() +{ + this->sat_.oil = 1.0 - this->sat_.water - this->sat_.gas; +} + +template +void PhaseSaturations::deriveGasSat() +{ + auto& sg = this->sat_.gas; + + const auto isIncr = true; // dPcgo/dSg >= 0 for all Sg. + const auto oilActive = this->evalPt_.ptable->oilActive(); + + if (this->isConstCapPress(this->gasPos())) { + // Sharp interface between phases. Can derive phase saturation + // directly from knowing where 'depth' of evaluation point is + // relative to depth of O/G contact. + const auto gas_contact = oilActive? this->evalPt_.region->zgoc() : this->evalPt_.region->zwoc(); + sg = this->fromDepthTable(gas_contact, + this->gasPos(), isIncr); + } + else { + // Capillary pressure curve is non-constant, meaning there is a + // transition zone between the gas and oil phases. Invert capillary + // pressure relation + // + // Pcgo(Sg) = Pg - Po + // + // Note that Pcgo is defined to be (Pg - Po), not (Po - Pg). + const auto pw = oilActive? this->press_.oil : this->press_.water; + const auto pcgo = this->press_.gas - pw; + sg = this->invertCapPress(pcgo, this->gasPos(), isIncr); + } +} + +template +void PhaseSaturations::deriveWaterSat() +{ + auto& sw = this->sat_.water; + + const auto oilActive = this->evalPt_.ptable->oilActive(); + if (!oilActive) { + // for 2p gas+water we set the water saturation to 1.0 - sg + sw = 1.0 - this->sat_.gas; + } + else { + const auto isIncr = false; // dPcow/dSw <= 0 for all Sw. + + if (this->isConstCapPress(this->waterPos())) { + // Sharp interface between phases. Can derive phase saturation + // directly from knowing where 'depth' of evaluation point is + // relative to depth of O/W contact. + sw = this->fromDepthTable(this->evalPt_.region->zwoc(), + this->waterPos(), isIncr); + } + else { + // Capillary pressure curve is non-constant, meaning there is a + // transition zone between the oil and water phases. Invert + // capillary pressure relation + // + // Pcow(Sw) = Po - Pw + // + // unless the model uses "SWATINIT". In the latter case, pick the + // saturation directly from the SWATINIT array of the pertinent + // cell. + const auto pcow = this->press_.oil - this->press_.water; + + if (this->swatInit_.empty()) { + sw = this->invertCapPress(pcow, this->waterPos(), isIncr); + } + else { + auto [swout, newSwatInit] = this->applySwatInit(pcow); + if (newSwatInit) + sw = this->invertCapPress(pcow, this->waterPos(), isIncr); + else { + sw = swout; + } + } + } + } +} + +template +void PhaseSaturations:: +fixUnphysicalTransition() +{ + auto& sg = this->sat_.gas; + auto& sw = this->sat_.water; + + // Overlapping gas/oil and oil/water transition zones can lead to + // unphysical phase saturations when individual saturations are derived + // directly from inverting O/G and O/W capillary pressure curves. + // + // Recalculate phase saturations using the implied gas/water capillary + // pressure: Pg - Pw. + const auto pcgw = this->press_.gas - this->press_.water; + if (! this->swatInit_.empty()) { + // Re-scale Pc to reflect imposed sw for vanishing oil phase. This + // seems consistent with ECLIPSE, but fails to honour SWATINIT in + // case of non-trivial gas/oil capillary pressure. + auto [swout, newSwatInit] = this->applySwatInit(pcgw, sw); + if (newSwatInit){ + const auto isIncr = false; // dPcow/dSw <= 0 for all Sw. + sw = this->invertCapPress(pcgw, this->waterPos(), isIncr); + } + else { + sw = swout; + } + } + + sw = satFromSumOfPcs + (this->matLawMgr_, this->waterPos(), this->gasPos(), + this->evalPt_.position->cell, pcgw); + sg = 1.0 - sw; + + this->fluidState_.setSaturation(this->oilPos(), 1.0 - sw - sg); + this->fluidState_.setSaturation(this->gasPos(), sg); + this->fluidState_.setSaturation(this->waterPos(), this->evalPt_ + .ptable->waterActive() ? sw : 0.0); + + // Pcgo = Pg - Po => Po = Pg - Pcgo + this->computeMaterialLawCapPress(); + this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil(); +} + +template +void PhaseSaturations:: +accountForScaledSaturations() +{ + const auto gasActive = this->evalPt_.ptable->gasActive(); + const auto watActive = this->evalPt_.ptable->waterActive(); + const auto oilActive = this->evalPt_.ptable->oilActive(); + + auto sg = gasActive? this->sat_.gas : 0.0; + auto sw = watActive? this->sat_.water : 0.0; + auto so = oilActive? this->sat_.oil : 0.0; + + this->fluidState_.setSaturation(this->waterPos(), sw); + this->fluidState_.setSaturation(this->oilPos(), so); + this->fluidState_.setSaturation(this->gasPos(), sg); + + const auto& scaledDrainageInfo = this->matLawMgr_ + .oilWaterScaledEpsInfoDrainage(this->evalPt_.position->cell); + + const auto thresholdSat = 1.0e-6; + if (watActive && ((sw + thresholdSat) > scaledDrainageInfo.Swu)) { + // Water saturation exceeds maximum possible value. Reset oil phase + // pressure to that which corresponds to maximum possible water + // saturation value. + this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swu); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swu); + } else if (gasActive) { + this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swu); + } + sw = scaledDrainageInfo.Swu; + this->computeMaterialLawCapPress(); + + if (oilActive) { + // Pcow = Po - Pw => Po = Pw + Pcow + this->press_.oil = this->press_.water + this->materialLawCapPressOilWater(); + } else { + // Pcgw = Pg - Pw => Pg = Pw + Pcgw + this->press_.gas = this->press_.water + this->materialLawCapPressGasWater(); + } + + } + if (gasActive && ((sg + thresholdSat) > scaledDrainageInfo.Sgu)) { + // Gas saturation exceeds maximum possible value. Reset oil phase + // pressure to that which corresponds to maximum possible gas + // saturation value. + this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgu); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgu); + } else if (watActive) { + this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgu); + } + sg = scaledDrainageInfo.Sgu; + this->computeMaterialLawCapPress(); + + if (oilActive) { + // Pcgo = Pg - Po => Po = Pg - Pcgo + this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil(); + } else { + // Pcgw = Pg - Pw => Pw = Pg - Pcgw + this->press_.water = this->press_.gas - this->materialLawCapPressGasWater(); + } + } + + if (watActive && ((sw - thresholdSat) < scaledDrainageInfo.Swl)) { + // Water saturation less than minimum possible value in cell. Reset + // water phase pressure to that which corresponds to minimum + // possible water saturation value. + this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swl); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swl); + } else if (gasActive) { + this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swl); + } + sw = scaledDrainageInfo.Swl; + this->computeMaterialLawCapPress(); + + if (oilActive) { + // Pcwo = Po - Pw => Pw = Po - Pcow + this->press_.water = this->press_.oil - this->materialLawCapPressOilWater(); + } else { + // Pcgw = Pg - Pw => Pw = Pg - Pcgw + this->press_.water = this->press_.gas - this->materialLawCapPressGasWater(); + } + } + + if (gasActive && ((sg - thresholdSat) < scaledDrainageInfo.Sgl)) { + // Gas saturation less than minimum possible value in cell. Reset + // gas phase pressure to that which corresponds to minimum possible + // gas saturation. + this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgl); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgl); + } else if (watActive) { + this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgl); + } + sg = scaledDrainageInfo.Sgl; + this->computeMaterialLawCapPress(); + + if (oilActive) { + // Pcgo = Pg - Po => Pg = Po + Pcgo + this->press_.gas = this->press_.oil + this->materialLawCapPressGasOil(); + } else { + // Pcgw = Pg - Pw => Pg = Pw + Pcgw + this->press_.gas = this->press_.water + this->materialLawCapPressGasWater(); + } + } +} + +template +std::pair +PhaseSaturations:: +applySwatInit(const Scalar pcow) +{ + return this->applySwatInit(pcow, this->swatInit_[this->evalPt_.position->cell]); +} + +template +std::pair +PhaseSaturations:: +applySwatInit(const Scalar pcow, const Scalar sw) +{ + return this->matLawMgr_.applySwatinit(this->evalPt_.position->cell, pcow, sw); +} + +template +void PhaseSaturations:: +computeMaterialLawCapPress() +{ + const auto& matParams = this->matLawMgr_ + .materialLawParams(this->evalPt_.position->cell); + + this->matLawCapPress_.fill(0.0); + MaterialLaw::capillaryPressures(this->matLawCapPress_, + matParams, this->fluidState_); +} + +template +typename FluidSystem::Scalar +PhaseSaturations:: +materialLawCapPressGasOil() const +{ + return this->matLawCapPress_[this->oilPos()] + + this->matLawCapPress_[this->gasPos()]; +} + +template +typename FluidSystem::Scalar +PhaseSaturations:: +materialLawCapPressOilWater() const +{ + return this->matLawCapPress_[this->oilPos()] + - this->matLawCapPress_[this->waterPos()]; +} + +template +typename FluidSystem::Scalar +PhaseSaturations:: +materialLawCapPressGasWater() const +{ + return this->matLawCapPress_[this->gasPos()] + - this->matLawCapPress_[this->waterPos()]; +} + +template +bool PhaseSaturations:: +isConstCapPress(const PhaseIdx phaseIdx) const +{ + return isConstPc + (this->matLawMgr_, phaseIdx, this->evalPt_.position->cell); +} + +template +bool PhaseSaturations:: +isOverlappingTransition() const +{ + return this->evalPt_.ptable->gasActive() + && this->evalPt_.ptable->waterActive() + && ((this->sat_.gas + this->sat_.water) > 1.0); +} + +template +typename FluidSystem::Scalar +PhaseSaturations:: +fromDepthTable(const Scalar contactdepth, + const PhaseIdx phasePos, + const bool isincr) const +{ + return satFromDepth + (this->matLawMgr_, this->evalPt_.position->depth, + contactdepth, static_cast(phasePos), + this->evalPt_.position->cell, isincr); +} + +template +typename FluidSystem::Scalar +PhaseSaturations:: +invertCapPress(const Scalar pc, + const PhaseIdx phasePos, + const bool isincr) const +{ + return satFromPc + (this->matLawMgr_, static_cast(phasePos), + this->evalPt_.position->cell, pc, isincr); +} + +template +PressureTable:: +PressureTable(const Scalar gravity, + const int samplePoints) + : gravity_(gravity) + , nsample_(samplePoints) +{ +} + +template +PressureTable:: +PressureTable(const PressureTable& rhs) + : gravity_(rhs.gravity_) + , nsample_(rhs.nsample_) +{ + this->copyInPointers(rhs); +} + +template +PressureTable:: +PressureTable(PressureTable&& rhs) + : gravity_(rhs.gravity_) + , nsample_(rhs.nsample_) + , oil_ (std::move(rhs.oil_)) + , gas_ (std::move(rhs.gas_)) + , wat_ (std::move(rhs.wat_)) +{ +} + +template +PressureTable& +PressureTable:: +operator=(const PressureTable& rhs) +{ + this->gravity_ = rhs.gravity_; + this->nsample_ = rhs.nsample_; + this->copyInPointers(rhs); + + return *this; +} + +template +PressureTable& +PressureTable:: +operator=(PressureTable&& rhs) +{ + this->gravity_ = rhs.gravity_; + this->nsample_ = rhs.nsample_; + + this->oil_ = std::move(rhs.oil_); + this->gas_ = std::move(rhs.gas_); + this->wat_ = std::move(rhs.wat_); + + return *this; +} + +template +void PressureTable:: +equilibrate(const Region& reg, + const VSpan& span) +{ + // One of the PressureTable::equil_*() member functions. + auto equil = this->selectEquilibrationStrategy(reg); + + (this->*equil)(reg, span); +} + +template +bool PressureTable:: +oilActive() const +{ + return FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); +} + +template +bool PressureTable:: +gasActive() const +{ + return FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); +} + +template +bool PressureTable:: +waterActive() const +{ + return FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); +} + +template +typename FluidSystem::Scalar +PressureTable:: +oil(const Scalar depth) const +{ + this->checkPtr(this->oil_.get(), "OIL"); + + return this->oil_->value(depth); +} + +template +typename FluidSystem::Scalar +PressureTable:: +gas(const Scalar depth) const +{ + this->checkPtr(this->gas_.get(), "GAS"); + + return this->gas_->value(depth); +} + + +template +typename FluidSystem::Scalar +PressureTable:: +water(const Scalar depth) const +{ + this->checkPtr(this->wat_.get(), "WATER"); + + return this->wat_->value(depth); +} + +template +void PressureTable:: +equil_WOG(const Region& reg, const VSpan& span) +{ + // Datum depth in water zone. Calculate phase pressure for water first, + // followed by oil and gas if applicable. + + if (! this->waterActive()) { + throw std::invalid_argument { + "Don't know how to interpret EQUIL datum depth in " + "WATER zone in model without active water phase" + }; + } + + { + const auto ic = typename WPress::InitCond { + reg.datum(), reg.pressure() + }; + + this->makeWatPressure(ic, reg, span); + } + + if (this->oilActive()) { + // Pcow = Po - Pw => Po = Pw + Pcow + const auto ic = typename OPress::InitCond { + reg.zwoc(), + this->water(reg.zwoc()) + reg.pcowWoc() + }; + + this->makeOilPressure(ic, reg, span); + } + + if (this->gasActive() && this->oilActive()) { + // Pcgo = Pg - Po => Pg = Po + Pcgo + const auto ic = typename GPress::InitCond { + reg.zgoc(), + this->oil(reg.zgoc()) + reg.pcgoGoc() + }; + + this->makeGasPressure(ic, reg, span); + } else if (this->gasActive() && !this->oilActive()) { + // No oil phase set Pg = Pw + Pcgw + const auto ic = typename GPress::InitCond { + reg.zwoc(), // The WOC is really the GWC for gas/water cases + this->water(reg.zwoc()) + reg.pcowWoc() // Pcow(WOC) is really Pcgw(GWC) for gas/water cases + }; + this->makeGasPressure(ic, reg, span); + } +} + +template +void PressureTable:: +equil_GOW(const Region& reg, const VSpan& span) +{ + // Datum depth in gas zone. Calculate phase pressure for gas first, + // followed by oil and water if applicable. + + if (! this->gasActive()) { + throw std::invalid_argument { + "Don't know how to interpret EQUIL datum depth in " + "GAS zone in model without active gas phase" + }; + } + + { + const auto ic = typename GPress::InitCond { + reg.datum(), reg.pressure() + }; + + this->makeGasPressure(ic, reg, span); + } + + if (this->oilActive()) { + // Pcgo = Pg - Po => Po = Pg - Pcgo + const auto ic = typename OPress::InitCond { + reg.zgoc(), + this->gas(reg.zgoc()) - reg.pcgoGoc() + }; + this->makeOilPressure(ic, reg, span); + } + + if (this->waterActive() && this->oilActive()) { + // Pcow = Po - Pw => Pw = Po - Pcow + const auto ic = typename WPress::InitCond { + reg.zwoc(), + this->oil(reg.zwoc()) - reg.pcowWoc() + }; + + this->makeWatPressure(ic, reg, span); + } else if (this->waterActive() && !this->oilActive()) { + // No oil phase set Pw = Pg - Pcgw + const auto ic = typename WPress::InitCond { + reg.zwoc(), // The WOC is really the GWC for gas/water cases + this->gas(reg.zwoc()) - reg.pcowWoc() // Pcow(WOC) is really Pcgw(GWC) for gas/water cases + }; + this->makeWatPressure(ic, reg, span); + } +} + +template +void PressureTable:: +equil_OWG(const Region& reg, const VSpan& span) +{ + // Datum depth in oil zone. Calculate phase pressure for oil first, + // followed by gas and water if applicable. + + if (! this->oilActive()) { + throw std::invalid_argument { + "Don't know how to interpret EQUIL datum depth in " + "OIL zone in model without active oil phase" + }; + } + + { + const auto ic = typename OPress::InitCond { + reg.datum(), reg.pressure() + }; + + this->makeOilPressure(ic, reg, span); + } + + if (this->waterActive()) { + // Pcow = Po - Pw => Pw = Po - Pcow + const auto ic = typename WPress::InitCond { + reg.zwoc(), + this->oil(reg.zwoc()) - reg.pcowWoc() + }; + + this->makeWatPressure(ic, reg, span); + } + + if (this->gasActive()) { + // Pcgo = Pg - Po => Pg = Po + Pcgo + const auto ic = typename GPress::InitCond { + reg.zgoc(), + this->oil(reg.zgoc()) + reg.pcgoGoc() + }; + this->makeGasPressure(ic, reg, span); + } +} + +template +void PressureTable:: +makeOilPressure(const typename OPress::InitCond& ic, + const Region& reg, + const VSpan& span) +{ + const auto drho = OilPressODE { + reg.tempVdTable(), reg.dissolutionCalculator(), + reg.pvtIdx(), this->gravity_ + }; + + this->oil_ = std::make_unique(drho, ic, this->nsample_, span); +} + +template +void PressureTable:: +makeGasPressure(const typename GPress::InitCond& ic, + const Region& reg, + const VSpan& span) +{ + const auto drho = GasPressODE { + reg.tempVdTable(), reg.evaporationCalculator(), reg.waterEvaporationCalculator(), + reg.pvtIdx(), this->gravity_ + }; + + this->gas_ = std::make_unique(drho, ic, this->nsample_, span); +} + +template +void PressureTable:: +makeWatPressure(const typename WPress::InitCond& ic, + const Region& reg, + const VSpan& span) +{ + const auto drho = WatPressODE { + reg.tempVdTable(), reg.saltVdTable(), reg.pvtIdx(), this->gravity_ + }; + + this->wat_ = std::make_unique(drho, ic, this->nsample_, span); +} + +} + +namespace DeckDependent { + +std::vector +getEquil(const EclipseState& state) +{ + const auto& init = state.getInitConfig(); + + if(!init.hasEquil()) { + throw std::domain_error("Deck does not provide equilibration data."); + } + + const auto& equil = init.getEquil(); + return { equil.begin(), equil.end() }; +} + +template +std::vector +equilnum(const EclipseState& eclipseState, + const GridView& gridview) +{ + std::vector eqlnum(gridview.size(0), 0); + + if (eclipseState.fieldProps().has_int("EQLNUM")) { + const auto& e = eclipseState.fieldProps().get_int("EQLNUM"); + std::ranges::transform(e, eqlnum.begin(), [](int n) { return n - 1; }); + } + OPM_BEGIN_PARALLEL_TRY_CATCH(); + const int num_regions = eclipseState.getTableManager().getEqldims().getNumEquilRegions(); + if (std::ranges::any_of(eqlnum, [num_regions](int n){return n >= num_regions;})) { + throw std::runtime_error("Values larger than maximum Equil regions " + + std::to_string(num_regions) + " provided in EQLNUM"); + } + if (std::ranges::any_of(eqlnum, [](int n){return n < 0;})) { + throw std::runtime_error("zero or negative values provided in EQLNUM"); + } + OPM_END_PARALLEL_TRY_CATCH("Invalied EQLNUM numbers: ", gridview.comm()); + + return eqlnum; +} + +template +InitialStateComputerBase:: +InitialStateComputerBase(const int numCells, + const Scalar defaultTemperature, + const CartesianIndexMapper& cartMapper, + const int num_pressure_points) + : temperature_(numCells, defaultTemperature) + , saltConcentration_(numCells) + , saltSaturation_(numCells) + , pp_(FluidSystem::numPhases, std::vector(numCells)) + , sat_(FluidSystem::numPhases, std::vector(numCells)) + , cartesianIndexMapper_(cartMapper) + , num_pressure_points_(num_pressure_points) +{ +} + +template +template +void InitialStateComputerBase:: +updateInitialTemperature_(const EclipseState& eclState, const RMap& reg) +{ + const int numEquilReg = regionPvtIdx_.size(); + tempVdTable_.resize(numEquilReg); + const auto& tables = eclState.getTableManager(); + if (!tables.hasTables("RTEMPVD")) { + 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); + } + } else { + const TableContainer& tempvdTables = tables.getRtempvdTables(); + for (std::size_t i = 0; i < tempvdTables.size(); ++i) { + const RtempvdTable& tempvdTable = tempvdTables.getTable(i); + tempVdTable_[i].setXYContainers(tempvdTable.getDepthColumn(), tempvdTable.getTemperatureColumn()); + const auto& cells = reg.cells(i); + for (const auto& cell : cells) { + const Scalar depth = cellCenterDepth_[cell]; + this->temperature_[cell] = tempVdTable_[i].eval(depth, /*extrapolate=*/true); + } + } + } +} + +template +template +void InitialStateComputerBase:: +updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg) +{ + const int numEquilReg = regionPvtIdx_.size(); + saltVdTable_.resize(numEquilReg); + const auto& tables = eclState.getTableManager(); + const TableContainer& saltvdTables = tables.getSaltvdTables(); + + // 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}; + for (auto& table : this->saltVdTable_) { + table.setXYContainers(x, y); + } + } else { + for (std::size_t i = 0; i < saltvdTables.size(); ++i) { + const SaltvdTable& saltvdTable = saltvdTables.getTable(i); + saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn()); + + const auto& cells = reg.cells(i); + for (const auto& cell : cells) { + const Scalar depth = cellCenterDepth_[cell]; + this->saltConcentration_[cell] = saltVdTable_[i].eval(depth, /*extrapolate=*/true); + } + } + } +} + +template +template +void InitialStateComputerBase:: +updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg) +{ + const int numEquilReg = regionPvtIdx_.size(); + saltpVdTable_.resize(numEquilReg); + const auto& tables = eclState.getTableManager(); + const TableContainer& saltpvdTables = tables.getSaltpvdTables(); + + for (std::size_t i = 0; i < saltpvdTables.size(); ++i) { + const SaltpvdTable& saltpvdTable = saltpvdTables.getTable(i); + saltpVdTable_[i].setXYContainers(saltpvdTable.getDepthColumn(), saltpvdTable.getSaltpColumn()); + + const auto& cells = reg.cells(i); + for (const auto& cell : cells) { + const Scalar depth = cellCenterDepth_[cell]; + this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth, /*extrapolate=*/true); + } + } +} + +template +void InitialStateComputerBase:: +updateCellProps_(const GridView& gridView, + const NumericalAquifers& aquifer) +{ + ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); + int numElements = gridView.size(/*codim=*/0); + cellCenterDepth_.resize(numElements); + cellCenterXY_.resize(numElements); + cellCorners_.resize(numElements); + cellZSpan_.resize(numElements); + cellZMinMax_.resize(numElements); + + auto elemIt = gridView.template begin(); + const auto& elemEndIt = gridView.template end(); + const auto num_aqu_cells = aquifer.allAquiferCells(); + for (; elemIt != elemEndIt; ++elemIt) { + const Element& element = *elemIt; + const unsigned int elemIdx = elemMapper.index(element); + cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element); + cellCenterXY_[elemIdx] = Details::cellCenterXY(element); + cellCorners_[elemIdx] = Details::getCellCornerXY(element); + const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); + 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 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; + cellZMinMax_[elemIdx].first += depth_change_num_aqu; + cellZMinMax_[elemIdx].second += depth_change_num_aqu; + } + } + } +} + +template +void InitialStateComputerBase:: +applyNumericalAquifers_(const GridView& gridView, + const NumericalAquifers& aquifer, + const bool co2store_or_h2store) +{ + const auto num_aqu_cells = aquifer.allAquiferCells(); + if (num_aqu_cells.empty()) return; + + // Check if water phase is active, or in the case of CO2STORE and H2STORE, water is modelled as oil phase + bool oil_as_brine = co2store_or_h2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + const auto watPos = oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx; + if (!FluidSystem::phaseIsActive(watPos)){ + throw std::logic_error { "Water phase has to be active for numerical aquifer case" }; + } + + ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); + auto elemIt = gridView.template begin(); + const auto& elemEndIt = gridView.template end(); + const auto oilPos = FluidSystem::oilPhaseIdx; + const auto gasPos = FluidSystem::gasPhaseIdx; + for (; elemIt != elemEndIt; ++elemIt) { + const Element& element = *elemIt; + const unsigned int elemIdx = elemMapper.index(element); + const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); + const auto search = num_aqu_cells.find(cartIx); + if (search != num_aqu_cells.end()) { + // numerical aquifer cells are filled with water initially + this->sat_[watPos][elemIdx] = 1.; + + if (!co2store_or_h2store && FluidSystem::phaseIsActive(oilPos)) { + this->sat_[oilPos][elemIdx] = 0.; + } + + if (FluidSystem::phaseIsActive(gasPos)) { + this->sat_[gasPos][elemIdx] = 0.; + } + const auto* aqu_cell = num_aqu_cells.at(cartIx); + const auto msg = fmt::format("FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL " + "AQUIFER {}, WATER SATURATION IS SET TO BE UNITY", + aqu_cell->I+1, aqu_cell->J+1, aqu_cell->K+1, aqu_cell->aquifer_id); + OpmLog::info(msg); + + // if pressure is specified for numerical aquifers, we use these pressure values + // for numerical aquifer cells + if (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; + } + if (FluidSystem::phaseIsActive(oilPos)) { + this->pp_[oilPos][elemIdx] = pres; + } + } + } + } +} + +template +template +void InitialStateComputerBase:: +setRegionPvtIdx(const EclipseState& eclState, const RMap& reg) +{ + const auto& pvtnumData = eclState.fieldProps().get_int("PVTNUM"); + + for (const auto& r : reg.activeRegions()) { + const auto& cells = reg.cells(r); + regionPvtIdx_[r] = pvtnumData[*cells.begin()] - 1; + } +} + +} +} // namespace EQUIL +} // namespace Opm + +#endif // OPM_INIT_STATE_EQUIL_BASE_IMPL_HPP diff --git a/opm/simulators/flow/equil/InitStateEquil_impl.hpp b/opm/simulators/flow/equil/InitStateEquil_impl.hpp index 92386a99623..2bab1086dcb 100644 --- a/opm/simulators/flow/equil/InitStateEquil_impl.hpp +++ b/opm/simulators/flow/equil/InitStateEquil_impl.hpp @@ -23,6 +23,9 @@ #ifndef OPM_INIT_STATE_EQUIL_IMPL_HPP #define OPM_INIT_STATE_EQUIL_IMPL_HPP +#include +#include + #include #include @@ -1423,46 +1426,10 @@ makeWatPressure(const typename WPress::InitCond& ic, } +======= +>>>>>>> 02e7fa473 (refactoring InitStateEquil) namespace DeckDependent { -std::vector -getEquil(const EclipseState& state) -{ - const auto& init = state.getInitConfig(); - - if(!init.hasEquil()) { - throw std::domain_error("Deck does not provide equilibration data."); - } - - const auto& equil = init.getEquil(); - return { equil.begin(), equil.end() }; -} - -template -std::vector -equilnum(const EclipseState& eclipseState, - const GridView& gridview) -{ - std::vector eqlnum(gridview.size(0), 0); - - if (eclipseState.fieldProps().has_int("EQLNUM")) { - const auto& e = eclipseState.fieldProps().get_int("EQLNUM"); - std::ranges::transform(e, eqlnum.begin(), [](int n) { return n - 1; }); - } - OPM_BEGIN_PARALLEL_TRY_CATCH(); - const int num_regions = eclipseState.getTableManager().getEqldims().getNumEquilRegions(); - if (std::ranges::any_of(eqlnum, [num_regions](int n){return n >= num_regions;})) { - throw std::runtime_error("Values larger than maximum Equil regions " + - std::to_string(num_regions) + " provided in EQLNUM"); - } - if (std::ranges::any_of(eqlnum, [](int n){return n < 0;})) { - throw std::runtime_error("zero or negative values provided in EQLNUM"); - } - OPM_END_PARALLEL_TRY_CATCH("Invalied EQLNUM numbers: ", gridview.comm()); - - return eqlnum; -} - template(grid.size(/*codim=*/0))), - sat_(FluidSystem::numPhases, - std::vector(grid.size(/*codim=*/0))), - rs_(grid.size(/*codim=*/0)), - rv_(grid.size(/*codim=*/0)), - rvw_(grid.size(/*codim=*/0)), - cartesianIndexMapper_(cartMapper), - num_pressure_points_(num_pressure_points) + : Base(grid.size(/*codim=*/0), + eclipseState.getTableManager().rtemp(), + cartMapper, + num_pressure_points) + , rs_(grid.size(/*codim=*/0)) + , rv_(grid.size(/*codim=*/0)) + , rvw_(grid.size(/*codim=*/0)) { //Check for presence of kw SWATINIT if (applySwatInit) { if (eclipseState.fieldProps().has_double("SWATINIT")) { if constexpr (std::is_same_v) { - swatInit_ = eclipseState.fieldProps().get_double("SWATINIT"); + this->swatInit_ = eclipseState.fieldProps().get_double("SWATINIT"); } else { const auto& input = eclipseState.fieldProps().get_double("SWATINIT"); - swatInit_.resize(input.size()); - std::ranges::copy(input, swatInit_.begin()); + this->swatInit_.resize(input.size()); + std::ranges::copy(input, this->swatInit_.begin()); } } } @@ -1511,7 +1473,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager, // Querry cell depth, cell top-bottom. // numerical aquifer cells might be specified with different depths. const auto& num_aquifers = eclipseState.aquifer().numericalAquifers(); - updateCellProps_(gridView, num_aquifers); + this->updateCellProps_(gridView, num_aquifers); // Get the equilibration records. const std::vector rec = getEquil(eclipseState); @@ -1519,8 +1481,8 @@ InitialStateComputer(MaterialLawManager& materialLawManager, // Create (inverse) region mapping. const RegionMapping<> eqlmap(equilnum(eclipseState, grid)); const int invalidRegion = -1; - regionPvtIdx_.resize(rec.size(), invalidRegion); - setRegionPvtIdx(eclipseState, eqlmap); + this->regionPvtIdx_.resize(rec.size(), invalidRegion); + this->setRegionPvtIdx(eclipseState, eqlmap); // Create Rs functions. rsFunc_.reserve(rec.size()); @@ -1543,7 +1505,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager, rsFunc_.push_back(std::shared_ptr>()); continue; } - const int pvtIdx = regionPvtIdx_[i]; + const int pvtIdx = this->regionPvtIdx_[i]; if (!rec[i].liveOilInitConstantRs()) { const TableContainer& rsvdTables = tables.getRsvdTables(); const TableContainer& pbvdTables = tables.getPbvdTables(); @@ -1619,7 +1581,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager, rvFunc_.push_back(std::shared_ptr>()); continue; } - const int pvtIdx = regionPvtIdx_[i]; + const int pvtIdx = this->regionPvtIdx_[i]; if (!rec[i].wetGasInitConstantRv()) { const TableContainer& rvvdTables = tables.getRvvdTables(); const TableContainer& pdvdTables = tables.getPdvdTables(); @@ -1666,7 +1628,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager, rvwFunc_.push_back(std::shared_ptr>()); continue; } - const int pvtIdx = regionPvtIdx_[i]; + const int pvtIdx = this->regionPvtIdx_[i]; if (!rec[i].humidGasInitConstantRvw()) { const TableContainer& rvwvdTables = tables.getRvwvdTables(); @@ -1723,269 +1685,27 @@ InitialStateComputer(MaterialLawManager& materialLawManager, } // EXTRACT the initial temperature - updateInitialTemperature_(eclipseState, eqlmap); + this->updateInitialTemperature_(eclipseState, eqlmap); // EXTRACT the initial salt concentration - updateInitialSaltConcentration_(eclipseState, eqlmap); + this->updateInitialSaltConcentration_(eclipseState, eqlmap); // EXTRACT the initial salt saturation - updateInitialSaltSaturation_(eclipseState, eqlmap); + this->updateInitialSaltSaturation_(eclipseState, eqlmap); // Compute pressures, saturations, rs and rv factors. const auto& comm = grid.comm(); calcPressSatRsRv(eqlmap, rec, materialLawManager, gridView, comm, grav); // modify the pressure and saturation for numerical aquifer cells - applyNumericalAquifers_(gridView, num_aquifers, - eclipseState.runspec().co2Storage() || - eclipseState.runspec().h2Storage()); + this->applyNumericalAquifers_(gridView, num_aquifers, + eclipseState.runspec().co2Storage() || + eclipseState.runspec().h2Storage()); // Modify oil pressure in no-oil regions so that the pressures of present phases can // be recovered from the oil pressure and capillary relations. } -template -template -void InitialStateComputer:: -updateInitialTemperature_(const EclipseState& eclState, const RMap& reg) -{ - const int numEquilReg = rsFunc_.size(); - tempVdTable_.resize(numEquilReg); - const auto& tables = eclState.getTableManager(); - if (!tables.hasTables("RTEMPVD")) { - 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); - } - } else { - const TableContainer& tempvdTables = tables.getRtempvdTables(); - for (std::size_t i = 0; i < tempvdTables.size(); ++i) { - const RtempvdTable& tempvdTable = tempvdTables.getTable(i); - tempVdTable_[i].setXYContainers(tempvdTable.getDepthColumn(), tempvdTable.getTemperatureColumn()); - const auto& cells = reg.cells(i); - for (const auto& cell : cells) { - const Scalar depth = cellCenterDepth_[cell]; - this->temperature_[cell] = tempVdTable_[i].eval(depth, /*extrapolate=*/true); - } - } - } -} - -template -template -void InitialStateComputer:: -updateInitialSaltConcentration_(const EclipseState& eclState, const RMap& reg) -{ - const int numEquilReg = rsFunc_.size(); - saltVdTable_.resize(numEquilReg); - const auto& tables = eclState.getTableManager(); - const TableContainer& saltvdTables = tables.getSaltvdTables(); - - // 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}; - for (auto& table : this->saltVdTable_) { - table.setXYContainers(x, y); - } - } else { - for (std::size_t i = 0; i < saltvdTables.size(); ++i) { - const SaltvdTable& saltvdTable = saltvdTables.getTable(i); - saltVdTable_[i].setXYContainers(saltvdTable.getDepthColumn(), saltvdTable.getSaltColumn()); - - const auto& cells = reg.cells(i); - for (const auto& cell : cells) { - const Scalar depth = cellCenterDepth_[cell]; - this->saltConcentration_[cell] = saltVdTable_[i].eval(depth, /*extrapolate=*/true); - } - } - } -} - -template -template -void InitialStateComputer:: -updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg) -{ - const int numEquilReg = rsFunc_.size(); - saltpVdTable_.resize(numEquilReg); - const auto& tables = eclState.getTableManager(); - const TableContainer& saltpvdTables = tables.getSaltpvdTables(); - - for (std::size_t i = 0; i < saltpvdTables.size(); ++i) { - const SaltpvdTable& saltpvdTable = saltpvdTables.getTable(i); - saltpVdTable_[i].setXYContainers(saltpvdTable.getDepthColumn(), saltpvdTable.getSaltpColumn()); - - const auto& cells = reg.cells(i); - for (const auto& cell : cells) { - const Scalar depth = cellCenterDepth_[cell]; - this->saltSaturation_[cell] = saltpVdTable_[i].eval(depth, /*extrapolate=*/true); - } - } -} - -template -void InitialStateComputer:: -updateCellProps_(const GridView& gridView, - const NumericalAquifers& aquifer) -{ - ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); - int numElements = gridView.size(/*codim=*/0); - cellCenterDepth_.resize(numElements); - cellCenterXY_.resize(numElements); - cellCorners_.resize(numElements); - cellZSpan_.resize(numElements); - cellZMinMax_.resize(numElements); - - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - const auto num_aqu_cells = aquifer.allAquiferCells(); - for (; elemIt != elemEndIt; ++elemIt) { - const Element& element = *elemIt; - const unsigned int elemIdx = elemMapper.index(element); - cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element); - cellCenterXY_[elemIdx] = Details::cellCenterXY(element); - cellCorners_[elemIdx] = Details::getCellCornerXY(element); - const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); - 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 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; - cellZMinMax_[elemIdx].first += depth_change_num_aqu; - cellZMinMax_[elemIdx].second += depth_change_num_aqu; - } - } - } -} - -template -void InitialStateComputer:: -applyNumericalAquifers_(const GridView& gridView, - const NumericalAquifers& aquifer, - const bool co2store_or_h2store) -{ - const auto num_aqu_cells = aquifer.allAquiferCells(); - if (num_aqu_cells.empty()) return; - - // Check if water phase is active, or in the case of CO2STORE and H2STORE, water is modelled as oil phase - bool oil_as_brine = co2store_or_h2store && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); - const auto watPos = oil_as_brine? FluidSystem::oilPhaseIdx : FluidSystem::waterPhaseIdx; - if (!FluidSystem::phaseIsActive(watPos)){ - throw std::logic_error { "Water phase has to be active for numerical aquifer case" }; - } - - ElementMapper elemMapper(gridView, Dune::mcmgElementLayout()); - auto elemIt = gridView.template begin(); - const auto& elemEndIt = gridView.template end(); - const auto oilPos = FluidSystem::oilPhaseIdx; - const auto gasPos = FluidSystem::gasPhaseIdx; - for (; elemIt != elemEndIt; ++elemIt) { - const Element& element = *elemIt; - const unsigned int elemIdx = elemMapper.index(element); - const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx); - const auto search = num_aqu_cells.find(cartIx); - if (search != num_aqu_cells.end()) { - // numerical aquifer cells are filled with water initially - this->sat_[watPos][elemIdx] = 1.; - - if (!co2store_or_h2store && FluidSystem::phaseIsActive(oilPos)) { - this->sat_[oilPos][elemIdx] = 0.; - } - - if (FluidSystem::phaseIsActive(gasPos)) { - this->sat_[gasPos][elemIdx] = 0.; - } - const auto* aqu_cell = num_aqu_cells.at(cartIx); - const auto msg = fmt::format("FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL " - "AQUIFER {}, WATER SATURATION IS SET TO BE UNITY", - aqu_cell->I+1, aqu_cell->J+1, aqu_cell->K+1, aqu_cell->aquifer_id); - OpmLog::info(msg); - - // if pressure is specified for numerical aquifers, we use these pressure values - // for numerical aquifer cells - if (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; - } - if (FluidSystem::phaseIsActive(oilPos)) { - this->pp_[oilPos][elemIdx] = pres; - } - } - } - } -} - -template -template -void InitialStateComputer:: -setRegionPvtIdx(const EclipseState& eclState, const RMap& reg) -{ - const auto& pvtnumData = eclState.fieldProps().get_int("PVTNUM"); - - for (const auto& r : reg.activeRegions()) { - const auto& cells = reg.cells(r); - regionPvtIdx_[r] = pvtnumData[*cells.begin()] - 1; - } -} - templatecellZMinMax_, comm, vspan); const auto acc = rec[r].initializationTargetAccuracy(); if (acc > 0) { @@ -2164,7 +1884,7 @@ equilibrateCellCentres(const CellRange& cells, Scalar& Rvw) -> void { const auto pos = CellPos { - cell, cellCenterDepth_[cell] + cell, this->cellCenterDepth_[cell] }; saturations = psat.deriveSaturations(pos, eqreg, ptable); @@ -2216,7 +1936,7 @@ equilibrateHorizontal(const CellRange& cells, saturations.reset(); Scalar totfrac = 0.0; - for (const auto& [depth, frac] : Details::horizontalSubdivision(cell, cellZSpan_[cell], acc)) { + for (const auto& [depth, frac] : Details::horizontalSubdivision(cell, this->cellZSpan_[cell], acc)) { const auto pos = CellPos { cell, depth }; saturations.axpy(psat.deriveSaturations(pos, eqreg, ptable), frac); @@ -2231,7 +1951,7 @@ equilibrateHorizontal(const CellRange& cells, } else { // Fall back to centre point method for zero-thickness cells. const auto pos = CellPos { - cell, cellCenterDepth_[cell] + cell, this->cellCenterDepth_[cell] }; saturations = psat.deriveSaturations(pos, eqreg, ptable); @@ -2239,7 +1959,7 @@ equilibrateHorizontal(const CellRange& cells, } const auto temp = this->temperature_[cell]; - const auto cz = cellCenterDepth_[cell]; + const auto cz = this->cellCenterDepth_[cell]; Rs = eqreg.dissolutionCalculator() (cz, pressures.oil, temp, saturations.gas); @@ -2279,19 +1999,19 @@ equilibrateTiltedFaultBlockSimple(const CellRange& cells, Scalar totalWeight = 0.0; // We assume grid blocks are treated as being tilted - const auto& [zmin, zmax] = cellZMinMax_[cell]; + const auto& [zmin, zmax] = this->cellZMinMax_[cell]; const Scalar cellThickness = zmax - zmin; const Scalar halfThickness = cellThickness / 2.0; // Calculate dip parameters from corner point geometry Scalar dipAngle, dipAzimuth; - Details::computeBlockDip(cellCorners_[cell], dipAngle, dipAzimuth); + Details::computeBlockDip(this->cellCorners_[cell], dipAngle, dipAzimuth); // Reference point for TVD calculations std::array referencePoint = { - cellCenterXY_[cell].first, - cellCenterXY_[cell].second, - cellCenterDepth_[cell] + this->cellCenterXY_[cell].first, + this->cellCenterXY_[cell].second, + this->cellCenterDepth_[cell] }; // We have acc levels within each half (upper and lower) of the block @@ -2323,7 +2043,7 @@ equilibrateTiltedFaultBlockSimple(const CellRange& cells, for (const auto& [depth, weight] : levels) { // Convert measured depth to True Vertical Depth for tilted blocks - const auto& [x, y] = cellCenterXY_[cell]; + const auto& [x, y] = this->cellCenterXY_[cell]; Scalar tvd = Details::calculateTrueVerticalDepth( depth, x, y, dipAngle, dipAzimuth, referencePoint); @@ -2344,9 +2064,9 @@ equilibrateTiltedFaultBlockSimple(const CellRange& cells, pressures /= totalWeight; } else { // Fallback to center point method using TVD - const auto& [x, y] = cellCenterXY_[cell]; + const auto& [x, y] = this->cellCenterXY_[cell]; Scalar tvdCenter = Details::calculateTrueVerticalDepth( - cellCenterDepth_[cell], x, y, dipAngle, dipAzimuth, referencePoint); + this->cellCenterDepth_[cell], x, y, dipAngle, dipAzimuth, referencePoint); const auto pos = CellPos{cell, tvdCenter}; saturations = psat.deriveSaturations(pos, eqreg, ptable); pressures = psat.correctedPhasePressures(); @@ -2354,9 +2074,9 @@ equilibrateTiltedFaultBlockSimple(const CellRange& cells, // Compute solution ratios at cell center TVD const auto temp = this->temperature_[cell]; - const auto& [x, y] = cellCenterXY_[cell]; + const auto& [x, y] = this->cellCenterXY_[cell]; Scalar tvdCenter = Details::calculateTrueVerticalDepth( - cellCenterDepth_[cell], x, y, dipAngle, dipAzimuth, referencePoint); + this->cellCenterDepth_[cell], x, y, dipAngle, dipAzimuth, referencePoint); Rs = eqreg.dissolutionCalculator()(tvdCenter, pressures.oil, temp, saturations.gas); Rv = eqreg.evaporationCalculator()(tvdCenter, pressures.gas, temp, saturations.oil); @@ -2492,7 +2212,7 @@ equilibrateTiltedFaultBlock(const CellRange& cells, std::array referencePoint = { this->cellCenterXY_[cell].first, this->cellCenterXY_[cell].second, - cellCenterDepth_[cell] + this->cellCenterDepth_[cell] }; // We have acc levels within each half (upper and lower) of the block From bc77a993b6392dcf4d1cc1f1e9ae539a7e84493c Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 14 Apr 2026 15:24:43 +0200 Subject: [PATCH 2/2] fixing the mistakes during rebasing --- .../flow/equil/InitStateEquil_impl.hpp | 1362 ----------------- 1 file changed, 1362 deletions(-) diff --git a/opm/simulators/flow/equil/InitStateEquil_impl.hpp b/opm/simulators/flow/equil/InitStateEquil_impl.hpp index 2bab1086dcb..006ba570e39 100644 --- a/opm/simulators/flow/equil/InitStateEquil_impl.hpp +++ b/opm/simulators/flow/equil/InitStateEquil_impl.hpp @@ -66,1368 +66,6 @@ namespace Opm { namespace EQUIL { -namespace Details { - -template -void verticalExtent(const CellRange& cells, - const std::vector>& cellZMinMax, - const Parallel::Communication& comm, - std::array& span) -{ - span[0] = std::numeric_limits::max(); - span[1] = std::numeric_limits::lowest(); - - // Define vertical span as - // - // [minimum(node depth(cells)), maximum(node depth(cells))] - // - // Note: The implementation of 'RK4IVP<>' implicitly - // imposes the requirement that cell centroids are all - // within this vertical span. That requirement is not - // checked. - for (const auto& cell : cells) { - if (cellZMinMax[cell].first < span[0]) { span[0] = cellZMinMax[cell].first; } - if (cellZMinMax[cell].second > span[1]) { span[1] = cellZMinMax[cell].second; } - } - span[0] = comm.min(span[0]); - span[1] = comm.max(span[1]); -} - -template -void subdivisionCentrePoints(const Scalar left, - const Scalar right, - const int numIntervals, - std::vector>& subdiv) -{ - const auto h = (right - left) / numIntervals; - - auto end = left; - for (auto i = 0*numIntervals; i < numIntervals; ++i) { - const auto start = end; - end = left + (i + 1)*h; - - subdiv.emplace_back((start + end) / 2, h); - } -} - -template -std::vector> -horizontalSubdivision(const CellID cell, - const std::pair topbot, - const int numIntervals) -{ - auto subdiv = std::vector>{}; - subdiv.reserve(2 * numIntervals); - - if (topbot.first > topbot.second) { - throw std::out_of_range { - "Negative thickness (inverted top/bottom faces) in cell " - + std::to_string(cell) - }; - } - - subdivisionCentrePoints(topbot.first, topbot.second, - 2*numIntervals, subdiv); - - return subdiv; -} - -template -Scalar cellCenterDepth(const Element& element) -{ - typedef typename Element::Geometry Geometry; - static constexpr int zCoord = Element::dimension - 1; - Scalar zz = 0.0; - - const Geometry& geometry = element.geometry(); - const int corners = geometry.corners(); - for (int i=0; i < corners; ++i) - zz += geometry.corner(i)[zCoord]; - - return zz/corners; -} - -template -std::pair cellCenterXY(const Element& element) -{ - typedef typename Element::Geometry Geometry; - static constexpr int xCoord = Element::dimension - 3; - static constexpr int yCoord = Element::dimension - 2; - Scalar yy = 0.0; - Scalar xx = 0.0; - - - const Geometry& geometry = element.geometry(); - const int corners = geometry.corners(); - for (int i=0; i < corners; ++i) { - xx += geometry.corner(i)[xCoord]; - yy += geometry.corner(i)[yCoord]; - } - return std::make_pair(xx/corners, yy/corners); -} - -template -std::pair cellZSpan(const Element& element) -{ - typedef typename Element::Geometry Geometry; - static constexpr int zCoord = Element::dimension - 1; - Scalar bot = 0.0; - Scalar top = 0.0; - - const Geometry& geometry = element.geometry(); - const int corners = geometry.corners(); - assert(corners == 8); - for (int i=0; i < 4; ++i) - bot += geometry.corner(i)[zCoord]; - for (int i=4; i < corners; ++i) - top += geometry.corner(i)[zCoord]; - - return std::make_pair(bot/4, top/4); -} - -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(); - - - for (int i=0; i < corners; ++i) { - 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 -void computeBlockDip(const CellCornerData& cellCorners, - Scalar& dipAngle, Scalar& dipAzimuth) -{ - const auto& Xc = cellCorners.X; - const auto& Yc = cellCorners.Y; - const auto& Zc = cellCorners.Z; - - Scalar v1x = Xc[1] - Xc[0]; - Scalar v1y = Yc[1] - Yc[0]; - Scalar v1z = Zc[1] - Zc[0]; - - Scalar v2x = Xc[2] - Xc[0]; - Scalar v2y = Yc[2] - Yc[0]; - Scalar v2z = Zc[2] - Zc[0]; - - // Cross product to get normal vector - Scalar nx = v1y * v2z - v1z * v2y; - Scalar ny = v1z * v2x - v1x * v2z; - Scalar nz = v1x * v2y - v1y * v2x; - - // Normalize the normal vector - Scalar norm = std::hypot(nx, ny, nz); - - if (norm > 1e-10) { - nx /= norm; - ny /= norm; - nz /= norm; - - // Dip angle is the angle between normal and vertical (0,0,1) - dipAngle = std::acos(std::abs(nz)); - - // Dip azimuth (direction of dip) - if (std::abs(nx) > 1e-10 || std::abs(ny) > 1e-10) { - dipAzimuth = std::atan2(ny, nx); - // Convert to 0-2π range - dipAzimuth = std::fmod(dipAzimuth + 2*std::numbers::pi_v, 2*std::numbers::pi_v); - } else { - dipAzimuth = 0.0; // Vertical cell - } - - // Clamp dip angle to reasonable values - const Scalar maxDip = std::numbers::pi_v/2 - static_cast(1e-6); - dipAngle = std::min(dipAngle, maxDip); - } else { - // Degenerate cell - assume horizontal - dipAngle = 0.0; - dipAzimuth = 0.0; - } -} - -template -CellCornerData getCellCornerXY(const Element& element) -{ - typedef typename Element::Geometry Geometry; - const Geometry& geometry = element.geometry(); - static constexpr int zCoord = Element::dimension - 1; - static constexpr int yCoord = Element::dimension - 2; - static constexpr int xCoord = Element::dimension - 3; - const int corners = geometry.corners(); - assert(corners == 8); - std::array X {}; - std::array Y {}; - std::array Z {}; - // Get all 8 corners of the hexahedral cell (maybe expensive) - for (int i = 0; i < corners; ++i) { - auto corner = geometry.corner(i); - X[i] = corner[xCoord]; - Y[i] = corner[yCoord]; - Z[i] = corner[zCoord]; - } - - return CellCornerData{X, Y, Z}; -} - -template -Scalar calculateTrueVerticalDepth(Scalar z, Scalar x, Scalar y, - Scalar dipAngle, Scalar dipAzimuth, - const std::array& referencePoint) -{ - // For True Vertical Depth calculation: - // TVD = reference_depth + (z - reference_z) * cos(dipAngle) - // + lateral_distance * sin(dipAngle) * cos(azimuth_difference) - - // Calculate lateral displacement from reference point - Scalar dx = x - referencePoint[0]; - Scalar dy = y - referencePoint[1]; - Scalar dz = z - referencePoint[2]; - - // If no dip, TVD is simply the depth - if (std::abs(dipAngle) < 1e-10) { - return referencePoint[2] + dz; - } - - // Calculate the direction from reference point to current point - Scalar pointAzimuth = std::atan2(dy, dx); - - // Calculate the angle between dip direction and point direction - Scalar azimuthDiff = pointAzimuth - dipAzimuth; - - // Calculate lateral distance - Scalar lateralDist = std::hypot(dx, dy); - - // Project lateral distance onto dip direction - Scalar lateralInDipDir = lateralDist * std::cos(azimuthDiff); - - // True Vertical Depth calculation - // TVD increases with depth (more negative z means deeper) - // For a dipping plane: TVD = vertical_component + dip_component - Scalar tvd = referencePoint[2] + dz * std::cos(dipAngle) + lateralInDipDir * std::sin(dipAngle); - - return tvd; -} - -template -RK4IVP::RK4IVP(const RHS& f, - const std::array& span, - const Scalar y0, - const int N) - : N_(N) - , span_(span) -{ - const Scalar h = stepsize(); - const Scalar h2 = h / 2; - const Scalar h6 = h / 6; - - y_.reserve(N + 1); - f_.reserve(N + 1); - - y_.push_back(y0); - f_.push_back(f(span_[0], y0)); - - for (int i = 0; i < N; ++i) { - const Scalar x = span_[0] + i*h; - const Scalar y = y_.back(); - - 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() == typename std::vector::size_type(N + 1)); -} - -template -Scalar RK4IVP:: -operator()(const Scalar x) const -{ - // Dense output (O(h**3)) according to Shampine - // (Hermite interpolation) - const Scalar h = stepsize(); - int i = (x - span_[0]) / 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 Scalar y0 = y_[i], y1 = y_[i + 1]; - const Scalar f0 = f_[i], f1 = f_[i + 1]; - - Scalar u = (1 - 2*t) * (y1 - y0); - u += h * ((t - 1)*f0 + t*f1); - u *= t * (t - 1); - u += (1 - t)*y0 + t*y1; - - return u; -} - -template -Scalar RK4IVP:: -stepsize() const -{ - return (span_[1] - span_[0]) / N_; -} - -namespace PhasePressODE { - -template -Water:: -Water(const TabulatedFunction& tempVdTable, - const TabulatedFunction& saltVdTable, - const int pvtRegionIdx, - const Scalar normGrav) - : tempVdTable_(tempVdTable) - , saltVdTable_(saltVdTable) - , pvtRegionIdx_(pvtRegionIdx) - , g_(normGrav) -{ -} - -template -typename Water::Scalar -Water:: -operator()(const Scalar depth, - const Scalar press) const -{ - return this->density(depth, press) * g_; -} - -template -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 - 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; -} - -template -Oil:: -Oil(const TabulatedFunction& tempVdTable, - const RS& rs, - const int pvtRegionIdx, - const Scalar normGrav) - : tempVdTable_(tempVdTable) - , rs_(rs) - , pvtRegionIdx_(pvtRegionIdx) - , g_(normGrav) -{ -} - -template -typename Oil::Scalar -Oil:: -operator()(const Scalar depth, - const Scalar press) const -{ - return this->density(depth, press) * g_; -} - -template -typename Oil::Scalar -Oil:: -density(const Scalar depth, - const Scalar press) const -{ - const Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true); - Scalar rs = 0.0; - if (FluidSystem::enableDissolvedGas() || FluidSystem::enableConstantRs()) - rs = rs_(depth, press, temp); - - 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); - } - Scalar rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); - if (FluidSystem::enableDissolvedGas() || FluidSystem::enableConstantRs()) { - rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - } - - return rho; -} - -template -Gas:: -Gas(const TabulatedFunction& tempVdTable, - const RV& rv, - const RVW& rvw, - const int pvtRegionIdx, - const Scalar normGrav) - : tempVdTable_(tempVdTable) - , rv_(rv) - , rvw_(rvw) - , pvtRegionIdx_(pvtRegionIdx) - , g_(normGrav) -{ -} - -template -typename Gas::Scalar -Gas:: -operator()(const Scalar depth, - const Scalar press) const -{ - return this->density(depth, press) * g_; -} - -template -typename Gas::Scalar -Gas:: -density(const Scalar depth, - const Scalar press) const -{ - const Scalar temp = tempVdTable_.eval(depth, /*extrapolate=*/true); - Scalar rv = 0.0; - if (FluidSystem::enableVaporizedOil()) - rv = rv_(depth, press, temp); - - Scalar rvw = 0.0; - if (FluidSystem::enableVaporizedWater()) - rvw = rvw_(depth, press, temp); - - Scalar bGas = 0.0; - - if (FluidSystem::enableVaporizedOil() && FluidSystem::enableVaporizedWater()) { - if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press) - && rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press)) - { - bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); - } else { - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp, press, rv, rvw); - } - 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; - } - - if (FluidSystem::enableVaporizedOil()){ - if (rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press)) { - bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); - } else { - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, - temp, - press, - rv, - Scalar{0.0}/*=rvw*/); - } - Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); - return rho; - } - - if (FluidSystem::enableVaporizedWater()){ - if (rvw >= FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press)) { - bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp, press); - } - else { - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, - temp, - press, - Scalar{0.0} /*=rv*/, - rvw); - } - 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, - Scalar{0.0} /*=rv*/, - Scalar{0.0} /*=rvw*/); - Scalar rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - - return rho; -} - -} - -template -template -PressureTable:: -PressureFunction::PressureFunction(const ODE& ode, - const InitCond& ic, - const int nsample, - const VSpan& span) - : initial_(ic) -{ - this->value_[Direction::Up] = std::make_unique - (ode, VSpan {{ ic.depth, span[0] }}, ic.pressure, nsample); - - this->value_[Direction::Down] = std::make_unique - (ode, VSpan {{ ic.depth, span[1] }}, ic.pressure, nsample); -} - -template -template -PressureTable:: -PressureFunction::PressureFunction(const PressureFunction& rhs) - : initial_(rhs.initial_) -{ - this->value_[Direction::Up] = - std::make_unique(*rhs.value_[Direction::Up]); - - this->value_[Direction::Down] = - std::make_unique(*rhs.value_[Direction::Down]); -} - -template -template -typename PressureTable::template PressureFunction& -PressureTable:: -PressureFunction:: -operator=(const PressureFunction& rhs) -{ - this->initial_ = rhs.initial_; - - this->value_[Direction::Up] = - std::make_unique(*rhs.value_[Direction::Up]); - - this->value_[Direction::Down] = - std::make_unique(*rhs.value_[Direction::Down]); - - return *this; -} - -template -template -typename PressureTable::template PressureFunction& -PressureTable:: -PressureFunction:: -operator=(PressureFunction&& rhs) -{ - this->initial_ = rhs.initial_; - this->value_ = std::move(rhs.value_); - - return *this; -} - -template -template -typename PressureTable::Scalar -PressureTable:: -PressureFunction:: -value(const Scalar depth) const -{ - if (depth < this->initial_.depth) { - // Value above initial condition depth. - return (*this->value_[Direction::Up])(depth); - } - else if (depth > this->initial_.depth) { - // Value below initial condition depth. - return (*this->value_[Direction::Down])(depth); - } - else { - // Value *at* initial condition depth. - return this->initial_.pressure; - } -} - - -template -template -void PressureTable:: -checkPtr(const PressFunc* phasePress, - const std::string& phaseName) const -{ - if (phasePress != nullptr) { return; } - - throw std::invalid_argument { - "Phase pressure function for \"" + phaseName - + "\" most not be null" - }; -} - -template -typename PressureTable::Strategy -PressureTable:: -selectEquilibrationStrategy(const Region& reg) const -{ - if (!this->oilActive()) { - if (reg.datum() > reg.zwoc()) { // Datum in water zone - return &PressureTable::equil_WOG; - } - return &PressureTable::equil_GOW; - } - - if (reg.datum() > reg.zwoc()) { // Datum in water zone - return &PressureTable::equil_WOG; - } - else if (reg.datum() < reg.zgoc()) { // Datum in gas zone - return &PressureTable::equil_GOW; - } - else { // Datum in oil zone - return &PressureTable::equil_OWG; - } -} - -template -void PressureTable:: -copyInPointers(const PressureTable& rhs) -{ - if (rhs.oil_ != nullptr) { - this->oil_ = std::make_unique(*rhs.oil_); - } - - if (rhs.gas_ != nullptr) { - this->gas_ = std::make_unique(*rhs.gas_); - } - - if (rhs.wat_ != nullptr) { - this->wat_ = std::make_unique(*rhs.wat_); - } -} - -template -PhaseSaturations:: -PhaseSaturations(MaterialLawManager& matLawMgr, - const std::vector& swatInit) - : matLawMgr_(matLawMgr) - , swatInit_ (swatInit) -{ -} - -template -PhaseSaturations:: -PhaseSaturations(const PhaseSaturations& rhs) - : matLawMgr_(rhs.matLawMgr_) - , swatInit_ (rhs.swatInit_) - , sat_ (rhs.sat_) - , press_ (rhs.press_) -{ - // Note: We don't need to do anything to the 'fluidState_' here. - this->setEvaluationPoint(*rhs.evalPt_.position, - *rhs.evalPt_.region, - *rhs.evalPt_.ptable); -} - -template -const PhaseQuantityValue& -PhaseSaturations:: -deriveSaturations(const Position& x, - const Region& reg, - const PTable& ptable) -{ - this->setEvaluationPoint(x, reg, ptable); - this->initializePhaseQuantities(); - - if (ptable.gasActive()) { this->deriveGasSat(); } - - if (ptable.waterActive()) { this->deriveWaterSat(); } - - - if (this->isOverlappingTransition()) { - this->fixUnphysicalTransition(); - } - - if (ptable.oilActive()) { this->deriveOilSat(); } - - this->accountForScaledSaturations(); - - return this->sat_; -} - -template -void PhaseSaturations:: -setEvaluationPoint(const Position& x, - const Region& reg, - const PTable& ptable) -{ - this->evalPt_.position = &x; - this->evalPt_.region = ® - this->evalPt_.ptable = &ptable; -} - -template -void PhaseSaturations:: -initializePhaseQuantities() -{ - this->sat_.reset(); - this->press_.reset(); - - const auto depth = this->evalPt_.position->depth; - const auto& ptable = *this->evalPt_.ptable; - - if (ptable.oilActive()) { - this->press_.oil = ptable.oil(depth); - } - - if (ptable.gasActive()) { - this->press_.gas = ptable.gas(depth); - } - - if (ptable.waterActive()) { - this->press_.water = ptable.water(depth); - } -} - -template -void PhaseSaturations::deriveOilSat() -{ - this->sat_.oil = 1.0 - this->sat_.water - this->sat_.gas; -} - -template -void PhaseSaturations::deriveGasSat() -{ - auto& sg = this->sat_.gas; - - const auto isIncr = true; // dPcgo/dSg >= 0 for all Sg. - const auto oilActive = this->evalPt_.ptable->oilActive(); - - if (this->isConstCapPress(this->gasPos())) { - // Sharp interface between phases. Can derive phase saturation - // directly from knowing where 'depth' of evaluation point is - // relative to depth of O/G contact. - const auto gas_contact = oilActive? this->evalPt_.region->zgoc() : this->evalPt_.region->zwoc(); - sg = this->fromDepthTable(gas_contact, - this->gasPos(), isIncr); - } - else { - // Capillary pressure curve is non-constant, meaning there is a - // transition zone between the gas and oil phases. Invert capillary - // pressure relation - // - // Pcgo(Sg) = Pg - Po - // - // Note that Pcgo is defined to be (Pg - Po), not (Po - Pg). - const auto pw = oilActive? this->press_.oil : this->press_.water; - const auto pcgo = this->press_.gas - pw; - sg = this->invertCapPress(pcgo, this->gasPos(), isIncr); - } -} - -template -void PhaseSaturations::deriveWaterSat() -{ - auto& sw = this->sat_.water; - - const auto oilActive = this->evalPt_.ptable->oilActive(); - if (!oilActive) { - // for 2p gas+water we set the water saturation to 1.0 - sg - sw = 1.0 - this->sat_.gas; - } - else { - const auto isIncr = false; // dPcow/dSw <= 0 for all Sw. - - if (this->isConstCapPress(this->waterPos())) { - // Sharp interface between phases. Can derive phase saturation - // directly from knowing where 'depth' of evaluation point is - // relative to depth of O/W contact. - sw = this->fromDepthTable(this->evalPt_.region->zwoc(), - this->waterPos(), isIncr); - } - else { - // Capillary pressure curve is non-constant, meaning there is a - // transition zone between the oil and water phases. Invert - // capillary pressure relation - // - // Pcow(Sw) = Po - Pw - // - // unless the model uses "SWATINIT". In the latter case, pick the - // saturation directly from the SWATINIT array of the pertinent - // cell. - const auto pcow = this->press_.oil - this->press_.water; - - if (this->swatInit_.empty()) { - sw = this->invertCapPress(pcow, this->waterPos(), isIncr); - } - else { - auto [swout, newSwatInit] = this->applySwatInit(pcow); - if (newSwatInit) - sw = this->invertCapPress(pcow, this->waterPos(), isIncr); - else { - sw = swout; - } - } - } - } -} - -template -void PhaseSaturations:: -fixUnphysicalTransition() -{ - auto& sg = this->sat_.gas; - auto& sw = this->sat_.water; - - // Overlapping gas/oil and oil/water transition zones can lead to - // unphysical phase saturations when individual saturations are derived - // directly from inverting O/G and O/W capillary pressure curves. - // - // Recalculate phase saturations using the implied gas/water capillary - // pressure: Pg - Pw. - const auto pcgw = this->press_.gas - this->press_.water; - if (! this->swatInit_.empty()) { - // Re-scale Pc to reflect imposed sw for vanishing oil phase. This - // seems consistent with ECLIPSE, but fails to honour SWATINIT in - // case of non-trivial gas/oil capillary pressure. - auto [swout, newSwatInit] = this->applySwatInit(pcgw, sw); - if (newSwatInit){ - const auto isIncr = false; // dPcow/dSw <= 0 for all Sw. - sw = this->invertCapPress(pcgw, this->waterPos(), isIncr); - } - else { - sw = swout; - } - } - - sw = satFromSumOfPcs - (this->matLawMgr_, this->waterPos(), this->gasPos(), - this->evalPt_.position->cell, pcgw); - sg = 1.0 - sw; - - this->fluidState_.setSaturation(this->oilPos(), 1.0 - sw - sg); - this->fluidState_.setSaturation(this->gasPos(), sg); - this->fluidState_.setSaturation(this->waterPos(), this->evalPt_ - .ptable->waterActive() ? sw : 0.0); - - // Pcgo = Pg - Po => Po = Pg - Pcgo - this->computeMaterialLawCapPress(); - this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil(); -} - -template -void PhaseSaturations:: -accountForScaledSaturations() -{ - const auto gasActive = this->evalPt_.ptable->gasActive(); - const auto watActive = this->evalPt_.ptable->waterActive(); - const auto oilActive = this->evalPt_.ptable->oilActive(); - - auto sg = gasActive? this->sat_.gas : 0.0; - auto sw = watActive? this->sat_.water : 0.0; - auto so = oilActive? this->sat_.oil : 0.0; - - this->fluidState_.setSaturation(this->waterPos(), sw); - this->fluidState_.setSaturation(this->oilPos(), so); - this->fluidState_.setSaturation(this->gasPos(), sg); - - const auto& scaledDrainageInfo = this->matLawMgr_ - .oilWaterScaledEpsInfoDrainage(this->evalPt_.position->cell); - - const auto thresholdSat = 1.0e-6; - if (watActive && ((sw + thresholdSat) > scaledDrainageInfo.Swu)) { - // Water saturation exceeds maximum possible value. Reset oil phase - // pressure to that which corresponds to maximum possible water - // saturation value. - this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swu); - if (oilActive) { - this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swu); - } else if (gasActive) { - this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swu); - } - sw = scaledDrainageInfo.Swu; - this->computeMaterialLawCapPress(); - - if (oilActive) { - // Pcow = Po - Pw => Po = Pw + Pcow - this->press_.oil = this->press_.water + this->materialLawCapPressOilWater(); - } else { - // Pcgw = Pg - Pw => Pg = Pw + Pcgw - this->press_.gas = this->press_.water + this->materialLawCapPressGasWater(); - } - - } - if (gasActive && ((sg + thresholdSat) > scaledDrainageInfo.Sgu)) { - // Gas saturation exceeds maximum possible value. Reset oil phase - // pressure to that which corresponds to maximum possible gas - // saturation value. - this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgu); - if (oilActive) { - this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgu); - } else if (watActive) { - this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgu); - } - sg = scaledDrainageInfo.Sgu; - this->computeMaterialLawCapPress(); - - if (oilActive) { - // Pcgo = Pg - Po => Po = Pg - Pcgo - this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil(); - } else { - // Pcgw = Pg - Pw => Pw = Pg - Pcgw - this->press_.water = this->press_.gas - this->materialLawCapPressGasWater(); - } - } - - if (watActive && ((sw - thresholdSat) < scaledDrainageInfo.Swl)) { - // Water saturation less than minimum possible value in cell. Reset - // water phase pressure to that which corresponds to minimum - // possible water saturation value. - this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swl); - if (oilActive) { - this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swl); - } else if (gasActive) { - this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swl); - } - sw = scaledDrainageInfo.Swl; - this->computeMaterialLawCapPress(); - - if (oilActive) { - // Pcwo = Po - Pw => Pw = Po - Pcow - this->press_.water = this->press_.oil - this->materialLawCapPressOilWater(); - } else { - // Pcgw = Pg - Pw => Pw = Pg - Pcgw - this->press_.water = this->press_.gas - this->materialLawCapPressGasWater(); - } - } - - if (gasActive && ((sg - thresholdSat) < scaledDrainageInfo.Sgl)) { - // Gas saturation less than minimum possible value in cell. Reset - // gas phase pressure to that which corresponds to minimum possible - // gas saturation. - this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgl); - if (oilActive) { - this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgl); - } else if (watActive) { - this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgl); - } - sg = scaledDrainageInfo.Sgl; - this->computeMaterialLawCapPress(); - - if (oilActive) { - // Pcgo = Pg - Po => Pg = Po + Pcgo - this->press_.gas = this->press_.oil + this->materialLawCapPressGasOil(); - } else { - // Pcgw = Pg - Pw => Pg = Pw + Pcgw - this->press_.gas = this->press_.water + this->materialLawCapPressGasWater(); - } - } -} - -template -std::pair -PhaseSaturations:: -applySwatInit(const Scalar pcow) -{ - return this->applySwatInit(pcow, this->swatInit_[this->evalPt_.position->cell]); -} - -template -std::pair -PhaseSaturations:: -applySwatInit(const Scalar pcow, const Scalar sw) -{ - return this->matLawMgr_.applySwatinit(this->evalPt_.position->cell, pcow, sw); -} - -template -void PhaseSaturations:: -computeMaterialLawCapPress() -{ - const auto& matParams = this->matLawMgr_ - .materialLawParams(this->evalPt_.position->cell); - - this->matLawCapPress_.fill(0.0); - MaterialLaw::capillaryPressures(this->matLawCapPress_, - matParams, this->fluidState_); -} - -template -typename FluidSystem::Scalar -PhaseSaturations:: -materialLawCapPressGasOil() const -{ - return this->matLawCapPress_[this->oilPos()] - + this->matLawCapPress_[this->gasPos()]; -} - -template -typename FluidSystem::Scalar -PhaseSaturations:: -materialLawCapPressOilWater() const -{ - return this->matLawCapPress_[this->oilPos()] - - this->matLawCapPress_[this->waterPos()]; -} - -template -typename FluidSystem::Scalar -PhaseSaturations:: -materialLawCapPressGasWater() const -{ - return this->matLawCapPress_[this->gasPos()] - - this->matLawCapPress_[this->waterPos()]; -} - -template -bool PhaseSaturations:: -isConstCapPress(const PhaseIdx phaseIdx) const -{ - return isConstPc - (this->matLawMgr_, phaseIdx, this->evalPt_.position->cell); -} - -template -bool PhaseSaturations:: -isOverlappingTransition() const -{ - return this->evalPt_.ptable->gasActive() - && this->evalPt_.ptable->waterActive() - && ((this->sat_.gas + this->sat_.water) > 1.0); -} - -template -typename FluidSystem::Scalar -PhaseSaturations:: -fromDepthTable(const Scalar contactdepth, - const PhaseIdx phasePos, - const bool isincr) const -{ - return satFromDepth - (this->matLawMgr_, this->evalPt_.position->depth, - contactdepth, static_cast(phasePos), - this->evalPt_.position->cell, isincr); -} - -template -typename FluidSystem::Scalar -PhaseSaturations:: -invertCapPress(const Scalar pc, - const PhaseIdx phasePos, - const bool isincr) const -{ - return satFromPc - (this->matLawMgr_, static_cast(phasePos), - this->evalPt_.position->cell, pc, isincr); -} - -template -PressureTable:: -PressureTable(const Scalar gravity, - const int samplePoints) - : gravity_(gravity) - , nsample_(samplePoints) -{ -} - -template -PressureTable:: -PressureTable(const PressureTable& rhs) - : gravity_(rhs.gravity_) - , nsample_(rhs.nsample_) -{ - this->copyInPointers(rhs); -} - -template -PressureTable:: -PressureTable(PressureTable&& rhs) - : gravity_(rhs.gravity_) - , nsample_(rhs.nsample_) - , oil_ (std::move(rhs.oil_)) - , gas_ (std::move(rhs.gas_)) - , wat_ (std::move(rhs.wat_)) -{ -} - -template -PressureTable& -PressureTable:: -operator=(const PressureTable& rhs) -{ - this->gravity_ = rhs.gravity_; - this->nsample_ = rhs.nsample_; - this->copyInPointers(rhs); - - return *this; -} - -template -PressureTable& -PressureTable:: -operator=(PressureTable&& rhs) -{ - this->gravity_ = rhs.gravity_; - this->nsample_ = rhs.nsample_; - - this->oil_ = std::move(rhs.oil_); - this->gas_ = std::move(rhs.gas_); - this->wat_ = std::move(rhs.wat_); - - return *this; -} - -template -void PressureTable:: -equilibrate(const Region& reg, - const VSpan& span) -{ - // One of the PressureTable::equil_*() member functions. - auto equil = this->selectEquilibrationStrategy(reg); - - (this->*equil)(reg, span); -} - -template -bool PressureTable:: -oilActive() const -{ - return FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); -} - -template -bool PressureTable:: -gasActive() const -{ - return FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); -} - -template -bool PressureTable:: -waterActive() const -{ - return FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); -} - -template -typename FluidSystem::Scalar -PressureTable:: -oil(const Scalar depth) const -{ - this->checkPtr(this->oil_.get(), "OIL"); - - return this->oil_->value(depth); -} - -template -typename FluidSystem::Scalar -PressureTable:: -gas(const Scalar depth) const -{ - this->checkPtr(this->gas_.get(), "GAS"); - - return this->gas_->value(depth); -} - - -template -typename FluidSystem::Scalar -PressureTable:: -water(const Scalar depth) const -{ - this->checkPtr(this->wat_.get(), "WATER"); - - return this->wat_->value(depth); -} - -template -void PressureTable:: -equil_WOG(const Region& reg, const VSpan& span) -{ - // Datum depth in water zone. Calculate phase pressure for water first, - // followed by oil and gas if applicable. - - if (! this->waterActive()) { - throw std::invalid_argument { - "Don't know how to interpret EQUIL datum depth in " - "WATER zone in model without active water phase" - }; - } - - { - const auto ic = typename WPress::InitCond { - reg.datum(), reg.pressure() - }; - - this->makeWatPressure(ic, reg, span); - } - - if (this->oilActive()) { - // Pcow = Po - Pw => Po = Pw + Pcow - const auto ic = typename OPress::InitCond { - reg.zwoc(), - this->water(reg.zwoc()) + reg.pcowWoc() - }; - - this->makeOilPressure(ic, reg, span); - } - - if (this->gasActive() && this->oilActive()) { - // Pcgo = Pg - Po => Pg = Po + Pcgo - const auto ic = typename GPress::InitCond { - reg.zgoc(), - this->oil(reg.zgoc()) + reg.pcgoGoc() - }; - - this->makeGasPressure(ic, reg, span); - } else if (this->gasActive() && !this->oilActive()) { - // No oil phase set Pg = Pw + Pcgw - const auto ic = typename GPress::InitCond { - reg.zwoc(), // The WOC is really the GWC for gas/water cases - this->water(reg.zwoc()) + reg.pcowWoc() // Pcow(WOC) is really Pcgw(GWC) for gas/water cases - }; - this->makeGasPressure(ic, reg, span); - } -} - -template -void PressureTable:: -equil_GOW(const Region& reg, const VSpan& span) -{ - // Datum depth in gas zone. Calculate phase pressure for gas first, - // followed by oil and water if applicable. - - if (! this->gasActive()) { - throw std::invalid_argument { - "Don't know how to interpret EQUIL datum depth in " - "GAS zone in model without active gas phase" - }; - } - - { - const auto ic = typename GPress::InitCond { - reg.datum(), reg.pressure() - }; - - this->makeGasPressure(ic, reg, span); - } - - if (this->oilActive()) { - // Pcgo = Pg - Po => Po = Pg - Pcgo - const auto ic = typename OPress::InitCond { - reg.zgoc(), - this->gas(reg.zgoc()) - reg.pcgoGoc() - }; - this->makeOilPressure(ic, reg, span); - } - - if (this->waterActive() && this->oilActive()) { - // Pcow = Po - Pw => Pw = Po - Pcow - const auto ic = typename WPress::InitCond { - reg.zwoc(), - this->oil(reg.zwoc()) - reg.pcowWoc() - }; - - this->makeWatPressure(ic, reg, span); - } else if (this->waterActive() && !this->oilActive()) { - // No oil phase set Pw = Pg - Pcgw - const auto ic = typename WPress::InitCond { - reg.zwoc(), // The WOC is really the GWC for gas/water cases - this->gas(reg.zwoc()) - reg.pcowWoc() // Pcow(WOC) is really Pcgw(GWC) for gas/water cases - }; - this->makeWatPressure(ic, reg, span); - } -} - -template -void PressureTable:: -equil_OWG(const Region& reg, const VSpan& span) -{ - // Datum depth in oil zone. Calculate phase pressure for oil first, - // followed by gas and water if applicable. - - if (! this->oilActive()) { - throw std::invalid_argument { - "Don't know how to interpret EQUIL datum depth in " - "OIL zone in model without active oil phase" - }; - } - - { - const auto ic = typename OPress::InitCond { - reg.datum(), reg.pressure() - }; - - this->makeOilPressure(ic, reg, span); - } - - if (this->waterActive()) { - // Pcow = Po - Pw => Pw = Po - Pcow - const auto ic = typename WPress::InitCond { - reg.zwoc(), - this->oil(reg.zwoc()) - reg.pcowWoc() - }; - - this->makeWatPressure(ic, reg, span); - } - - if (this->gasActive()) { - // Pcgo = Pg - Po => Pg = Po + Pcgo - const auto ic = typename GPress::InitCond { - reg.zgoc(), - this->oil(reg.zgoc()) + reg.pcgoGoc() - }; - this->makeGasPressure(ic, reg, span); - } -} - -template -void PressureTable:: -makeOilPressure(const typename OPress::InitCond& ic, - const Region& reg, - const VSpan& span) -{ - const auto drho = OilPressODE { - reg.tempVdTable(), reg.dissolutionCalculator(), - reg.pvtIdx(), this->gravity_ - }; - - this->oil_ = std::make_unique(drho, ic, this->nsample_, span); -} - -template -void PressureTable:: -makeGasPressure(const typename GPress::InitCond& ic, - const Region& reg, - const VSpan& span) -{ - const auto drho = GasPressODE { - reg.tempVdTable(), reg.evaporationCalculator(), reg.waterEvaporationCalculator(), - reg.pvtIdx(), this->gravity_ - }; - - this->gas_ = std::make_unique(drho, ic, this->nsample_, span); -} - -template -void PressureTable:: -makeWatPressure(const typename WPress::InitCond& ic, - const Region& reg, - const VSpan& span) -{ - const auto drho = WatPressODE { - reg.tempVdTable(), reg.saltVdTable(), reg.pvtIdx(), this->gravity_ - }; - - this->wat_ = std::make_unique(drho, ic, this->nsample_, span); -} - -} - -======= ->>>>>>> 02e7fa473 (refactoring InitStateEquil) namespace DeckDependent { template