OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Remove BlackoilProps from equil initalization code

Browse Source 
Use FluidSystem and materialLaw from opm-material
directly not via the BlackoilProps in opm-core
This commit is contained in:
Tor Harald Sandve 2017-11-16 13:04:26 +01:00
parent d049125205
commit 8def11ff93
5 changed files with 715 additions and 538 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
examples/compute_initial_state.cpp
View File

@ -29,11 +29,16 @@
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/compressedToCartesian.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <boost/filesystem.hpp>
#include <fstream>
@ -97,10 +102,24 @@ try
BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param);
warnIfUnusedParams(param);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
// Forward declaring the MaterialLawManager template.
typedef Opm::ThreePhaseMaterialTraits<double,
/*wettingPhaseIdx=*/Opm::BlackoilPhases::Aqua,
/*nonWettingPhaseIdx=*/Opm::BlackoilPhases::Liquid,
/*gasPhaseIdx=*/Opm::BlackoilPhases::Vapour> MaterialTraits;
typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
// Initialisation.
//initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
BlackoilState state( UgGridHelpers::numCells(grid) , UgGridHelpers::numFaces(grid), 3);
initStateEquil(grid, props, deck, eclipseState, grav, state);
initStateEquil(grid, materialLawManager, deck, eclipseState, grav, state);
// Output.
const std::string output_dir = param.getDefault<std::string>("output_dir", "output");

368
opm/core/simulator/EquilibrationHelpers.hpp
View File

@ -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

131
opm/core/simulator/initStateEquil.hpp
View File

@ -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_);

305
opm/core/simulator/initStateEquil_impl.hpp
View File

@ -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);
}

428
tests/test_equil.cpp
View File

@ -22,10 +22,15 @@
#include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/utility/compressedToCartesian.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
@ -49,6 +54,15 @@
#include <string>
#include <vector>
// Forward declaring the MaterialLawManager template.
typedef Opm::ThreePhaseMaterialTraits<double,
/*wettingPhaseIdx=*/Opm::BlackoilPhases::Aqua,
/*nonWettingPhaseIdx=*/Opm::BlackoilPhases::Liquid,
/*gasPhaseIdx=*/Opm::BlackoilPhases::Vapour> MaterialTraits;
typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
#define CHECK(value, expected, reltol) \
{ \
if (std::fabs((expected)) < 1.e-14) \
@ -57,6 +71,87 @@
BOOST_CHECK_CLOSE((value), (expected), (reltol)); \
}
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
Opm::PhaseUsage initDefaultFluidSystem() {
std::vector<std::pair<double, double> > Bo = {
{ 101353, 1. },
{ 6.21542e+07, 1 }
};
std::vector<std::pair<double, double> > muo = {
{ 101353, 1. },
{ 6.21542e+07, 1 }
};
std::vector<std::pair<double, double> > Bg = {
{ 101353, 1. },
{ 6.21542e+07, 1 }
};
std::vector<std::pair<double, double> > mug = {
{ 101353, 1. },
{ 6.21542e+07, 1 }
};
double rhoRefO = 700; // [kg]
double rhoRefG = 1000; // [kg]
double rhoRefW = 1000; // [kg]
FluidSystem::initBegin(/*numPvtRegions=*/1);
FluidSystem::setEnableDissolvedGas(false);
FluidSystem::setEnableVaporizedOil(false);
FluidSystem::setReferenceDensities(rhoRefO, rhoRefW, rhoRefG, /*regionIdx=*/0);
Opm::GasPvtMultiplexer<double> *gasPvt = new Opm::GasPvtMultiplexer<double>;
gasPvt->setApproach(Opm::GasPvtMultiplexer<double>::DryGasPvt);
auto& dryGasPvt = gasPvt->template getRealPvt<Opm::GasPvtMultiplexer<double>::DryGasPvt>();
dryGasPvt.setNumRegions(/*numPvtRegion=*/1);
dryGasPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW);
dryGasPvt.setGasFormationVolumeFactor(/*regionIdx=*/0, Bg);
dryGasPvt.setGasViscosity(/*regionIdx=*/0, mug);
Opm::OilPvtMultiplexer<double> *oilPvt = new Opm::OilPvtMultiplexer<double>;
oilPvt->setApproach(Opm::OilPvtMultiplexer<double>::DeadOilPvt);
auto& deadOilPvt = oilPvt->template getRealPvt<Opm::OilPvtMultiplexer<double>::DeadOilPvt>();
deadOilPvt.setNumRegions(/*numPvtRegion=*/1);
deadOilPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW);
deadOilPvt.setInverseOilFormationVolumeFactor(/*regionIdx=*/0, Bo);
deadOilPvt.setOilViscosity(/*regionIdx=*/0, muo);
Opm::WaterPvtMultiplexer<double> *waterPvt = new Opm::WaterPvtMultiplexer<double>;
waterPvt->setApproach(Opm::WaterPvtMultiplexer<double>::ConstantCompressibilityWaterPvt);
auto& ccWaterPvt = waterPvt->template getRealPvt<Opm::WaterPvtMultiplexer<double>::ConstantCompressibilityWaterPvt>();
ccWaterPvt.setNumRegions(/*numPvtRegions=*/1);
ccWaterPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW);
ccWaterPvt.setViscosity(/*regionIdx=*/0, 1);
ccWaterPvt.setCompressibility(/*regionIdx=*/0, 0);
gasPvt->initEnd();
oilPvt->initEnd();
waterPvt->initEnd();
typedef std::shared_ptr<Opm::GasPvtMultiplexer<double> > GasPvtSharedPtr;
FluidSystem::setGasPvt(GasPvtSharedPtr(gasPvt));
typedef std::shared_ptr<Opm::OilPvtMultiplexer<double> > OilPvtSharedPtr;
FluidSystem::setOilPvt(OilPvtSharedPtr(oilPvt));
typedef std::shared_ptr<Opm::WaterPvtMultiplexer<double> > WaterPvtSharedPtr;
FluidSystem::setWaterPvt(WaterPvtSharedPtr(waterPvt));
FluidSystem::initEnd();
Opm::PhaseUsage pu;
pu.num_phases = 2;
// Might just as well assume water-oil.
pu.phase_used[Opm::BlackoilPhases::Aqua] = true;
pu.phase_used[Opm::BlackoilPhases::Liquid] = true;
pu.phase_used[Opm::BlackoilPhases::Vapour] = false;
pu.phase_pos[Opm::BlackoilPhases::Aqua] = 0;
pu.phase_pos[Opm::BlackoilPhases::Liquid] = 1;
pu.phase_pos[Opm::BlackoilPhases::Vapour] = 1; // Unused.
return pu;
}
BOOST_AUTO_TEST_SUITE ()
static Opm::EquilRecord mkEquilRecord( double datd, double datp,
@ -123,35 +218,22 @@ BOOST_AUTO_TEST_CASE (PhasePressure)
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::ParameterGroup param;
{
using Opm::unit::kilogram;
using Opm::unit::meter;
using Opm::unit::cubic;
std::stringstream dens; dens << 700*kilogram/cubic(meter);
param.insertParameter("rho2", dens.str());
}
typedef Opm::BlackoilPropertiesBasic Props;
Props props(param, G->dimensions, G->number_of_cells);
typedef Opm::EQUIL::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
RhoCalc calc(props, 0);
auto record = mkEquilRecord( 0, 1e5, 5, 0, 0, 0 );
Opm::EQUIL::EquilReg<RhoCalc>
region(record, calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage());
Opm::PhaseUsage pu = initDefaultFluidSystem();
Opm::EQUIL::EquilReg
region(record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu);
std::vector<int> cells(G->number_of_cells);
std::iota(cells.begin(), cells.end(), 0);
const double grav = 10;
const PPress ppress = Opm::EQUIL::phasePressures(*G, region, cells, grav);
const PPress ppress = Opm::EQUIL::phasePressures<FluidSystem>(*G, region, cells, grav);
const int first = 0, last = G->number_of_cells - 1;
const double reltol = 1.0e-8;
@ -169,47 +251,38 @@ BOOST_AUTO_TEST_CASE (CellSubset)
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::ParameterGroup param;
{
using Opm::unit::kilogram;
using Opm::unit::meter;
using Opm::unit::cubic;
Opm::PhaseUsage pu = initDefaultFluidSystem();
std::stringstream dens; dens << 700*kilogram/cubic(meter);
param.insertParameter("rho2", dens.str());
}
typedef Opm::BlackoilPropertiesBasic Props;
Props props(param, G->dimensions, G->number_of_cells);
typedef Opm::EQUIL::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
RhoCalc calc(props, 0);
Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
Opm::EQUIL::EquilReg<RhoCalc> region[] =
{
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
};
Opm::EQUIL::EquilReg region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
};
const int cdim[] = { 2, 1, 2 };
int ncoarse = cdim[0];
@ -239,7 +312,7 @@ BOOST_AUTO_TEST_CASE (CellSubset)
const int rno = int(r - cells.begin());
const double grav = 10;
const PPress p =
Opm::EQUIL::phasePressures(*G, region[rno], *r, grav);
Opm::EQUIL::phasePressures<FluidSystem>(*G, region[rno], *r, grav);
PVal::size_type i = 0;
for (std::vector<int>::const_iterator
@ -272,46 +345,36 @@ BOOST_AUTO_TEST_CASE (RegMapping)
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::ParameterGroup param;
{
using Opm::unit::kilogram;
using Opm::unit::meter;
using Opm::unit::cubic;
std::stringstream dens; dens << 700*kilogram/cubic(meter);
param.insertParameter("rho2", dens.str());
}
typedef Opm::BlackoilPropertiesBasic Props;
Props props(param, G->dimensions, G->number_of_cells);
typedef Opm::EQUIL::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
RhoCalc calc(props, 0);
Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
Opm::EQUIL::EquilReg<RhoCalc> region[] =
{
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
Opm::PhaseUsage pu = initDefaultFluidSystem();
Opm::EQUIL::EquilReg region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
};
std::vector<int> eqlnum(G->number_of_cells);
@ -336,7 +399,7 @@ BOOST_AUTO_TEST_CASE (RegMapping)
const int rno = r;
const double grav = 10;
const PPress p =
Opm::EQUIL::phasePressures(*G, region[rno], rng, grav);
Opm::EQUIL::phasePressures<FluidSystem>(*G, region[rno], rng, grav);
PVal::size_type i = 0;
for (const auto& c : rng) {
@ -366,9 +429,16 @@ BOOST_AUTO_TEST_CASE (DeckAllDead)
Opm::ParseContext parseContext;
Opm::Parser parser;
Opm::Deck deck = parser.parseFile("deadfluids.DATA" , parseContext);
Opm::EclipseState eclipseState(deck, parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, *grid, false);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, *grid, 10.0);
Opm::EclipseState eclipseState(deck, parseContext); // Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid->number_of_cells, grid->global_cell);
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, *grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells);
@ -395,8 +465,20 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
typedef MaterialLawManager::MaterialLaw MaterialLaw;
if (!FluidSystem::isInitialized()) {
// make sure that we don't initialize the fluid system twice
FluidSystem::initFromDeck(deck, eclipseState);
}
// Test the capillary inversion for oil-water.
const int cell = 0;
const double reltol = 1.0e-7;
@ -407,7 +489,7 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
const std::vector<double> s = { 0.2, 0.2, 0.2, 0.466666666666, 0.733333333333, 1.0, 1.0, 1.0, 1.0 };
BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromPc(props, phase, cell, pc[i], increasing);
const double s_computed = Opm::EQUIL::satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, phase, cell, pc[i], increasing);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
@ -420,7 +502,7 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
const std::vector<double> s = { 0.8, 0.8, 0.8, 0.533333333333, 0.266666666666, 0.0, 0.0, 0.0, 0.0 };
BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromPc(props, phase, cell, pc[i], increasing);
const double s_computed = Opm::EQUIL::satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, phase, cell, pc[i], increasing);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
@ -433,7 +515,7 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
const std::vector<double> s = { 0.2, 0.333333333333, 0.6, 0.866666666666, 1.0 };
BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromSumOfPcs(props, water, gas, cell, pc[i]);
const double s_computed = Opm::EQUIL::satFromSumOfPcs<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, water, gas, cell, pc[i]);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
@ -449,9 +531,16 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 10.0);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -490,9 +579,15 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary_overlap.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -553,9 +648,15 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("equil_liveoil.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -633,9 +734,15 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("equil_livegas.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -715,9 +822,15 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
Opm::Parser parser;
Opm::Deck deck = parser.parseFile("equil_rsvd_and_rvvd.DATA", parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -817,8 +930,15 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
Opm::EclipseState eclipseState(deck , parseContext);
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
Opm::BlackoilPropertiesFromDeck propsScaled(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
auto materialLawManagerScaled = std::make_shared<MaterialLawManager>();
materialLawManagerScaled->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::BlackoilState state( Opm::UgGridHelpers::numCells( grid ) , Opm::UgGridHelpers::numFaces( grid ) , 3);
@ -848,12 +968,43 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
}
}
// 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 MaterialLawManager::MaterialLaw MaterialLaw;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
// reference pcs
const int numCells = Opm::UgGridHelpers::numCells(grid);
std::vector<int> cells(numCells);
for (int c = 0; c < numCells; ++c) { cells[c] = c; }
std::vector<double> pc_original = state.saturation();
props.capPress(numCells, sats.data(), cells.data(), pc_original.data(), nullptr);
const int numCells = Opm::UgGridHelpers::numCells(grid);
for (int c = 0; c < numCells; ++c) {
std::vector<double> pc = {0,0,0};
double sw = sats[3*c + 0];
double so = sats[3*c + 1];
double sg = sats[3*c + 2];
fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg);
const auto& matParams = materialLawManager->materialLawParams(c);
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
pc_original[3*c + 0] = pc[FluidSystem::oilPhaseIdx] - pc[FluidSystem::waterPhaseIdx];
pc_original[3*c + 1] = 0.0;
pc_original[3*c + 2] = pc[FluidSystem::oilPhaseIdx] + pc[FluidSystem::gasPhaseIdx];
}
std::vector<double> pc_scaled_truth = pc_original;
@ -876,17 +1027,46 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
// apply swatinit
Opm::BlackoilState state_scaled = state;
initStateEquil(grid, propsScaled, deck, eclipseState, 9.81, state_scaled, true);
initStateEquil(grid, materialLawManagerScaled, deck, eclipseState, 9.81, state_scaled, true);
// don't apply swatinit
Opm::BlackoilState state_unscaled = state;
initStateEquil(grid, props, deck, eclipseState, 9.81, state_unscaled, false);
initStateEquil(grid, materialLawManager, deck, eclipseState, 9.81, state_unscaled, false);
// compute pc
std::vector<double> pc_scaled= state.saturation();
propsScaled.capPress(numCells, state_scaled.saturation().data(), cells.data(), pc_scaled.data(), nullptr);
for (int c = 0; c < numCells; ++c) {
std::vector<double> pc = {0,0,0};
double sw = state_scaled.saturation().data()[3*c + 0];
double so = state_scaled.saturation().data()[3*c + 1];
double sg = state_scaled.saturation().data()[3*c + 2];
fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg);
const auto& matParams = materialLawManagerScaled->materialLawParams(c);
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
pc_scaled[3*c + 0] = pc[FluidSystem::oilPhaseIdx] - pc[FluidSystem::waterPhaseIdx];
pc_scaled[3*c + 1] = 0.0;
pc_scaled[3*c + 2] = pc[FluidSystem::oilPhaseIdx] + pc[FluidSystem::gasPhaseIdx];
}
std::vector<double> pc_unscaled= state.saturation();
props.capPress(numCells, state_unscaled.saturation().data(), cells.data(), pc_unscaled.data(), nullptr);
for (int c = 0; c < numCells; ++c) {
std::vector<double> pc = {0,0,0};
double sw = state_unscaled.saturation().data()[3*c + 0];
double so = state_unscaled.saturation().data()[3*c + 1];
double sg = state_unscaled.saturation().data()[3*c + 2];
fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg);
const auto& matParams = materialLawManager->materialLawParams(c);
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
pc_unscaled[3*c + 0] = pc[FluidSystem::oilPhaseIdx] - pc[FluidSystem::waterPhaseIdx];
pc_unscaled[3*c + 1] = 0.0;
pc_unscaled[3*c + 2] = pc[FluidSystem::oilPhaseIdx] + pc[FluidSystem::gasPhaseIdx];
}
// test
const double reltol = 1.0e-3;