mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Moved equilibration utilities to separate file.
This commit is contained in:
Atgeirr Flø Rasmussen
Andreas Lauser
parent
e2913a178c
commit
e1dac5269b
631
opm/core/simulator/EquilibrationHelpers.hpp
Normal file
631
opm/core/simulator/EquilibrationHelpers.hpp
Normal file
@ -0,0 +1,631 @@
|
||||
/*
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
|
||||
#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/utility/linearInterpolation.hpp>
|
||||
#include <opm/core/utility/RegionMapping.hpp>
|
||||
#include <opm/core/utility/RootFinders.hpp>
|
||||
|
||||
|
||||
/*
|
||||
---- synopsis of EquilibrationHelpers.hpp ----
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
namespace Equil {
|
||||
|
||||
template <class Props>
|
||||
class DensityCalculator;
|
||||
|
||||
template <>
|
||||
class DensityCalculator< BlackoilPropertiesInterface >;
|
||||
|
||||
namespace Miscibility {
|
||||
struct NoMixing;
|
||||
class RsVD;
|
||||
class RsSatAtContact;
|
||||
}
|
||||
|
||||
struct EquilRecord;
|
||||
|
||||
template <class DensCalc,
|
||||
class RS = Miscibility::NoMixing,
|
||||
class RV = Miscibility::NoMixing>
|
||||
class EquilReg;
|
||||
|
||||
|
||||
struct PcEq;
|
||||
|
||||
inline double satFromPc(const BlackoilPropertiesInterface& props,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double target_pc,
|
||||
const bool increasing = false);
|
||||
struct PcEqSum
|
||||
inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
|
||||
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 {
|
||||
|
||||
|
||||
template <class Props>
|
||||
class DensityCalculator;
|
||||
|
||||
/**
|
||||
* Facility for calculating phase densities based on the
|
||||
* BlackoilPropertiesInterface.
|
||||
*
|
||||
* Implements the crucial <CODE>operator()(p,svol)</CODE>
|
||||
* function that is expected by class EquilReg.
|
||||
*/
|
||||
template <>
|
||||
class DensityCalculator< BlackoilPropertiesInterface > {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] props Implementation of the
|
||||
* BlackoilPropertiesInterface.
|
||||
*
|
||||
* \param[in] c Single cell used as a representative cell
|
||||
* in a PVT region.
|
||||
*/
|
||||
DensityCalculator(const BlackoilPropertiesInterface& props,
|
||||
const int c)
|
||||
: props_(props)
|
||||
, c_(1, c)
|
||||
{
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute phase densities of all phases at phase point
|
||||
* given by (pressure, surface volume) tuple.
|
||||
*
|
||||
* \param[in] p Fluid pressure.
|
||||
*
|
||||
* \param[in] z Surface volumes of all phases.
|
||||
*
|
||||
* \return Phase densities at phase point.
|
||||
*/
|
||||
std::vector<double>
|
||||
operator()(const double p,
|
||||
const std::vector<double>& z) const
|
||||
{
|
||||
const int np = props_.numPhases();
|
||||
std::vector<double> A(np * np, 0);
|
||||
|
||||
assert (z.size() == std::vector<double>::size_type(np));
|
||||
|
||||
double* dAdp = 0;
|
||||
props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
|
||||
|
||||
std::vector<double> rho(np, 0.0);
|
||||
props_.density(1, &A[0], &rho[0]);
|
||||
|
||||
return rho;
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const std::vector<int> c_;
|
||||
};
|
||||
|
||||
/**
|
||||
* Types and routines relating to phase mixing in
|
||||
* equilibration calculations.
|
||||
*/
|
||||
namespace Miscibility {
|
||||
/**
|
||||
* Type that implements "no phase mixing" policy.
|
||||
*/
|
||||
struct NoMixing {
|
||||
/**
|
||||
* Function call.
|
||||
*
|
||||
* \param[in] depth Depth at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \param[in] press Pressure 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
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
};
|
||||
|
||||
/**
|
||||
* Type that implements "dissolved gas-oil ratio"
|
||||
* tabulated as a function of depth policy. Data
|
||||
* typically taken from keyword 'RSVD'.
|
||||
*/
|
||||
class RsVD {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] depth Depth nodes.
|
||||
* \param[in] rs Dissolved gas-oil ratio at @c depth.
|
||||
*/
|
||||
RsVD(const std::vector<double>& depth,
|
||||
const std::vector<double>& rs)
|
||||
: 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.
|
||||
*
|
||||
* \return Dissolved gas-oil ratio (RS) at depth @c
|
||||
* depth and pressure @c press.
|
||||
*/
|
||||
double
|
||||
operator()(const double depth,
|
||||
const double /* press */) const
|
||||
{
|
||||
return linearInterpolation(depth_, rs_, depth);
|
||||
}
|
||||
|
||||
private:
|
||||
std::vector<double> depth_; /**< Depth nodes */
|
||||
std::vector<double> rs_; /**< Dissolved gas-oil ratio */
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* 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.
|
||||
*/
|
||||
class RsSatAtContact {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] props property object
|
||||
* \param[in] cell any cell in the pvt region
|
||||
* \param[in] p_contact oil pressure at the contact
|
||||
*/
|
||||
RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
|
||||
: props_(props), cell_(cell)
|
||||
{
|
||||
auto pu = props_.phaseUsage();
|
||||
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
|
||||
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
|
||||
rs_sat_contact_ = satRs(p_contact);
|
||||
}
|
||||
|
||||
/**
|
||||
* Function call.
|
||||
*
|
||||
* \param[in] depth Depth at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \param[in] press Pressure 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
|
||||
{
|
||||
return std::max(satRs(press), rs_sat_contact_);
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int cell_;
|
||||
double z_[BlackoilPhases::MaxNumPhases];
|
||||
double rs_sat_contact_;
|
||||
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
|
||||
|
||||
double satRs(const double press) const
|
||||
{
|
||||
props_.matrix(1, &press, z_, &cell_, A_, 0);
|
||||
// Rs/Bo is in the gas row and oil column of A_.
|
||||
// 1/Bo is in the oil row and column.
|
||||
// Recall also that it is stored in column-major order.
|
||||
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
|
||||
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
|
||||
const int np = props_.numPhases();
|
||||
return A_[np*opos + gpos] / A_[np*opos + opos];
|
||||
}
|
||||
};
|
||||
} // namespace Miscibility
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Equilibration record.
|
||||
*
|
||||
* Layout and contents inspired by first six items of
|
||||
* ECLIPSE's 'EQUIL' records. This is the minimum amount of
|
||||
* input data needed to define phase pressures in an
|
||||
* equilibration region.
|
||||
*
|
||||
* Data consists of three pairs of depth and pressure values:
|
||||
* 1. main
|
||||
* - @c depth Main datum depth.
|
||||
* - @c press Pressure at datum depth.
|
||||
*
|
||||
* 2. woc
|
||||
* - @c depth Depth of water-oil contact
|
||||
* - @c press water-oil capillary pressure at water-oil contact.
|
||||
* Capillary pressure defined as "P_oil - P_water".
|
||||
*
|
||||
* 3. goc
|
||||
* - @c depth Depth of gas-oil contact
|
||||
* - @c press Gas-oil capillary pressure at gas-oil contact.
|
||||
* Capillary pressure defined as "P_gas - P_oil".
|
||||
*/
|
||||
struct EquilRecord {
|
||||
struct {
|
||||
double depth;
|
||||
double press;
|
||||
} main, woc, goc;
|
||||
};
|
||||
|
||||
/**
|
||||
* 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
|
||||
* <CODE>
|
||||
* std::vector<double>
|
||||
* operator()(const double press,
|
||||
* const std::vector<double>& svol )
|
||||
* </CODE>
|
||||
* that calculates the phase densities of all phases in @c
|
||||
* svol at fluid pressure @c press.
|
||||
*
|
||||
* \tparam RS Type that provides access to a calculator for
|
||||
* (initial) dissolved gas-oil ratios as a function of depth
|
||||
* and (oil) pressure. Must implement an operator() declared
|
||||
* as
|
||||
* <CODE>
|
||||
* double
|
||||
* operator()(const double depth,
|
||||
* const double press)
|
||||
* </CODE>
|
||||
* that calculates the dissolved gas-oil ratio at depth @c
|
||||
* depth and (oil) pressure @c press.
|
||||
*
|
||||
* \tparam RV Type that provides access to a calculator for
|
||||
* (initial) vapourised oil-gas ratios as a function of depth
|
||||
* and (gas) pressure. Must implement an operator() declared
|
||||
* as
|
||||
* <CODE>
|
||||
* double
|
||||
* operator()(const double depth,
|
||||
* const double press)
|
||||
* </CODE>
|
||||
* that calculates the vapourised oil-gas ratio at depth @c
|
||||
* depth and (gas) pressure @c press.
|
||||
*/
|
||||
template <class DensCalc,
|
||||
class RS = Miscibility::NoMixing,
|
||||
class RV = Miscibility::NoMixing>
|
||||
class EquilReg {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] rec Equilibration data of current region.
|
||||
* \param[in] density Density calculator of current region.
|
||||
* \param[in] rs Calculator of dissolved gas-oil ratio.
|
||||
* \param[in] rv Calculator of vapourised oil-gas ratio.
|
||||
* \param[in] pu Summary of current active phases.
|
||||
*/
|
||||
EquilReg(const EquilRecord& rec,
|
||||
const DensCalc& density,
|
||||
const RS& rs,
|
||||
const RV& rv,
|
||||
const PhaseUsage& pu)
|
||||
: rec_ (rec)
|
||||
, density_(density)
|
||||
, rs_ (rs)
|
||||
, rv_ (rv)
|
||||
, pu_ (pu)
|
||||
{
|
||||
}
|
||||
|
||||
/**
|
||||
* Type of density calculator.
|
||||
*/
|
||||
typedef DensCalc CalcDensity;
|
||||
|
||||
/**
|
||||
* Type of dissolved gas-oil ratio calculator.
|
||||
*/
|
||||
typedef RS CalcDissolution;
|
||||
|
||||
/**
|
||||
* Type of vapourised oil-gas ratio calculator.
|
||||
*/
|
||||
typedef RV CalcEvaporation;
|
||||
|
||||
/**
|
||||
* Datum depth in current region
|
||||
*/
|
||||
double datum() const { return this->rec_.main.depth; }
|
||||
|
||||
/**
|
||||
* Pressure at datum depth in current region.
|
||||
*/
|
||||
double pressure() const { return this->rec_.main.press; }
|
||||
|
||||
/**
|
||||
* Depth of water-oil contact.
|
||||
*/
|
||||
double zwoc() const { return this->rec_.woc .depth; }
|
||||
|
||||
/**
|
||||
* water-oil capillary pressure at water-oil contact.
|
||||
*
|
||||
* \return P_o - P_w at WOC.
|
||||
*/
|
||||
double pcow_woc() const { return this->rec_.woc .press; }
|
||||
|
||||
/**
|
||||
* Depth of gas-oil contact.
|
||||
*/
|
||||
double zgoc() const { return this->rec_.goc .depth; }
|
||||
|
||||
/**
|
||||
* Gas-oil capillary pressure at gas-oil contact.
|
||||
*
|
||||
* \return P_g - P_o at GOC.
|
||||
*/
|
||||
double pcgo_goc() const { return this->rec_.goc .press; }
|
||||
|
||||
/**
|
||||
* Retrieve phase density calculator of current region.
|
||||
*/
|
||||
const CalcDensity&
|
||||
densityCalculator() const { return this->density_; }
|
||||
|
||||
/**
|
||||
* 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 active fluid phase summary.
|
||||
*/
|
||||
const PhaseUsage&
|
||||
phaseUsage() const { return this->pu_; }
|
||||
|
||||
private:
|
||||
EquilRecord rec_; /**< Equilibration data */
|
||||
DensCalc density_; /**< Density calculator */
|
||||
RS rs_; /**< RS calculator */
|
||||
RV rv_; /**< RV calculator */
|
||||
PhaseUsage pu_; /**< Active phase summary */
|
||||
};
|
||||
|
||||
|
||||
|
||||
/// Functor for inverting capillary pressure function.
|
||||
/// Function represented is
|
||||
/// f(s) = pc(s) - target_pc
|
||||
struct PcEq
|
||||
{
|
||||
PcEq(const BlackoilPropertiesInterface& props,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
: props_(props),
|
||||
phase_(phase),
|
||||
cell_(cell),
|
||||
target_pc_(target_pc)
|
||||
{
|
||||
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
}
|
||||
double operator()(double s) const
|
||||
{
|
||||
s_[phase_] = s;
|
||||
props_.capPress(1, s_, &cell_, pc_, 0);
|
||||
return pc_[phase_] - target_pc_;
|
||||
}
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int phase_;
|
||||
const int cell_;
|
||||
const int target_pc_;
|
||||
mutable double s_[BlackoilPhases::MaxNumPhases];
|
||||
mutable double pc_[BlackoilPhases::MaxNumPhases];
|
||||
};
|
||||
|
||||
|
||||
|
||||
/// Compute saturation of some phase corresponding to a given
|
||||
/// capillary pressure.
|
||||
inline double satFromPc(const BlackoilPropertiesInterface& props,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double target_pc,
|
||||
const bool increasing = false)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double sminarr[BlackoilPhases::MaxNumPhases];
|
||||
double smaxarr[BlackoilPhases::MaxNumPhases];
|
||||
props.satRange(1, &cell, sminarr, smaxarr);
|
||||
const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
|
||||
const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
|
||||
|
||||
// Create the equation f(s) = pc(s) - target_pc
|
||||
const PcEq f(props, 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 = 30;
|
||||
const double tol = 1e-6;
|
||||
int iter_used = -1;
|
||||
typedef RegulaFalsi<ThrowOnError> 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
|
||||
struct PcEqSum
|
||||
{
|
||||
PcEqSum(const BlackoilPropertiesInterface& props,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
: props_(props),
|
||||
phase1_(phase1),
|
||||
phase2_(phase2),
|
||||
cell_(cell),
|
||||
target_pc_(target_pc)
|
||||
{
|
||||
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
}
|
||||
double operator()(double s) const
|
||||
{
|
||||
s_[phase1_] = s;
|
||||
s_[phase2_] = 1.0 - s;
|
||||
props_.capPress(1, s_, &cell_, pc_, 0);
|
||||
return pc_[phase1_] + pc_[phase2_] - target_pc_;
|
||||
}
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int phase1_;
|
||||
const int phase2_;
|
||||
const int cell_;
|
||||
const int target_pc_;
|
||||
mutable double s_[BlackoilPhases::MaxNumPhases];
|
||||
mutable double pc_[BlackoilPhases::MaxNumPhases];
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
/// 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.
|
||||
inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double sminarr[BlackoilPhases::MaxNumPhases];
|
||||
double smaxarr[BlackoilPhases::MaxNumPhases];
|
||||
props.satRange(1, &cell, sminarr, smaxarr);
|
||||
const double smin = sminarr[phase1];
|
||||
const double smax = smaxarr[phase1];
|
||||
|
||||
// Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
|
||||
const PcEqSum f(props, 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<ThrowOnError> ScalarSolver;
|
||||
const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used);
|
||||
return sol;
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Equil
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
#endif // OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
|
||||
@ -20,6 +20,7 @@
|
||||
#ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
|
||||
#define OPM_INITSTATEEQUIL_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/simulator/EquilibrationHelpers.hpp>
|
||||
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
|
||||
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
@ -50,545 +51,6 @@ namespace Opm
|
||||
* with other parts of OPM.
|
||||
*/
|
||||
namespace Equil {
|
||||
template <class Props>
|
||||
class DensityCalculator;
|
||||
|
||||
/**
|
||||
* Facility for calculating phase densities based on the
|
||||
* BlackoilPropertiesInterface.
|
||||
*
|
||||
* Implements the crucial <CODE>operator()(p,svol)</CODE>
|
||||
* function that is expected by class EquilReg.
|
||||
*/
|
||||
template <>
|
||||
class DensityCalculator< BlackoilPropertiesInterface > {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] props Implementation of the
|
||||
* BlackoilPropertiesInterface.
|
||||
*
|
||||
* \param[in] c Single cell used as a representative cell
|
||||
* in a PVT region.
|
||||
*/
|
||||
DensityCalculator(const BlackoilPropertiesInterface& props,
|
||||
const int c)
|
||||
: props_(props)
|
||||
, c_(1, c)
|
||||
{
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute phase densities of all phases at phase point
|
||||
* given by (pressure, surface volume) tuple.
|
||||
*
|
||||
* \param[in] p Fluid pressure.
|
||||
*
|
||||
* \param[in] z Surface volumes of all phases.
|
||||
*
|
||||
* \return Phase densities at phase point.
|
||||
*/
|
||||
std::vector<double>
|
||||
operator()(const double p,
|
||||
const std::vector<double>& z) const
|
||||
{
|
||||
const int np = props_.numPhases();
|
||||
std::vector<double> A(np * np, 0);
|
||||
|
||||
assert (z.size() == std::vector<double>::size_type(np));
|
||||
|
||||
double* dAdp = 0;
|
||||
props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
|
||||
|
||||
std::vector<double> rho(np, 0.0);
|
||||
props_.density(1, &A[0], &rho[0]);
|
||||
|
||||
return rho;
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const std::vector<int> c_;
|
||||
};
|
||||
|
||||
/**
|
||||
* Types and routines relating to phase mixing in
|
||||
* equilibration calculations.
|
||||
*/
|
||||
namespace Miscibility {
|
||||
/**
|
||||
* Type that implements "no phase mixing" policy.
|
||||
*/
|
||||
struct NoMixing {
|
||||
/**
|
||||
* Function call.
|
||||
*
|
||||
* \param[in] depth Depth at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \param[in] press Pressure 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
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
};
|
||||
|
||||
/**
|
||||
* Type that implements "dissolved gas-oil ratio"
|
||||
* tabulated as a function of depth policy. Data
|
||||
* typically taken from keyword 'RSVD'.
|
||||
*/
|
||||
class RsVD {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] depth Depth nodes.
|
||||
* \param[in] rs Dissolved gas-oil ratio at @c depth.
|
||||
*/
|
||||
RsVD(const std::vector<double>& depth,
|
||||
const std::vector<double>& rs)
|
||||
: 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.
|
||||
*
|
||||
* \return Dissolved gas-oil ratio (RS) at depth @c
|
||||
* depth and pressure @c press.
|
||||
*/
|
||||
double
|
||||
operator()(const double depth,
|
||||
const double /* press */) const
|
||||
{
|
||||
return linearInterpolation(depth_, rs_, depth);
|
||||
}
|
||||
|
||||
private:
|
||||
std::vector<double> depth_; /**< Depth nodes */
|
||||
std::vector<double> rs_; /**< Dissolved gas-oil ratio */
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* 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.
|
||||
*/
|
||||
class RsSatAtContact {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] props property object
|
||||
* \param[in] cell any cell in the pvt region
|
||||
* \param[in] p_contact oil pressure at the contact
|
||||
*/
|
||||
RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
|
||||
: props_(props), cell_(cell)
|
||||
{
|
||||
auto pu = props_.phaseUsage();
|
||||
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
|
||||
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
|
||||
rs_sat_contact_ = satRs(p_contact);
|
||||
}
|
||||
|
||||
/**
|
||||
* Function call.
|
||||
*
|
||||
* \param[in] depth Depth at which to calculate RS
|
||||
* value.
|
||||
*
|
||||
* \param[in] press Pressure 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
|
||||
{
|
||||
return std::max(satRs(press), rs_sat_contact_);
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int cell_;
|
||||
double z_[BlackoilPhases::MaxNumPhases];
|
||||
double rs_sat_contact_;
|
||||
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
|
||||
|
||||
double satRs(const double press) const
|
||||
{
|
||||
props_.matrix(1, &press, z_, &cell_, A_, 0);
|
||||
// Rs/Bo is in the gas row and oil column of A_.
|
||||
// 1/Bo is in the oil row and column.
|
||||
// Recall also that it is stored in column-major order.
|
||||
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
|
||||
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
|
||||
const int np = props_.numPhases();
|
||||
return A_[np*opos + gpos] / A_[np*opos + opos];
|
||||
}
|
||||
};
|
||||
} // namespace Miscibility
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Equilibration record.
|
||||
*
|
||||
* Layout and contents inspired by first six items of
|
||||
* ECLIPSE's 'EQUIL' records. This is the minimum amount of
|
||||
* input data needed to define phase pressures in an
|
||||
* equilibration region.
|
||||
*
|
||||
* Data consists of three pairs of depth and pressure values:
|
||||
* 1. main
|
||||
* - @c depth Main datum depth.
|
||||
* - @c press Pressure at datum depth.
|
||||
*
|
||||
* 2. woc
|
||||
* - @c depth Depth of water-oil contact
|
||||
* - @c press water-oil capillary pressure at water-oil contact.
|
||||
* Capillary pressure defined as "P_oil - P_water".
|
||||
*
|
||||
* 3. goc
|
||||
* - @c depth Depth of gas-oil contact
|
||||
* - @c press Gas-oil capillary pressure at gas-oil contact.
|
||||
* Capillary pressure defined as "P_gas - P_oil".
|
||||
*/
|
||||
struct EquilRecord {
|
||||
struct {
|
||||
double depth;
|
||||
double press;
|
||||
} main, woc, goc;
|
||||
};
|
||||
|
||||
/**
|
||||
* 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
|
||||
* <CODE>
|
||||
* std::vector<double>
|
||||
* operator()(const double press,
|
||||
* const std::vector<double>& svol )
|
||||
* </CODE>
|
||||
* that calculates the phase densities of all phases in @c
|
||||
* svol at fluid pressure @c press.
|
||||
*
|
||||
* \tparam RS Type that provides access to a calculator for
|
||||
* (initial) dissolved gas-oil ratios as a function of depth
|
||||
* and (oil) pressure. Must implement an operator() declared
|
||||
* as
|
||||
* <CODE>
|
||||
* double
|
||||
* operator()(const double depth,
|
||||
* const double press)
|
||||
* </CODE>
|
||||
* that calculates the dissolved gas-oil ratio at depth @c
|
||||
* depth and (oil) pressure @c press.
|
||||
*
|
||||
* \tparam RV Type that provides access to a calculator for
|
||||
* (initial) vapourised oil-gas ratios as a function of depth
|
||||
* and (gas) pressure. Must implement an operator() declared
|
||||
* as
|
||||
* <CODE>
|
||||
* double
|
||||
* operator()(const double depth,
|
||||
* const double press)
|
||||
* </CODE>
|
||||
* that calculates the vapourised oil-gas ratio at depth @c
|
||||
* depth and (gas) pressure @c press.
|
||||
*/
|
||||
template <class DensCalc,
|
||||
class RS = Miscibility::NoMixing,
|
||||
class RV = Miscibility::NoMixing>
|
||||
class EquilReg {
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* \param[in] rec Equilibration data of current region.
|
||||
* \param[in] density Density calculator of current region.
|
||||
* \param[in] rs Calculator of dissolved gas-oil ratio.
|
||||
* \param[in] rv Calculator of vapourised oil-gas ratio.
|
||||
* \param[in] pu Summary of current active phases.
|
||||
*/
|
||||
EquilReg(const EquilRecord& rec,
|
||||
const DensCalc& density,
|
||||
const RS& rs,
|
||||
const RV& rv,
|
||||
const PhaseUsage& pu)
|
||||
: rec_ (rec)
|
||||
, density_(density)
|
||||
, rs_ (rs)
|
||||
, rv_ (rv)
|
||||
, pu_ (pu)
|
||||
{
|
||||
}
|
||||
|
||||
/**
|
||||
* Type of density calculator.
|
||||
*/
|
||||
typedef DensCalc CalcDensity;
|
||||
|
||||
/**
|
||||
* Type of dissolved gas-oil ratio calculator.
|
||||
*/
|
||||
typedef RS CalcDissolution;
|
||||
|
||||
/**
|
||||
* Type of vapourised oil-gas ratio calculator.
|
||||
*/
|
||||
typedef RV CalcEvaporation;
|
||||
|
||||
/**
|
||||
* Datum depth in current region
|
||||
*/
|
||||
double datum() const { return this->rec_.main.depth; }
|
||||
|
||||
/**
|
||||
* Pressure at datum depth in current region.
|
||||
*/
|
||||
double pressure() const { return this->rec_.main.press; }
|
||||
|
||||
/**
|
||||
* Depth of water-oil contact.
|
||||
*/
|
||||
double zwoc() const { return this->rec_.woc .depth; }
|
||||
|
||||
/**
|
||||
* water-oil capillary pressure at water-oil contact.
|
||||
*
|
||||
* \return P_o - P_w at WOC.
|
||||
*/
|
||||
double pcow_woc() const { return this->rec_.woc .press; }
|
||||
|
||||
/**
|
||||
* Depth of gas-oil contact.
|
||||
*/
|
||||
double zgoc() const { return this->rec_.goc .depth; }
|
||||
|
||||
/**
|
||||
* Gas-oil capillary pressure at gas-oil contact.
|
||||
*
|
||||
* \return P_g - P_o at GOC.
|
||||
*/
|
||||
double pcgo_goc() const { return this->rec_.goc .press; }
|
||||
|
||||
/**
|
||||
* Retrieve phase density calculator of current region.
|
||||
*/
|
||||
const CalcDensity&
|
||||
densityCalculator() const { return this->density_; }
|
||||
|
||||
/**
|
||||
* 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 active fluid phase summary.
|
||||
*/
|
||||
const PhaseUsage&
|
||||
phaseUsage() const { return this->pu_; }
|
||||
|
||||
private:
|
||||
EquilRecord rec_; /**< Equilibration data */
|
||||
DensCalc density_; /**< Density calculator */
|
||||
RS rs_; /**< RS calculator */
|
||||
RV rv_; /**< RV calculator */
|
||||
PhaseUsage pu_; /**< Active phase summary */
|
||||
};
|
||||
|
||||
|
||||
|
||||
/// Functor for inverting capillary pressure function.
|
||||
/// Function represented is
|
||||
/// f(s) = pc(s) - target_pc
|
||||
struct PcEq
|
||||
{
|
||||
PcEq(const BlackoilPropertiesInterface& props,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
: props_(props),
|
||||
phase_(phase),
|
||||
cell_(cell),
|
||||
target_pc_(target_pc)
|
||||
{
|
||||
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
}
|
||||
double operator()(double s) const
|
||||
{
|
||||
s_[phase_] = s;
|
||||
props_.capPress(1, s_, &cell_, pc_, 0);
|
||||
return pc_[phase_] - target_pc_;
|
||||
}
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int phase_;
|
||||
const int cell_;
|
||||
const int target_pc_;
|
||||
mutable double s_[BlackoilPhases::MaxNumPhases];
|
||||
mutable double pc_[BlackoilPhases::MaxNumPhases];
|
||||
};
|
||||
|
||||
|
||||
|
||||
/// Compute saturation of some phase corresponding to a given
|
||||
/// capillary pressure.
|
||||
inline double satFromPc(const BlackoilPropertiesInterface& props,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double target_pc,
|
||||
const bool increasing = false)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double sminarr[BlackoilPhases::MaxNumPhases];
|
||||
double smaxarr[BlackoilPhases::MaxNumPhases];
|
||||
props.satRange(1, &cell, sminarr, smaxarr);
|
||||
const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
|
||||
const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
|
||||
|
||||
// Create the equation f(s) = pc(s) - target_pc
|
||||
const PcEq f(props, 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 = 30;
|
||||
const double tol = 1e-6;
|
||||
int iter_used = -1;
|
||||
typedef RegulaFalsi<ThrowOnError> 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
|
||||
struct PcEqSum
|
||||
{
|
||||
PcEqSum(const BlackoilPropertiesInterface& props,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
: props_(props),
|
||||
phase1_(phase1),
|
||||
phase2_(phase2),
|
||||
cell_(cell),
|
||||
target_pc_(target_pc)
|
||||
{
|
||||
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
}
|
||||
double operator()(double s) const
|
||||
{
|
||||
s_[phase1_] = s;
|
||||
s_[phase2_] = 1.0 - s;
|
||||
props_.capPress(1, s_, &cell_, pc_, 0);
|
||||
return pc_[phase1_] + pc_[phase2_] - target_pc_;
|
||||
}
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int phase1_;
|
||||
const int phase2_;
|
||||
const int cell_;
|
||||
const int target_pc_;
|
||||
mutable double s_[BlackoilPhases::MaxNumPhases];
|
||||
mutable double pc_[BlackoilPhases::MaxNumPhases];
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
/// 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.
|
||||
inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double sminarr[BlackoilPhases::MaxNumPhases];
|
||||
double smaxarr[BlackoilPhases::MaxNumPhases];
|
||||
props.satRange(1, &cell, sminarr, smaxarr);
|
||||
const double smin = sminarr[phase1];
|
||||
const double smax = smaxarr[phase1];
|
||||
|
||||
// Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
|
||||
const PcEqSum f(props, 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<ThrowOnError> ScalarSolver;
|
||||
const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used);
|
||||
return sol;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Compute initial phase pressures by means of equilibration.
|
||||
|
||||
Loading…
Reference in New Issue
Block a user