OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Remove 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
committed by Andreas Lauser
parent 8efc60957a
commit 4f915d116f
3 changed files with 75 additions and 227 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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