From c2e9a7a518dd2bed42ff75af6ab89f0bf1e084af Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 2 Jan 2018 12:43:56 +0100 Subject: [PATCH] move the hydrostatic equilibrium code to its proper location and make it compile this just moves the hydrostatic equilibrium code from its historc location at opm/core to ebos/equil and adds minimal changes to make it compile. this allows to clean up that code without disturbing the legacy simulators. --- opm/core/simulator/EquilibrationHelpers.hpp | 799 -------------------- opm/core/simulator/initStateEquil.hpp | 446 ----------- opm/core/simulator/initStateEquil_impl.hpp | 791 ------------------- opm/core/utility/RegionMapping.hpp | 181 ----- 4 files changed, 2217 deletions(-) delete mode 100644 opm/core/simulator/EquilibrationHelpers.hpp delete mode 100644 opm/core/simulator/initStateEquil.hpp delete mode 100644 opm/core/simulator/initStateEquil_impl.hpp delete mode 100644 opm/core/utility/RegionMapping.hpp diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp deleted file mode 100644 index 05e841db2..000000000 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ /dev/null @@ -1,799 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - Copyright 2017 IRIS - - 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 . -*/ - -#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED -#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED - -#include -#include -#include - -#include - -#include -#include -#include - -#include - - -/* ----- synopsis of EquilibrationHelpers.hpp ---- - -namespace Opm -{ - namespace EQUIL { - - namespace Miscibility { - class RsFunction; - class NoMixing; - template - class RsVD; - template - class RsSatAtContact; - } - - class EquilReg; - - - template - struct PcEq; - - template - inline double satFromPc(const MaterialLawManager& materialLawManager, - const int phase, - const int cell, - const double target_pc, - const bool increasing = false) - template - inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager, - const int phase1, - const int phase2, - const int cell, - const double target_pc) - } // namespace Equil -} // namespace Opm - ----- end of synopsis of EquilibrationHelpers.hpp ---- -*/ - - -namespace Opm -{ - /** - * 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 { - - - typedef Opm::FluidSystems::BlackOil FluidSystemSimple; - - // Adjust oil pressure according to gas saturation and cap pressure - typedef Opm::SimpleModularFluidState SatOnlyFluidState; - - /** - * Types and routines relating to phase mixing in - * equilibration calculations. - */ - namespace Miscibility { - - /** - * Base class for phase mixing functions. - */ - class RsFunction - { - public: - /** - * Function call operator. - * - * \param[in] depth Depth at which to calculate RS - * value. - * - * \param[in] press Pressure at which to calculate RS - * value. - * - * \param[in] temp Temperature at which to calculate RS - * value. - * - * \return Dissolved gas-oil ratio (RS) at depth @c - * depth and pressure @c press. - */ - virtual double operator()(const double depth, - const double press, - const double temp, - const double sat = 0.0) const = 0; - }; - - - /** - * Type that implements "no phase mixing" policy. - */ - class NoMixing : public RsFunction { - public: - /** - * Function call. - * - * \param[in] depth Depth at which to calculate RS - * value. - * - * \param[in] press Pressure at which to calculate RS - * value. - * - * \param[in] temp Temperature at which to calculate RS - * value. - * - * \return Dissolved gas-oil ratio (RS) at depth @c - * depth and pressure @c press. In "no mixing - * policy", this is identically zero. - */ - double - operator()(const double /* depth */, - const double /* press */, - const double /* temp */, - const double /* sat */ = 0.0) const - { - return 0.0; - } - }; - - - /** - * Type that implements "dissolved gas-oil ratio" - * tabulated as a function of depth policy. Data - * typically taken from keyword 'RSVD'. - */ - template - class RsVD : public RsFunction { - public: - /** - * Constructor. - * - * \param[in] pvtRegionIdx The pvt region index - * \param[in] depth Depth nodes. - * \param[in] rs Dissolved gas-oil ratio at @c depth. - */ - RsVD(const int pvtRegionIdx, - const std::vector& depth, - const std::vector& rs) - : pvtRegionIdx_(pvtRegionIdx) - , depth_(depth) - , rs_(rs) - { - } - - /** - * Function call. - * - * \param[in] depth Depth at which to calculate RS - * value. - * - * \param[in] press Pressure at which to calculate RS - * value. - * - * \param[in] temp Temperature at which to calculate RS - * value. - * - * \return Dissolved gas-oil ratio (RS) at depth @c - * depth and pressure @c press. - */ - double - operator()(const double depth, - const double press, - const double temp, - const double sat_gas = 0.0) const - { - if (sat_gas > 0.0) { - return satRs(press, temp); - } else { - return std::min(satRs(press, temp), linearInterpolationNoExtrapolation(depth_, rs_, depth)); - } - } - - private: - const int pvtRegionIdx_; - std::vector depth_; /**< Depth nodes */ - std::vector rs_; /**< Dissolved gas-oil ratio */ - - double satRs(const double press, const double temp) const - { - return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); - } - }; - - - /** - * Type that implements "vaporized oil-gas ratio" - * tabulated as a function of depth policy. Data - * typically taken from keyword 'RVVD'. - */ - template - class RvVD : public RsFunction { - public: - /** - * Constructor. - * - * \param[in] pvtRegionIdx The pvt region index - * \param[in] depth Depth nodes. - * \param[in] rv Dissolved gas-oil ratio at @c depth. - */ - RvVD(const int pvtRegionIdx, - const std::vector& depth, - const std::vector& rv) - : pvtRegionIdx_(pvtRegionIdx) - , depth_(depth) - , rv_(rv) - { - } - - /** - * Function call. - * - * \param[in] depth Depth at which to calculate RV - * value. - * - * \param[in] press Pressure at which to calculate RV - * value. - * - * \param[in] temp Temperature at which to calculate RV - * value. - * - * \return Vaporized oil-gas ratio (RV) at depth @c - * depth and pressure @c press. - */ - double - operator()(const double depth, - const double press, - const double temp, - const double sat_oil = 0.0 ) const - { - if (std::abs(sat_oil) > 1e-16) { - return satRv(press, temp); - } else { - return std::min(satRv(press, temp), linearInterpolationNoExtrapolation(depth_, rv_, depth)); - } - } - - private: - const int pvtRegionIdx_; - std::vector depth_; /**< Depth nodes */ - std::vector rv_; /**< Vaporized oil-gas ratio */ - - double satRv(const double press, const double temp) const - { - return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press); - } - }; - - - /** - * Class that implements "dissolved gas-oil ratio" (Rs) - * as function of depth and pressure as follows: - * - * 1. The Rs at the gas-oil contact is equal to the - * saturated Rs value, Rs_sat_contact. - * - * 2. The Rs elsewhere is equal to Rs_sat_contact, but - * constrained to the saturated value as given by the - * local pressure. - * - * This should yield Rs-values that are constant below the - * contact, and decreasing above the contact. - */ - template - class RsSatAtContact : public RsFunction { - public: - /** - * Constructor. - * - * \param[in] pvtRegionIdx The pvt region index - * \param[in] p_contact oil pressure at the contact - * \param[in] T_contact temperature at the contact - */ - RsSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact) - : pvtRegionIdx_(pvtRegionIdx) - { - rs_sat_contact_ = satRs(p_contact, T_contact); - } - - /** - * Function call. - * - * \param[in] depth Depth at which to calculate RS - * value. - * - * \param[in] press Pressure at which to calculate RS - * value. - * - * \param[in] temp Temperature at which to calculate RS - * value. - * - * \return Dissolved gas-oil ratio (RS) at depth @c - * depth and pressure @c press. - */ - double - operator()(const double /* depth */, - const double press, - const double temp, - const double sat_gas = 0.0) const - { - if (sat_gas > 0.0) { - return satRs(press, temp); - } else { - return std::min(satRs(press, temp), rs_sat_contact_); - } - } - - private: - const int pvtRegionIdx_; - double rs_sat_contact_; - - double satRs(const double press, const double temp) const - { - return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press); - } - }; - - - /** - * Class that implements "vaporized oil-gas ratio" (Rv) - * as function of depth and pressure as follows: - * - * 1. The Rv at the gas-oil contact is equal to the - * saturated Rv value, Rv_sat_contact. - * - * 2. The Rv elsewhere is equal to Rv_sat_contact, but - * constrained to the saturated value as given by the - * local pressure. - * - * This should yield Rv-values that are constant below the - * contact, and decreasing above the contact. - */ - template - class RvSatAtContact : public RsFunction { - public: - /** - * Constructor. - * - * \param[in] pvtRegionIdx The pvt region index - * \param[in] p_contact oil pressure at the contact - * \param[in] T_contact temperature at the contact - */ - RvSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact) - :pvtRegionIdx_(pvtRegionIdx) - { - rv_sat_contact_ = satRv(p_contact, T_contact); - } - - /** - * Function call. - * - * \param[in] depth Depth at which to calculate RV - * value. - * - * \param[in] press Pressure at which to calculate RV - * value. - * - * \param[in] temp Temperature at which to calculate RV - * value. - * - * \return Dissolved oil-gas ratio (RV) at depth @c - * depth and pressure @c press. - */ - double - operator()(const double /*depth*/, - const double press, - const double temp, - const double sat_oil = 0.0) const - { - if (sat_oil > 0.0) { - return satRv(press, temp); - } else { - return std::min(satRv(press, temp), rv_sat_contact_); - } - } - - private: - const int pvtRegionIdx_; - double rv_sat_contact_; - - double satRv(const double press, const double temp) const - { - return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);; - } - }; - - } // namespace Miscibility - - /** - * Aggregate information base of an equilibration region. - * - * Provides inquiry methods for retrieving depths of contacs - * and pressure values as well as a means of calculating fluid - * densities, dissolved gas-oil ratio and vapourised oil-gas - * ratios. - * - * \tparam DensCalc Type that provides access to a phase - * density calculation facility. Must implement an operator() - * declared as - * - * std::vector - * operator()(const double press, - * const std::vector& svol ) - * - * that calculates the phase densities of all phases in @c - * svol at fluid pressure @c press. - */ - class EquilReg { - public: - /** - * Constructor. - * - * \param[in] rec Equilibration data of current region. - * \param[in] rs Calculator of dissolved gas-oil ratio. - * \param[in] rv Calculator of vapourised oil-gas ratio. - * \param[in] pvtRegionIdx The pvt region index - */ - EquilReg(const EquilRecord& rec, - std::shared_ptr rs, - std::shared_ptr rv, - const int pvtIdx) - : rec_ (rec) - , rs_ (rs) - , rv_ (rv) - , pvtIdx_ (pvtIdx) - { - } - - /** - * Type of dissolved gas-oil ratio calculator. - */ - typedef Miscibility::RsFunction CalcDissolution; - - /** - * Type of vapourised oil-gas ratio calculator. - */ - typedef Miscibility::RsFunction CalcEvaporation; - - /** - * Datum depth in current region - */ - double datum() const { return this->rec_.datumDepth(); } - - /** - * Pressure at datum depth in current region. - */ - double pressure() const { return this->rec_.datumDepthPressure(); } - - /** - * Depth of water-oil contact. - */ - double zwoc() const { return this->rec_.waterOilContactDepth(); } - - /** - * water-oil capillary pressure at water-oil contact. - * - * \return P_o - P_w at WOC. - */ - double pcow_woc() const { return this->rec_.waterOilContactCapillaryPressure(); } - - /** - * Depth of gas-oil contact. - */ - double zgoc() const { return this->rec_.gasOilContactDepth(); } - - /** - * Gas-oil capillary pressure at gas-oil contact. - * - * \return P_g - P_o at GOC. - */ - double pcgo_goc() const { return this->rec_.gasOilContactCapillaryPressure(); } - - - /** - * Retrieve dissolved gas-oil ratio calculator of current - * region. - */ - const CalcDissolution& - dissolutionCalculator() const { return *this->rs_; } - - /** - * Retrieve vapourised oil-gas ratio calculator of current - * region. - */ - const CalcEvaporation& - evaporationCalculator() const { return *this->rv_; } - - /** - * Retrieve pvtIdx of the region. - */ - int pvtIdx() const { return this->pvtIdx_; } - - - private: - EquilRecord rec_; /**< Equilibration data */ - std::shared_ptr rs_; /**< RS calculator */ - std::shared_ptr rv_; /**< RV calculator */ - const int pvtIdx_; - }; - - - - /// Functor for inverting capillary pressure function. - /// Function represented is - /// f(s) = pc(s) - target_pc - template - struct PcEq - { - PcEq(const MaterialLawManager& materialLawManager, - const int phase, - const int cell, - const double target_pc) - : materialLawManager_(materialLawManager), - phase_(phase), - cell_(cell), - target_pc_(target_pc) - { - - } - double operator()(double s) const - { - const auto& matParams = materialLawManager_.materialLawParams(cell_); - SatOnlyFluidState fluidState; - fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); - fluidState.setSaturation(phase_, s); - - double pc[FluidSystem::numPhases]; - std::fill(pc, pc + FluidSystem::numPhases, 0.0); - MaterialLaw::capillaryPressures(pc, matParams, fluidState); - double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; - double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_]; - return pcPhase - target_pc_; - } - private: - const MaterialLawManager& materialLawManager_; - const int phase_; - const int cell_; - const double target_pc_; - }; - - template - double minSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) { - const auto& scaledDrainageInfo = - materialLawManager.oilWaterScaledEpsInfoDrainage(cell); - - // Find minimum and maximum saturations. - switch(phase) { - case FluidSystem::waterPhaseIdx : - { - return scaledDrainageInfo.Swl; - break; - } - case FluidSystem::gasPhaseIdx : - { - return scaledDrainageInfo.Sgl; - break; - } - case FluidSystem::oilPhaseIdx : - { - OPM_THROW(std::runtime_error, "Min saturation not implemented for oil phase."); - break; - } - default: OPM_THROW(std::runtime_error, "Unknown phaseIdx ."); - } - return -1.0; - } - - template - double maxSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) { - const auto& scaledDrainageInfo = - materialLawManager.oilWaterScaledEpsInfoDrainage(cell); - - // Find minimum and maximum saturations. - switch(phase) { - case FluidSystem::waterPhaseIdx : - { - return scaledDrainageInfo.Swu; - break; - } - case FluidSystem::gasPhaseIdx : - { - return scaledDrainageInfo.Sgu; - break; - } - case FluidSystem::oilPhaseIdx : - { - OPM_THROW(std::runtime_error, "Max saturation not implemented for oil phase."); - break; - } - default: OPM_THROW(std::runtime_error, "Unknown phaseIdx ."); - } - return -1.0; - } - - - /// Compute saturation of some phase corresponding to a given - /// capillary pressure. - template - inline double satFromPc(const MaterialLawManager& materialLawManager, - const int phase, - const int cell, - const double target_pc, - const bool increasing = false) - { - // Find minimum and maximum saturations. - const double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); - const double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); - - // Create the equation f(s) = pc(s) - target_pc - const PcEq f(materialLawManager, phase, cell, target_pc); - const double f0 = f(s0); - const double f1 = f(s1); - - if (f0 <= 0.0) { - return s0; - } else if (f1 > 0.0) { - return s1; - } else { - const int max_iter = 60; - const double tol = 1e-6; - int iter_used = -1; - typedef RegulaFalsi ScalarSolver; - const double sol = ScalarSolver::solve(f, std::min(s0, s1), std::max(s0, s1), max_iter, tol, iter_used); - return sol; - } - } - - - /// Functor for inverting a sum of capillary pressure functions. - /// Function represented is - /// f(s) = pc1(s) + pc2(1 - s) - target_pc - template - struct PcEqSum - { - PcEqSum(const MaterialLawManager& materialLawManager, - const int phase1, - const int phase2, - const int cell, - const double target_pc) - : materialLawManager_(materialLawManager), - phase1_(phase1), - phase2_(phase2), - cell_(cell), - target_pc_(target_pc) - { - } - double operator()(double s) const - { - const auto& matParams = materialLawManager_.materialLawParams(cell_); - SatOnlyFluidState fluidState; - fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); - fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); - fluidState.setSaturation(phase1_, s); - fluidState.setSaturation(phase2_, 1.0 - s); - - double pc[FluidSystem::numPhases]; - std::fill(pc, pc + FluidSystem::numPhases, 0.0); - - MaterialLaw::capillaryPressures(pc, matParams, fluidState); - double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; - double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_]; - double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0; - double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_]; - return pc1 + pc2 - target_pc_; - } - private: - const MaterialLawManager& materialLawManager_; - const int phase1_; - const int phase2_; - const int cell_; - const double target_pc_; - }; - - - - - /// Compute saturation of some phase corresponding to a given - /// capillary pressure, where the capillary pressure function - /// is given as a sum of two other functions. - template - inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager, - const int phase1, - const int phase2, - const int cell, - const double target_pc) - { - // Find minimum and maximum saturations. - const double smin = minSaturations(materialLawManager, phase1, cell); - const double smax = maxSaturations(materialLawManager, phase1, cell); - - // Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc - const PcEqSum f(materialLawManager, phase1, phase2, cell, target_pc); - const double f0 = f(smin); - const double f1 = f(smax); - if (f0 <= 0.0) { - return smin; - } else if (f1 > 0.0) { - return smax; - } else { - const int max_iter = 30; - const double tol = 1e-6; - int iter_used = -1; - typedef RegulaFalsi ScalarSolver; - const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used); - return sol; - } - } - - /// Compute saturation from depth. Used for constant capillary pressure function - template - inline double satFromDepth(const MaterialLawManager& materialLawManager, - const double cellDepth, - const double contactDepth, - const int phase, - const int cell, - const bool increasing = false) - { - const double s0 = increasing ? maxSaturations(materialLawManager, phase, cell) : minSaturations(materialLawManager, phase, cell); - const double s1 = increasing ? minSaturations(materialLawManager, phase, cell) : maxSaturations(materialLawManager, phase, cell); - - if (cellDepth < contactDepth){ - return s0; - } else { - return s1; - } - - } - - /// Return true if capillary pressure function is constant - template - inline bool isConstPc(const MaterialLawManager& materialLawManager, - const int phase, - const int cell) - { - // Create the equation f(s) = pc(s); - const PcEq f(materialLawManager, phase, cell, 0); - const double f0 = f(minSaturations(materialLawManager, phase, cell)); - const double f1 = f(maxSaturations(materialLawManager, phase, cell)); - return std::abs(f0 - f1) < std::numeric_limits::epsilon(); - } - - } // namespace Equil -} // namespace Opm - - -#endif // OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp deleted file mode 100644 index 218a2ddb0..000000000 --- a/opm/core/simulator/initStateEquil.hpp +++ /dev/null @@ -1,446 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services - Copyright 2015 NTNU - Copyright 2017 IRIS - - 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 . -*/ - -#ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED -#define OPM_INITSTATEEQUIL_HEADER_INCLUDED - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include - -/** - * \file - * Facilities for an ECLIPSE-style equilibration-based - * initialisation scheme (keyword 'EQUIL'). - */ -struct UnstructuredGrid; - -namespace Opm -{ - - - /** - * 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 { - - /** - * Compute initial phase pressures by means of equilibration. - * - * This function uses the information contained in an - * equilibration record (i.e., depths and pressurs) as well as - * a density calculator and related data to vertically - * integrate the phase pressure ODE - * \f[ - * \frac{\mathrm{d}p_{\alpha}}{\mathrm{d}z} = - * \rho_{\alpha}(z,p_{\alpha})\cdot g - * \f] - * in which \f$\rho_{\alpha}$ denotes the fluid density of - * fluid phase \f$\alpha\f$, \f$p_{\alpha}\f$ is the - * corresponding phase pressure, \f$z\f$ is the depth and - * \f$g\f$ is the acceleration due to gravity (assumed - * directed downwords, in the positive \f$z\f$ direction). - * - * \tparam Region Type of an equilibration region information - * base. Typically an instance of the EquilReg - * class template. - * - * \tparam CellRange Type of cell range that demarcates the - * cells pertaining to the current - * equilibration region. Must implement - * methods begin() and end() to bound the range - * as well as provide an inner type, - * const_iterator, to traverse the range. - * - * \param[in] G Grid. - * \param[in] reg Current equilibration region. - * \param[in] cells Range that spans the cells of the current - * equilibration region. - * \param[in] grav Acceleration of gravity. - * - * \return Phase pressures, one vector for each active phase, - * of pressure values in each cell in the current - * equilibration region. - */ - template - std::vector< std::vector > - phasePressures(const Grid& G, - const Region& reg, - const CellRange& cells, - const double grav = unit::gravity); - - - - /** - * Compute initial phase saturations by means of equilibration. - * - * \tparam FluidSystem The FluidSystem from opm-material - * Must be initialized before used. - * - * \tparam Grid Type of the grid - * - * \tparam Region Type of an equilibration region information - * base. Typically an instance of the EquilReg - * class template. - * - * \tparam CellRange Type of cell range that demarcates the - * cells pertaining to the current - * equilibration region. Must implement - * methods begin() and end() to bound the range - * as well as provide an inner type, - * const_iterator, to traverse the range. - * - * \tparam MaterialLawManager The MaterialLawManager from opm-material - * - * \param[in] G Grid. - * \param[in] reg Current equilibration region. - * \param[in] cells Range that spans the cells of the current - * equilibration region. - * \param[in] materialLawManager The MaterialLawManager from opm-material - * \param[in] swat_init A vector of initial water saturations. - * The capillary pressure is scaled to fit these values - * \param[in] phase_pressures Phase pressures, one vector for each active phase, - * of pressure values in each cell in the current - * equilibration region. - * \return Phase saturations, one vector for each phase, each containing - * one saturation value per cell in the region. - */ - template - std::vector< std::vector > - phaseSaturations(const Grid& grid, - const Region& reg, - const CellRange& cells, - MaterialLawManager& materialLawManager, - const std::vector swat_init, - std::vector< std::vector >& phase_pressures); - - - - /** - * Compute initial Rs values. - * - * \tparam CellRangeType Type of cell range that demarcates the - * cells pertaining to the current - * equilibration region. Must implement - * methods begin() and end() to bound the range - * as well as provide an inner type, - * const_iterator, to traverse the range. - * - * \param[in] grid Grid. - * \param[in] cells Range that spans the cells of the current - * equilibration region. - * \param[in] oil_pressure Oil pressure for each cell in range. - * \param[in] temperature Temperature for each cell in range. - * \param[in] rs_func Rs as function of pressure and depth. - * \return Rs values, one for each cell in the 'cells' range. - */ - template - std::vector computeRs(const Grid& grid, - const CellRangeType& cells, - const std::vector oil_pressure, - const std::vector& temperature, - const Miscibility::RsFunction& rs_func, - const std::vector gas_saturation); - - namespace DeckDependent { - inline - std::vector - getEquil(const Opm::EclipseState& state) - { - const auto& init = state.getInitConfig(); - - if( !init.hasEquil() ) { - OPM_THROW(std::domain_error, "Deck does not provide equilibration data."); - } - - const auto& equil = init.getEquil(); - return { equil.begin(), equil.end() }; - } - - template - inline - std::vector - equilnum(const Opm::EclipseState& eclipseState, - const Grid& G ) - { - std::vector eqlnum; - if (eclipseState.get3DProperties().hasDeckIntGridProperty("EQLNUM")) { - const int nc = UgGridHelpers::numCells(G); - eqlnum.resize(nc); - const std::vector& e = - eclipseState.get3DProperties().getIntGridProperty("EQLNUM").getData(); - const int* gc = UgGridHelpers::globalCell(G); - for (int cell = 0; cell < nc; ++cell) { - const int deck_pos = (gc == NULL) ? cell : gc[cell]; - eqlnum[cell] = e[deck_pos] - 1; - } - } - else { - // No explicit equilibration region. - // All cells in region zero. - eqlnum.assign(UgGridHelpers::numCells(G), 0); - } - - return eqlnum; - } - - template - class InitialStateComputer { - public: - template - InitialStateComputer(MaterialLawManager& materialLawManager, - const Opm::EclipseState& eclipseState, - const Grid& G , - const double grav = unit::gravity, - const bool applySwatInit = true - ) - : pp_(FluidSystem::numPhases, - std::vector(UgGridHelpers::numCells(G))), - sat_(FluidSystem::numPhases, - std::vector(UgGridHelpers::numCells(G))), - rs_(UgGridHelpers::numCells(G)), - rv_(UgGridHelpers::numCells(G)) - { - //Check for presence of kw SWATINIT - if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatInit) { - const std::vector& swat_init_ecl = eclipseState. - get3DProperties().getDoubleGridProperty("SWATINIT").getData(); - const int nc = UgGridHelpers::numCells(G); - swat_init_.resize(nc); - const int* gc = UgGridHelpers::globalCell(G); - for (int c = 0; c < nc; ++c) { - const int deck_pos = (gc == NULL) ? c : gc[c]; - swat_init_[c] = swat_init_ecl[deck_pos]; - } - } - // Get the equilibration records. - const std::vector rec = getEquil(eclipseState); - const auto& tables = eclipseState.getTableManager(); - // Create (inverse) region mapping. - const RegionMapping<> eqlmap(equilnum(eclipseState, G)); - const int invalidRegion = -1; - regionPvtIdx_.resize(rec.size(), invalidRegion); - setRegionPvtIdx(G, eclipseState, eqlmap); - - // Create Rs functions. - rs_func_.reserve(rec.size()); - if (FluidSystem::enableDissolvedGas()) { - const TableContainer& rsvdTables = tables.getRsvdTables(); - for (size_t i = 0; i < rec.size(); ++i) { - if (eqlmap.cells(i).empty()) - { - rs_func_.push_back(std::shared_ptr>()); - continue; - } - const int pvtIdx = regionPvtIdx_[i]; - if (!rec[i].liveOilInitConstantRs()) { - if (rsvdTables.size() <= 0 ) { - OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table not available."); - } - const RsvdTable& rsvdTable = rsvdTables.getTable(i); - std::vector depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy(); - std::vector rsColumn = rsvdTable.getColumn("RS").vectorCopy(); - rs_func_.push_back(std::make_shared>(pvtIdx, - depthColumn , rsColumn)); - } else { - if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { - OPM_THROW(std::runtime_error, - "Cannot initialise: when no explicit RSVD table is given, \n" - "datum depth must be at the gas-oil-contact. " - "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold."); - } - const double p_contact = rec[i].datumDepthPressure(); - const double T_contact = 273.15 + 20; // standard temperature for now - rs_func_.push_back(std::make_shared>(pvtIdx, p_contact, T_contact)); - } - } - } else { - for (size_t i = 0; i < rec.size(); ++i) { - rs_func_.push_back(std::make_shared()); - } - } - - rv_func_.reserve(rec.size()); - if (FluidSystem::enableVaporizedOil()) { - const TableContainer& rvvdTables = tables.getRvvdTables(); - for (size_t i = 0; i < rec.size(); ++i) { - if (eqlmap.cells(i).empty()) - { - rv_func_.push_back(std::shared_ptr>()); - continue; - } - const int pvtIdx = regionPvtIdx_[i]; - if (!rec[i].wetGasInitConstantRv()) { - if (rvvdTables.size() <= 0) { - OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table not available."); - } - - const RvvdTable& rvvdTable = rvvdTables.getTable(i); - std::vector depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy(); - std::vector rvColumn = rvvdTable.getColumn("RV").vectorCopy(); - rv_func_.push_back(std::make_shared>(pvtIdx, - depthColumn , rvColumn)); - - } else { - if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { - OPM_THROW(std::runtime_error, - "Cannot initialise: when no explicit RVVD table is given, \n" - "datum depth must be at the gas-oil-contact. " - "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold."); - } - const double p_contact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); - const double T_contact = 273.15 + 20; // standard temperature for now - rv_func_.push_back(std::make_shared>(pvtIdx ,p_contact, T_contact)); - } - } - } else { - for (size_t i = 0; i < rec.size(); ++i) { - rv_func_.push_back(std::make_shared()); - } - } - - // Compute pressures, saturations, rs and rv factors. - calcPressSatRsRv(eqlmap, rec, materialLawManager, G, grav); - - // Modify oil pressure in no-oil regions so that the pressures of present phases can - // be recovered from the oil pressure and capillary relations. - } - - typedef std::vector Vec; - typedef std::vector PVec; // One per phase. - - const PVec& press() const { return pp_; } - const PVec& saturation() const { return sat_; } - const Vec& rs() const { return rs_; } - const Vec& rv() const { return rv_; } - - private: - typedef EquilReg EqReg; - std::vector< std::shared_ptr > rs_func_; - std::vector< std::shared_ptr > rv_func_; - std::vector regionPvtIdx_; - PVec pp_; - PVec sat_; - Vec rs_; - Vec rv_; - Vec swat_init_; - - template - void setRegionPvtIdx(const Grid& G, const Opm::EclipseState& eclipseState, const RMap& reg) { - - std::vector cellPvtRegionIdx; - extractPvtTableIndex(cellPvtRegionIdx, eclipseState, UgGridHelpers::numCells(G), UgGridHelpers::globalCell(G)); - for (const auto& r : reg.activeRegions()) { - const auto& cells = reg.cells(r); - const int cell = *(cells.begin()); - regionPvtIdx_[r] = cellPvtRegionIdx[cell]; - } - } - - template - void - calcPressSatRsRv(const RMap& reg , - const std::vector< EquilRecord >& rec , - MaterialLawManager& materialLawManager, - const Grid& G , - const double grav) - { - for (const auto& r : reg.activeRegions()) { - const auto& cells = reg.cells(r); - if (cells.empty()) - { - OpmLog::warning("Equilibration region " + std::to_string(r + 1) - + " has no active cells"); - continue; - } - - const EqReg eqreg(rec[r], rs_func_[r], rv_func_[r], regionPvtIdx_[r]); - - PVec pressures = phasePressures(G, eqreg, cells, grav); - const std::vector& temp = temperature(G, eqreg, cells); - const PVec sat = phaseSaturations(G, eqreg, cells, materialLawManager, swat_init_, pressures); - - const int np = FluidSystem::numPhases; - for (int p = 0; p < np; ++p) { - copyFromRegion(pressures[p], cells, pp_[p]); - copyFromRegion(sat[p], cells, sat_[p]); - } - const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); - const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - if (oil && gas) { - const int oilpos = FluidSystem::oilPhaseIdx; - const int gaspos = FluidSystem::gasPhaseIdx; - const Vec rs_vals = computeRs(G, cells, pressures[oilpos], temp, *(rs_func_[r]), sat[gaspos]); - const Vec rv_vals = computeRs(G, cells, pressures[gaspos], temp, *(rv_func_[r]), sat[oilpos]); - copyFromRegion(rs_vals, cells, rs_); - copyFromRegion(rv_vals, cells, rv_); - } - } - } - - template - void copyFromRegion(const Vec& source, - const CellRangeType& cells, - Vec& destination) - { - auto s = source.begin(); - auto c = cells.begin(); - const auto e = cells.end(); - for (; c != e; ++c, ++s) { - destination[*c] = *s; - } - } - - }; - } // namespace DeckDependent - } // namespace EQUIL -} // namespace Opm - -#include - -#endif // OPM_INITSTATEEQUIL_HEADER_INCLUDED diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp deleted file mode 100644 index ef5d69962..000000000 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ /dev/null @@ -1,791 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services - Copyright 2015 NTNU - Copyright 2017 IRIS - - 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 . -*/ - -#ifndef OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED -#define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED - -#include -#include - -#include - -#include -#include -#include -#include - -namespace Opm -{ - namespace Details { - template - class RK4IVP { - public: - RK4IVP(const RHS& f , - const std::array& span, - const double y0 , - const int N ) - : N_(N) - , span_(span) - { - const double h = stepsize(); - const double h2 = h / 2; - const double 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 double x = span_[0] + i*h; - const double y = y_.back(); - - const double k1 = f_[i]; - const double k2 = f(x + h2, y + h2*k1); - const double k3 = f(x + h2, y + h2*k2); - const double k4 = f(x + h , y + h*k3); - - y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4)); - f_.push_back(f(x + h, y_.back())); - } - - assert (y_.size() == std::vector::size_type(N + 1)); - } - - double - operator()(const double x) const - { - // Dense output (O(h**3)) according to Shampine - // (Hermite interpolation) - const double h = stepsize(); - int i = (x - span_[0]) / h; - const double t = (x - (span_[0] + i*h)) / h; - - // Crude handling of evaluation point outside "span_"; - if (i < 0) { i = 0; } - if (N_ <= i) { i = N_ - 1; } - - const double y0 = y_[i], y1 = y_[i + 1]; - const double f0 = f_[i], f1 = f_[i + 1]; - - double 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; - } - - private: - int N_; - std::array span_; - std::vector y_; - std::vector f_; - - double - stepsize() const { return (span_[1] - span_[0]) / N_; } - }; - - namespace PhasePressODE { - template - class Water { - public: - Water(const double temp, - const int pvtRegionIdx, - const double norm_grav) - : temp_(temp) - , pvtRegionIdx_(pvtRegionIdx) - , g_(norm_grav) - { - } - - double - operator()(const double /* depth */, - const double press) const - { - return this->density(press) * g_; - } - - private: - const double temp_; - const int pvtRegionIdx_; - const double g_; - - double - density(const double press) const - { - double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); - return rho; - } - }; - - template - class Oil { - public: - Oil(const double temp, - const RS& rs, - const int pvtRegionIdx, - const double norm_grav) - : temp_(temp) - , rs_(rs) - , pvtRegionIdx_(pvtRegionIdx) - , g_(norm_grav) - { - } - - double - operator()(const double depth, - const double press) const - { - return this->density(depth, press) * g_; - } - - private: - const double temp_; - const RS& rs_; - const int pvtRegionIdx_; - const double g_; - - double - density(const double depth, - const double press) const - { - double rs = rs_(depth, press, temp_); - double bOil = 0.0; - if ( !FluidSystem::enableDissolvedGas() || rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press) ) { - bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - } else { - bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rs); - } - double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); - if (FluidSystem::enableDissolvedGas()) { - rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - } - - return rho; - } - }; - - template - class Gas { - public: - Gas(const double temp, - const RV& rv, - const int pvtRegionIdx, - const double norm_grav) - : temp_(temp) - , rv_(rv) - , pvtRegionIdx_(pvtRegionIdx) - , g_(norm_grav) - { - } - - double - operator()(const double depth, - const double press) const - { - return this->density(depth, press) * g_; - } - - private: - const double temp_; - const RV& rv_; - const int pvtRegionIdx_; - const double g_; - - double - density(const double depth, - const double press) const - { - double rv = rv_(depth, press, temp_); - double bGas = 0.0; - if ( !FluidSystem::enableVaporizedOil() || rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press) ) { - bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - } else { - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv); - } - double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - if (FluidSystem::enableVaporizedOil()) { - rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); - } - - return rho; - } - }; - } // namespace PhasePressODE - - - namespace PhasePressure { - template - void - assign(const Grid& G , - const std::array& f , - const double split, - const CellRange& cells, - std::vector& p ) - { - - enum { up = 0, down = 1 }; - - std::vector::size_type c = 0; - for (typename CellRange::const_iterator - ci = cells.begin(), ce = cells.end(); - ci != ce; ++ci, ++c) - { - assert (c < p.size()); - - const double z = UgGridHelpers::cellCenterDepth(G, *ci); - p[c] = (z < split) ? f[up](z) : f[down](z); - } - } - - template < class FluidSystem, - class Grid, - class Region, - class CellRange> - void - water(const Grid& G , - const Region& reg , - const std::array& span , - const double grav , - double& po_woc, - const CellRange& cells , - std::vector& press ) - { - using PhasePressODE::Water; - typedef Water ODE; - - const double T = 273.15 + 20; // standard temperature for now - ODE drho(T, reg.pvtIdx() , grav); - - double z0; - double p0; - if (reg.datum() > reg.zwoc()) {//Datum in water zone - z0 = reg.datum(); - p0 = reg.pressure(); - } else { - z0 = reg.zwoc(); - p0 = po_woc - reg.pcow_woc(); // Water pressure at contact - } - - std::array up = {{ z0, span[0] }}; - std::array down = {{ z0, span[1] }}; - - typedef Details::RK4IVP WPress; - std::array wpress = { - { - WPress(drho, up , p0, 2000) - , - WPress(drho, down, p0, 2000) - } - }; - - assign(G, wpress, z0, cells, press); - - if (reg.datum() > reg.zwoc()) { - // Return oil pressure at contact - po_woc = wpress[0](reg.zwoc()) + reg.pcow_woc(); - } - } - - template - void - oil(const Grid& G , - const Region& reg , - const std::array& span , - const double grav , - const CellRange& cells , - std::vector& press , - double& po_woc, - double& po_goc) - { - using PhasePressODE::Oil; - typedef Oil ODE; - - const double T = 273.15 + 20; // standard temperature for now - ODE drho(T, reg.dissolutionCalculator(), - reg.pvtIdx(), grav); - - double z0; - double p0; - if (reg.datum() > reg.zwoc()) {//Datum in water zone, po_woc given - z0 = reg.zwoc(); - p0 = po_woc; - } else if (reg.datum() < reg.zgoc()) {//Datum in gas zone, po_goc given - z0 = reg.zgoc(); - p0 = po_goc; - } else { //Datum in oil zone - z0 = reg.datum(); - p0 = reg.pressure(); - } - - std::array up = {{ z0, span[0] }}; - std::array down = {{ z0, span[1] }}; - - typedef Details::RK4IVP OPress; - std::array opress = { - { - OPress(drho, up , p0, 2000) - , - OPress(drho, down, p0, 2000) - } - }; - - assign(G, opress, z0, cells, press); - - const double woc = reg.zwoc(); - if (z0 > woc) { po_woc = opress[0](woc); } // WOC above datum - else if (z0 < woc) { po_woc = opress[1](woc); } // WOC below datum - else { po_woc = p0; } // WOC *at* datum - - const double goc = reg.zgoc(); - if (z0 > goc) { po_goc = opress[0](goc); } // GOC above datum - else if (z0 < goc) { po_goc = opress[1](goc); } // GOC below datum - else { po_goc = p0; } // GOC *at* datum - } - - template - void - gas(const Grid& G , - const Region& reg , - const std::array& span , - const double grav , - double& po_goc, - const CellRange& cells , - std::vector& press ) - { - using PhasePressODE::Gas; - typedef Gas ODE; - - const double T = 273.15 + 20; // standard temperature for now - ODE drho(T, reg.evaporationCalculator(), - reg.pvtIdx(), grav); - - double z0; - double p0; - if (reg.datum() < reg.zgoc()) {//Datum in gas zone - z0 = reg.datum(); - p0 = reg.pressure(); - } else { - z0 = reg.zgoc(); - p0 = po_goc + reg.pcgo_goc(); // Gas pressure at contact - } - - std::array up = {{ z0, span[0] }}; - std::array down = {{ z0, span[1] }}; - - typedef Details::RK4IVP GPress; - std::array gpress = { - { - GPress(drho, up , p0, 2000) - , - GPress(drho, down, p0, 2000) - } - }; - - assign(G, gpress, z0, cells, press); - - if (reg.datum() < reg.zgoc()) { - // Return oil pressure at contact - po_goc = gpress[1](reg.zgoc()) - reg.pcgo_goc(); - } - } - } // namespace PhasePressure - - template - void - equilibrateOWG(const Grid& G, - const Region& reg, - const double grav, - const std::array& span, - const CellRange& cells, - std::vector< std::vector >& press) - { - const bool water = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); - const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); - const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - const int oilpos = FluidSystem::oilPhaseIdx; - const int waterpos = FluidSystem::waterPhaseIdx; - const int gaspos = FluidSystem::gasPhaseIdx; - - if (reg.datum() > reg.zwoc()) { // Datum in water zone - double po_woc = -1; - double po_goc = -1; - - if (water) { - PhasePressure::water(G, reg, span, grav, po_woc, - cells, press[ waterpos ]); - } - - if (oil) { - PhasePressure::oil(G, reg, span, grav, cells, - press[ oilpos ], po_woc, po_goc); - } - - if (gas) { - PhasePressure::gas(G, reg, span, grav, po_goc, - cells, press[ gaspos ]); - } - } else if (reg.datum() < reg.zgoc()) { // Datum in gas zone - double po_woc = -1; - double po_goc = -1; - - if (gas) { - PhasePressure::gas(G, reg, span, grav, po_goc, - cells, press[ gaspos ]); - } - - if (oil) { - PhasePressure::oil(G, reg, span, grav, cells, - press[ oilpos ], po_woc, po_goc); - } - - if (water) { - PhasePressure::water(G, reg, span, grav, po_woc, - cells, press[ waterpos ]); - } - } else { // Datum in oil zone - double po_woc = -1; - double po_goc = -1; - - if (oil) { - PhasePressure::oil(G, reg, span, grav, cells, - press[ oilpos ], po_woc, po_goc); - } - - if (water) { - PhasePressure::water(G, reg, span, grav, po_woc, - cells, press[ waterpos ]); - } - - if (gas) { - PhasePressure::gas(G, reg, span, grav, po_goc, - cells, press[ gaspos ]); - } - } - } - } // namespace Details - - - namespace EQUIL { - - - template - std::vector< std::vector > - phasePressures(const Grid& G, - const Region& reg, - const CellRange& cells, - const double grav) - { - std::array span = - {{ std::numeric_limits::max() , - -std::numeric_limits::max() }}; // Symm. about 0. - - int ncell = 0; - { - // This code is only supported in three space dimensions - assert (UgGridHelpers::dimensions(G) == 3); - - const int nd = UgGridHelpers::dimensions(G); - - // Define vertical span as - // - // [minimum(node depth(cells)), maximum(node depth(cells))] - // - // Note: We use a sledgehammer approach--looping all - // the nodes of all the faces of all the 'cells'--to - // compute those bounds. This necessarily entails - // visiting some nodes (and faces) multiple times. - // - // Note: The implementation of 'RK4IVP<>' implicitly - // imposes the requirement that cell centroids are all - // within this vertical span. That requirement is not - // checked. - auto cell2Faces = UgGridHelpers::cell2Faces(G); - auto faceVertices = UgGridHelpers::face2Vertices(G); - - for (typename CellRange::const_iterator - ci = cells.begin(), ce = cells.end(); - ci != ce; ++ci, ++ncell) - { - for (auto fi=cell2Faces[*ci].begin(), - fe=cell2Faces[*ci].end(); - fi != fe; - ++fi) - { - for (auto i = faceVertices[*fi].begin(), e = faceVertices[*fi].end(); - i != e; ++i) - { - const double z = UgGridHelpers::vertexCoordinates(G, *i)[nd-1]; - - if (z < span[0]) { span[0] = z; } - if (z > span[1]) { span[1] = z; } - } - } - } - } - const int np = FluidSystem::numPhases; //reg.phaseUsage().num_phases; - - typedef std::vector pval; - std::vector press(np, pval(ncell, 0.0)); - - const double zwoc = reg.zwoc (); - const double zgoc = reg.zgoc (); - - // make sure goc and woc is within the span for the phase pressure calculation - span[0] = std::min(span[0],zgoc); - span[1] = std::max(span[1],zwoc); - - Details::equilibrateOWG(G, reg, grav, span, cells, press); - - return press; - } - - template - std::vector - temperature(const Grid& /* G */, - const Region& /* reg */, - const CellRange& cells) - { - // use the standard temperature for everything for now - return std::vector(cells.size(), 273.15 + 20.0); - } - - template - std::vector< std::vector > - phaseSaturations(const Grid& G, - const Region& reg, - const CellRange& cells, - MaterialLawManager& materialLawManager, - const std::vector swat_init, - std::vector< std::vector >& phase_pressures) - { - if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - OPM_THROW(std::runtime_error, "Cannot initialise: not handling water-gas cases."); - } - - std::vector< std::vector > phase_saturations = phase_pressures; // Just to get the right size. - - // Adjust oil pressure according to gas saturation and cap pressure - typedef Opm::SimpleModularFluidState SatOnlyFluidState; - - SatOnlyFluidState fluidState; - typedef typename MaterialLawManager::MaterialLaw MaterialLaw; - - const bool water = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); - const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - const int oilpos = FluidSystem::oilPhaseIdx; - const int waterpos = FluidSystem::waterPhaseIdx; - const int gaspos = FluidSystem::gasPhaseIdx; - std::vector::size_type local_index = 0; - for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) { - const int cell = *ci; - const auto& scaledDrainageInfo = - materialLawManager.oilWaterScaledEpsInfoDrainage(cell); - const auto& matParams = materialLawManager.materialLawParams(cell); - - // Find saturations from pressure differences by - // inverting capillary pressure functions. - double sw = 0.0; - if (water) { - if (isConstPc(materialLawManager,FluidSystem::waterPhaseIdx, cell)){ - const double cellDepth = UgGridHelpers::cellCenterDepth(G, - cell); - sw = satFromDepth(materialLawManager,cellDepth,reg.zwoc(),waterpos,cell,false); - phase_saturations[waterpos][local_index] = sw; - } - else{ - const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index]; - if (swat_init.empty()) { // Invert Pc to find sw - sw = satFromPc(materialLawManager, waterpos, cell, pcov); - phase_saturations[waterpos][local_index] = sw; - } else { // Scale Pc to reflect imposed sw - sw = swat_init[cell]; - sw = materialLawManager.applySwatinit(cell, pcov, sw); - phase_saturations[waterpos][local_index] = sw; - } - } - } - double sg = 0.0; - if (gas) { - if (isConstPc(materialLawManager,FluidSystem::gasPhaseIdx,cell)){ - const double cellDepth = UgGridHelpers::cellCenterDepth(G, - cell); - sg = satFromDepth(materialLawManager,cellDepth,reg.zgoc(),gaspos,cell,true); - phase_saturations[gaspos][local_index] = sg; - } - else{ - // Note that pcog is defined to be (pg - po), not (po - pg). - const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index]; - const double increasing = true; // pcog(sg) expected to be increasing function - sg = satFromPc(materialLawManager, gaspos, cell, pcog, increasing); - phase_saturations[gaspos][local_index] = sg; - } - } - if (gas && water && (sg + sw > 1.0)) { - // Overlapping gas-oil and oil-water transition - // zones can lead to unphysical saturations when - // treated as above. Must recalculate using gas-water - // capillary pressure. - const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index]; - if (! swat_init.empty()) { - // Re-scale Pc to reflect imposed sw for vanishing oil phase. - // This seems consistent with ecl, and fails to honour - // swat_init in case of non-trivial gas-oil cap pressure. - sw = materialLawManager.applySwatinit(cell, pcgw, sw); - } - sw = satFromSumOfPcs(materialLawManager, waterpos, gaspos, cell, pcgw); - sg = 1.0 - sw; - phase_saturations[waterpos][local_index] = sw; - phase_saturations[gaspos][local_index] = sg; - if ( water ) { - fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw); - } - else { - fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); - } - fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0 - sw - sg); - fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg); - - double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 }; - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx]; - phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pcGas; - } - phase_saturations[oilpos][local_index] = 1.0 - sw - sg; - - // Adjust phase pressures for max and min saturation ... - double threshold_sat = 1.0e-6; - - double so = 1.0; - double pC[FluidSystem::numPhases] = { 0.0, 0.0, 0.0 }; - if (water) { - double swu = scaledDrainageInfo.Swu; - fluidState.setSaturation(FluidSystem::waterPhaseIdx, swu); - so -= swu; - } - if (gas) { - double sgu = scaledDrainageInfo.Sgu; - fluidState.setSaturation(FluidSystem::gasPhaseIdx, sgu); - so-= sgu; - } - fluidState.setSaturation(FluidSystem::oilPhaseIdx, so); - - if (water && sw > scaledDrainageInfo.Swu-threshold_sat ) { - fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swu); - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx]; - phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pcWat; - } else if (gas && sg > scaledDrainageInfo.Sgu-threshold_sat) { - fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgu); - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx]; - phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pcGas; - } - if (gas && sg < scaledDrainageInfo.Sgl+threshold_sat) { - fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgl); - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx]; - phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pcGas; - } - if (water && sw < scaledDrainageInfo.Swl+threshold_sat) { - fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swl); - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx]; - phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pcWat; - } - } - return phase_saturations; - } - - - /** - * Compute initial Rs values. - * - * \tparam CellRangeType Type of cell range that demarcates the - * cells pertaining to the current - * equilibration region. Must implement - * methods begin() and end() to bound the range - * as well as provide an inner type, - * const_iterator, to traverse the range. - * - * \param[in] grid Grid. - * \param[in] cells Range that spans the cells of the current - * equilibration region. - * \param[in] oil_pressure Oil pressure for each cell in range. - * \param[in] temperature Temperature for each cell in range. - * \param[in] rs_func Rs as function of pressure and depth. - * \return Rs values, one for each cell in the 'cells' range. - */ - template - std::vector computeRs(const Grid& grid, - const CellRangeType& cells, - const std::vector oil_pressure, - const std::vector& temperature, - const Miscibility::RsFunction& rs_func, - const std::vector gas_saturation) - { - assert(UgGridHelpers::dimensions(grid) == 3); - std::vector rs(cells.size()); - int count = 0; - for (auto it = cells.begin(); it != cells.end(); ++it, ++count) { - const double depth = UgGridHelpers::cellCenterDepth(grid, *it); - rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]); - } - return rs; - } - - } // namespace Equil - - -} // namespace Opm - -#endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED diff --git a/opm/core/utility/RegionMapping.hpp b/opm/core/utility/RegionMapping.hpp deleted file mode 100644 index 63d412b0f..000000000 --- a/opm/core/utility/RegionMapping.hpp +++ /dev/null @@ -1,181 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - - 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 . -*/ - -#ifndef OPM_REGIONMAPPING_HEADER_INCLUDED -#define OPM_REGIONMAPPING_HEADER_INCLUDED - -#include - -#include -#include - -namespace Opm -{ - - /** - * Forward and reverse mappings between cells and - * regions/partitions (e.g., the ECLIPSE-style 'SATNUM', - * 'PVTNUM', or 'EQUILNUM' arrays). - * - * \tparam Region Type of a forward region mapping. Expected - * to provide indexed access through - * operator[]() as well as inner types - * 'value_type', 'size_type', and - * 'const_iterator'. - */ - template < class Region = std::vector > - class RegionMapping { - public: - /** - * Constructor. - * - * \param[in] reg Forward region mapping, restricted to - * active cells only. - */ - explicit - RegionMapping(const Region& reg) - : reg_(reg) - { - rev_.init(reg_); - } - - /** - * Type of forward (cell-to-region) mapping result. - * Expected to be an integer. - */ - typedef typename Region::value_type RegionId; - - /** - * Type of reverse (region-to-cell) mapping (element) - * result. - */ - typedef typename Region::size_type CellId; - - /** - * Type of reverse region-to-cell range bounds and - * iterators. - */ - typedef typename std::vector::const_iterator CellIter; - - typedef boost::iterator_range Range; - - /** - * Compute region number of given active cell. - * - * \param[in] c Active cell - * \return Region to which @c c belongs. - */ - RegionId - region(const CellId c) const { return reg_[c]; } - - const std::vector& - activeRegions() const - { - return rev_.active; - } - - /** - * Extract active cells in particular region. - * - * \param[in] r Region number - * - * \return Range of active cells in region @c r. Empty if @c r is - * not an active region. - */ - Range - cells(const RegionId r) const { - const auto id = rev_.binid.find(r); - - if (id == rev_.binid.end()) { - // Region 'r' not an active region. Return empty. - return Range(rev_.c.end(), rev_.c.end()); - } - - const auto i = id->second; - - return Range(rev_.c.begin() + rev_.p[i + 0], - rev_.c.begin() + rev_.p[i + 1]); - } - - private: - /** - * Copy of forward region mapping (cell-to-region). - */ - Region reg_; - - /** - * Reverse mapping (region-to-cell). - */ - struct { - typedef typename std::vector::size_type Pos; - - std::unordered_map binid; - std::vector active; - - std::vector p; /**< Region start pointers */ - std::vector c; /**< Region cells */ - - /** - * Compute reverse mapping. Standard linear insertion - * sort algorithm. - */ - void - init(const Region& reg) - { - binid.clear(); - for (const auto& r : reg) { - ++binid[r]; - } - - p .clear(); p.emplace_back(0); - active.clear(); - { - Pos n = 0; - for (auto& id : binid) { - active.push_back(id.first); - p .push_back(id.second); - - id.second = n++; - } - } - - for (decltype(p.size()) i = 1, sz = p.size(); i < sz; ++i) { - p[0] += p[i]; - p[i] = p[0] - p[i]; - } - - assert (p[0] == static_cast(reg.size())); - - c.resize(reg.size()); - { - CellId i = 0; - for (const auto& r : reg) { - auto& pos = p[ binid[r] + 1 ]; - c[ pos++ ] = i++; - } - } - - p[0] = 0; - } - } rev_; /**< Reverse mapping instance */ - }; - -} // namespace Opm - -#endif // OPM_REGIONMAPPING_HEADER_INCLUDED