OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-07-07 04:53:03 -05:00
Code Issues Packages Projects Releases Wiki Activity

Remove blackoilPhases and phaseUsage from the initialization code

Browse Source 
Note 1: The initialization code now always consider 3 phases.
For 2-phase cases a trivial (0) state is returned.
Note 2: The initialization code does not compute a BlackoilStats,
but instead pass the initialization object with the initial state.
This commit is contained in:
Tor Harald Sandve 2017-11-20 08:47:41 +01:00
parent 36e8f1bfec
commit 51f48fcd13
5 changed files with 205 additions and 317 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
examples/compute_initial_state.cpp
View File

@ -28,6 +28,8 @@
#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>
@ -75,12 +77,45 @@ 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)
@ -119,7 +154,32 @@ try
// Initialisation.
//initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
BlackoilState state( UgGridHelpers::numCells(grid) , UgGridHelpers::numFaces(grid), 3);
initStateEquil(grid, materialLawManager, deck, eclipseState, grav, state);
typedef FluidSystems::BlackOil<double> FluidSystem;
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);
if (state.hasCellData(std::string("GASOILRATIO"))) {
std::vector<double>& rs = state.getCellData(std::string("GASOILRATIO"));
rs = isc.rs();
}
if (state.hasCellData(std::string("RV"))){
std::vector<double>& rv = state.getCellData(std::string("RV"));
rv = isc.rv();
}
// Output.
const std::string output_dir = param.getDefault<std::string>("output_dir", "output");
@ -132,3 +192,5 @@ catch (const std::exception& e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

20
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,7 +21,6 @@
#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/utility/RootFinders.hpp>
@ -478,13 +478,11 @@ namespace Opm
EquilReg(const EquilRecord& rec,
std::shared_ptr<Miscibility::RsFunction> rs,
std::shared_ptr<Miscibility::RsFunction> rv,
const int pvtIdx,
const PhaseUsage& pu)
const int pvtIdx)
: rec_ (rec)
, rs_ (rs)
, rv_ (rv)
, pvtIdx_ (pvtIdx)
, pu_ (pu)
{
}
@ -554,18 +552,12 @@ namespace Opm
const int
pvtIdx() const { return this->pvtIdx_; }
/**
* Retrieve active fluid phase summary.
*/
const PhaseUsage&
phaseUsage() const { return this->pu_; }
private:
EquilRecord rec_; /**< Equilibration data */
std::shared_ptr<Miscibility::RsFunction> rs_; /**< RS calculator */
std::shared_ptr<Miscibility::RsFunction> rv_; /**< RV calculator */
const int pvtIdx_;
PhaseUsage pu_; /**< Active phase summary */
};
@ -588,7 +580,7 @@ namespace Opm
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_ + FluidSystem::numPhases, 0.0);
}
@ -609,7 +601,7 @@ namespace Opm
const int cell_;
const double target_pc_;
mutable SatOnlyFluidState fluidState_;
mutable double pc_[BlackoilPhases::MaxNumPhases];
mutable double pc_[FluidSystem::numPhases];
};
template <class FluidSystem, class MaterialLawManager>
@ -722,7 +714,7 @@ namespace Opm
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_ + FluidSystem::numPhases, 0.0);
}
double operator()(double s) const
{
@ -745,7 +737,7 @@ namespace Opm
const int cell_;
const double target_pc_;
mutable SatOnlyFluidState fluidState_;
mutable double pc_[BlackoilPhases::MaxNumPhases];
mutable double pc_[FluidSystem::numPhases];
};

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

@ -24,9 +24,6 @@
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/simulator/EquilibrationHelpers.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/utility/extractPvtTableIndex.hpp>
@ -40,6 +37,8 @@
#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>
@ -198,12 +197,12 @@ 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 =
@ -223,45 +222,48 @@ namespace Opm
return eqlnum;
}
template<class FluidSystem>
class InitialStateComputer {
public:
template<class MaterialLawManager, class Grid>
InitialStateComputer(std::shared_ptr<MaterialLawManager> materialLawManager,
const PhaseUsage& phaseUsage,
const Opm::Deck& deck,
const Opm::EclipseState& eclipseState,
const Grid& G ,
const double grav = unit::gravity,
const std::vector<double>& swat_init = {}
const bool applySwatInit = true
)
: pp_(phaseUsage.num_phases,
: pp_(FluidSystem::numPhases,
std::vector<double>(UgGridHelpers::numCells(G))),
sat_(phaseUsage.num_phases,
sat_(FluidSystem::numPhases,
std::vector<double>(UgGridHelpers::numCells(G))),
rs_(UgGridHelpers::numCells(G)),
rv_(UgGridHelpers::numCells(G)),
swat_init_(swat_init),
phaseUsage_(phaseUsage)
rv_(UgGridHelpers::numCells(G))
{
typedef FluidSystems::BlackOil<double> FluidSystem;
// Initialize the fluid system
FluidSystem::initFromDeck(deck, eclipseState);
//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())
@ -299,7 +301,7 @@ 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())
@ -338,7 +340,7 @@ namespace Opm
}
// Compute pressures, saturations, rs and rv factors.
calcPressSatRsRv<FluidSystem>(eqlmap, rec, materialLawManager, 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.
@ -366,7 +368,6 @@ namespace Opm
Vec rs_;
Vec rv_;
Vec swat_init_;
PhaseUsage phaseUsage_;
template<class Grid, class RMap>
void setRegionPvtIdx(const Grid& G, const Opm::EclipseState& eclipseState, const RMap& reg) {
@ -381,7 +382,7 @@ namespace Opm
}
}
template <class FluidSystem, class RMap, class MaterialLawManager, class Grid>
template <class RMap, class MaterialLawManager, class Grid>
void
calcPressSatRsRv(const RMap& reg ,
const std::vector< EquilRecord >& rec ,
@ -398,23 +399,24 @@ namespace Opm
continue;
}
const EqReg eqreg(rec[r],
rs_func_[r], rv_func_[r], regionPvtIdx_[r],
phaseUsage_);
const EqReg eqreg(rec[r], rs_func_[r], rv_func_[r], regionPvtIdx_[r]);
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 int np = phaseUsage_.num_phases;
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 (phaseUsage_.phase_used[BlackoilPhases::Liquid]
&& phaseUsage_.phase_used[BlackoilPhases::Vapour]) {
const int oilpos = phaseUsage_.phase_pos[BlackoilPhases::Liquid];
const int gaspos = phaseUsage_.phase_pos[BlackoilPhases::Vapour];
const 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_);

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

@ -25,7 +25,7 @@
#include <opm/core/grid.h>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <cassert>
#include <cmath>
@ -237,60 +237,6 @@ namespace Opm
};
} // 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,
@ -491,70 +437,66 @@ namespace Opm
const CellRange& cells,
std::vector< std::vector<double> >& press)
{
const PhaseUsage& pu = reg.phaseUsage();
const bool water = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
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);
if (water) {
PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
cells, press[ wix ]);
cells, press[ waterpos ]);
}
if (PhaseUsed::oil(pu)) {
const int oix = PhaseIndex::oil(pu);
if (oil) {
PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
press[ oix ], po_woc, po_goc);
press[ oilpos ], po_woc, po_goc);
}
if (PhaseUsed::gas(pu)) {
const int gix = PhaseIndex::gas(pu);
if (gas) {
PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
cells, press[ gix ]);
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);
if (gas) {
PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
cells, press[ gix ]);
cells, press[ gaspos ]);
}
if (PhaseUsed::oil(pu)) {
const int oix = PhaseIndex::oil(pu);
if (oil) {
PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
press[ oix ], po_woc, po_goc);
press[ oilpos ], po_woc, po_goc);
}
if (PhaseUsed::water(pu)) {
const int wix = PhaseIndex::water(pu);
if (water) {
PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
cells, press[ wix ]);
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);
if (oil) {
PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
press[ oix ], po_woc, po_goc);
press[ oilpos ], po_woc, po_goc);
}
if (PhaseUsed::water(pu)) {
const int wix = PhaseIndex::water(pu);
if (water) {
PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
cells, press[ wix ]);
cells, press[ waterpos ]);
}
if (PhaseUsed::gas(pu)) {
const int gix = PhaseIndex::gas(pu);
if (gas) {
PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
cells, press[ gix ]);
cells, press[ gaspos ]);
}
}
}
@ -621,7 +563,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));
@ -659,7 +601,7 @@ namespace Opm
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.");
}
@ -682,11 +624,11 @@ namespace Opm
SatOnlyFluidState fluidState;
typedef typename MaterialLawManager::MaterialLaw MaterialLaw;
const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
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];
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;
@ -769,7 +711,7 @@ namespace Opm
double threshold_sat = 1.0e-6;
double so = 1.0;
double pC[/*numPhases=*/BlackoilPhases::MaxNumPhases] = { 0.0, 0.0, 0.0 };
double pC[FluidSystem::numPhases] = { 0.0, 0.0, 0.0 };
if (water) {
double swu = scaledDrainageInfo.Swu;
@ -850,94 +792,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 MaterialLawManager, class Grid>
void initStateEquil(const Grid& grid,
std::shared_ptr<MaterialLawManager> materialLawManager,
const Opm::Deck& deck,
const Opm::EclipseState& eclipseState,
const double gravity,
BlackoilState& state,
bool applySwatinit = true)
{
typedef EQUIL::DeckDependent::InitialStateComputer ISC;
PhaseUsage pu = phaseUsageFromDeck(deck);
//Check for presence of kw SWATINIT
std::vector<double> swat_init = {};
if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatinit) {
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(materialLawManager, pu, deck, eclipseState, grid, gravity, swat_init);
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
? pu.phase_pos[BlackoilPhases::Liquid]
: pu.phase_pos[BlackoilPhases::Aqua];
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

156
tests/test_equil.cpp
View File

@ -23,13 +23,11 @@
#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/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
@ -42,7 +40,6 @@
#include <opm/core/pressure/msmfem/partition.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <array>
@ -55,11 +52,12 @@
#include <vector>
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
// Forward declaring the MaterialLawManager template.
typedef Opm::ThreePhaseMaterialTraits<double,
/*wettingPhaseIdx=*/Opm::BlackoilPhases::Aqua,
/*nonWettingPhaseIdx=*/Opm::BlackoilPhases::Liquid,
/*gasPhaseIdx=*/Opm::BlackoilPhases::Vapour> MaterialTraits;
/*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
/*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
/*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> MaterialTraits;
typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
@ -71,8 +69,8 @@ typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
BOOST_CHECK_CLOSE((value), (expected), (reltol)); \
}
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
Opm::PhaseUsage initDefaultFluidSystem() {
void initDefaultFluidSystem() {
std::vector<std::pair<double, double> > Bo = {
{ 101353, 1. },
{ 6.21542e+07, 1 }
@ -139,17 +137,6 @@ Opm::PhaseUsage initDefaultFluidSystem() {
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 ()
@ -220,14 +207,13 @@ BOOST_AUTO_TEST_CASE (PhasePressure)
auto record = mkEquilRecord( 0, 1e5, 5, 0, 0, 0 );
Opm::PhaseUsage pu = initDefaultFluidSystem();
initDefaultFluidSystem();
Opm::EQUIL::EquilReg
region(record,
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu);
0);
std::vector<int> cells(G->number_of_cells);
std::iota(cells.begin(), cells.end(), 0);
@ -251,8 +237,7 @@ BOOST_AUTO_TEST_CASE (CellSubset)
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::PhaseUsage pu = initDefaultFluidSystem();
initDefaultFluidSystem();
Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
@ -262,26 +247,22 @@ BOOST_AUTO_TEST_CASE (CellSubset)
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
0)
};
const int cdim[] = { 2, 1, 2 };
@ -348,33 +329,29 @@ BOOST_AUTO_TEST_CASE (RegMapping)
Opm::EquilRecord record[] = { mkEquilRecord( 0, 1e5, 2.5, -0.075e5, 0, 0 ),
mkEquilRecord( 5, 1.35e5, 7.5, -0.225e5, 5, 0 ) };
Opm::PhaseUsage pu = initDefaultFluidSystem();
initDefaultFluidSystem();
Opm::EQUIL::EquilReg region[] =
{
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
0)
,
Opm::EQUIL::EquilReg(record[0],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
0)
,
Opm::EQUIL::EquilReg(record[1],
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
std::make_shared<Opm::EQUIL::Miscibility::NoMixing>(),
0,
pu)
0)
};
std::vector<int> eqlnum(G->number_of_cells);
@ -436,9 +413,12 @@ BOOST_AUTO_TEST_CASE (DeckAllDead)
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, *grid, 10.0);
// 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);
@ -538,9 +518,13 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary)
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
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);
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
@ -585,9 +569,13 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
// 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);
@ -654,9 +642,11 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
// 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);
@ -740,9 +730,12 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
// 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);
@ -828,9 +821,12 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::PhaseUsage pu = phaseUsageFromDeck(deck);
typedef Opm::FluidSystems::BlackOil<double> FluidSystem;
Opm::EQUIL::DeckDependent::InitialStateComputer comp(materialLawManager, pu, deck, eclipseState, grid, 9.80665);
// 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);
@ -940,8 +936,6 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
auto materialLawManagerScaled = std::make_shared<MaterialLawManager>();
materialLawManagerScaled->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
Opm::BlackoilState state( Opm::UgGridHelpers::numCells( grid ) , Opm::UgGridHelpers::numFaces( grid ) , 3);
// reference saturations
const std::vector<double> s[3]{
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.42528761746004229, 0.77462669821009045, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
@ -955,18 +949,6 @@ 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,
@ -989,13 +971,13 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
FluidSystem::initFromDeck(deck, eclipseState);
// reference pcs
std::vector<double> pc_original = state.saturation();
const int numCells = Opm::UgGridHelpers::numCells(grid);
std::vector<double> pc_original(numCells * FluidSystem::numPhases);
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];
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);
@ -1024,22 +1006,18 @@ 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, materialLawManagerScaled, 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, materialLawManager, 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();
std::vector<double> pc_scaled(numCells * FluidSystem::numPhases);
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];
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);
@ -1050,12 +1028,12 @@ BOOST_AUTO_TEST_CASE (DeckWithSwatinit)
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(numCells * FluidSystem::numPhases);
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];
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);
@ -1079,8 +1057,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);
}
}
}