mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-11-30 13:03:49 -06:00
Remove BlackoilProps from equil initalization code
Use FluidSystem and materialLaw from opm-material directly not via the BlackoilProps in opm-core
This commit is contained in:
Tor Harald Sandve
Andreas Lauser
parent
546759f775
commit
8efc60957a
@ -20,7 +20,6 @@
|
||||
#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>
|
||||
@ -28,6 +27,10 @@
|
||||
|
||||
#include <opm/parser/eclipse/EclipseState/InitConfig/Equil.hpp>
|
||||
|
||||
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
|
||||
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
|
||||
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
|
||||
|
||||
#include <memory>
|
||||
|
||||
|
||||
@ -87,71 +90,21 @@ namespace 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] T Temperature.
|
||||
*
|
||||
* \param[in] z Surface volumes of all phases.
|
||||
*
|
||||
* \return Phase densities at phase point.
|
||||
*/
|
||||
std::vector<double>
|
||||
operator()(const double p,
|
||||
const double T,
|
||||
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, &T, &z[0], &c_[0], &A[0], dAdp);
|
||||
|
||||
std::vector<double> rho(np, 0.0);
|
||||
props_.density(1, &A[0], &c_[0], &rho[0]);
|
||||
|
||||
return rho;
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const std::vector<int> c_;
|
||||
};
|
||||
typedef Opm::FluidSystems::BlackOil<double> FluidSystemSimple;
|
||||
|
||||
// Adjust oil pressure according to gas saturation and cap pressure
|
||||
typedef Opm::SimpleModularFluidState<double,
|
||||
/*numPhases=*/3,
|
||||
/*numComponents=*/3,
|
||||
FluidSystemSimple,
|
||||
/*storePressure=*/false,
|
||||
/*storeTemperature=*/false,
|
||||
/*storeComposition=*/false,
|
||||
/*storeFugacity=*/false,
|
||||
/*storeSaturation=*/true,
|
||||
/*storeDensity=*/false,
|
||||
/*storeViscosity=*/false,
|
||||
/*storeEnthalpy=*/false> SatOnlyFluidState;
|
||||
|
||||
/**
|
||||
* Types and routines relating to phase mixing in
|
||||
@ -224,6 +177,7 @@ namespace Opm
|
||||
* tabulated as a function of depth policy. Data
|
||||
* typically taken from keyword 'RSVD'.
|
||||
*/
|
||||
template <class FluidSystem>
|
||||
class RsVD : public RsFunction {
|
||||
public:
|
||||
/**
|
||||
@ -234,19 +188,14 @@ namespace Opm
|
||||
* \param[in] depth Depth nodes.
|
||||
* \param[in] rs Dissolved gas-oil ratio at @c depth.
|
||||
*/
|
||||
RsVD(const BlackoilPropertiesInterface& props,
|
||||
const int cell,
|
||||
RsVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& rs)
|
||||
: props_(props)
|
||||
, cell_(cell)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, depth_(depth)
|
||||
, rs_(rs)
|
||||
{
|
||||
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;
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
@ -278,23 +227,13 @@ namespace Opm
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int cell_;
|
||||
const int pvtRegionIdx_;
|
||||
std::vector<double> depth_; /**< Depth nodes */
|
||||
std::vector<double> rs_; /**< Dissolved gas-oil ratio */
|
||||
double z_[BlackoilPhases::MaxNumPhases];
|
||||
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
|
||||
|
||||
double satRs(const double press, const double temp) const
|
||||
{
|
||||
props_.matrix(1, &press, &temp, 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];
|
||||
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
};
|
||||
|
||||
@ -304,6 +243,7 @@ namespace Opm
|
||||
* tabulated as a function of depth policy. Data
|
||||
* typically taken from keyword 'RVVD'.
|
||||
*/
|
||||
template <class FluidSystem>
|
||||
class RvVD : public RsFunction {
|
||||
public:
|
||||
/**
|
||||
@ -314,19 +254,14 @@ namespace Opm
|
||||
* \param[in] depth Depth nodes.
|
||||
* \param[in] rv Dissolved gas-oil ratio at @c depth.
|
||||
*/
|
||||
RvVD(const BlackoilPropertiesInterface& props,
|
||||
const int cell,
|
||||
RvVD(const int pvtRegionIdx,
|
||||
const std::vector<double>& depth,
|
||||
const std::vector<double>& rv)
|
||||
: props_(props)
|
||||
, cell_(cell)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
, depth_(depth)
|
||||
, rv_(rv)
|
||||
{
|
||||
auto pu = props_.phaseUsage();
|
||||
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
|
||||
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
@ -358,23 +293,13 @@ namespace Opm
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int cell_;
|
||||
const int pvtRegionIdx_;
|
||||
std::vector<double> depth_; /**< Depth nodes */
|
||||
std::vector<double> rv_; /**< Vaporized oil-gas ratio */
|
||||
double z_[BlackoilPhases::MaxNumPhases];
|
||||
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
|
||||
|
||||
double satRv(const double press, const double temp) const
|
||||
{
|
||||
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
|
||||
// Rv/Bg is in the oil row and gas column of A_.
|
||||
// 1/Bg is in the gas 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*gpos + opos] / A_[np*gpos + gpos];
|
||||
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
};
|
||||
|
||||
@ -393,6 +318,7 @@ namespace Opm
|
||||
* This should yield Rs-values that are constant below the
|
||||
* contact, and decreasing above the contact.
|
||||
*/
|
||||
template <class FluidSystem>
|
||||
class RsSatAtContact : public RsFunction {
|
||||
public:
|
||||
/**
|
||||
@ -403,13 +329,9 @@ namespace Opm
|
||||
* \param[in] p_contact oil pressure at the contact
|
||||
* \param[in] T_contact temperature at the contact
|
||||
*/
|
||||
RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
|
||||
: props_(props), cell_(cell)
|
||||
RsSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact)
|
||||
: pvtRegionIdx_(pvtRegionIdx)
|
||||
{
|
||||
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, T_contact);
|
||||
}
|
||||
|
||||
@ -442,22 +364,12 @@ namespace Opm
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int cell_;
|
||||
double z_[BlackoilPhases::MaxNumPhases];
|
||||
const int pvtRegionIdx_;
|
||||
double rs_sat_contact_;
|
||||
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
|
||||
|
||||
double satRs(const double press, const double temp) const
|
||||
{
|
||||
props_.matrix(1, &press, &temp, 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];
|
||||
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
|
||||
}
|
||||
};
|
||||
|
||||
@ -476,6 +388,7 @@ namespace Opm
|
||||
* This should yield Rv-values that are constant below the
|
||||
* contact, and decreasing above the contact.
|
||||
*/
|
||||
template <class FluidSystem>
|
||||
class RvSatAtContact : public RsFunction {
|
||||
public:
|
||||
/**
|
||||
@ -486,13 +399,9 @@ namespace Opm
|
||||
* \param[in] p_contact oil pressure at the contact
|
||||
* \param[in] T_contact temperature at the contact
|
||||
*/
|
||||
RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
|
||||
: props_(props), cell_(cell)
|
||||
RvSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact)
|
||||
:pvtRegionIdx_(pvtRegionIdx)
|
||||
{
|
||||
auto pu = props_.phaseUsage();
|
||||
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
|
||||
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
|
||||
rv_sat_contact_ = satRv(p_contact, T_contact);
|
||||
}
|
||||
|
||||
@ -525,22 +434,12 @@ namespace Opm
|
||||
}
|
||||
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
const int cell_;
|
||||
double z_[BlackoilPhases::MaxNumPhases];
|
||||
const int pvtRegionIdx_;
|
||||
double rv_sat_contact_;
|
||||
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
|
||||
|
||||
double satRv(const double press, const double temp) const
|
||||
{
|
||||
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
|
||||
// Rv/Bg is in the oil row and gas column of A_.
|
||||
// 1/Bg is in the gas 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*gpos + opos] / A_[np*gpos + gpos];
|
||||
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
|
||||
}
|
||||
};
|
||||
|
||||
@ -565,7 +464,6 @@ namespace Opm
|
||||
* that calculates the phase densities of all phases in @c
|
||||
* svol at fluid pressure @c press.
|
||||
*/
|
||||
template <class DensCalc>
|
||||
class EquilReg {
|
||||
public:
|
||||
/**
|
||||
@ -578,23 +476,19 @@ namespace Opm
|
||||
* \param[in] pu Summary of current active phases.
|
||||
*/
|
||||
EquilReg(const EquilRecord& rec,
|
||||
const DensCalc& density,
|
||||
std::shared_ptr<Miscibility::RsFunction> rs,
|
||||
std::shared_ptr<Miscibility::RsFunction> rv,
|
||||
const int pvtIdx,
|
||||
const PhaseUsage& pu)
|
||||
: rec_ (rec)
|
||||
, density_(density)
|
||||
, rs_ (rs)
|
||||
, rv_ (rv)
|
||||
, pvtIdx_ (pvtIdx)
|
||||
, pu_ (pu)
|
||||
|
||||
{
|
||||
}
|
||||
|
||||
/**
|
||||
* Type of density calculator.
|
||||
*/
|
||||
typedef DensCalc CalcDensity;
|
||||
|
||||
/**
|
||||
* Type of dissolved gas-oil ratio calculator.
|
||||
*/
|
||||
@ -639,11 +533,6 @@ namespace Opm
|
||||
*/
|
||||
double pcgo_goc() const { return this->rec_.gasOilContactCapillaryPressure(); }
|
||||
|
||||
/**
|
||||
* Retrieve phase density calculator of current region.
|
||||
*/
|
||||
const CalcDensity&
|
||||
densityCalculator() const { return this->density_; }
|
||||
|
||||
/**
|
||||
* Retrieve dissolved gas-oil ratio calculator of current
|
||||
@ -659,6 +548,12 @@ namespace Opm
|
||||
const CalcEvaporation&
|
||||
evaporationCalculator() const { return *this->rv_; }
|
||||
|
||||
/**
|
||||
* Retrieve active fluid phase summary.
|
||||
*/
|
||||
const int
|
||||
pvtIdx() const { return this->pvtIdx_; }
|
||||
|
||||
/**
|
||||
* Retrieve active fluid phase summary.
|
||||
*/
|
||||
@ -667,9 +562,9 @@ namespace Opm
|
||||
|
||||
private:
|
||||
EquilRecord rec_; /**< Equilibration data */
|
||||
DensCalc density_; /**< Density calculator */
|
||||
std::shared_ptr<Miscibility::RsFunction> rs_; /**< RS calculator */
|
||||
std::shared_ptr<Miscibility::RsFunction> rv_; /**< RV calculator */
|
||||
const int pvtIdx_;
|
||||
PhaseUsage pu_; /**< Active phase summary */
|
||||
};
|
||||
|
||||
@ -678,54 +573,117 @@ namespace Opm
|
||||
/// Functor for inverting capillary pressure function.
|
||||
/// Function represented is
|
||||
/// f(s) = pc(s) - target_pc
|
||||
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
|
||||
struct PcEq
|
||||
{
|
||||
PcEq(const BlackoilPropertiesInterface& props,
|
||||
PcEq(std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const int phase,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
: props_(props),
|
||||
: materialLawManager_(materialLawManager),
|
||||
phase_(phase),
|
||||
cell_(cell),
|
||||
target_pc_(target_pc)
|
||||
{
|
||||
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
fluidState_.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
|
||||
fluidState_.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
|
||||
fluidState_.setSaturation(FluidSystem::gasPhaseIdx, 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_;
|
||||
const auto& matParams = materialLawManager_->materialLawParams(cell_);
|
||||
fluidState_.setSaturation(phase_, s);
|
||||
MaterialLaw::capillaryPressures(pc_, matParams, fluidState_);
|
||||
double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
|
||||
double pc = pc_[FluidSystem::oilPhaseIdx] + sign * pc_[phase_];
|
||||
return pc - target_pc_;
|
||||
}
|
||||
private:
|
||||
const BlackoilPropertiesInterface& props_;
|
||||
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager_;
|
||||
const int phase_;
|
||||
const int cell_;
|
||||
const double target_pc_;
|
||||
mutable double s_[BlackoilPhases::MaxNumPhases];
|
||||
mutable SatOnlyFluidState fluidState_;
|
||||
mutable double pc_[BlackoilPhases::MaxNumPhases];
|
||||
};
|
||||
|
||||
template <class FluidSystem, class MaterialLawManager>
|
||||
double minSaturations(std::shared_ptr<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 <class FluidSystem, class MaterialLawManager>
|
||||
double maxSaturations(std::shared_ptr<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.
|
||||
inline double satFromPc(const BlackoilPropertiesInterface& props,
|
||||
template <class FluidSystem, class MaterialLaw, class MaterialLawManager >
|
||||
inline double satFromPc(std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
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];
|
||||
// Find minimum and maximum saturations.
|
||||
const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
|
||||
// Create the equation f(s) = pc(s) - target_pc
|
||||
const PcEq f(props, phase, cell, target_pc);
|
||||
const PcEq<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase, cell, target_pc);
|
||||
const double f0 = f(s0);
|
||||
const double f1 = f(s1);
|
||||
|
||||
@ -747,36 +705,46 @@ namespace Opm
|
||||
/// Functor for inverting a sum of capillary pressure functions.
|
||||
/// Function represented is
|
||||
/// f(s) = pc1(s) + pc2(1 - s) - target_pc
|
||||
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
|
||||
struct PcEqSum
|
||||
{
|
||||
PcEqSum(const BlackoilPropertiesInterface& props,
|
||||
PcEqSum(std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const int phase1,
|
||||
const int phase2,
|
||||
const int cell,
|
||||
const double target_pc)
|
||||
: props_(props),
|
||||
: materialLawManager_(materialLawManager),
|
||||
phase1_(phase1),
|
||||
phase2_(phase2),
|
||||
cell_(cell),
|
||||
target_pc_(target_pc)
|
||||
{
|
||||
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
|
||||
fluidState_.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
|
||||
fluidState_.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
|
||||
fluidState_.setSaturation(FluidSystem::gasPhaseIdx, 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_;
|
||||
const auto& matParams = materialLawManager_->materialLawParams(cell_);
|
||||
fluidState_.setSaturation(phase1_, s);
|
||||
fluidState_.setSaturation(phase2_, 1.0 - s);
|
||||
|
||||
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 BlackoilPropertiesInterface& props_;
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager_;
|
||||
const int phase1_;
|
||||
const int phase2_;
|
||||
const int cell_;
|
||||
const double target_pc_;
|
||||
mutable double s_[BlackoilPhases::MaxNumPhases];
|
||||
mutable SatOnlyFluidState fluidState_;
|
||||
mutable double pc_[BlackoilPhases::MaxNumPhases];
|
||||
};
|
||||
|
||||
@ -786,21 +754,20 @@ namespace Opm
|
||||
/// 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,
|
||||
|
||||
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
|
||||
inline double satFromSumOfPcs(std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
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];
|
||||
const double smin = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
|
||||
const double smax = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
|
||||
|
||||
// Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
|
||||
const PcEqSum f(props, phase1, phase2, cell, target_pc);
|
||||
const PcEqSum<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase1, phase2, cell, target_pc);
|
||||
const double f0 = f(smin);
|
||||
const double f1 = f(smax);
|
||||
if (f0 <= 0.0) {
|
||||
@ -818,19 +785,16 @@ namespace Opm
|
||||
}
|
||||
|
||||
/// Compute saturation from depth. Used for constant capillary pressure function
|
||||
inline double satFromDepth(const BlackoilPropertiesInterface& props,
|
||||
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
|
||||
inline double satFromDepth(std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const double cellDepth,
|
||||
const double contactDepth,
|
||||
const int phase,
|
||||
const int cell,
|
||||
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];
|
||||
const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
|
||||
|
||||
if (cellDepth < contactDepth){
|
||||
return s0;
|
||||
@ -841,22 +805,20 @@ namespace Opm
|
||||
}
|
||||
|
||||
/// Return true if capillary pressure function is constant
|
||||
inline bool isConstPc(const BlackoilPropertiesInterface& props,
|
||||
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
|
||||
inline bool isConstPc(std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const int phase,
|
||||
const int cell)
|
||||
{
|
||||
// Find minimum and maximum saturations.
|
||||
double sminarr[BlackoilPhases::MaxNumPhases];
|
||||
double smaxarr[BlackoilPhases::MaxNumPhases];
|
||||
props.satRange(1, &cell, sminarr, smaxarr);
|
||||
|
||||
// Create the equation f(s) = pc(s);
|
||||
const PcEq f(props, phase, cell, 0);
|
||||
const double f0 = f(sminarr[phase]);
|
||||
const double f1 = f(smaxarr[phase]);
|
||||
const PcEq<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase, cell, 0);
|
||||
const double f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
|
||||
const double f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
|
||||
return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
|
||||
}
|
||||
|
||||
|
||||
|
||||
} // namespace Equil
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
@ -25,9 +25,11 @@
|
||||
#include <opm/core/grid/GridHelpers.hpp>
|
||||
#include <opm/core/simulator/EquilibrationHelpers.hpp>
|
||||
#include <opm/core/simulator/BlackoilState.hpp>
|
||||
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/props/phaseUsageFromDeck.hpp>
|
||||
#include <opm/core/utility/RegionMapping.hpp>
|
||||
#include <opm/core/utility/extractPvtTableIndex.hpp>
|
||||
|
||||
#include <opm/parser/eclipse/Units/Units.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
|
||||
#include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
|
||||
@ -39,6 +41,11 @@
|
||||
#include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
|
||||
#include <opm/common/OpmLog/OpmLog.hpp>
|
||||
|
||||
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
|
||||
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
|
||||
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
|
||||
|
||||
|
||||
#include <array>
|
||||
#include <cassert>
|
||||
#include <utility>
|
||||
@ -54,32 +61,6 @@ struct UnstructuredGrid;
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
/**
|
||||
* Compute initial state by an equilibration procedure.
|
||||
*
|
||||
* The following state fields are modified:
|
||||
* pressure(),
|
||||
* saturation(),
|
||||
* surfacevol(),
|
||||
* gasoilratio(),
|
||||
* rv().
|
||||
*
|
||||
* \param[in] grid Grid.
|
||||
* \param[in] props Property object, pvt and capillary properties are used.
|
||||
* \param[in] deck Simulation deck, used to obtain EQUIL and related data.
|
||||
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
|
||||
* \param[in] applySwatInit Make it possible to not apply SWATINIT even if it
|
||||
* is present in the deck
|
||||
*/
|
||||
template<class Grid>
|
||||
void initStateEquil(const Grid& grid,
|
||||
const BlackoilPropertiesInterface& props,
|
||||
const Opm::Deck& deck,
|
||||
const Opm::EclipseState& eclipseState,
|
||||
const double gravity,
|
||||
BlackoilState& state,
|
||||
bool applySwatInit = true);
|
||||
|
||||
|
||||
/**
|
||||
* Types and routines that collectively implement a basic
|
||||
@ -90,6 +71,7 @@ namespace Opm
|
||||
*/
|
||||
namespace EQUIL {
|
||||
|
||||
|
||||
/**
|
||||
* Compute initial phase pressures by means of equilibration.
|
||||
*
|
||||
@ -128,7 +110,7 @@ namespace Opm
|
||||
* of pressure values in each cell in the current
|
||||
* equilibration region.
|
||||
*/
|
||||
template <class Grid, class Region, class CellRange>
|
||||
template <class Grid, class Region, class CellRange, class FluidSystem>
|
||||
std::vector< std::vector<double> >
|
||||
phasePressures(const Grid& G,
|
||||
const Region& reg,
|
||||
@ -161,12 +143,12 @@ namespace Opm
|
||||
* \return Phase saturations, one vector for each phase, each containing
|
||||
* one saturation value per cell in the region.
|
||||
*/
|
||||
template <class Grid, class Region, class CellRange>
|
||||
template <class FluidSystem, class Grid, class Region, class CellRange, class MaterialLawManager>
|
||||
std::vector< std::vector<double> >
|
||||
phaseSaturations(const Grid& grid,
|
||||
const Region& reg,
|
||||
const CellRange& cells,
|
||||
BlackoilPropertiesFromDeck& props,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const std::vector<double> swat_init,
|
||||
std::vector< std::vector<double> >& phase_pressures);
|
||||
|
||||
@ -244,29 +226,39 @@ namespace Opm
|
||||
|
||||
class InitialStateComputer {
|
||||
public:
|
||||
template<class Grid>
|
||||
InitialStateComputer(BlackoilPropertiesInterface& props,
|
||||
template<class MaterialLawManager, class Grid>
|
||||
InitialStateComputer(std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const PhaseUsage& phaseUsage,
|
||||
const Opm::Deck& deck,
|
||||
const Opm::EclipseState& eclipseState,
|
||||
const Grid& G ,
|
||||
const double grav = unit::gravity,
|
||||
const std::vector<double>& swat_init = {}
|
||||
)
|
||||
: pp_(props.numPhases(),
|
||||
: pp_(phaseUsage.num_phases,
|
||||
std::vector<double>(UgGridHelpers::numCells(G))),
|
||||
sat_(props.numPhases(),
|
||||
sat_(phaseUsage.num_phases,
|
||||
std::vector<double>(UgGridHelpers::numCells(G))),
|
||||
rs_(UgGridHelpers::numCells(G)),
|
||||
rv_(UgGridHelpers::numCells(G)),
|
||||
swat_init_(swat_init)
|
||||
swat_init_(swat_init),
|
||||
phaseUsage_(phaseUsage)
|
||||
|
||||
{
|
||||
|
||||
typedef FluidSystems::BlackOil<double> FluidSystem;
|
||||
|
||||
// Initialize the fluid system
|
||||
FluidSystem::initFromDeck(deck, eclipseState);
|
||||
|
||||
// Get the equilibration records.
|
||||
const std::vector<EquilRecord> rec = getEquil(eclipseState);
|
||||
const auto& tables = eclipseState.getTableManager();
|
||||
// Create (inverse) region mapping.
|
||||
const RegionMapping<> eqlmap(equilnum(deck, eclipseState, G));
|
||||
|
||||
setRegionPvtIdx(G, eclipseState, eqlmap);
|
||||
|
||||
// Create Rs functions.
|
||||
rs_func_.reserve(rec.size());
|
||||
if (deck.hasKeyword("DISGAS")) {
|
||||
@ -274,10 +266,11 @@ namespace Opm
|
||||
for (size_t i = 0; i < rec.size(); ++i) {
|
||||
if (eqlmap.cells(i).empty())
|
||||
{
|
||||
rs_func_.push_back(std::shared_ptr<Miscibility::RsVD>());
|
||||
rs_func_.push_back(std::shared_ptr<Miscibility::RsVD<FluidSystem>>());
|
||||
continue;
|
||||
}
|
||||
const int cell = *(eqlmap.cells(i).begin());
|
||||
const int pvtIdx = regionPvtIdx_[i];
|
||||
|
||||
if (!rec[i].liveOilInitConstantRs()) {
|
||||
if (rsvdTables.size() <= 0 ) {
|
||||
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table not available.");
|
||||
@ -285,8 +278,7 @@ namespace Opm
|
||||
const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i);
|
||||
std::vector<double> depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy();
|
||||
std::vector<double> rsColumn = rsvdTable.getColumn("RS").vectorCopy();
|
||||
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props,
|
||||
cell,
|
||||
rs_func_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
|
||||
depthColumn , rsColumn));
|
||||
} else {
|
||||
if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
|
||||
@ -297,7 +289,7 @@ namespace Opm
|
||||
}
|
||||
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<Miscibility::RsSatAtContact>(props, cell, p_contact, T_contact));
|
||||
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact<FluidSystem>>(pvtIdx, p_contact, T_contact));
|
||||
}
|
||||
}
|
||||
} else {
|
||||
@ -312,10 +304,10 @@ namespace Opm
|
||||
for (size_t i = 0; i < rec.size(); ++i) {
|
||||
if (eqlmap.cells(i).empty())
|
||||
{
|
||||
rv_func_.push_back(std::shared_ptr<Miscibility::RvVD>());
|
||||
rv_func_.push_back(std::shared_ptr<Miscibility::RvVD<FluidSystem>>());
|
||||
continue;
|
||||
}
|
||||
const int cell = *(eqlmap.cells(i).begin());
|
||||
const int pvtIdx = regionPvtIdx_[i];
|
||||
if (!rec[i].wetGasInitConstantRv()) {
|
||||
if (rvvdTables.size() <= 0) {
|
||||
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table not available.");
|
||||
@ -324,9 +316,7 @@ namespace Opm
|
||||
const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i);
|
||||
std::vector<double> depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy();
|
||||
std::vector<double> rvColumn = rvvdTable.getColumn("RV").vectorCopy();
|
||||
|
||||
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props,
|
||||
cell,
|
||||
rv_func_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
|
||||
depthColumn , rvColumn));
|
||||
|
||||
} else {
|
||||
@ -338,7 +328,7 @@ namespace Opm
|
||||
}
|
||||
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<Miscibility::RvSatAtContact>(props, cell, p_contact, T_contact));
|
||||
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact<FluidSystem>>(pvtIdx ,p_contact, T_contact));
|
||||
}
|
||||
}
|
||||
} else {
|
||||
@ -348,7 +338,7 @@ namespace Opm
|
||||
}
|
||||
|
||||
// Compute pressures, saturations, rs and rv factors.
|
||||
calcPressSatRsRv(eqlmap, rec, props, G, grav);
|
||||
calcPressSatRsRv<FluidSystem>(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.
|
||||
@ -363,23 +353,39 @@ namespace Opm
|
||||
const Vec& rv() const { return rv_; }
|
||||
|
||||
private:
|
||||
typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
|
||||
typedef EquilReg<RhoCalc> EqReg;
|
||||
|
||||
typedef EquilReg EqReg;
|
||||
|
||||
std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
|
||||
std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
|
||||
|
||||
std::vector<int> regionPvtIdx_;
|
||||
|
||||
PVec pp_;
|
||||
PVec sat_;
|
||||
Vec rs_;
|
||||
Vec rv_;
|
||||
Vec swat_init_;
|
||||
PhaseUsage phaseUsage_;
|
||||
|
||||
template <class RMap, class Grid>
|
||||
template<class Grid, class RMap>
|
||||
void setRegionPvtIdx(const Grid& G, const Opm::EclipseState& eclipseState, const RMap& reg) {
|
||||
|
||||
regionPvtIdx_.reserve(reg.activeRegions().size());
|
||||
std::vector<int> 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 <class FluidSystem, class RMap, class MaterialLawManager, class Grid>
|
||||
void
|
||||
calcPressSatRsRv(const RMap& reg ,
|
||||
const std::vector< EquilRecord >& rec ,
|
||||
Opm::BlackoilPropertiesInterface& props,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const Grid& G ,
|
||||
const double grav)
|
||||
{
|
||||
@ -391,27 +397,24 @@ namespace Opm
|
||||
+ " has no active cells");
|
||||
continue;
|
||||
}
|
||||
const int repcell = *cells.begin();
|
||||
|
||||
const RhoCalc calc(props, repcell);
|
||||
const EqReg eqreg(rec[r], calc,
|
||||
rs_func_[r], rv_func_[r],
|
||||
props.phaseUsage());
|
||||
const EqReg eqreg(rec[r],
|
||||
rs_func_[r], rv_func_[r], regionPvtIdx_[r],
|
||||
phaseUsage_);
|
||||
|
||||
PVec pressures = phasePressures(G, eqreg, cells, grav);
|
||||
PVec pressures = phasePressures<FluidSystem>(G, eqreg, cells, grav);
|
||||
const std::vector<double>& temp = temperature(G, eqreg, cells);
|
||||
const PVec sat = phaseSaturations<FluidSystem>(G, eqreg, cells, materialLawManager, swat_init_, pressures);
|
||||
|
||||
const PVec sat = phaseSaturations(G, eqreg, cells, props, swat_init_, pressures);
|
||||
|
||||
const int np = props.numPhases();
|
||||
const int np = phaseUsage_.num_phases;
|
||||
for (int p = 0; p < np; ++p) {
|
||||
copyFromRegion(pressures[p], cells, pp_[p]);
|
||||
copyFromRegion(sat[p], cells, sat_[p]);
|
||||
}
|
||||
if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
|
||||
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
|
||||
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
|
||||
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
|
||||
if (phaseUsage_.phase_used[BlackoilPhases::Liquid]
|
||||
&& phaseUsage_.phase_used[BlackoilPhases::Vapour]) {
|
||||
const int oilpos = phaseUsage_.phase_pos[BlackoilPhases::Liquid];
|
||||
const int gaspos = phaseUsage_.phase_pos[BlackoilPhases::Vapour];
|
||||
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_);
|
||||
|
||||
@ -24,8 +24,8 @@
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/grid/GridHelpers.hpp>
|
||||
|
||||
#include <opm/core/props/BlackoilPhases.hpp>
|
||||
#include <opm/core/simulator/initState.hpp>
|
||||
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
@ -35,6 +35,9 @@
|
||||
namespace Opm
|
||||
{
|
||||
namespace Details {
|
||||
|
||||
|
||||
|
||||
template <class RHS>
|
||||
class RK4IVP {
|
||||
public:
|
||||
@ -105,22 +108,17 @@ namespace Opm
|
||||
stepsize() const { return (span_[1] - span_[0]) / N_; }
|
||||
};
|
||||
|
||||
namespace PhasePressODE {
|
||||
template <class Density>
|
||||
namespace PhasePressODE {
|
||||
template <class FluidSystem>
|
||||
class Water {
|
||||
public:
|
||||
Water(const double temp,
|
||||
const Density& rho,
|
||||
const int np,
|
||||
const int ix,
|
||||
const int pvtRegionIdx,
|
||||
const double norm_grav)
|
||||
: temp_(temp)
|
||||
, rho_(rho)
|
||||
, svol_(np, 0)
|
||||
, ix_(ix)
|
||||
, pvtRegionIdx_(pvtRegionIdx)
|
||||
, g_(norm_grav)
|
||||
{
|
||||
svol_[ix_] = 1.0;
|
||||
}
|
||||
|
||||
double
|
||||
@ -132,39 +130,30 @@ namespace Opm
|
||||
|
||||
private:
|
||||
const double temp_;
|
||||
const Density& rho_;
|
||||
std::vector<double> svol_;
|
||||
const int ix_;
|
||||
const int pvtRegionIdx_;
|
||||
const double g_;
|
||||
|
||||
double
|
||||
density(const double press) const
|
||||
{
|
||||
const std::vector<double>& rho = rho_(press, temp_, svol_);
|
||||
|
||||
return rho[ix_];
|
||||
double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
|
||||
rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
|
||||
return rho;
|
||||
}
|
||||
};
|
||||
|
||||
template <class Density, class RS>
|
||||
template <class FluidSystem, class RS>
|
||||
class Oil {
|
||||
public:
|
||||
Oil(const double temp,
|
||||
const Density& rho,
|
||||
const RS& rs,
|
||||
const int np,
|
||||
const int oix,
|
||||
const int gix,
|
||||
const int pvtRegionIdx,
|
||||
const double norm_grav)
|
||||
: temp_(temp)
|
||||
, rho_(rho)
|
||||
, rs_(rs)
|
||||
, svol_(np, 0)
|
||||
, oix_(oix)
|
||||
, gix_(gix)
|
||||
, pvtRegionIdx_(pvtRegionIdx)
|
||||
, g_(norm_grav)
|
||||
{
|
||||
svol_[oix_] = 1.0;
|
||||
}
|
||||
|
||||
double
|
||||
@ -176,45 +165,42 @@ namespace Opm
|
||||
|
||||
private:
|
||||
const double temp_;
|
||||
const Density& rho_;
|
||||
const RS& rs_;
|
||||
mutable std::vector<double> svol_;
|
||||
const int oix_;
|
||||
const int gix_;
|
||||
const int pvtRegionIdx_;
|
||||
const double g_;
|
||||
|
||||
double
|
||||
density(const double depth,
|
||||
const double press) const
|
||||
{
|
||||
if (gix_ >= 0) {
|
||||
svol_[gix_] = rs_(depth, press, temp_);
|
||||
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_);
|
||||
}
|
||||
|
||||
const std::vector<double>& rho = rho_(press, temp_, svol_);
|
||||
return rho[oix_];
|
||||
return rho;
|
||||
}
|
||||
};
|
||||
|
||||
template <class Density, class RV>
|
||||
template <class FluidSystem, class RV>
|
||||
class Gas {
|
||||
public:
|
||||
Gas(const double temp,
|
||||
const Density& rho,
|
||||
const RV& rv,
|
||||
const int np,
|
||||
const int gix,
|
||||
const int oix,
|
||||
const int pvtRegionIdx,
|
||||
const double norm_grav)
|
||||
: temp_(temp)
|
||||
, rho_(rho)
|
||||
, rv_(rv)
|
||||
, svol_(np, 0)
|
||||
, gix_(gix)
|
||||
, oix_(oix)
|
||||
, pvtRegionIdx_(pvtRegionIdx)
|
||||
, g_(norm_grav)
|
||||
{
|
||||
svol_[gix_] = 1.0;
|
||||
}
|
||||
|
||||
double
|
||||
@ -226,23 +212,27 @@ namespace Opm
|
||||
|
||||
private:
|
||||
const double temp_;
|
||||
const Density& rho_;
|
||||
const RV& rv_;
|
||||
mutable std::vector<double> svol_;
|
||||
const int gix_;
|
||||
const int oix_;
|
||||
const int pvtRegionIdx_;
|
||||
const double g_;
|
||||
|
||||
double
|
||||
density(const double depth,
|
||||
const double press) const
|
||||
{
|
||||
if (oix_ >= 0) {
|
||||
svol_[oix_] = rv_(depth, press, temp_);
|
||||
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_);
|
||||
}
|
||||
|
||||
const std::vector<double>& rho = rho_(press, temp_, svol_);
|
||||
return rho[gix_];
|
||||
return rho;
|
||||
}
|
||||
};
|
||||
} // namespace PhasePressODE
|
||||
@ -328,9 +318,10 @@ namespace Opm
|
||||
}
|
||||
}
|
||||
|
||||
template <class Grid,
|
||||
class Region,
|
||||
class CellRange>
|
||||
template < class FluidSystem,
|
||||
class Grid,
|
||||
class Region,
|
||||
class CellRange>
|
||||
void
|
||||
water(const Grid& G ,
|
||||
const Region& reg ,
|
||||
@ -341,13 +332,10 @@ namespace Opm
|
||||
std::vector<double>& press )
|
||||
{
|
||||
using PhasePressODE::Water;
|
||||
typedef Water<typename Region::CalcDensity> ODE;
|
||||
typedef Water<FluidSystem> ODE;
|
||||
|
||||
const PhaseUsage& pu = reg.phaseUsage();
|
||||
|
||||
const int wix = PhaseIndex::water(pu);
|
||||
const double T = 273.15 + 20; // standard temperature for now
|
||||
ODE drho(T, reg.densityCalculator(), pu.num_phases, wix, grav);
|
||||
ODE drho(T, reg.pvtIdx() , grav);
|
||||
|
||||
double z0;
|
||||
double p0;
|
||||
@ -379,9 +367,11 @@ namespace Opm
|
||||
}
|
||||
}
|
||||
|
||||
template <class Grid,
|
||||
template <class FluidSystem,
|
||||
class Grid,
|
||||
class Region,
|
||||
class CellRange>
|
||||
|
||||
void
|
||||
oil(const Grid& G ,
|
||||
const Region& reg ,
|
||||
@ -393,18 +383,12 @@ namespace Opm
|
||||
double& po_goc)
|
||||
{
|
||||
using PhasePressODE::Oil;
|
||||
typedef Oil<typename Region::CalcDensity,
|
||||
typename Region::CalcDissolution> ODE;
|
||||
typedef Oil<FluidSystem, typename Region::CalcDissolution> ODE;
|
||||
|
||||
const PhaseUsage& pu = reg.phaseUsage();
|
||||
|
||||
const int oix = PhaseIndex::oil(pu);
|
||||
const int gix = PhaseIndex::gas(pu);
|
||||
const double T = 273.15 + 20; // standard temperature for now
|
||||
ODE drho(T,
|
||||
reg.densityCalculator(),
|
||||
reg.dissolutionCalculator(),
|
||||
pu.num_phases, oix, gix, grav);
|
||||
|
||||
ODE drho(T, reg.dissolutionCalculator(),
|
||||
reg.pvtIdx(), grav);
|
||||
|
||||
double z0;
|
||||
double p0;
|
||||
@ -444,7 +428,8 @@ namespace Opm
|
||||
else { po_goc = p0; } // GOC *at* datum
|
||||
}
|
||||
|
||||
template <class Grid,
|
||||
template <class FluidSystem,
|
||||
class Grid,
|
||||
class Region,
|
||||
class CellRange>
|
||||
void
|
||||
@ -457,19 +442,11 @@ namespace Opm
|
||||
std::vector<double>& press )
|
||||
{
|
||||
using PhasePressODE::Gas;
|
||||
typedef Gas<typename Region::CalcDensity,
|
||||
typename Region::CalcEvaporation> ODE;
|
||||
|
||||
const PhaseUsage& pu = reg.phaseUsage();
|
||||
|
||||
const int gix = PhaseIndex::gas(pu);
|
||||
const int oix = PhaseIndex::oil(pu);
|
||||
typedef Gas<FluidSystem, typename Region::CalcEvaporation> ODE;
|
||||
|
||||
const double T = 273.15 + 20; // standard temperature for now
|
||||
ODE drho(T,
|
||||
reg.densityCalculator(),
|
||||
reg.evaporationCalculator(),
|
||||
pu.num_phases, gix, oix, grav);
|
||||
ODE drho(T, reg.evaporationCalculator(),
|
||||
reg.pvtIdx(), grav);
|
||||
|
||||
double z0;
|
||||
double p0;
|
||||
@ -502,7 +479,8 @@ namespace Opm
|
||||
}
|
||||
} // namespace PhasePressure
|
||||
|
||||
template <class Grid,
|
||||
template <class FluidSystem,
|
||||
class Grid,
|
||||
class Region,
|
||||
class CellRange>
|
||||
void
|
||||
@ -521,19 +499,19 @@ namespace Opm
|
||||
|
||||
if (PhaseUsed::water(pu)) {
|
||||
const int wix = PhaseIndex::water(pu);
|
||||
PhasePressure::water(G, reg, span, grav, po_woc,
|
||||
PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
|
||||
cells, press[ wix ]);
|
||||
}
|
||||
|
||||
if (PhaseUsed::oil(pu)) {
|
||||
const int oix = PhaseIndex::oil(pu);
|
||||
PhasePressure::oil(G, reg, span, grav, cells,
|
||||
PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
|
||||
press[ oix ], po_woc, po_goc);
|
||||
}
|
||||
|
||||
if (PhaseUsed::gas(pu)) {
|
||||
const int gix = PhaseIndex::gas(pu);
|
||||
PhasePressure::gas(G, reg, span, grav, po_goc,
|
||||
PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
|
||||
cells, press[ gix ]);
|
||||
}
|
||||
} else if (reg.datum() < reg.zgoc()) { // Datum in gas zone
|
||||
@ -542,19 +520,19 @@ namespace Opm
|
||||
|
||||
if (PhaseUsed::gas(pu)) {
|
||||
const int gix = PhaseIndex::gas(pu);
|
||||
PhasePressure::gas(G, reg, span, grav, po_goc,
|
||||
PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
|
||||
cells, press[ gix ]);
|
||||
}
|
||||
|
||||
if (PhaseUsed::oil(pu)) {
|
||||
const int oix = PhaseIndex::oil(pu);
|
||||
PhasePressure::oil(G, reg, span, grav, cells,
|
||||
PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
|
||||
press[ oix ], po_woc, po_goc);
|
||||
}
|
||||
|
||||
if (PhaseUsed::water(pu)) {
|
||||
const int wix = PhaseIndex::water(pu);
|
||||
PhasePressure::water(G, reg, span, grav, po_woc,
|
||||
PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
|
||||
cells, press[ wix ]);
|
||||
}
|
||||
} else { // Datum in oil zone
|
||||
@ -563,19 +541,19 @@ namespace Opm
|
||||
|
||||
if (PhaseUsed::oil(pu)) {
|
||||
const int oix = PhaseIndex::oil(pu);
|
||||
PhasePressure::oil(G, reg, span, grav, cells,
|
||||
press[ oix ], po_woc, po_goc);
|
||||
PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
|
||||
press[ oix ], po_woc, po_goc);
|
||||
}
|
||||
|
||||
if (PhaseUsed::water(pu)) {
|
||||
const int wix = PhaseIndex::water(pu);
|
||||
PhasePressure::water(G, reg, span, grav, po_woc,
|
||||
PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
|
||||
cells, press[ wix ]);
|
||||
}
|
||||
|
||||
if (PhaseUsed::gas(pu)) {
|
||||
const int gix = PhaseIndex::gas(pu);
|
||||
PhasePressure::gas(G, reg, span, grav, po_goc,
|
||||
PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
|
||||
cells, press[ gix ]);
|
||||
}
|
||||
}
|
||||
@ -586,7 +564,8 @@ namespace Opm
|
||||
namespace EQUIL {
|
||||
|
||||
|
||||
template <class Grid,
|
||||
template <class FluidSystem,
|
||||
class Grid,
|
||||
class Region,
|
||||
class CellRange>
|
||||
std::vector< std::vector<double> >
|
||||
@ -654,7 +633,7 @@ namespace Opm
|
||||
span[0] = std::min(span[0],zgoc);
|
||||
span[1] = std::max(span[1],zwoc);
|
||||
|
||||
Details::equilibrateOWG(G, reg, grav, span, cells, press);
|
||||
Details::equilibrateOWG<FluidSystem>(G, reg, grav, span, cells, press);
|
||||
|
||||
return press;
|
||||
}
|
||||
@ -671,12 +650,12 @@ namespace Opm
|
||||
return std::vector<double>(cells.size(), 273.15 + 20.0);
|
||||
}
|
||||
|
||||
template <class Grid, class Region, class CellRange>
|
||||
template <class FluidSystem, class Grid, class Region, class CellRange, class MaterialLawManager>
|
||||
std::vector< std::vector<double> >
|
||||
phaseSaturations(const Grid& G,
|
||||
const Region& reg,
|
||||
const CellRange& cells,
|
||||
BlackoilPropertiesInterface& props,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const std::vector<double> swat_init,
|
||||
std::vector< std::vector<double> >& phase_pressures)
|
||||
{
|
||||
@ -685,8 +664,23 @@ namespace Opm
|
||||
}
|
||||
|
||||
std::vector< std::vector<double> > phase_saturations = phase_pressures; // Just to get the right size.
|
||||
double smin[BlackoilPhases::MaxNumPhases] = { 0.0 };
|
||||
double smax[BlackoilPhases::MaxNumPhases] = { 0.0 };
|
||||
|
||||
// Adjust oil pressure according to gas saturation and cap pressure
|
||||
typedef Opm::SimpleModularFluidState<double,
|
||||
/*numPhases=*/3,
|
||||
/*numComponents=*/3,
|
||||
FluidSystem,
|
||||
/*storePressure=*/false,
|
||||
/*storeTemperature=*/false,
|
||||
/*storeComposition=*/false,
|
||||
/*storeFugacity=*/false,
|
||||
/*storeSaturation=*/true,
|
||||
/*storeDensity=*/false,
|
||||
/*storeViscosity=*/false,
|
||||
/*storeEnthalpy=*/false> SatOnlyFluidState;
|
||||
|
||||
SatOnlyFluidState fluidState;
|
||||
typedef typename MaterialLawManager::MaterialLaw MaterialLaw;
|
||||
|
||||
const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
|
||||
const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
|
||||
@ -696,42 +690,45 @@ namespace Opm
|
||||
std::vector<double>::size_type local_index = 0;
|
||||
for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) {
|
||||
const int cell = *ci;
|
||||
props.satRange(1, &cell, smin, smax);
|
||||
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(props,waterpos,cell)){
|
||||
if (isConstPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,FluidSystem::waterPhaseIdx, cell)){
|
||||
const double cellDepth = UgGridHelpers::cellCenterDepth(G,
|
||||
cell);
|
||||
sw = satFromDepth(props,cellDepth,reg.zwoc(),waterpos,cell,false);
|
||||
sw = satFromDepth<FluidSystem, MaterialLaw, MaterialLawManager>(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(props, waterpos, cell, pcov);
|
||||
sw = satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, waterpos, cell, pcov);
|
||||
phase_saturations[waterpos][local_index] = sw;
|
||||
} else { // Scale Pc to reflect imposed sw
|
||||
sw = swat_init[cell];
|
||||
props.swatInitScaling(cell, pcov, sw);
|
||||
sw = materialLawManager->applySwatinit(cell, pcov, sw);
|
||||
phase_saturations[waterpos][local_index] = sw;
|
||||
}
|
||||
}
|
||||
}
|
||||
double sg = 0.0;
|
||||
if (gas) {
|
||||
if (isConstPc(props,gaspos,cell)){
|
||||
if (isConstPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,FluidSystem::gasPhaseIdx,cell)){
|
||||
const double cellDepth = UgGridHelpers::cellCenterDepth(G,
|
||||
cell);
|
||||
sg = satFromDepth(props,cellDepth,reg.zgoc(),gaspos,cell,true);
|
||||
sg = satFromDepth<FluidSystem, MaterialLaw, MaterialLawManager>(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(props, gaspos, cell, pcog, increasing);
|
||||
sg = satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, gaspos, cell, pcog, increasing);
|
||||
phase_saturations[gaspos][local_index] = sg;
|
||||
}
|
||||
}
|
||||
@ -745,55 +742,69 @@ namespace Opm
|
||||
// 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.
|
||||
props.swatInitScaling(cell, pcgw, sw);
|
||||
sw = materialLawManager->applySwatinit(cell, pcgw, sw);
|
||||
}
|
||||
sw = satFromSumOfPcs(props, waterpos, gaspos, cell, pcgw);
|
||||
sw = satFromSumOfPcs<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, waterpos, gaspos, cell, pcgw);
|
||||
sg = 1.0 - sw;
|
||||
phase_saturations[waterpos][local_index] = sw;
|
||||
phase_saturations[gaspos][local_index] = sg;
|
||||
// Adjust oil pressure according to gas saturation and cap pressure
|
||||
double pc[BlackoilPhases::MaxNumPhases];
|
||||
double sat[BlackoilPhases::MaxNumPhases];
|
||||
sat[waterpos] = sw;
|
||||
sat[gaspos] = sg;
|
||||
sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
|
||||
|
||||
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 pc[BlackoilPhases::MaxNumPhases];
|
||||
double sat[BlackoilPhases::MaxNumPhases];
|
||||
double threshold_sat = 1.0e-6;
|
||||
|
||||
sat[oilpos] = 1.0;
|
||||
double so = 1.0;
|
||||
double pC[/*numPhases=*/BlackoilPhases::MaxNumPhases] = { 0.0, 0.0, 0.0 };
|
||||
|
||||
if (water) {
|
||||
sat[waterpos] = smax[waterpos];
|
||||
sat[oilpos] -= sat[waterpos];
|
||||
double swu = scaledDrainageInfo.Swu;
|
||||
fluidState.setSaturation(FluidSystem::waterPhaseIdx, swu);
|
||||
so -= swu;
|
||||
}
|
||||
if (gas) {
|
||||
sat[gaspos] = smax[gaspos];
|
||||
sat[oilpos] -= sat[gaspos];
|
||||
double sgu = scaledDrainageInfo.Sgu;
|
||||
fluidState.setSaturation(FluidSystem::gasPhaseIdx, sgu);
|
||||
so-= sgu;
|
||||
}
|
||||
if (water && sw > smax[waterpos]-threshold_sat ) {
|
||||
sat[waterpos] = smax[waterpos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
|
||||
} else if (gas && sg > smax[gaspos]-threshold_sat) {
|
||||
sat[gaspos] = smax[gaspos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
|
||||
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 < smin[gaspos]+threshold_sat) {
|
||||
sat[gaspos] = smin[gaspos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
|
||||
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 < smin[waterpos]+threshold_sat) {
|
||||
sat[waterpos] = smin[waterpos];
|
||||
props.capPress(1, sat, &cell, pc, 0);
|
||||
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
|
||||
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;
|
||||
@ -885,9 +896,9 @@ namespace Opm
|
||||
* \param[in] applySwatInit Make it possible to not apply SWATINIT even if it
|
||||
* is present in the deck
|
||||
*/
|
||||
template<class Grid>
|
||||
template<class MaterialLawManager, class Grid>
|
||||
void initStateEquil(const Grid& grid,
|
||||
BlackoilPropertiesFromDeck& props,
|
||||
std::shared_ptr<MaterialLawManager> materialLawManager,
|
||||
const Opm::Deck& deck,
|
||||
const Opm::EclipseState& eclipseState,
|
||||
const double gravity,
|
||||
@ -896,6 +907,9 @@ namespace Opm
|
||||
{
|
||||
|
||||
typedef EQUIL::DeckDependent::InitialStateComputer ISC;
|
||||
|
||||
PhaseUsage pu = phaseUsageFromDeck(deck);
|
||||
|
||||
//Check for presence of kw SWATINIT
|
||||
std::vector<double> swat_init = {};
|
||||
if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatinit) {
|
||||
@ -910,8 +924,7 @@ namespace Opm
|
||||
}
|
||||
}
|
||||
|
||||
ISC isc(props, deck, eclipseState, grid, gravity, swat_init);
|
||||
const auto pu = props.phaseUsage();
|
||||
ISC isc(materialLawManager, pu, deck, eclipseState, grid, gravity, swat_init);
|
||||
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
|
||||
? pu.phase_pos[BlackoilPhases::Liquid]
|
||||
: pu.phase_pos[BlackoilPhases::Aqua];
|
||||
@ -920,7 +933,7 @@ namespace Opm
|
||||
state.gasoilratio() = isc.rs();
|
||||
state.rv() = isc.rv();
|
||||
|
||||
initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
|
||||
//initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
|
||||
}
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user