OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues 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 committed by Andreas Lauser
parent 9e95d662bd
commit b9b75a109c
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/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp> #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/simulator/BlackoilState.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/ParseContext.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp> #include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp> #include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <boost/filesystem.hpp> #include <boost/filesystem.hpp>
#include <fstream> #include <fstream>
@ -97,10 +102,24 @@ try
BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param); BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param);
warnIfUnusedParams(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. // Initialisation.
//initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state); //initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
BlackoilState state( UgGridHelpers::numCells(grid) , UgGridHelpers::numFaces(grid), 3); BlackoilState state( UgGridHelpers::numCells(grid) , UgGridHelpers::numFaces(grid), 3);
initStateEquil(grid, props, deck, eclipseState, grav, state); initStateEquil(grid, materialLawManager, deck, eclipseState, grav, state);
// Output. // Output.
const std::string output_dir = param.getDefault<std::string>("output_dir", "output"); const std::string output_dir = param.getDefault<std::string>("output_dir", "output");

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

@ -20,7 +20,6 @@
#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED #ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED #define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp> #include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/linearInterpolation.hpp> #include <opm/core/utility/linearInterpolation.hpp>
#include <opm/core/utility/RegionMapping.hpp> #include <opm/core/utility/RegionMapping.hpp>
@ -28,6 +27,10 @@
#include <opm/parser/eclipse/EclipseState/InitConfig/Equil.hpp> #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> #include <memory>
@ -87,71 +90,21 @@ namespace Opm
namespace EQUIL { namespace EQUIL {
template <class Props> typedef Opm::FluidSystems::BlackOil<double> FluidSystemSimple;
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_;
};
// 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 * Types and routines relating to phase mixing in
@ -224,6 +177,7 @@ namespace Opm
* tabulated as a function of depth policy. Data * tabulated as a function of depth policy. Data
* typically taken from keyword 'RSVD'. * typically taken from keyword 'RSVD'.
*/ */
template <class FluidSystem>
class RsVD : public RsFunction { class RsVD : public RsFunction {
public: public:
/** /**
@ -234,19 +188,14 @@ namespace Opm
* \param[in] depth Depth nodes. * \param[in] depth Depth nodes.
* \param[in] rs Dissolved gas-oil ratio at @c depth. * \param[in] rs Dissolved gas-oil ratio at @c depth.
*/ */
RsVD(const BlackoilPropertiesInterface& props, RsVD(const int pvtRegionIdx,
const int cell,
const std::vector<double>& depth, const std::vector<double>& depth,
const std::vector<double>& rs) const std::vector<double>& rs)
: props_(props) : pvtRegionIdx_(pvtRegionIdx)
, cell_(cell)
, depth_(depth) , depth_(depth)
, rs_(rs) , 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: private:
const BlackoilPropertiesInterface& props_; const int pvtRegionIdx_;
const int cell_;
std::vector<double> depth_; /**< Depth nodes */ std::vector<double> depth_; /**< Depth nodes */
std::vector<double> rs_; /**< Dissolved gas-oil ratio */ 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 double satRs(const double press, const double temp) const
{ {
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0); return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
// 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];
} }
}; };
@ -304,6 +243,7 @@ namespace Opm
* tabulated as a function of depth policy. Data * tabulated as a function of depth policy. Data
* typically taken from keyword 'RVVD'. * typically taken from keyword 'RVVD'.
*/ */
template <class FluidSystem>
class RvVD : public RsFunction { class RvVD : public RsFunction {
public: public:
/** /**
@ -314,19 +254,14 @@ namespace Opm
* \param[in] depth Depth nodes. * \param[in] depth Depth nodes.
* \param[in] rv Dissolved gas-oil ratio at @c depth. * \param[in] rv Dissolved gas-oil ratio at @c depth.
*/ */
RvVD(const BlackoilPropertiesInterface& props, RvVD(const int pvtRegionIdx,
const int cell,
const std::vector<double>& depth, const std::vector<double>& depth,
const std::vector<double>& rv) const std::vector<double>& rv)
: props_(props) : pvtRegionIdx_(pvtRegionIdx)
, cell_(cell)
, depth_(depth) , depth_(depth)
, rv_(rv) , 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: private:
const BlackoilPropertiesInterface& props_; const int pvtRegionIdx_;
const int cell_;
std::vector<double> depth_; /**< Depth nodes */ std::vector<double> depth_; /**< Depth nodes */
std::vector<double> rv_; /**< Vaporized oil-gas ratio */ 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 double satRv(const double press, const double temp) const
{ {
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0); return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
// 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];
} }
}; };
@ -393,6 +318,7 @@ namespace Opm
* This should yield Rs-values that are constant below the * This should yield Rs-values that are constant below the
* contact, and decreasing above the contact. * contact, and decreasing above the contact.
*/ */
template <class FluidSystem>
class RsSatAtContact : public RsFunction { class RsSatAtContact : public RsFunction {
public: public:
/** /**
@ -403,13 +329,9 @@ namespace Opm
* \param[in] p_contact oil pressure at the contact * \param[in] p_contact oil pressure at the contact
* \param[in] T_contact temperature 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) RsSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact)
: props_(props), cell_(cell) : 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); rs_sat_contact_ = satRs(p_contact, T_contact);
} }
@ -442,22 +364,12 @@ namespace Opm
} }
private: private:
const BlackoilPropertiesInterface& props_; const int pvtRegionIdx_;
const int cell_;
double z_[BlackoilPhases::MaxNumPhases];
double rs_sat_contact_; double rs_sat_contact_;
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRs(const double press, const double temp) const double satRs(const double press, const double temp) const
{ {
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0); return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
// 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];
} }
}; };
@ -476,6 +388,7 @@ namespace Opm
* This should yield Rv-values that are constant below the * This should yield Rv-values that are constant below the
* contact, and decreasing above the contact. * contact, and decreasing above the contact.
*/ */
template <class FluidSystem>
class RvSatAtContact : public RsFunction { class RvSatAtContact : public RsFunction {
public: public:
/** /**
@ -486,13 +399,9 @@ namespace Opm
* \param[in] p_contact oil pressure at the contact * \param[in] p_contact oil pressure at the contact
* \param[in] T_contact temperature 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) RvSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact)
: props_(props), cell_(cell) :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); rv_sat_contact_ = satRv(p_contact, T_contact);
} }
@ -525,22 +434,12 @@ namespace Opm
} }
private: private:
const BlackoilPropertiesInterface& props_; const int pvtRegionIdx_;
const int cell_;
double z_[BlackoilPhases::MaxNumPhases];
double rv_sat_contact_; double rv_sat_contact_;
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRv(const double press, const double temp) const double satRv(const double press, const double temp) const
{ {
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0); return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
// 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];
} }
}; };
@ -565,7 +464,6 @@ namespace Opm
* that calculates the phase densities of all phases in @c * that calculates the phase densities of all phases in @c
* svol at fluid pressure @c press. * svol at fluid pressure @c press.
*/ */
template <class DensCalc>
class EquilReg { class EquilReg {
public: public:
/** /**
@ -578,23 +476,19 @@ namespace Opm
* \param[in] pu Summary of current active phases. * \param[in] pu Summary of current active phases.
*/ */
EquilReg(const EquilRecord& rec, EquilReg(const EquilRecord& rec,
const DensCalc& density,
std::shared_ptr<Miscibility::RsFunction> rs, std::shared_ptr<Miscibility::RsFunction> rs,
std::shared_ptr<Miscibility::RsFunction> rv, std::shared_ptr<Miscibility::RsFunction> rv,
const int pvtIdx,
const PhaseUsage& pu) const PhaseUsage& pu)
: rec_ (rec) : rec_ (rec)
, density_(density)
, rs_ (rs) , rs_ (rs)
, rv_ (rv) , rv_ (rv)
, pvtIdx_ (pvtIdx)
, pu_ (pu) , pu_ (pu)
{ {
} }
/**
* Type of density calculator.
*/
typedef DensCalc CalcDensity;
/** /**
* Type of dissolved gas-oil ratio calculator. * Type of dissolved gas-oil ratio calculator.
*/ */
@ -639,11 +533,6 @@ namespace Opm
*/ */
double pcgo_goc() const { return this->rec_.gasOilContactCapillaryPressure(); } 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 * Retrieve dissolved gas-oil ratio calculator of current
@ -659,6 +548,12 @@ namespace Opm
const CalcEvaporation& const CalcEvaporation&
evaporationCalculator() const { return *this->rv_; } evaporationCalculator() const { return *this->rv_; }
/**
* Retrieve active fluid phase summary.
*/
const int
pvtIdx() const { return this->pvtIdx_; }
/** /**
* Retrieve active fluid phase summary. * Retrieve active fluid phase summary.
*/ */
@ -667,9 +562,9 @@ namespace Opm
private: private:
EquilRecord rec_; /**< Equilibration data */ EquilRecord rec_; /**< Equilibration data */
DensCalc density_; /**< Density calculator */
std::shared_ptr<Miscibility::RsFunction> rs_; /**< RS calculator */ std::shared_ptr<Miscibility::RsFunction> rs_; /**< RS calculator */
std::shared_ptr<Miscibility::RsFunction> rv_; /**< RV calculator */ std::shared_ptr<Miscibility::RsFunction> rv_; /**< RV calculator */
const int pvtIdx_;
PhaseUsage pu_; /**< Active phase summary */ PhaseUsage pu_; /**< Active phase summary */
}; };
@ -678,54 +573,117 @@ namespace Opm
/// Functor for inverting capillary pressure function. /// Functor for inverting capillary pressure function.
/// Function represented is /// Function represented is
/// f(s) = pc(s) - target_pc /// f(s) = pc(s) - target_pc
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
struct PcEq struct PcEq
{ {
PcEq(const BlackoilPropertiesInterface& props, PcEq(std::shared_ptr<MaterialLawManager> materialLawManager,
const int phase, const int phase,
const int cell, const int cell,
const double target_pc) const double target_pc)
: props_(props), : materialLawManager_(materialLawManager),
phase_(phase), phase_(phase),
cell_(cell), cell_(cell),
target_pc_(target_pc) 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); std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
} }
double operator()(double s) const double operator()(double s) const
{ {
s_[phase_] = s; const auto& matParams = materialLawManager_->materialLawParams(cell_);
props_.capPress(1, s_, &cell_, pc_, 0); fluidState_.setSaturation(phase_, s);
return pc_[phase_] - target_pc_; 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: private:
const BlackoilPropertiesInterface& props_;
std::shared_ptr<MaterialLawManager> materialLawManager_;
const int phase_; const int phase_;
const int cell_; const int cell_;
const double target_pc_; const double target_pc_;
mutable double s_[BlackoilPhases::MaxNumPhases]; mutable SatOnlyFluidState fluidState_;
mutable double pc_[BlackoilPhases::MaxNumPhases]; 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 /// Compute saturation of some phase corresponding to a given
/// capillary pressure. /// 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 phase,
const int cell, const int cell,
const double target_pc, const double target_pc,
const bool increasing = false) const bool increasing = false)
{ {
// Find minimum and maximum saturations. // Find minimum and maximum saturations.
double sminarr[BlackoilPhases::MaxNumPhases]; const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
double smaxarr[BlackoilPhases::MaxNumPhases]; const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
props.satRange(1, &cell, sminarr, smaxarr);
const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
// Create the equation f(s) = pc(s) - target_pc // 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 f0 = f(s0);
const double f1 = f(s1); const double f1 = f(s1);
@ -747,36 +705,46 @@ namespace Opm
/// Functor for inverting a sum of capillary pressure functions. /// Functor for inverting a sum of capillary pressure functions.
/// Function represented is /// Function represented is
/// f(s) = pc1(s) + pc2(1 - s) - target_pc /// f(s) = pc1(s) + pc2(1 - s) - target_pc
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
struct PcEqSum struct PcEqSum
{ {
PcEqSum(const BlackoilPropertiesInterface& props, PcEqSum(std::shared_ptr<MaterialLawManager> materialLawManager,
const int phase1, const int phase1,
const int phase2, const int phase2,
const int cell, const int cell,
const double target_pc) const double target_pc)
: props_(props), : materialLawManager_(materialLawManager),
phase1_(phase1), phase1_(phase1),
phase2_(phase2), phase2_(phase2),
cell_(cell), cell_(cell),
target_pc_(target_pc) 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); std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
} }
double operator()(double s) const double operator()(double s) const
{ {
s_[phase1_] = s; const auto& matParams = materialLawManager_->materialLawParams(cell_);
s_[phase2_] = 1.0 - s; fluidState_.setSaturation(phase1_, s);
props_.capPress(1, s_, &cell_, pc_, 0); fluidState_.setSaturation(phase2_, 1.0 - s);
return pc_[phase1_] + pc_[phase2_] - target_pc_;
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: private:
const BlackoilPropertiesInterface& props_; std::shared_ptr<MaterialLawManager> materialLawManager_;
const int phase1_; const int phase1_;
const int phase2_; const int phase2_;
const int cell_; const int cell_;
const double target_pc_; const double target_pc_;
mutable double s_[BlackoilPhases::MaxNumPhases]; mutable SatOnlyFluidState fluidState_;
mutable double pc_[BlackoilPhases::MaxNumPhases]; mutable double pc_[BlackoilPhases::MaxNumPhases];
}; };
@ -786,21 +754,20 @@ namespace Opm
/// Compute saturation of some phase corresponding to a given /// Compute saturation of some phase corresponding to a given
/// capillary pressure, where the capillary pressure function /// capillary pressure, where the capillary pressure function
/// is given as a sum of two other functions. /// 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 phase1,
const int phase2, const int phase2,
const int cell, const int cell,
const double target_pc) const double target_pc)
{ {
// Find minimum and maximum saturations. // Find minimum and maximum saturations.
double sminarr[BlackoilPhases::MaxNumPhases]; const double smin = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
double smaxarr[BlackoilPhases::MaxNumPhases]; const double smax = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
props.satRange(1, &cell, sminarr, smaxarr);
const double smin = sminarr[phase1];
const double smax = smaxarr[phase1];
// Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc // 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 f0 = f(smin);
const double f1 = f(smax); const double f1 = f(smax);
if (f0 <= 0.0) { if (f0 <= 0.0) {
@ -818,19 +785,16 @@ namespace Opm
} }
/// Compute saturation from depth. Used for constant capillary pressure function /// 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 cellDepth,
const double contactDepth, const double contactDepth,
const int phase, const int phase,
const int cell, const int cell,
const bool increasing = false) const bool increasing = false)
{ {
// Find minimum and maximum saturations. const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
double sminarr[BlackoilPhases::MaxNumPhases]; const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
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];
if (cellDepth < contactDepth){ if (cellDepth < contactDepth){
return s0; return s0;
@ -841,22 +805,20 @@ namespace Opm
} }
/// Return true if capillary pressure function is constant /// 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 phase,
const int cell) 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); // Create the equation f(s) = pc(s);
const PcEq f(props, phase, cell, 0); const PcEq<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase, cell, 0);
const double f0 = f(sminarr[phase]); const double f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
const double f1 = f(smaxarr[phase]); const double f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon(); return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
} }
} // namespace Equil } // namespace Equil
} // namespace Opm } // namespace Opm

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

@ -25,9 +25,11 @@
#include <opm/core/grid/GridHelpers.hpp> #include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/simulator/EquilibrationHelpers.hpp> #include <opm/core/simulator/EquilibrationHelpers.hpp>
#include <opm/core/simulator/BlackoilState.hpp> #include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp> #include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/RegionMapping.hpp> #include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/utility/extractPvtTableIndex.hpp>
#include <opm/parser/eclipse/Units/Units.hpp> #include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp> #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp> #include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
@ -39,6 +41,11 @@
#include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp> #include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
#include <opm/common/OpmLog/OpmLog.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 <array>
#include <cassert> #include <cassert>
#include <utility> #include <utility>
@ -54,32 +61,6 @@ struct UnstructuredGrid;
namespace Opm 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 * Types and routines that collectively implement a basic
@ -90,6 +71,7 @@ namespace Opm
*/ */
namespace EQUIL { namespace EQUIL {
/** /**
* Compute initial phase pressures by means of equilibration. * Compute initial phase pressures by means of equilibration.
* *
@ -128,7 +110,7 @@ namespace Opm
* of pressure values in each cell in the current * of pressure values in each cell in the current
* equilibration region. * equilibration region.
*/ */
template <class Grid, class Region, class CellRange> template <class Grid, class Region, class CellRange, class FluidSystem>
std::vector< std::vector<double> > std::vector< std::vector<double> >
phasePressures(const Grid& G, phasePressures(const Grid& G,
const Region& reg, const Region& reg,
@ -161,12 +143,12 @@ namespace Opm
* \return Phase saturations, one vector for each phase, each containing * \return Phase saturations, one vector for each phase, each containing
* one saturation value per cell in the region. * 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> > std::vector< std::vector<double> >
phaseSaturations(const Grid& grid, phaseSaturations(const Grid& grid,
const Region& reg, const Region& reg,
const CellRange& cells, const CellRange& cells,
BlackoilPropertiesFromDeck& props, std::shared_ptr<MaterialLawManager> materialLawManager,
const std::vector<double> swat_init, const std::vector<double> swat_init,
std::vector< std::vector<double> >& phase_pressures); std::vector< std::vector<double> >& phase_pressures);
@ -244,29 +226,39 @@ namespace Opm
class InitialStateComputer { class InitialStateComputer {
public: public:
template<class Grid> template<class MaterialLawManager, class Grid>
InitialStateComputer(BlackoilPropertiesInterface& props, InitialStateComputer(std::shared_ptr<MaterialLawManager> materialLawManager,
const PhaseUsage& phaseUsage,
const Opm::Deck& deck, const Opm::Deck& deck,
const Opm::EclipseState& eclipseState, const Opm::EclipseState& eclipseState,
const Grid& G , const Grid& G ,
const double grav = unit::gravity, const double grav = unit::gravity,
const std::vector<double>& swat_init = {} const std::vector<double>& swat_init = {}
) )
: pp_(props.numPhases(), : pp_(phaseUsage.num_phases,
std::vector<double>(UgGridHelpers::numCells(G))), std::vector<double>(UgGridHelpers::numCells(G))),
sat_(props.numPhases(), sat_(phaseUsage.num_phases,
std::vector<double>(UgGridHelpers::numCells(G))), std::vector<double>(UgGridHelpers::numCells(G))),
rs_(UgGridHelpers::numCells(G)), rs_(UgGridHelpers::numCells(G)),
rv_(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. // Get the equilibration records.
const std::vector<EquilRecord> rec = getEquil(eclipseState); const std::vector<EquilRecord> rec = getEquil(eclipseState);
const auto& tables = eclipseState.getTableManager(); const auto& tables = eclipseState.getTableManager();
// Create (inverse) region mapping. // Create (inverse) region mapping.
const RegionMapping<> eqlmap(equilnum(deck, eclipseState, G)); const RegionMapping<> eqlmap(equilnum(deck, eclipseState, G));
setRegionPvtIdx(G, eclipseState, eqlmap);
// Create Rs functions. // Create Rs functions.
rs_func_.reserve(rec.size()); rs_func_.reserve(rec.size());
if (deck.hasKeyword("DISGAS")) { if (deck.hasKeyword("DISGAS")) {
@ -274,10 +266,11 @@ namespace Opm
for (size_t i = 0; i < rec.size(); ++i) { for (size_t i = 0; i < rec.size(); ++i) {
if (eqlmap.cells(i).empty()) 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; continue;
} }
const int cell = *(eqlmap.cells(i).begin()); const int pvtIdx = regionPvtIdx_[i];
if (!rec[i].liveOilInitConstantRs()) { if (!rec[i].liveOilInitConstantRs()) {
if (rsvdTables.size() <= 0 ) { if (rsvdTables.size() <= 0 ) {
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table not available."); OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table not available.");
@ -285,8 +278,7 @@ namespace Opm
const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i); const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i);
std::vector<double> depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy(); std::vector<double> depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> rsColumn = rsvdTable.getColumn("RS").vectorCopy(); std::vector<double> rsColumn = rsvdTable.getColumn("RS").vectorCopy();
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, rs_func_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
cell,
depthColumn , rsColumn)); depthColumn , rsColumn));
} else { } else {
if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
@ -297,7 +289,7 @@ namespace Opm
} }
const double p_contact = rec[i].datumDepthPressure(); const double p_contact = rec[i].datumDepthPressure();
const double T_contact = 273.15 + 20; // standard temperature for now 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 { } else {
@ -312,10 +304,10 @@ namespace Opm
for (size_t i = 0; i < rec.size(); ++i) { for (size_t i = 0; i < rec.size(); ++i) {
if (eqlmap.cells(i).empty()) 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; continue;
} }
const int cell = *(eqlmap.cells(i).begin()); const int pvtIdx = regionPvtIdx_[i];
if (!rec[i].wetGasInitConstantRv()) { if (!rec[i].wetGasInitConstantRv()) {
if (rvvdTables.size() <= 0) { if (rvvdTables.size() <= 0) {
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table not available."); OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table not available.");
@ -324,9 +316,7 @@ namespace Opm
const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i); const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i);
std::vector<double> depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy(); std::vector<double> depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> rvColumn = rvvdTable.getColumn("RV").vectorCopy(); std::vector<double> rvColumn = rvvdTable.getColumn("RV").vectorCopy();
rv_func_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props,
cell,
depthColumn , rvColumn)); depthColumn , rvColumn));
} else { } else {
@ -338,7 +328,7 @@ namespace Opm
} }
const double p_contact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); const double p_contact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
const double T_contact = 273.15 + 20; // standard temperature for now 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 { } else {
@ -348,7 +338,7 @@ namespace Opm
} }
// Compute pressures, saturations, rs and rv factors. // 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 // Modify oil pressure in no-oil regions so that the pressures of present phases can
// be recovered from the oil pressure and capillary relations. // be recovered from the oil pressure and capillary relations.
@ -363,23 +353,39 @@ namespace Opm
const Vec& rv() const { return rv_; } const Vec& rv() const { return rv_; }
private: 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> > rs_func_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_; std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
std::vector<int> regionPvtIdx_;
PVec pp_; PVec pp_;
PVec sat_; PVec sat_;
Vec rs_; Vec rs_;
Vec rv_; Vec rv_;
Vec swat_init_; 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 void
calcPressSatRsRv(const RMap& reg , calcPressSatRsRv(const RMap& reg ,
const std::vector< EquilRecord >& rec , const std::vector< EquilRecord >& rec ,
Opm::BlackoilPropertiesInterface& props, std::shared_ptr<MaterialLawManager> materialLawManager,
const Grid& G , const Grid& G ,
const double grav) const double grav)
{ {
@ -391,27 +397,24 @@ namespace Opm
+ " has no active cells"); + " has no active cells");
continue; continue;
} }
const int repcell = *cells.begin();
const RhoCalc calc(props, repcell); const EqReg eqreg(rec[r],
const EqReg eqreg(rec[r], calc, rs_func_[r], rv_func_[r], regionPvtIdx_[r],
rs_func_[r], rv_func_[r], phaseUsage_);
props.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 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 = phaseUsage_.num_phases;
const int np = props.numPhases();
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
copyFromRegion(pressures[p], cells, pp_[p]); copyFromRegion(pressures[p], cells, pp_[p]);
copyFromRegion(sat[p], cells, sat_[p]); copyFromRegion(sat[p], cells, sat_[p]);
} }
if (props.phaseUsage().phase_used[BlackoilPhases::Liquid] if (phaseUsage_.phase_used[BlackoilPhases::Liquid]
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) { && phaseUsage_.phase_used[BlackoilPhases::Vapour]) {
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid]; const int oilpos = phaseUsage_.phase_pos[BlackoilPhases::Liquid];
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour]; 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 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]); const Vec rv_vals = computeRs(G, cells, pressures[gaspos], temp, *(rv_func_[r]), sat[oilpos]);
copyFromRegion(rs_vals, cells, rs_); copyFromRegion(rs_vals, cells, rs_);

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

@ -24,8 +24,8 @@
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/grid/GridHelpers.hpp> #include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/props/BlackoilPhases.hpp> #include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/simulator/initState.hpp>
#include <cassert> #include <cassert>
#include <cmath> #include <cmath>
@ -35,6 +35,9 @@
namespace Opm namespace Opm
{ {
namespace Details { namespace Details {
template <class RHS> template <class RHS>
class RK4IVP { class RK4IVP {
public: public:
@ -106,21 +109,16 @@ namespace Opm
}; };
namespace PhasePressODE { namespace PhasePressODE {
template <class Density> template <class FluidSystem>
class Water { class Water {
public: public:
Water(const double temp, Water(const double temp,
const Density& rho, const int pvtRegionIdx,
const int np,
const int ix,
const double norm_grav) const double norm_grav)
: temp_(temp) : temp_(temp)
, rho_(rho) , pvtRegionIdx_(pvtRegionIdx)
, svol_(np, 0)
, ix_(ix)
, g_(norm_grav) , g_(norm_grav)
{ {
svol_[ix_] = 1.0;
} }
double double
@ -132,39 +130,30 @@ namespace Opm
private: private:
const double temp_; const double temp_;
const Density& rho_; const int pvtRegionIdx_;
std::vector<double> svol_;
const int ix_;
const double g_; const double g_;
double double
density(const double press) const density(const double press) const
{ {
const std::vector<double>& rho = rho_(press, temp_, svol_); double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
return rho[ix_]; return rho;
} }
}; };
template <class Density, class RS> template <class FluidSystem, class RS>
class Oil { class Oil {
public: public:
Oil(const double temp, Oil(const double temp,
const Density& rho,
const RS& rs, const RS& rs,
const int np, const int pvtRegionIdx,
const int oix,
const int gix,
const double norm_grav) const double norm_grav)
: temp_(temp) : temp_(temp)
, rho_(rho)
, rs_(rs) , rs_(rs)
, svol_(np, 0) , pvtRegionIdx_(pvtRegionIdx)
, oix_(oix)
, gix_(gix)
, g_(norm_grav) , g_(norm_grav)
{ {
svol_[oix_] = 1.0;
} }
double double
@ -176,45 +165,42 @@ namespace Opm
private: private:
const double temp_; const double temp_;
const Density& rho_;
const RS& rs_; const RS& rs_;
mutable std::vector<double> svol_; const int pvtRegionIdx_;
const int oix_;
const int gix_;
const double g_; const double g_;
double double
density(const double depth, density(const double depth,
const double press) const const double press) const
{ {
if (gix_ >= 0) { double rs = rs_(depth, press, temp_);
svol_[gix_] = 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;
return rho[oix_];
} }
}; };
template <class Density, class RV> template <class FluidSystem, class RV>
class Gas { class Gas {
public: public:
Gas(const double temp, Gas(const double temp,
const Density& rho,
const RV& rv, const RV& rv,
const int np, const int pvtRegionIdx,
const int gix,
const int oix,
const double norm_grav) const double norm_grav)
: temp_(temp) : temp_(temp)
, rho_(rho)
, rv_(rv) , rv_(rv)
, svol_(np, 0) , pvtRegionIdx_(pvtRegionIdx)
, gix_(gix)
, oix_(oix)
, g_(norm_grav) , g_(norm_grav)
{ {
svol_[gix_] = 1.0;
} }
double double
@ -226,23 +212,27 @@ namespace Opm
private: private:
const double temp_; const double temp_;
const Density& rho_;
const RV& rv_; const RV& rv_;
mutable std::vector<double> svol_; const int pvtRegionIdx_;
const int gix_;
const int oix_;
const double g_; const double g_;
double double
density(const double depth, density(const double depth,
const double press) const const double press) const
{ {
if (oix_ >= 0) { double rv = rv_(depth, press, temp_);
svol_[oix_] = 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;
return rho[gix_];
} }
}; };
} // namespace PhasePressODE } // namespace PhasePressODE
@ -328,7 +318,8 @@ namespace Opm
} }
} }
template <class Grid, template < class FluidSystem,
class Grid,
class Region, class Region,
class CellRange> class CellRange>
void void
@ -341,13 +332,10 @@ namespace Opm
std::vector<double>& press ) std::vector<double>& press )
{ {
using PhasePressODE::Water; 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 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 z0;
double p0; double p0;
@ -379,9 +367,11 @@ namespace Opm
} }
} }
template <class Grid, template <class FluidSystem,
class Grid,
class Region, class Region,
class CellRange> class CellRange>
void void
oil(const Grid& G , oil(const Grid& G ,
const Region& reg , const Region& reg ,
@ -393,18 +383,12 @@ namespace Opm
double& po_goc) double& po_goc)
{ {
using PhasePressODE::Oil; using PhasePressODE::Oil;
typedef Oil<typename Region::CalcDensity, typedef Oil<FluidSystem, typename Region::CalcDissolution> ODE;
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 const double T = 273.15 + 20; // standard temperature for now
ODE drho(T,
reg.densityCalculator(), ODE drho(T, reg.dissolutionCalculator(),
reg.dissolutionCalculator(), reg.pvtIdx(), grav);
pu.num_phases, oix, gix, grav);
double z0; double z0;
double p0; double p0;
@ -444,7 +428,8 @@ namespace Opm
else { po_goc = p0; } // GOC *at* datum else { po_goc = p0; } // GOC *at* datum
} }
template <class Grid, template <class FluidSystem,
class Grid,
class Region, class Region,
class CellRange> class CellRange>
void void
@ -457,19 +442,11 @@ namespace Opm
std::vector<double>& press ) std::vector<double>& press )
{ {
using PhasePressODE::Gas; using PhasePressODE::Gas;
typedef Gas<typename Region::CalcDensity, typedef Gas<FluidSystem, typename Region::CalcEvaporation> ODE;
typename Region::CalcEvaporation> ODE;
const PhaseUsage& pu = reg.phaseUsage();
const int gix = PhaseIndex::gas(pu);
const int oix = PhaseIndex::oil(pu);
const double T = 273.15 + 20; // standard temperature for now const double T = 273.15 + 20; // standard temperature for now
ODE drho(T, ODE drho(T, reg.evaporationCalculator(),
reg.densityCalculator(), reg.pvtIdx(), grav);
reg.evaporationCalculator(),
pu.num_phases, gix, oix, grav);
double z0; double z0;
double p0; double p0;
@ -502,7 +479,8 @@ namespace Opm
} }
} // namespace PhasePressure } // namespace PhasePressure
template <class Grid, template <class FluidSystem,
class Grid,
class Region, class Region,
class CellRange> class CellRange>
void void
@ -521,19 +499,19 @@ namespace Opm
if (PhaseUsed::water(pu)) { if (PhaseUsed::water(pu)) {
const int wix = PhaseIndex::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 ]); cells, press[ wix ]);
} }
if (PhaseUsed::oil(pu)) { if (PhaseUsed::oil(pu)) {
const int oix = PhaseIndex::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); press[ oix ], po_woc, po_goc);
} }
if (PhaseUsed::gas(pu)) { if (PhaseUsed::gas(pu)) {
const int gix = PhaseIndex::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 ]); cells, press[ gix ]);
} }
} else if (reg.datum() < reg.zgoc()) { // Datum in gas zone } else if (reg.datum() < reg.zgoc()) { // Datum in gas zone
@ -542,19 +520,19 @@ namespace Opm
if (PhaseUsed::gas(pu)) { if (PhaseUsed::gas(pu)) {
const int gix = PhaseIndex::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 ]); cells, press[ gix ]);
} }
if (PhaseUsed::oil(pu)) { if (PhaseUsed::oil(pu)) {
const int oix = PhaseIndex::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); press[ oix ], po_woc, po_goc);
} }
if (PhaseUsed::water(pu)) { if (PhaseUsed::water(pu)) {
const int wix = PhaseIndex::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 ]); cells, press[ wix ]);
} }
} else { // Datum in oil zone } else { // Datum in oil zone
@ -563,19 +541,19 @@ namespace Opm
if (PhaseUsed::oil(pu)) { if (PhaseUsed::oil(pu)) {
const int oix = PhaseIndex::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); press[ oix ], po_woc, po_goc);
} }
if (PhaseUsed::water(pu)) { if (PhaseUsed::water(pu)) {
const int wix = PhaseIndex::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 ]); cells, press[ wix ]);
} }
if (PhaseUsed::gas(pu)) { if (PhaseUsed::gas(pu)) {
const int gix = PhaseIndex::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 ]); cells, press[ gix ]);
} }
} }
@ -586,7 +564,8 @@ namespace Opm
namespace EQUIL { namespace EQUIL {
template <class Grid, template <class FluidSystem,
class Grid,
class Region, class Region,
class CellRange> class CellRange>
std::vector< std::vector<double> > std::vector< std::vector<double> >
@ -654,7 +633,7 @@ namespace Opm
span[0] = std::min(span[0],zgoc); span[0] = std::min(span[0],zgoc);
span[1] = std::max(span[1],zwoc); 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; return press;
} }
@ -671,12 +650,12 @@ namespace Opm
return std::vector<double>(cells.size(), 273.15 + 20.0); 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> > std::vector< std::vector<double> >
phaseSaturations(const Grid& G, phaseSaturations(const Grid& G,
const Region& reg, const Region& reg,
const CellRange& cells, const CellRange& cells,
BlackoilPropertiesInterface& props, std::shared_ptr<MaterialLawManager> materialLawManager,
const std::vector<double> swat_init, const std::vector<double> swat_init,
std::vector< std::vector<double> >& phase_pressures) 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. 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 water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour]; const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
@ -696,42 +690,45 @@ namespace Opm
std::vector<double>::size_type local_index = 0; std::vector<double>::size_type local_index = 0;
for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) { for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) {
const int cell = *ci; 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 // Find saturations from pressure differences by
// inverting capillary pressure functions. // inverting capillary pressure functions.
double sw = 0.0; double sw = 0.0;
if (water) { if (water) {
if (isConstPc(props,waterpos,cell)){ if (isConstPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,FluidSystem::waterPhaseIdx, cell)){
const double cellDepth = UgGridHelpers::cellCenterDepth(G, const double cellDepth = UgGridHelpers::cellCenterDepth(G,
cell); 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; phase_saturations[waterpos][local_index] = sw;
} }
else{ else{
const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index]; const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
if (swat_init.empty()) { // Invert Pc to find sw 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; phase_saturations[waterpos][local_index] = sw;
} else { // Scale Pc to reflect imposed sw } else { // Scale Pc to reflect imposed sw
sw = swat_init[cell]; sw = swat_init[cell];
props.swatInitScaling(cell, pcov, sw); sw = materialLawManager->applySwatinit(cell, pcov, sw);
phase_saturations[waterpos][local_index] = sw; phase_saturations[waterpos][local_index] = sw;
} }
} }
} }
double sg = 0.0; double sg = 0.0;
if (gas) { if (gas) {
if (isConstPc(props,gaspos,cell)){ if (isConstPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,FluidSystem::gasPhaseIdx,cell)){
const double cellDepth = UgGridHelpers::cellCenterDepth(G, const double cellDepth = UgGridHelpers::cellCenterDepth(G,
cell); 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; phase_saturations[gaspos][local_index] = sg;
} }
else{ else{
// Note that pcog is defined to be (pg - po), not (po - pg). // 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 pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index];
const double increasing = true; // pcog(sg) expected to be increasing function 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; phase_saturations[gaspos][local_index] = sg;
} }
} }
@ -745,55 +742,69 @@ namespace Opm
// Re-scale Pc to reflect imposed sw for vanishing oil phase. // Re-scale Pc to reflect imposed sw for vanishing oil phase.
// This seems consistent with ecl, and fails to honour // This seems consistent with ecl, and fails to honour
// swat_init in case of non-trivial gas-oil cap pressure. // 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; sg = 1.0 - sw;
phase_saturations[waterpos][local_index] = sw; phase_saturations[waterpos][local_index] = sw;
phase_saturations[gaspos][local_index] = sg; phase_saturations[gaspos][local_index] = sg;
// Adjust oil pressure according to gas saturation and cap pressure
double pc[BlackoilPhases::MaxNumPhases]; if ( water ) {
double sat[BlackoilPhases::MaxNumPhases]; fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
sat[waterpos] = sw; }
sat[gaspos] = sg; else {
sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos]; fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
props.capPress(1, sat, &cell, pc, 0); }
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos]; 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; phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
// Adjust phase pressures for max and min saturation ... // Adjust phase pressures for max and min saturation ...
double pc[BlackoilPhases::MaxNumPhases];
double sat[BlackoilPhases::MaxNumPhases];
double threshold_sat = 1.0e-6; 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) { if (water) {
sat[waterpos] = smax[waterpos]; double swu = scaledDrainageInfo.Swu;
sat[oilpos] -= sat[waterpos]; fluidState.setSaturation(FluidSystem::waterPhaseIdx, swu);
so -= swu;
} }
if (gas) { if (gas) {
sat[gaspos] = smax[gaspos]; double sgu = scaledDrainageInfo.Sgu;
sat[oilpos] -= sat[gaspos]; fluidState.setSaturation(FluidSystem::gasPhaseIdx, sgu);
so-= sgu;
} }
if (water && sw > smax[waterpos]-threshold_sat ) { fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
sat[waterpos] = smax[waterpos];
props.capPress(1, sat, &cell, pc, 0); if (water && sw > scaledDrainageInfo.Swu-threshold_sat ) {
phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos]; fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swu);
} else if (gas && sg > smax[gaspos]-threshold_sat) { MaterialLaw::capillaryPressures(pC, matParams, fluidState);
sat[gaspos] = smax[gaspos]; double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx];
props.capPress(1, sat, &cell, pc, 0); phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pcWat;
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos]; } 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) { if (gas && sg < scaledDrainageInfo.Sgl+threshold_sat) {
sat[gaspos] = smin[gaspos]; fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgl);
props.capPress(1, sat, &cell, pc, 0); MaterialLaw::capillaryPressures(pC, matParams, fluidState);
phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos]; 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) { if (water && sw < scaledDrainageInfo.Swl+threshold_sat) {
sat[waterpos] = smin[waterpos]; fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swl);
props.capPress(1, sat, &cell, pc, 0); MaterialLaw::capillaryPressures(pC, matParams, fluidState);
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos]; double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx];
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pcWat;
} }
} }
return phase_saturations; return phase_saturations;
@ -885,9 +896,9 @@ namespace Opm
* \param[in] applySwatInit Make it possible to not apply SWATINIT even if it * \param[in] applySwatInit Make it possible to not apply SWATINIT even if it
* is present in the deck * is present in the deck
*/ */
template<class Grid> template<class MaterialLawManager, class Grid>
void initStateEquil(const Grid& grid, void initStateEquil(const Grid& grid,
BlackoilPropertiesFromDeck& props, std::shared_ptr<MaterialLawManager> materialLawManager,
const Opm::Deck& deck, const Opm::Deck& deck,
const Opm::EclipseState& eclipseState, const Opm::EclipseState& eclipseState,
const double gravity, const double gravity,
@ -896,6 +907,9 @@ namespace Opm
{ {
typedef EQUIL::DeckDependent::InitialStateComputer ISC; typedef EQUIL::DeckDependent::InitialStateComputer ISC;
PhaseUsage pu = phaseUsageFromDeck(deck);
//Check for presence of kw SWATINIT //Check for presence of kw SWATINIT
std::vector<double> swat_init = {}; std::vector<double> swat_init = {};
if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatinit) { if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatinit) {
@ -910,8 +924,7 @@ namespace Opm
} }
} }
ISC isc(props, deck, eclipseState, grid, gravity, swat_init); ISC isc(materialLawManager, pu, deck, eclipseState, grid, gravity, swat_init);
const auto pu = props.phaseUsage();
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid] const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
? pu.phase_pos[BlackoilPhases::Liquid] ? pu.phase_pos[BlackoilPhases::Liquid]
: pu.phase_pos[BlackoilPhases::Aqua]; : pu.phase_pos[BlackoilPhases::Aqua];
@ -920,7 +933,7 @@ namespace Opm
state.gasoilratio() = isc.rs(); state.gasoilratio() = isc.rs();
state.rv() = isc.rv(); state.rv() = isc.rv();
initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state); //initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
} }

374
tests/test_equil.cpp
View File

@ -22,10 +22,15 @@
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h> #include <opm/core/grid/cart_grid.h>
#include <opm/core/grid/GridManager.hpp> #include <opm/core/grid/GridManager.hpp>
#include <opm/core/utility/compressedToCartesian.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp> #include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp> #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.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/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp> #include <opm/parser/eclipse/Parser/ParseContext.hpp>
@ -49,6 +54,15 @@
#include <string> #include <string>
#include <vector> #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) \ #define CHECK(value, expected, reltol) \
{ \ { \
if (std::fabs((expected)) < 1.e-14) \ if (std::fabs((expected)) < 1.e-14) \
@ -57,6 +71,87 @@
BOOST_CHECK_CLOSE((value), (expected), (reltol)); \ 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 () BOOST_AUTO_TEST_SUITE ()
static Opm::EquilRecord mkEquilRecord( double datd, double datp, static Opm::EquilRecord mkEquilRecord( double datd, double datp,
@ -123,35 +218,22 @@ BOOST_AUTO_TEST_CASE (PhasePressure)
std::shared_ptr<UnstructuredGrid> std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid); 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 ); auto record = mkEquilRecord( 0, 1e5, 5, 0, 0, 0 );
Opm::EQUIL::EquilReg<RhoCalc> Opm::PhaseUsage pu = initDefaultFluidSystem();
region(record, calc,
Opm::EQUIL::EquilReg
region(record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage()); 0,
pu);
std::vector<int> cells(G->number_of_cells); std::vector<int> cells(G->number_of_cells);
std::iota(cells.begin(), cells.end(), 0); std::iota(cells.begin(), cells.end(), 0);
const double grav = 10; 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 int first = 0, last = G->number_of_cells - 1;
const double reltol = 1.0e-8; const double reltol = 1.0e-8;
@ -169,46 +251,37 @@ BOOST_AUTO_TEST_CASE (CellSubset)
std::shared_ptr<UnstructuredGrid> std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid); G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::ParameterGroup param; Opm::PhaseUsage pu = initDefaultFluidSystem();
{
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 ), Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) }; mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
Opm::EQUIL::EquilReg<RhoCalc> region[] = Opm::EQUIL::EquilReg region[] =
{ {
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc, Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage()) 0,
pu)
, ,
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc, Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage()) 0,
pu)
, ,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc, Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage()) 0,
pu)
, ,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc, Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage()) 0,
pu)
}; };
const int cdim[] = { 2, 1, 2 }; const int cdim[] = { 2, 1, 2 };
@ -239,7 +312,7 @@ BOOST_AUTO_TEST_CASE (CellSubset)
const int rno = int(r - cells.begin()); const int rno = int(r - cells.begin());
const double grav = 10; const double grav = 10;
const PPress p = 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; PVal::size_type i = 0;
for (std::vector<int>::const_iterator for (std::vector<int>::const_iterator
@ -272,46 +345,36 @@ BOOST_AUTO_TEST_CASE (RegMapping)
std::shared_ptr<UnstructuredGrid> std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid); 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 ), Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) }; mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
Opm::EQUIL::EquilReg<RhoCalc> region[] = Opm::PhaseUsage pu = initDefaultFluidSystem();
Opm::EQUIL::EquilReg region[] =
{ {
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc, Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage()) 0,
pu)
, ,
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc, Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage()) 0,
pu)
, ,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc, Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage()) 0,
pu)
, ,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc, Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(), std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage()) 0,
pu)
}; };
std::vector<int> eqlnum(G->number_of_cells); std::vector<int> eqlnum(G->number_of_cells);
@ -336,7 +399,7 @@ BOOST_AUTO_TEST_CASE (RegMapping)
const int rno = r; const int rno = r;
const double grav = 10; const double grav = 10;
const PPress p = 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; PVal::size_type i = 0;
for (const auto& c : rng) { for (const auto& c : rng) {
@ -366,9 +429,16 @@ BOOST_AUTO_TEST_CASE (DeckAllDead)
Opm::ParseContext parseContext; Opm::ParseContext parseContext;
Opm::Parser parser; Opm::Parser parser;
Opm::Deck deck = parser.parseFile("deadfluids.DATA" , parseContext); Opm::Deck deck = parser.parseFile("deadfluids.DATA" , parseContext);
Opm::EclipseState eclipseState(deck, parseContext); Opm::EclipseState eclipseState(deck, parseContext); // Create material law manager.
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, *grid, false);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, *grid, 10.0); 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(); const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells); BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells);
@ -395,8 +465,20 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
Opm::ParseContext parseContext; Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary.DATA" , parseContext); Opm::Deck deck = parser.parseFile("capillary.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , 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. // Test the capillary inversion for oil-water.
const int cell = 0; const int cell = 0;
const double reltol = 1.0e-7; 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 }; 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()); BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) { 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); 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 }; 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()); BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) { 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); 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 }; const std::vector<double> s = { 0.2, 0.333333333333, 0.6, 0.866666666666, 1.0 };
BOOST_REQUIRE(pc.size() == s.size()); BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) { 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); BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
} }
} }
@ -449,9 +531,16 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
Opm::ParseContext parseContext; Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary.DATA" , parseContext); Opm::Deck deck = parser.parseFile("capillary.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , 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(); const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -490,9 +579,15 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
Opm::ParseContext parseContext; Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary_overlap.DATA" , parseContext); Opm::Deck deck = parser.parseFile("capillary_overlap.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , 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(); const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -553,9 +648,15 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
Opm::ParseContext parseContext; Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("equil_liveoil.DATA" , parseContext); Opm::Deck deck = parser.parseFile("equil_liveoil.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , 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(); const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -633,9 +734,15 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
Opm::ParseContext parseContext; Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("equil_livegas.DATA" , parseContext); Opm::Deck deck = parser.parseFile("equil_livegas.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , 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(); const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -715,9 +822,15 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
Opm::Parser parser; Opm::Parser parser;
Opm::Deck deck = parser.parseFile("equil_rsvd_and_rvvd.DATA", parseContext); Opm::Deck deck = parser.parseFile("equil_rsvd_and_rvvd.DATA", parseContext);
Opm::EclipseState eclipseState(deck , 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(); const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); 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::EclipseState eclipseState(deck , parseContext);
Opm::GridManager gm(eclipseState.getInputGrid()); Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid()); 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); 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 // 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(); 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; std::vector<double> pc_scaled_truth = pc_original;
@ -876,17 +1027,46 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
// apply swatinit // apply swatinit
Opm::BlackoilState state_scaled = state; 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 // don't apply swatinit
Opm::BlackoilState state_unscaled = state; 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 // compute pc
std::vector<double> pc_scaled= state.saturation(); 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(); 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 // test
const double reltol = 1.0e-3; const double reltol = 1.0e-3;