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

Merge pull request #1187 from totto82/removeBlackOilPropsFromInit

Browse Source 
Remove black oil props from init
This commit is contained in:
Atgeirr Flø Rasmussen
2017-11-27 07:26:45 +01:00
committed by GitHub
parent 37421c62da 729d6cf30e
commit 2a88805822
5 changed files with 822 additions and 806 deletions
Download Patch File Download Diff File Expand all files Collapse all files

71
examples/compute_initial_state.cpp
View File

@@ -1,5 +1,6 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2017 IRIS
This file is part of the Open Porous Media project (OPM).
@@ -27,13 +28,17 @@
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/compressedToCartesian.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <boost/filesystem.hpp>
#include <fstream>
@@ -70,12 +75,36 @@ namespace
std::copy(data.begin(), data.end(), std::ostream_iterator<double>(file, "\n"));
}
/// Convert saturations from a vector of individual phase saturation vectors
/// to an interleaved format where all values for a given cell come before all
/// values for the next cell, all in a single vector.
template <class FluidSystem>
void convertSats(std::vector<double>& sat_interleaved, const std::vector< std::vector<double> >& sat, const Opm::PhaseUsage& pu)
{
assert(sat.size() == 3);
const auto nc = sat[0].size();
const auto np = sat_interleaved.size() / nc;
for (size_t c = 0; c < nc; ++c) {
if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int opos = pu.phase_pos[Opm::BlackoilPhases::Liquid];
const std::vector<double>& sat_p = sat[ FluidSystem::oilPhaseIdx];
sat_interleaved[np*c + opos] = sat_p[c];
}
if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int wpos = pu.phase_pos[Opm::BlackoilPhases::Aqua];
const std::vector<double>& sat_p = sat[ FluidSystem::waterPhaseIdx];
sat_interleaved[np*c + wpos] = sat_p[c];
}
if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gpos = pu.phase_pos[Opm::BlackoilPhases::Vapour];
const std::vector<double>& sat_p = sat[ FluidSystem::gasPhaseIdx];
sat_interleaved[np*c + gpos] = sat_p[c];
}
}
}
} // anon namespace
// ----------------- Main program -----------------
int
main(int argc, char** argv)
@@ -94,13 +123,43 @@ try
const double grav = param.getDefault("gravity", unit::gravity);
GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *gm.c_grid();
BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param);
warnIfUnusedParams(param);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
typedef FluidSystems::BlackOil<double> FluidSystem;
// Forward declaring the MaterialLawManager template.
typedef Opm::ThreePhaseMaterialTraits<double,
/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> MaterialTraits;
typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
// Initialisation.
//initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
BlackoilState state( UgGridHelpers::numCells(grid) , UgGridHelpers::numFaces(grid), 3);
initStateEquil(grid, props, deck, eclipseState, grav, state);
FluidSystem::initFromDeck(deck, eclipseState);
PhaseUsage pu = phaseUsageFromDeck(deck);
typedef EQUIL::DeckDependent::InitialStateComputer<FluidSystem> ISC;
ISC isc(materialLawManager, eclipseState, grid, grav);
const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
const int oilpos = FluidSystem::oilPhaseIdx;
const int waterpos = FluidSystem::waterPhaseIdx;
const int ref_phase = oil ? oilpos : waterpos;
state.pressure() = isc.press()[ref_phase];
convertSats<FluidSystem>(state.saturation(), isc.saturation(), pu);
state.gasoilratio() = isc.rs();
state.rv() = isc.rv();
// Output.
const std::string output_dir = param.getDefault<std::string>("output_dir", "output");

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

@@ -1,5 +1,6 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2017 IRIS
This file is part of the Open Porous Media project (OPM).
@@ -20,14 +21,16 @@
#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/utility/RootFinders.hpp>
#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>
@@ -38,36 +41,33 @@ namespace Opm
{
namespace EQUIL {
template <class Props>
class DensityCalculator;
template <>
class DensityCalculator< BlackoilPropertiesInterface >;
namespace Miscibility {
class RsFunction;
class NoMixing;
template <class FluidSystem>
class RsVD;
template <class FluidSystem>
class RsSatAtContact;
}
template <class DensCalc>
class EquilReg;
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
struct PcEq;
inline double satFromPc(const BlackoilPropertiesInterface& props,
template <class FluidSystem, class MaterialLaw, class MaterialLawManager >
inline double satFromPc(const MaterialLawManager& materialLawManager,
const int phase,
const int cell,
const double target_pc,
const bool increasing = false);
struct PcEqSum
inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
const bool increasing = false)
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
const int phase1,
const int phase2,
const int cell,
const double target_pc);
const double target_pc)
} // namespace Equil
} // namespace Opm
@@ -87,71 +87,21 @@ namespace Opm
namespace EQUIL {
template <class Props>
class DensityCalculator;
/**
* Facility for calculating phase densities based on the
* BlackoilPropertiesInterface.
*
* Implements the crucial <CODE>operator()(p,svol)</CODE>
* function that is expected by class EquilReg.
*/
template <>
class DensityCalculator< BlackoilPropertiesInterface > {
public:
/**
* Constructor.
*
* \param[in] props Implementation of the
* BlackoilPropertiesInterface.
*
* \param[in] c Single cell used as a representative cell
* in a PVT region.
*/
DensityCalculator(const BlackoilPropertiesInterface& props,
const int c)
: props_(props)
, c_(1, c)
{
}
/**
* Compute phase densities of all phases at phase point
* given by (pressure, surface volume) tuple.
*
* \param[in] p Fluid pressure.
*
* \param[in] T Temperature.
*
* \param[in] z Surface volumes of all phases.
*
* \return Phase densities at phase point.
*/
std::vector<double>
operator()(const double p,
const double T,
const std::vector<double>& z) const
{
const int np = props_.numPhases();
std::vector<double> A(np * np, 0);
assert (z.size() == std::vector<double>::size_type(np));
double* dAdp = 0;
props_.matrix(1, &p, &T, &z[0], &c_[0], &A[0], dAdp);
std::vector<double> rho(np, 0.0);
props_.density(1, &A[0], &c_[0], &rho[0]);
return rho;
}
private:
const BlackoilPropertiesInterface& props_;
const std::vector<int> c_;
};
typedef Opm::FluidSystems::BlackOil<double> FluidSystemSimple;
// Adjust oil pressure according to gas saturation and cap pressure
typedef Opm::SimpleModularFluidState<double,
/*numPhases=*/3,
/*numComponents=*/3,
FluidSystemSimple,
/*storePressure=*/false,
/*storeTemperature=*/false,
/*storeComposition=*/false,
/*storeFugacity=*/false,
/*storeSaturation=*/true,
/*storeDensity=*/false,
/*storeViscosity=*/false,
/*storeEnthalpy=*/false> SatOnlyFluidState;
/**
* Types and routines relating to phase mixing in
@@ -224,29 +174,23 @@ namespace Opm
* tabulated as a function of depth policy. Data
* typically taken from keyword 'RSVD'.
*/
template <class FluidSystem>
class RsVD : public RsFunction {
public:
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] rs Dissolved gas-oil ratio at @c depth.
*/
RsVD(const BlackoilPropertiesInterface& props,
const int cell,
RsVD(const int pvtRegionIdx,
const std::vector<double>& depth,
const std::vector<double>& rs)
: props_(props)
, cell_(cell)
: pvtRegionIdx_(pvtRegionIdx)
, depth_(depth)
, rs_(rs)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
}
/**
@@ -278,23 +222,13 @@ namespace Opm
}
private:
const BlackoilPropertiesInterface& props_;
const int cell_;
const int pvtRegionIdx_;
std::vector<double> depth_; /**< Depth nodes */
std::vector<double> rs_; /**< Dissolved gas-oil ratio */
double z_[BlackoilPhases::MaxNumPhases];
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRs(const double press, const double temp) const
{
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rs/Bo is in the gas row and oil column of A_.
// 1/Bo is in the oil row and column.
// Recall also that it is stored in column-major order.
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const int np = props_.numPhases();
return A_[np*opos + gpos] / A_[np*opos + opos];
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
}
};
@@ -304,29 +238,23 @@ namespace Opm
* tabulated as a function of depth policy. Data
* typically taken from keyword 'RVVD'.
*/
template <class FluidSystem>
class RvVD : public RsFunction {
public:
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] pvtRegionIdx The pvt region index
* \param[in] depth Depth nodes.
* \param[in] rv Dissolved gas-oil ratio at @c depth.
*/
RvVD(const BlackoilPropertiesInterface& props,
const int cell,
RvVD(const int pvtRegionIdx,
const std::vector<double>& depth,
const std::vector<double>& rv)
: props_(props)
, cell_(cell)
: pvtRegionIdx_(pvtRegionIdx)
, depth_(depth)
, rv_(rv)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
}
/**
@@ -358,23 +286,13 @@ namespace Opm
}
private:
const BlackoilPropertiesInterface& props_;
const int cell_;
const int pvtRegionIdx_;
std::vector<double> depth_; /**< Depth nodes */
std::vector<double> rv_; /**< Vaporized oil-gas ratio */
double z_[BlackoilPhases::MaxNumPhases];
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRv(const double press, const double temp) const
{
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rv/Bg is in the oil row and gas column of A_.
// 1/Bg is in the gas row and column.
// Recall also that it is stored in column-major order.
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const int np = props_.numPhases();
return A_[np*gpos + opos] / A_[np*gpos + gpos];
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
}
};
@@ -393,23 +311,19 @@ namespace Opm
* This should yield Rs-values that are constant below the
* contact, and decreasing above the contact.
*/
template <class FluidSystem>
class RsSatAtContact : public RsFunction {
public:
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] pvtRegionIdx The pvt region index
* \param[in] p_contact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/
RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
: props_(props), cell_(cell)
RsSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact)
: pvtRegionIdx_(pvtRegionIdx)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
rs_sat_contact_ = satRs(p_contact, T_contact);
}
@@ -442,22 +356,12 @@ namespace Opm
}
private:
const BlackoilPropertiesInterface& props_;
const int cell_;
double z_[BlackoilPhases::MaxNumPhases];
const int pvtRegionIdx_;
double rs_sat_contact_;
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRs(const double press, const double temp) const
{
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rs/Bo is in the gas row and oil column of A_.
// 1/Bo is in the oil row and column.
// Recall also that it is stored in column-major order.
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const int np = props_.numPhases();
return A_[np*opos + gpos] / A_[np*opos + opos];
return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
}
};
@@ -476,23 +380,19 @@ namespace Opm
* This should yield Rv-values that are constant below the
* contact, and decreasing above the contact.
*/
template <class FluidSystem>
class RvSatAtContact : public RsFunction {
public:
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] pvtRegionIdx The pvt region index
* \param[in] p_contact oil pressure at the contact
* \param[in] T_contact temperature at the contact
*/
RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
: props_(props), cell_(cell)
RvSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact)
:pvtRegionIdx_(pvtRegionIdx)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
rv_sat_contact_ = satRv(p_contact, T_contact);
}
@@ -525,22 +425,12 @@ namespace Opm
}
private:
const BlackoilPropertiesInterface& props_;
const int cell_;
double z_[BlackoilPhases::MaxNumPhases];
const int pvtRegionIdx_;
double rv_sat_contact_;
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRv(const double press, const double temp) const
{
props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
// Rv/Bg is in the oil row and gas column of A_.
// 1/Bg is in the gas row and column.
// Recall also that it is stored in column-major order.
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const int np = props_.numPhases();
return A_[np*gpos + opos] / A_[np*gpos + gpos];
return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
}
};
@@ -565,36 +455,27 @@ namespace Opm
* that calculates the phase densities of all phases in @c
* svol at fluid pressure @c press.
*/
template <class DensCalc>
class EquilReg {
public:
/**
* Constructor.
*
* \param[in] rec Equilibration data of current region.
* \param[in] density Density calculator of current region.
* \param[in] rs Calculator of dissolved gas-oil ratio.
* \param[in] rv Calculator of vapourised oil-gas ratio.
* \param[in] pu Summary of current active phases.
* \param[in] pvtRegionIdx The pvt region index
*/
EquilReg(const EquilRecord& rec,
const DensCalc& density,
std::shared_ptr<Miscibility::RsFunction> rs,
std::shared_ptr<Miscibility::RsFunction> rv,
const PhaseUsage& pu)
const int pvtIdx)
: rec_ (rec)
, density_(density)
, rs_ (rs)
, rv_ (rv)
, pu_ (pu)
, pvtIdx_ (pvtIdx)
{
}
/**
* Type of density calculator.
*/
typedef DensCalc CalcDensity;
/**
* Type of dissolved gas-oil ratio calculator.
*/
@@ -639,11 +520,6 @@ namespace Opm
*/
double pcgo_goc() const { return this->rec_.gasOilContactCapillaryPressure(); }
/**
* Retrieve phase density calculator of current region.
*/
const CalcDensity&
densityCalculator() const { return this->density_; }
/**
* Retrieve dissolved gas-oil ratio calculator of current
@@ -660,17 +536,16 @@ namespace Opm
evaporationCalculator() const { return *this->rv_; }
/**
* Retrieve active fluid phase summary.
* Retrieve pvtIdx of the region.
*/
const PhaseUsage&
phaseUsage() const { return this->pu_; }
int pvtIdx() const { return this->pvtIdx_; }
private:
EquilRecord rec_; /**< Equilibration data */
DensCalc density_; /**< Density calculator */
std::shared_ptr<Miscibility::RsFunction> rs_; /**< RS calculator */
std::shared_ptr<Miscibility::RsFunction> rv_; /**< RV calculator */
PhaseUsage pu_; /**< Active phase summary */
const int pvtIdx_;
};
@@ -678,54 +553,113 @@ namespace Opm
/// Functor for inverting capillary pressure function.
/// Function represented is
/// f(s) = pc(s) - target_pc
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
struct PcEq
{
PcEq(const BlackoilPropertiesInterface& props,
PcEq(const MaterialLawManager& materialLawManager,
const int phase,
const int cell,
const double target_pc)
: props_(props),
: materialLawManager_(materialLawManager),
phase_(phase),
cell_(cell),
target_pc_(target_pc)
{
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
}
double operator()(double s) const
{
s_[phase_] = s;
props_.capPress(1, s_, &cell_, pc_, 0);
return pc_[phase_] - target_pc_;
const auto& matParams = materialLawManager_.materialLawParams(cell_);
SatOnlyFluidState fluidState;
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
fluidState.setSaturation(phase_, s);
double pc[FluidSystem::numPhases];
std::fill(pc, pc + FluidSystem::numPhases, 0.0);
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_];
return pcPhase - target_pc_;
}
private:
const BlackoilPropertiesInterface& props_;
const MaterialLawManager& materialLawManager_;
const int phase_;
const int cell_;
const double target_pc_;
mutable double s_[BlackoilPhases::MaxNumPhases];
mutable double pc_[BlackoilPhases::MaxNumPhases];
};
template <class FluidSystem, class MaterialLawManager>
double minSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) {
const auto& scaledDrainageInfo =
materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
// Find minimum and maximum saturations.
switch(phase) {
case FluidSystem::waterPhaseIdx :
{
return scaledDrainageInfo.Swl;
break;
}
case FluidSystem::gasPhaseIdx :
{
return scaledDrainageInfo.Sgl;
break;
}
case FluidSystem::oilPhaseIdx :
{
OPM_THROW(std::runtime_error, "Min saturation not implemented for oil phase.");
break;
}
default: OPM_THROW(std::runtime_error, "Unknown phaseIdx .");
}
return -1.0;
}
template <class FluidSystem, class MaterialLawManager>
double maxSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) {
const auto& scaledDrainageInfo =
materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
// Find minimum and maximum saturations.
switch(phase) {
case FluidSystem::waterPhaseIdx :
{
return scaledDrainageInfo.Swu;
break;
}
case FluidSystem::gasPhaseIdx :
{
return scaledDrainageInfo.Sgu;
break;
}
case FluidSystem::oilPhaseIdx :
{
OPM_THROW(std::runtime_error, "Max saturation not implemented for oil phase.");
break;
}
default: OPM_THROW(std::runtime_error, "Unknown phaseIdx .");
}
return -1.0;
}
/// Compute saturation of some phase corresponding to a given
/// capillary pressure.
inline double satFromPc(const BlackoilPropertiesInterface& props,
template <class FluidSystem, class MaterialLaw, class MaterialLawManager >
inline double satFromPc(const MaterialLawManager& materialLawManager,
const int phase,
const int cell,
const double target_pc,
const bool increasing = false)
{
// Find minimum and maximum saturations.
double sminarr[BlackoilPhases::MaxNumPhases];
double smaxarr[BlackoilPhases::MaxNumPhases];
props.satRange(1, &cell, sminarr, smaxarr);
const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
// Find minimum and maximum saturations.
const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
// Create the equation f(s) = pc(s) - target_pc
const PcEq f(props, phase, cell, target_pc);
const PcEq<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase, cell, target_pc);
const double f0 = f(s0);
const double f1 = f(s1);
@@ -747,37 +681,47 @@ namespace Opm
/// Functor for inverting a sum of capillary pressure functions.
/// Function represented is
/// f(s) = pc1(s) + pc2(1 - s) - target_pc
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
struct PcEqSum
{
PcEqSum(const BlackoilPropertiesInterface& props,
PcEqSum(const MaterialLawManager& materialLawManager,
const int phase1,
const int phase2,
const int cell,
const double target_pc)
: props_(props),
: materialLawManager_(materialLawManager),
phase1_(phase1),
phase2_(phase2),
cell_(cell),
target_pc_(target_pc)
{
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
}
double operator()(double s) const
{
s_[phase1_] = s;
s_[phase2_] = 1.0 - s;
props_.capPress(1, s_, &cell_, pc_, 0);
return pc_[phase1_] + pc_[phase2_] - target_pc_;
const auto& matParams = materialLawManager_.materialLawParams(cell_);
SatOnlyFluidState fluidState;
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
fluidState.setSaturation(phase1_, s);
fluidState.setSaturation(phase2_, 1.0 - s);
double pc[FluidSystem::numPhases];
std::fill(pc, pc + FluidSystem::numPhases, 0.0);
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_];
double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_];
return pc1 + pc2 - target_pc_;
}
private:
const BlackoilPropertiesInterface& props_;
const MaterialLawManager& materialLawManager_;
const int phase1_;
const int phase2_;
const int cell_;
const double target_pc_;
mutable double s_[BlackoilPhases::MaxNumPhases];
mutable double pc_[BlackoilPhases::MaxNumPhases];
};
@@ -786,21 +730,19 @@ namespace Opm
/// Compute saturation of some phase corresponding to a given
/// capillary pressure, where the capillary pressure function
/// is given as a sum of two other functions.
inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
const int phase1,
const int phase2,
const int cell,
const double target_pc)
{
// Find minimum and maximum saturations.
double sminarr[BlackoilPhases::MaxNumPhases];
double smaxarr[BlackoilPhases::MaxNumPhases];
props.satRange(1, &cell, sminarr, smaxarr);
const double smin = sminarr[phase1];
const double smax = smaxarr[phase1];
const double smin = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
const double smax = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
// Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
const PcEqSum f(props, phase1, phase2, cell, target_pc);
const PcEqSum<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase1, phase2, cell, target_pc);
const double f0 = f(smin);
const double f1 = f(smax);
if (f0 <= 0.0) {
@@ -818,19 +760,16 @@ namespace Opm
}
/// Compute saturation from depth. Used for constant capillary pressure function
inline double satFromDepth(const BlackoilPropertiesInterface& props,
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
inline double satFromDepth(const MaterialLawManager& materialLawManager,
const double cellDepth,
const double contactDepth,
const int phase,
const int cell,
const bool increasing = false)
{
// Find minimum and maximum saturations.
double sminarr[BlackoilPhases::MaxNumPhases];
double smaxarr[BlackoilPhases::MaxNumPhases];
props.satRange(1, &cell, sminarr, smaxarr);
const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
if (cellDepth < contactDepth){
return s0;
@@ -841,19 +780,15 @@ namespace Opm
}
/// Return true if capillary pressure function is constant
inline bool isConstPc(const BlackoilPropertiesInterface& props,
template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
inline bool isConstPc(const MaterialLawManager& materialLawManager,
const int phase,
const int cell)
{
// Find minimum and maximum saturations.
double sminarr[BlackoilPhases::MaxNumPhases];
double smaxarr[BlackoilPhases::MaxNumPhases];
props.satRange(1, &cell, sminarr, smaxarr);
// Create the equation f(s) = pc(s);
const PcEq f(props, phase, cell, 0);
const double f0 = f(sminarr[phase]);
const double f1 = f(smaxarr[phase]);
const PcEq<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase, cell, 0);
const double f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
const double f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
}

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

@@ -2,6 +2,7 @@
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
Copyright 2017 IRIS
This file is part of the Open Porous Media project (OPM).
@@ -24,10 +25,9 @@
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/simulator/EquilibrationHelpers.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/utility/extractPvtTableIndex.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
@@ -38,6 +38,11 @@
#include <opm/parser/eclipse/EclipseState/Tables/RsvdTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <array>
#include <cassert>
@@ -54,32 +59,6 @@ struct UnstructuredGrid;
namespace Opm
{
/**
* Compute initial state by an equilibration procedure.
*
* The following state fields are modified:
* pressure(),
* saturation(),
* surfacevol(),
* gasoilratio(),
* rv().
*
* \param[in] grid Grid.
* \param[in] props Property object, pvt and capillary properties are used.
* \param[in] deck Simulation deck, used to obtain EQUIL and related data.
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
* \param[in] applySwatInit Make it possible to not apply SWATINIT even if it
* is present in the deck
*/
template<class Grid>
void initStateEquil(const Grid& grid,
const BlackoilPropertiesInterface& props,
const Opm::Deck& deck,
const Opm::EclipseState& eclipseState,
const double gravity,
BlackoilState& state,
bool applySwatInit = true);
/**
* Types and routines that collectively implement a basic
@@ -128,7 +107,7 @@ namespace Opm
* of pressure values in each cell in the current
* equilibration region.
*/
template <class Grid, class Region, class CellRange>
template <class Grid, class Region, class CellRange, class FluidSystem>
std::vector< std::vector<double> >
phasePressures(const Grid& G,
const Region& reg,
@@ -140,6 +119,11 @@ namespace Opm
/**
* Compute initial phase saturations by means of equilibration.
*
* \tparam FluidSystem The FluidSystem from opm-material
* Must be initialized before used.
*
* \tparam Grid Type of the grid
*
* \tparam Region Type of an equilibration region information
* base. Typically an instance of the EquilReg
* class template.
@@ -151,22 +135,27 @@ namespace Opm
* as well as provide an inner type,
* const_iterator, to traverse the range.
*
* \tparam MaterialLawManager The MaterialLawManager from opm-material
*
* \param[in] G Grid.
* \param[in] reg Current equilibration region.
* \param[in] cells Range that spans the cells of the current
* equilibration region.
* \param[in] props Property object, needed for capillary functions.
* \param[in] materialLawManager The MaterialLawManager from opm-material
* \param[in] swat_init A vector of initial water saturations.
* The capillary pressure is scaled to fit these values
* \param[in] phase_pressures Phase pressures, one vector for each active phase,
* of pressure values in each cell in the current
* equilibration region.
* \return Phase saturations, one vector for each phase, each containing
* one saturation value per cell in the region.
*/
template <class Grid, class Region, class CellRange>
template <class FluidSystem, class Grid, class Region, class CellRange, class MaterialLawManager>
std::vector< std::vector<double> >
phaseSaturations(const Grid& grid,
const Region& reg,
const CellRange& cells,
BlackoilPropertiesFromDeck& props,
MaterialLawManager& materialLawManager,
const std::vector<double> swat_init,
std::vector< std::vector<double> >& phase_pressures);
@@ -216,12 +205,11 @@ namespace Opm
template<class Grid>
inline
std::vector<int>
equilnum(const Opm::Deck& deck,
const Opm::EclipseState& eclipseState,
equilnum(const Opm::EclipseState& eclipseState,
const Grid& G )
{
std::vector<int> eqlnum;
if (deck.hasKeyword("EQLNUM")) {
if (eclipseState.get3DProperties().hasDeckIntGridProperty("EQLNUM")) {
const int nc = UgGridHelpers::numCells(G);
eqlnum.resize(nc);
const std::vector<int>& e =
@@ -241,43 +229,53 @@ namespace Opm
return eqlnum;
}
template<class FluidSystem>
class InitialStateComputer {
public:
template<class Grid>
InitialStateComputer(BlackoilPropertiesInterface& props,
const Opm::Deck& deck,
template<class MaterialLawManager, class Grid>
InitialStateComputer(MaterialLawManager& materialLawManager,
const Opm::EclipseState& eclipseState,
const Grid& G ,
const double grav = unit::gravity,
const std::vector<double>& swat_init = {}
const bool applySwatInit = true
)
: pp_(props.numPhases(),
: pp_(FluidSystem::numPhases,
std::vector<double>(UgGridHelpers::numCells(G))),
sat_(props.numPhases(),
sat_(FluidSystem::numPhases,
std::vector<double>(UgGridHelpers::numCells(G))),
rs_(UgGridHelpers::numCells(G)),
rv_(UgGridHelpers::numCells(G)),
swat_init_(swat_init)
rv_(UgGridHelpers::numCells(G))
{
//Check for presence of kw SWATINIT
if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatInit) {
const std::vector<double>& swat_init_ecl = eclipseState.
get3DProperties().getDoubleGridProperty("SWATINIT").getData();
const int nc = UgGridHelpers::numCells(G);
swat_init_.resize(nc);
const int* gc = UgGridHelpers::globalCell(G);
for (int c = 0; c < nc; ++c) {
const int deck_pos = (gc == NULL) ? c : gc[c];
swat_init_[c] = swat_init_ecl[deck_pos];
}
}
// Get the equilibration records.
const std::vector<EquilRecord> rec = getEquil(eclipseState);
const auto& tables = eclipseState.getTableManager();
// Create (inverse) region mapping.
const RegionMapping<> eqlmap(equilnum(deck, eclipseState, G));
const RegionMapping<> eqlmap(equilnum(eclipseState, G));
setRegionPvtIdx(G, eclipseState, eqlmap);
// Create Rs functions.
rs_func_.reserve(rec.size());
if (deck.hasKeyword("DISGAS")) {
if (FluidSystem::enableDissolvedGas()) {
const TableContainer& rsvdTables = tables.getRsvdTables();
for (size_t i = 0; i < rec.size(); ++i) {
if (eqlmap.cells(i).empty())
{
rs_func_.push_back(std::shared_ptr<Miscibility::RsVD>());
rs_func_.push_back(std::shared_ptr<Miscibility::RsVD<FluidSystem>>());
continue;
}
const int cell = *(eqlmap.cells(i).begin());
const int pvtIdx = regionPvtIdx_[i];
if (!rec[i].liveOilInitConstantRs()) {
if (rsvdTables.size() <= 0 ) {
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table not available.");
@@ -285,8 +283,7 @@ namespace Opm
const RsvdTable& rsvdTable = rsvdTables.getTable<RsvdTable>(i);
std::vector<double> depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> rsColumn = rsvdTable.getColumn("RS").vectorCopy();
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props,
cell,
rs_func_.push_back(std::make_shared<Miscibility::RsVD<FluidSystem>>(pvtIdx,
depthColumn , rsColumn));
} else {
if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) {
@@ -297,7 +294,7 @@ namespace Opm
}
const double p_contact = rec[i].datumDepthPressure();
const double T_contact = 273.15 + 20; // standard temperature for now
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact, T_contact));
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact<FluidSystem>>(pvtIdx, p_contact, T_contact));
}
}
} else {
@@ -307,15 +304,15 @@ namespace Opm
}
rv_func_.reserve(rec.size());
if (deck.hasKeyword("VAPOIL")) {
if (FluidSystem::enableVaporizedOil()) {
const TableContainer& rvvdTables = tables.getRvvdTables();
for (size_t i = 0; i < rec.size(); ++i) {
if (eqlmap.cells(i).empty())
{
rv_func_.push_back(std::shared_ptr<Miscibility::RvVD>());
rv_func_.push_back(std::shared_ptr<Miscibility::RvVD<FluidSystem>>());
continue;
}
const int cell = *(eqlmap.cells(i).begin());
const int pvtIdx = regionPvtIdx_[i];
if (!rec[i].wetGasInitConstantRv()) {
if (rvvdTables.size() <= 0) {
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table not available.");
@@ -324,9 +321,7 @@ namespace Opm
const RvvdTable& rvvdTable = rvvdTables.getTable<RvvdTable>(i);
std::vector<double> depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy();
std::vector<double> rvColumn = rvvdTable.getColumn("RV").vectorCopy();
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props,
cell,
rv_func_.push_back(std::make_shared<Miscibility::RvVD<FluidSystem>>(pvtIdx,
depthColumn , rvColumn));
} else {
@@ -338,7 +333,7 @@ namespace Opm
}
const double p_contact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure();
const double T_contact = 273.15 + 20; // standard temperature for now
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact, T_contact));
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact<FluidSystem>>(pvtIdx ,p_contact, T_contact));
}
}
} else {
@@ -348,7 +343,7 @@ namespace Opm
}
// Compute pressures, saturations, rs and rv factors.
calcPressSatRsRv(eqlmap, rec, props, G, grav);
calcPressSatRsRv(eqlmap, rec, materialLawManager, G, grav);
// Modify oil pressure in no-oil regions so that the pressures of present phases can
// be recovered from the oil pressure and capillary relations.
@@ -363,23 +358,34 @@ namespace Opm
const Vec& rv() const { return rv_; }
private:
typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
typedef EquilReg<RhoCalc> EqReg;
typedef EquilReg EqReg;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
std::vector<int> regionPvtIdx_;
PVec pp_;
PVec sat_;
Vec rs_;
Vec rv_;
Vec swat_init_;
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 RMap, class MaterialLawManager, class Grid>
void
calcPressSatRsRv(const RMap& reg ,
const std::vector< EquilRecord >& rec ,
Opm::BlackoilPropertiesInterface& props,
MaterialLawManager& materialLawManager,
const Grid& G ,
const double grav)
{
@@ -391,27 +397,23 @@ namespace Opm
+ " has no active cells");
continue;
}
const int repcell = *cells.begin();
const RhoCalc calc(props, repcell);
const EqReg eqreg(rec[r], calc,
rs_func_[r], rv_func_[r],
props.phaseUsage());
const EqReg eqreg(rec[r], rs_func_[r], rv_func_[r], regionPvtIdx_[r]);
PVec pressures = phasePressures(G, eqreg, cells, grav);
PVec pressures = phasePressures<FluidSystem>(G, eqreg, cells, grav);
const std::vector<double>& temp = temperature(G, eqreg, cells);
const PVec sat = phaseSaturations<FluidSystem>(G, eqreg, cells, materialLawManager, swat_init_, pressures);
const PVec sat = phaseSaturations(G, eqreg, cells, props, swat_init_, pressures);
const int np = props.numPhases();
const int np = FluidSystem::numPhases;
for (int p = 0; p < np; ++p) {
copyFromRegion(pressures[p], cells, pp_[p]);
copyFromRegion(sat[p], cells, sat_[p]);
}
if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
if (oil && gas) {
const int oilpos = FluidSystem::oilPhaseIdx;
const int gaspos = FluidSystem::gasPhaseIdx;
const Vec rs_vals = computeRs(G, cells, pressures[oilpos], temp, *(rs_func_[r]), sat[gaspos]);
const Vec rv_vals = computeRs(G, cells, pressures[gaspos], temp, *(rv_func_[r]), sat[oilpos]);
copyFromRegion(rs_vals, cells, rs_);

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

@@ -2,6 +2,7 @@
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 NTNU
Copyright 2017 IRIS
This file is part of the Open Porous Media project (OPM).
@@ -24,8 +25,8 @@
#include <opm/core/grid.h>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <cassert>
#include <cmath>
@@ -106,21 +107,16 @@ namespace Opm
};
namespace PhasePressODE {
template <class Density>
template <class FluidSystem>
class Water {
public:
Water(const double temp,
const Density& rho,
const int np,
const int ix,
const int pvtRegionIdx,
const double norm_grav)
: temp_(temp)
, rho_(rho)
, svol_(np, 0)
, ix_(ix)
, pvtRegionIdx_(pvtRegionIdx)
, g_(norm_grav)
{
svol_[ix_] = 1.0;
}
double
@@ -132,39 +128,30 @@ namespace Opm
private:
const double temp_;
const Density& rho_;
std::vector<double> svol_;
const int ix_;
const int pvtRegionIdx_;
const double g_;
double
density(const double press) const
{
const std::vector<double>& rho = rho_(press, temp_, svol_);
return rho[ix_];
double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
return rho;
}
};
template <class Density, class RS>
template <class FluidSystem, class RS>
class Oil {
public:
Oil(const double temp,
const Density& rho,
const RS& rs,
const int np,
const int oix,
const int gix,
const int pvtRegionIdx,
const double norm_grav)
: temp_(temp)
, rho_(rho)
, rs_(rs)
, svol_(np, 0)
, oix_(oix)
, gix_(gix)
, pvtRegionIdx_(pvtRegionIdx)
, g_(norm_grav)
{
svol_[oix_] = 1.0;
}
double
@@ -176,45 +163,42 @@ namespace Opm
private:
const double temp_;
const Density& rho_;
const RS& rs_;
mutable std::vector<double> svol_;
const int oix_;
const int gix_;
const int pvtRegionIdx_;
const double g_;
double
density(const double depth,
const double press) const
{
if (gix_ >= 0) {
svol_[gix_] = rs_(depth, press, temp_);
double rs = rs_(depth, press, temp_);
double bOil = 0.0;
if ( !FluidSystem::enableDissolvedGas() || rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press) ) {
bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
} else {
bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rs);
}
double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
if (FluidSystem::enableDissolvedGas()) {
rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
}
const std::vector<double>& rho = rho_(press, temp_, svol_);
return rho[oix_];
return rho;
}
};
template <class Density, class RV>
template <class FluidSystem, class RV>
class Gas {
public:
Gas(const double temp,
const Density& rho,
const RV& rv,
const int np,
const int gix,
const int oix,
const int pvtRegionIdx,
const double norm_grav)
: temp_(temp)
, rho_(rho)
, rv_(rv)
, svol_(np, 0)
, gix_(gix)
, oix_(oix)
, pvtRegionIdx_(pvtRegionIdx)
, g_(norm_grav)
{
svol_[gix_] = 1.0;
}
double
@@ -226,81 +210,31 @@ namespace Opm
private:
const double temp_;
const Density& rho_;
const RV& rv_;
mutable std::vector<double> svol_;
const int gix_;
const int oix_;
const int pvtRegionIdx_;
const double g_;
double
density(const double depth,
const double press) const
{
if (oix_ >= 0) {
svol_[oix_] = rv_(depth, press, temp_);
double rv = rv_(depth, press, temp_);
double bGas = 0.0;
if ( !FluidSystem::enableVaporizedOil() || rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press) ) {
bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
} else {
bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv);
}
double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
if (FluidSystem::enableVaporizedOil()) {
rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
}
const std::vector<double>& rho = rho_(press, temp_, svol_);
return rho[gix_];
return rho;
}
};
} // namespace PhasePressODE
namespace PhaseUsed {
inline bool
water(const PhaseUsage& pu)
{
return bool(pu.phase_used[ Opm::BlackoilPhases::Aqua ]);
}
inline bool
oil(const PhaseUsage& pu)
{
return bool(pu.phase_used[ Opm::BlackoilPhases::Liquid ]);
}
inline bool
gas(const PhaseUsage& pu)
{
return bool(pu.phase_used[ Opm::BlackoilPhases::Vapour ]);
}
} // namespace PhaseUsed
namespace PhaseIndex {
inline int
water(const PhaseUsage& pu)
{
int i = -1;
if (PhaseUsed::water(pu)) {
i = pu.phase_pos[ Opm::BlackoilPhases::Aqua ];
}
return i;
}
inline int
oil(const PhaseUsage& pu)
{
int i = -1;
if (PhaseUsed::oil(pu)) {
i = pu.phase_pos[ Opm::BlackoilPhases::Liquid ];
}
return i;
}
inline int
gas(const PhaseUsage& pu)
{
int i = -1;
if (PhaseUsed::gas(pu)) {
i = pu.phase_pos[ Opm::BlackoilPhases::Vapour ];
}
return i;
}
} // namespace PhaseIndex
namespace PhasePressure {
template <class Grid,
@@ -328,9 +262,10 @@ namespace Opm
}
}
template <class Grid,
class Region,
class CellRange>
template < class FluidSystem,
class Grid,
class Region,
class CellRange>
void
water(const Grid& G ,
const Region& reg ,
@@ -341,13 +276,10 @@ namespace Opm
std::vector<double>& press )
{
using PhasePressODE::Water;
typedef Water<typename Region::CalcDensity> ODE;
typedef Water<FluidSystem> ODE;
const PhaseUsage& pu = reg.phaseUsage();
const int wix = PhaseIndex::water(pu);
const double T = 273.15 + 20; // standard temperature for now
ODE drho(T, reg.densityCalculator(), pu.num_phases, wix, grav);
ODE drho(T, reg.pvtIdx() , grav);
double z0;
double p0;
@@ -379,7 +311,8 @@ namespace Opm
}
}
template <class Grid,
template <class FluidSystem,
class Grid,
class Region,
class CellRange>
void
@@ -393,18 +326,11 @@ namespace Opm
double& po_goc)
{
using PhasePressODE::Oil;
typedef Oil<typename Region::CalcDensity,
typename Region::CalcDissolution> ODE;
typedef Oil<FluidSystem, typename Region::CalcDissolution> ODE;
const PhaseUsage& pu = reg.phaseUsage();
const int oix = PhaseIndex::oil(pu);
const int gix = PhaseIndex::gas(pu);
const double T = 273.15 + 20; // standard temperature for now
ODE drho(T,
reg.densityCalculator(),
reg.dissolutionCalculator(),
pu.num_phases, oix, gix, grav);
ODE drho(T, reg.dissolutionCalculator(),
reg.pvtIdx(), grav);
double z0;
double p0;
@@ -444,7 +370,8 @@ namespace Opm
else { po_goc = p0; } // GOC *at* datum
}
template <class Grid,
template <class FluidSystem,
class Grid,
class Region,
class CellRange>
void
@@ -457,19 +384,11 @@ namespace Opm
std::vector<double>& press )
{
using PhasePressODE::Gas;
typedef Gas<typename Region::CalcDensity,
typename Region::CalcEvaporation> ODE;
const PhaseUsage& pu = reg.phaseUsage();
const int gix = PhaseIndex::gas(pu);
const int oix = PhaseIndex::oil(pu);
typedef Gas<FluidSystem, typename Region::CalcEvaporation> ODE;
const double T = 273.15 + 20; // standard temperature for now
ODE drho(T,
reg.densityCalculator(),
reg.evaporationCalculator(),
pu.num_phases, gix, oix, grav);
ODE drho(T, reg.evaporationCalculator(),
reg.pvtIdx(), grav);
double z0;
double p0;
@@ -502,7 +421,8 @@ namespace Opm
}
} // namespace PhasePressure
template <class Grid,
template <class FluidSystem,
class Grid,
class Region,
class CellRange>
void
@@ -513,70 +433,66 @@ namespace Opm
const CellRange& cells,
std::vector< std::vector<double> >& press)
{
const PhaseUsage& pu = reg.phaseUsage();
const bool water = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
const int oilpos = FluidSystem::oilPhaseIdx;
const int waterpos = FluidSystem::waterPhaseIdx;
const int gaspos = FluidSystem::gasPhaseIdx;
if (reg.datum() > reg.zwoc()) { // Datum in water zone
double po_woc = -1;
double po_goc = -1;
if (PhaseUsed::water(pu)) {
const int wix = PhaseIndex::water(pu);
PhasePressure::water(G, reg, span, grav, po_woc,
cells, press[ wix ]);
if (water) {
PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
cells, press[ waterpos ]);
}
if (PhaseUsed::oil(pu)) {
const int oix = PhaseIndex::oil(pu);
PhasePressure::oil(G, reg, span, grav, cells,
press[ oix ], po_woc, po_goc);
if (oil) {
PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
press[ oilpos ], po_woc, po_goc);
}
if (PhaseUsed::gas(pu)) {
const int gix = PhaseIndex::gas(pu);
PhasePressure::gas(G, reg, span, grav, po_goc,
cells, press[ gix ]);
if (gas) {
PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
cells, press[ gaspos ]);
}
} else if (reg.datum() < reg.zgoc()) { // Datum in gas zone
double po_woc = -1;
double po_goc = -1;
if (PhaseUsed::gas(pu)) {
const int gix = PhaseIndex::gas(pu);
PhasePressure::gas(G, reg, span, grav, po_goc,
cells, press[ gix ]);
if (gas) {
PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
cells, press[ gaspos ]);
}
if (PhaseUsed::oil(pu)) {
const int oix = PhaseIndex::oil(pu);
PhasePressure::oil(G, reg, span, grav, cells,
press[ oix ], po_woc, po_goc);
if (oil) {
PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
press[ oilpos ], po_woc, po_goc);
}
if (PhaseUsed::water(pu)) {
const int wix = PhaseIndex::water(pu);
PhasePressure::water(G, reg, span, grav, po_woc,
cells, press[ wix ]);
if (water) {
PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
cells, press[ waterpos ]);
}
} else { // Datum in oil zone
double po_woc = -1;
double po_goc = -1;
if (PhaseUsed::oil(pu)) {
const int oix = PhaseIndex::oil(pu);
PhasePressure::oil(G, reg, span, grav, cells,
press[ oix ], po_woc, po_goc);
if (oil) {
PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
press[ oilpos ], po_woc, po_goc);
}
if (PhaseUsed::water(pu)) {
const int wix = PhaseIndex::water(pu);
PhasePressure::water(G, reg, span, grav, po_woc,
cells, press[ wix ]);
if (water) {
PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
cells, press[ waterpos ]);
}
if (PhaseUsed::gas(pu)) {
const int gix = PhaseIndex::gas(pu);
PhasePressure::gas(G, reg, span, grav, po_goc,
cells, press[ gix ]);
if (gas) {
PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
cells, press[ gaspos ]);
}
}
}
@@ -586,7 +502,8 @@ namespace Opm
namespace EQUIL {
template <class Grid,
template <class FluidSystem,
class Grid,
class Region,
class CellRange>
std::vector< std::vector<double> >
@@ -642,7 +559,7 @@ namespace Opm
}
}
}
const int np = reg.phaseUsage().num_phases;
const int np = FluidSystem::numPhases; //reg.phaseUsage().num_phases;
typedef std::vector<double> pval;
std::vector<pval> press(np, pval(ncell, 0.0));
@@ -654,7 +571,7 @@ namespace Opm
span[0] = std::min(span[0],zgoc);
span[1] = std::max(span[1],zwoc);
Details::equilibrateOWG(G, reg, grav, span, cells, press);
Details::equilibrateOWG<FluidSystem>(G, reg, grav, span, cells, press);
return press;
}
@@ -671,67 +588,85 @@ namespace Opm
return std::vector<double>(cells.size(), 273.15 + 20.0);
}
template <class Grid, class Region, class CellRange>
template <class FluidSystem, class Grid, class Region, class CellRange, class MaterialLawManager>
std::vector< std::vector<double> >
phaseSaturations(const Grid& G,
const Region& reg,
const CellRange& cells,
BlackoilPropertiesInterface& props,
MaterialLawManager& materialLawManager,
const std::vector<double> swat_init,
std::vector< std::vector<double> >& phase_pressures)
{
if (!reg.phaseUsage().phase_used[BlackoilPhases::Liquid]) {
if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
OPM_THROW(std::runtime_error, "Cannot initialise: not handling water-gas cases.");
}
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 };
const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua];
const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour];
// 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 = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
const int oilpos = FluidSystem::oilPhaseIdx;
const int waterpos = FluidSystem::waterPhaseIdx;
const int gaspos = FluidSystem::gasPhaseIdx;
std::vector<double>::size_type local_index = 0;
for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) {
const int cell = *ci;
props.satRange(1, &cell, smin, smax);
const auto& scaledDrainageInfo =
materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
const auto& matParams = materialLawManager.materialLawParams(cell);
// Find saturations from pressure differences by
// inverting capillary pressure functions.
double sw = 0.0;
if (water) {
if (isConstPc(props,waterpos,cell)){
if (isConstPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,FluidSystem::waterPhaseIdx, cell)){
const double cellDepth = UgGridHelpers::cellCenterDepth(G,
cell);
sw = satFromDepth(props,cellDepth,reg.zwoc(),waterpos,cell,false);
sw = satFromDepth<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,cellDepth,reg.zwoc(),waterpos,cell,false);
phase_saturations[waterpos][local_index] = sw;
}
else{
const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
if (swat_init.empty()) { // Invert Pc to find sw
sw = satFromPc(props, waterpos, cell, pcov);
sw = satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, waterpos, cell, pcov);
phase_saturations[waterpos][local_index] = sw;
} else { // Scale Pc to reflect imposed sw
sw = swat_init[cell];
props.swatInitScaling(cell, pcov, sw);
sw = materialLawManager.applySwatinit(cell, pcov, sw);
phase_saturations[waterpos][local_index] = sw;
}
}
}
double sg = 0.0;
if (gas) {
if (isConstPc(props,gaspos,cell)){
if (isConstPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,FluidSystem::gasPhaseIdx,cell)){
const double cellDepth = UgGridHelpers::cellCenterDepth(G,
cell);
sg = satFromDepth(props,cellDepth,reg.zgoc(),gaspos,cell,true);
sg = satFromDepth<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,cellDepth,reg.zgoc(),gaspos,cell,true);
phase_saturations[gaspos][local_index] = sg;
}
else{
// Note that pcog is defined to be (pg - po), not (po - pg).
const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index];
const double increasing = true; // pcog(sg) expected to be increasing function
sg = satFromPc(props, gaspos, cell, pcog, increasing);
sg = satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, gaspos, cell, pcog, increasing);
phase_saturations[gaspos][local_index] = sg;
}
}
@@ -745,56 +680,68 @@ namespace Opm
// Re-scale Pc to reflect imposed sw for vanishing oil phase.
// This seems consistent with ecl, and fails to honour
// swat_init in case of non-trivial gas-oil cap pressure.
props.swatInitScaling(cell, pcgw, sw);
sw = materialLawManager.applySwatinit(cell, pcgw, sw);
}
sw = satFromSumOfPcs(props, waterpos, gaspos, cell, pcgw);
sw = satFromSumOfPcs<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, waterpos, gaspos, cell, pcgw);
sg = 1.0 - sw;
phase_saturations[waterpos][local_index] = sw;
phase_saturations[gaspos][local_index] = sg;
// Adjust oil pressure according to gas saturation and cap pressure
double pc[BlackoilPhases::MaxNumPhases];
double sat[BlackoilPhases::MaxNumPhases];
sat[waterpos] = sw;
sat[gaspos] = sg;
sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos];
props.capPress(1, sat, &cell, pc, 0);
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
if ( water ) {
fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
}
else {
fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
}
fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0 - sw - sg);
fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg);
double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 };
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx];
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pcGas;
}
phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
// Adjust phase pressures for max and min saturation ...
double pc[BlackoilPhases::MaxNumPhases];
double sat[BlackoilPhases::MaxNumPhases];
double threshold_sat = 1.0e-6;
sat[oilpos] = 1.0;
if (water) {
sat[waterpos] = smax[waterpos];
sat[oilpos] -= sat[waterpos];
}
if (gas) {
sat[gaspos] = smax[gaspos];
sat[oilpos] -= sat[gaspos];
}
if (water && sw > smax[waterpos]-threshold_sat ) {
sat[waterpos] = smax[waterpos];
props.capPress(1, sat, &cell, pc, 0);
phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
} else if (gas && sg > smax[gaspos]-threshold_sat) {
sat[gaspos] = smax[gaspos];
props.capPress(1, sat, &cell, pc, 0);
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
}
if (gas && sg < smin[gaspos]+threshold_sat) {
sat[gaspos] = smin[gaspos];
props.capPress(1, sat, &cell, pc, 0);
phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
}
if (water && sw < smin[waterpos]+threshold_sat) {
sat[waterpos] = smin[waterpos];
props.capPress(1, sat, &cell, pc, 0);
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
}
double so = 1.0;
double pC[FluidSystem::numPhases] = { 0.0, 0.0, 0.0 };
if (water) {
double swu = scaledDrainageInfo.Swu;
fluidState.setSaturation(FluidSystem::waterPhaseIdx, swu);
so -= swu;
}
if (gas) {
double sgu = scaledDrainageInfo.Sgu;
fluidState.setSaturation(FluidSystem::gasPhaseIdx, sgu);
so-= sgu;
}
fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
if (water && sw > scaledDrainageInfo.Swu-threshold_sat ) {
fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swu);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx];
phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pcWat;
} else if (gas && sg > scaledDrainageInfo.Sgu-threshold_sat) {
fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgu);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx];
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pcGas;
}
if (gas && sg < scaledDrainageInfo.Sgl+threshold_sat) {
fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgl);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx];
phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pcGas;
}
if (water && sw < scaledDrainageInfo.Swl+threshold_sat) {
fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swl);
MaterialLaw::capillaryPressures(pC, matParams, fluidState);
double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx];
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pcWat;
}
}
return phase_saturations;
}
@@ -839,92 +786,6 @@ namespace Opm
} // namespace Equil
namespace Details
{
/// Convert saturations from a vector of individual phase saturation vectors
/// to an interleaved format where all values for a given cell come before all
/// values for the next cell, all in a single vector.
inline std::vector<double>
convertSats(const std::vector< std::vector<double> >& sat)
{
const auto np = sat.size();
const auto nc = sat[0].size();
std::vector<double> s(np * nc);
for (decltype(sat.size()) p = 0; p < np; ++p) {
const auto& sat_p = sat[p];
double* sp = & s[0*nc + p];
for (decltype(sat[0].size()) c = 0;
c < nc; ++c, sp += np)
{
*sp = sat_p[c];
}
}
return s;
}
} // namespace Details
/**
* 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,
BlackoilPropertiesFromDeck& props,
const Opm::Deck& deck,
const Opm::EclipseState& eclipseState,
const double gravity,
BlackoilState& state,
bool applySwatinit = true)
{
typedef EQUIL::DeckDependent::InitialStateComputer ISC;
//Check for presence of kw SWATINIT
std::vector<double> swat_init = {};
if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatinit) {
const std::vector<double>& swat_init_ecl = eclipseState.
get3DProperties().getDoubleGridProperty("SWATINIT").getData();
const int nc = UgGridHelpers::numCells(grid);
swat_init.resize(nc);
const int* gc = UgGridHelpers::globalCell(grid);
for (int c = 0; c < nc; ++c) {
const int deck_pos = (gc == NULL) ? c : gc[c];
swat_init[c] = swat_init_ecl[deck_pos];
}
}
ISC isc(props, deck, eclipseState, grid, gravity, swat_init);
const auto pu = props.phaseUsage();
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
? pu.phase_pos[BlackoilPhases::Liquid]
: pu.phase_pos[BlackoilPhases::Aqua];
state.pressure() = isc.press()[ref_phase];
state.saturation() = Details::convertSats(isc.saturation());
state.gasoilratio() = isc.rs();
state.rv() = isc.rv();
initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
}
} // namespace Opm
#endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED

483
tests/test_equil.cpp
View File

@@ -22,10 +22,13 @@
#include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/utility/compressedToCartesian.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
@@ -35,9 +38,6 @@
#include <opm/parser/eclipse/EclipseState/InitConfig/Equil.hpp>
#include <opm/parser/eclipse/Units/Dimension.hpp>
#include <opm/core/pressure/msmfem/partition.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <array>
@@ -49,6 +49,16 @@
#include <string>
#include <vector>
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Forward declaring the MaterialLawManager template.
typedef Opm::ThreePhaseMaterialTraits<double,
/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> MaterialTraits;
typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
#define CHECK(value, expected, reltol) \
{ \
if (std::fabs((expected)) < 1.e-14) \
@@ -57,6 +67,71 @@
BOOST_CHECK_CLOSE((value), (expected), (reltol)); \
}
namespace
{
static void 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/m3]
double rhoRefG = 1000; // [kg/m3]
double rhoRefW = 1000; // [kg/m3]
FluidSystem::initBegin(/*numPvtRegions=*/1);
FluidSystem::setEnableDissolvedGas(false);
FluidSystem::setEnableVaporizedOil(false);
FluidSystem::setReferenceDensities(rhoRefO, rhoRefW, rhoRefG, /*regionIdx=*/0);
auto gasPvt = std::make_shared<Opm::GasPvtMultiplexer<double>>();
gasPvt->setApproach(Opm::GasPvtMultiplexer<double>::DryGasPvt);
auto& dryGasPvt = gasPvt->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);
auto oilPvt = std::make_shared<Opm::OilPvtMultiplexer<double>>();
oilPvt->setApproach(Opm::OilPvtMultiplexer<double>::DeadOilPvt);
auto& deadOilPvt = oilPvt->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);
auto waterPvt = std::make_shared<Opm::WaterPvtMultiplexer<double>>();
waterPvt->setApproach(Opm::WaterPvtMultiplexer<double>::ConstantCompressibilityWaterPvt);
auto& ccWaterPvt = waterPvt->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();
FluidSystem::setGasPvt(std::move(gasPvt));
FluidSystem::setOilPvt(std::move(oilPvt));
FluidSystem::setWaterPvt(std::move(waterPvt));
FluidSystem::initEnd();
}
}
BOOST_AUTO_TEST_SUITE ()
static Opm::EquilRecord mkEquilRecord( double datd, double datp,
@@ -123,35 +198,21 @@ BOOST_AUTO_TEST_CASE (PhasePressure)
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::ParameterGroup param;
{
using Opm::unit::kilogram;
using Opm::unit::meter;
using Opm::unit::cubic;
std::stringstream dens; dens << 700*kilogram/cubic(meter);
param.insertParameter("rho2", dens.str());
}
typedef Opm::BlackoilPropertiesBasic Props;
Props props(param, G->dimensions, G->number_of_cells);
typedef Opm::EQUIL::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
RhoCalc calc(props, 0);
auto record = mkEquilRecord( 0, 1e5, 5, 0, 0, 0 );
Opm::EQUIL::EquilReg<RhoCalc>
region(record, calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage());
initDefaultFluidSystem();
Opm::EQUIL::EquilReg
region(record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0);
std::vector<int> cells(G->number_of_cells);
std::iota(cells.begin(), cells.end(), 0);
const double grav = 10;
const PPress ppress = Opm::EQUIL::phasePressures(*G, region, cells, grav);
const PPress ppress = Opm::EQUIL::phasePressures<FluidSystem>(*G, region, cells, grav);
const int first = 0, last = G->number_of_cells - 1;
const double reltol = 1.0e-8;
@@ -169,47 +230,33 @@ BOOST_AUTO_TEST_CASE (CellSubset)
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::ParameterGroup param;
{
using Opm::unit::kilogram;
using Opm::unit::meter;
using Opm::unit::cubic;
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);
initDefaultFluidSystem();
Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
Opm::EQUIL::EquilReg<RhoCalc> region[] =
{
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
};
Opm::EQUIL::EquilReg region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0)
};
const int cdim[] = { 2, 1, 2 };
int ncoarse = cdim[0];
@@ -239,7 +286,7 @@ BOOST_AUTO_TEST_CASE (CellSubset)
const int rno = int(r - cells.begin());
const double grav = 10;
const PPress p =
Opm::EQUIL::phasePressures(*G, region[rno], *r, grav);
Opm::EQUIL::phasePressures<FluidSystem>(*G, region[rno], *r, grav);
PVal::size_type i = 0;
for (std::vector<int>::const_iterator
@@ -272,61 +319,55 @@ BOOST_AUTO_TEST_CASE (RegMapping)
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::ParameterGroup param;
{
using Opm::unit::kilogram;
using Opm::unit::meter;
using Opm::unit::cubic;
std::stringstream dens; dens << 700*kilogram/cubic(meter);
param.insertParameter("rho2", dens.str());
}
typedef Opm::BlackoilPropertiesBasic Props;
Props props(param, G->dimensions, G->number_of_cells);
typedef Opm::EQUIL::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
RhoCalc calc(props, 0);
Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
Opm::EQUIL::EquilReg<RhoCalc> region[] =
{
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::EQUIL::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
props.phaseUsage())
initDefaultFluidSystem();
Opm::EQUIL::EquilReg region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0)
};
std::vector<int> eqlnum(G->number_of_cells);
// [ 0 1; 2 3]
{
std::vector<int> cells(G->number_of_cells);
std::iota(cells.begin(), cells.end(), 0);
const int cdim[] = { 2, 1, 2 };
int ncoarse = cdim[0];
for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; }
partition_unif_idx(G->dimensions, G->number_of_cells,
G->cartdims, cdim,
&cells[0], &eqlnum[0]);
for (int i = 0; i < 5; ++i) {
for (int j = 0; j < 5; ++j) {
eqlnum[i*10 + j] = 0;
}
for (int j = 5; j < 10; ++j) {
eqlnum[i*10 + j] = 1;
}
}
for (int i = 5; i < 10; ++i) {
for (int j = 0; j < 5; ++j) {
eqlnum[i*10 + j] = 2;
}
for (int j = 5; j < 10; ++j) {
eqlnum[i*10 + j] = 3;
}
}
}
Opm::RegionMapping<> eqlmap(eqlnum);
PPress ppress(2, PVal(G->number_of_cells, 0));
@@ -336,7 +377,7 @@ BOOST_AUTO_TEST_CASE (RegMapping)
const int rno = r;
const double grav = 10;
const PPress p =
Opm::EQUIL::phasePressures(*G, region[rno], rng, grav);
Opm::EQUIL::phasePressures<FluidSystem>(*G, region[rno], rng, grav);
PVal::size_type i = 0;
for (const auto& c : rng) {
@@ -366,9 +407,19 @@ BOOST_AUTO_TEST_CASE (DeckAllDead)
Opm::ParseContext parseContext;
Opm::Parser parser;
Opm::Deck deck = parser.parseFile("deadfluids.DATA" , parseContext);
Opm::EclipseState eclipseState(deck, parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, *grid, false);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, *grid, 10.0);
Opm::EclipseState eclipseState(deck, parseContext); // Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid->number_of_cells, grid->global_cell);
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, *grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells);
@@ -395,8 +446,20 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
typedef MaterialLawManager::MaterialLaw MaterialLaw;
if (!FluidSystem::isInitialized()) {
// make sure that we don't initialize the fluid system twice
FluidSystem::initFromDeck(deck, eclipseState);
}
// Test the capillary inversion for oil-water.
const int cell = 0;
const double reltol = 1.0e-7;
@@ -407,7 +470,7 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
const std::vector<double> s = { 0.2, 0.2, 0.2, 0.466666666666, 0.733333333333, 1.0, 1.0, 1.0, 1.0 };
BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromPc(props, phase, cell, pc[i], increasing);
const double s_computed = Opm::EQUIL::satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, phase, cell, pc[i], increasing);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
@@ -420,7 +483,7 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
const std::vector<double> s = { 0.8, 0.8, 0.8, 0.533333333333, 0.266666666666, 0.0, 0.0, 0.0, 0.0 };
BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromPc(props, phase, cell, pc[i], increasing);
const double s_computed = Opm::EQUIL::satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, phase, cell, pc[i], increasing);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
@@ -433,7 +496,7 @@ BOOST_AUTO_TEST_CASE (CapillaryInversion)
const std::vector<double> s = { 0.2, 0.333333333333, 0.6, 0.866666666666, 1.0 };
BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::EQUIL::satFromSumOfPcs(props, water, gas, cell, pc[i]);
const double s_computed = Opm::EQUIL::satFromSumOfPcs<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, water, gas, cell, pc[i]);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
@@ -449,9 +512,20 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 10.0);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -490,9 +564,19 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("capillary_overlap.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -553,9 +637,17 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("equil_liveoil.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -633,9 +725,18 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
Opm::ParseContext parseContext;
Opm::Deck deck = parser.parseFile("equil_livegas.DATA" , parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -715,9 +816,18 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
Opm::Parser parser;
Opm::Deck deck = parser.parseFile("equil_rsvd_and_rvvd.DATA", parseContext);
Opm::EclipseState eclipseState(deck , parseContext);
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> comp(materialLawManager, eclipseState, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@@ -817,10 +927,15 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
Opm::EclipseState eclipseState(deck , parseContext);
Opm::GridManager gm(eclipseState.getInputGrid());
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false);
Opm::BlackoilPropertiesFromDeck propsScaled(deck, eclipseState, grid, false);
Opm::BlackoilState state( Opm::UgGridHelpers::numCells( grid ) , Opm::UgGridHelpers::numFaces( grid ) , 3);
// Create material law manager.
std::vector<int> compressedToCartesianIdx
= Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
MaterialLawManager materialLawManager = MaterialLawManager();
materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
MaterialLawManager materialLawManagerScaled = MaterialLawManager();
materialLawManagerScaled.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
// reference saturations
const std::vector<double> s[3]{
@@ -835,25 +950,44 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
{ 0, 0, 0, 0.014813991154779993, 0.78525420807446045, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.78518600884522005, 0.014745791925539575, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
std::vector<double> sats = state.saturation();
for (int phase = 0; phase < 3; ++phase) {
for (size_t i = 0; i < 20; ++i) {
sats[3*i + phase] = s[phase][i];
}
}
std::vector<double> sats_swatinit = state.saturation();
for (int phase = 0; phase < 3; ++phase) {
for (size_t i = 0; i < 20; ++i) {
sats_swatinit[3*i + phase] = swatinit[phase][i];
}
}
// Adjust oil pressure according to gas saturation and cap pressure
typedef Opm::SimpleModularFluidState<double,
/*numPhases=*/3,
/*numComponents=*/3,
FluidSystem,
/*storePressure=*/false,
/*storeTemperature=*/false,
/*storeComposition=*/false,
/*storeFugacity=*/false,
/*storeSaturation=*/true,
/*storeDensity=*/false,
/*storeViscosity=*/false,
/*storeEnthalpy=*/false> SatOnlyFluidState;
SatOnlyFluidState fluidState;
typedef MaterialLawManager::MaterialLaw MaterialLaw;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
// reference pcs
const int numCells = Opm::UgGridHelpers::numCells(grid);
std::vector<int> cells(numCells);
for (int c = 0; c < numCells; ++c) { cells[c] = c; }
std::vector<double> pc_original = state.saturation();
props.capPress(numCells, sats.data(), cells.data(), pc_original.data(), nullptr);
std::vector<double> pc_original(numCells * FluidSystem::numPhases);
for (int c = 0; c < numCells; ++c) {
std::vector<double> pc = {0,0,0};
double sw = s[0][c];
double so = s[1][c];
double sg = s[2][c];
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;
@@ -873,20 +1007,45 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
pc_scaled_truth[3*11 + 0] = 5364.1;
// compute the initial state
// apply swatinit
Opm::BlackoilState state_scaled = state;
initStateEquil(grid, propsScaled, deck, eclipseState, 9.81, state_scaled, true);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> compScaled(materialLawManagerScaled, eclipseState, grid, 9.81, true);
// don't apply swatinit
Opm::BlackoilState state_unscaled = state;
initStateEquil(grid, props, deck, eclipseState, 9.81, state_unscaled, false);
Opm::EQUIL::DeckDependent::InitialStateComputer<FluidSystem> compUnscaled(materialLawManager, eclipseState, grid, 9.81, false);
// compute pc
std::vector<double> pc_scaled= state.saturation();
propsScaled.capPress(numCells, state_scaled.saturation().data(), cells.data(), pc_scaled.data(), nullptr);
std::vector<double> pc_unscaled= state.saturation();
props.capPress(numCells, state_unscaled.saturation().data(), cells.data(), pc_unscaled.data(), nullptr);
std::vector<double> pc_scaled(numCells * FluidSystem::numPhases);
for (int c = 0; c < numCells; ++c) {
std::vector<double> pc = {0,0,0};
double sw = compScaled.saturation().data()[0][c];
double so = compScaled.saturation().data()[1][c];
double sg = compScaled.saturation().data()[2][c];
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(numCells * FluidSystem::numPhases);
for (int c = 0; c < numCells; ++c) {
std::vector<double> pc = {0,0,0};
double sw = compUnscaled.saturation().data()[0][c];
double so = compUnscaled.saturation().data()[1][c];
double sg = compUnscaled.saturation().data()[2][c];
fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg);
const auto& matParams = materialLawManager.materialLawParams(c);
MaterialLaw::capillaryPressures(pc, matParams, fluidState);
pc_unscaled[3*c + 0] = pc[FluidSystem::oilPhaseIdx] - pc[FluidSystem::waterPhaseIdx];
pc_unscaled[3*c + 1] = 0.0;
pc_unscaled[3*c + 2] = pc[FluidSystem::oilPhaseIdx] + pc[FluidSystem::gasPhaseIdx];
}
// test
const double reltol = 1.0e-3;
@@ -899,8 +1058,8 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
for (int phase = 0; phase < 3; ++phase) {
for (size_t i = 0; i < 20; ++i) {
CHECK(state_unscaled.saturation()[3*i + phase], s[phase][i], reltol);
CHECK(state_scaled.saturation()[3*i + phase], swatinit[phase][i], reltol);
CHECK(compUnscaled.saturation()[phase][i], s[phase][i], reltol);
CHECK(compScaled.saturation()[phase][i], swatinit[phase][i], reltol);
}
}
}