From 10f5b07915e4c2ad076853877829d5861221b75e Mon Sep 17 00:00:00 2001 From: osae Date: Tue, 10 Jun 2014 14:02:22 +0200 Subject: [PATCH] SWATINIT: Initialisation and capillary pressure scaling. --- opm/core/simulator/initStateEquil.hpp | 33 ++++++++++++++++------ opm/core/simulator/initStateEquil_impl.hpp | 15 +++++++--- 2 files changed, 36 insertions(+), 12 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index a4ec1ee4f..273271f2d 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -150,8 +150,9 @@ namespace Opm std::vector< std::vector > phaseSaturations(const Region& reg, const CellRange& cells, - const BlackoilPropertiesInterface& props, - const std::vector< std::vector >& phase_pressures); + BlackoilPropertiesInterface& props, + const std::vector swat_init, + std::vector< std::vector >& phase_pressures); @@ -254,7 +255,7 @@ namespace Opm class InitialStateComputer { public: - InitialStateComputer(const BlackoilPropertiesInterface& props, + InitialStateComputer(BlackoilPropertiesInterface& props, const Opm::DeckConstPtr deck, const UnstructuredGrid& G , const double grav = unit::gravity) @@ -333,6 +334,19 @@ namespace Opm rv_func_.push_back(std::make_shared()); } } + + + // Check for presence of kw SWATINIT + if (deck->hasKeyword("SWATINIT")) { + const std::vector& swat_init = deck->getKeyword("SWATINIT")->getSIDoubleData(); + swat_init_.resize(G.number_of_cells); + const int* gc = G.global_cell; + for (int c = 0; c < G.number_of_cells; ++c) { + const int deck_pos = (gc == NULL) ? c : gc[c]; + swat_init_[c] = swat_init[deck_pos]; + } + + } // Compute pressures, saturations, rs and rv factors. calcPressSatRsRv(eqlmap, rec, props, G, grav); @@ -360,15 +374,18 @@ namespace Opm PVec sat_; Vec rs_; Vec rv_; + Vec swat_init_; template void - calcPressSatRsRv(const RMap& reg , - const std::vector< EquilRecord >& rec , - const Opm::BlackoilPropertiesInterface& props, - const UnstructuredGrid& G , + calcPressSatRsRv(const RMap& reg , + const std::vector< EquilRecord >& rec , + Opm::BlackoilPropertiesInterface& props, + const UnstructuredGrid& G , const double grav) { + typedef Miscibility::NoMixing NoMix; + for (typename RMap::RegionId r = 0, nr = reg.numRegions(); r < nr; ++r) @@ -383,7 +400,7 @@ namespace Opm PVec press = phasePressures(G, eqreg, cells, grav); - const PVec sat = phaseSaturations(eqreg, cells, props, press); + const PVec sat = phaseSaturations(eqreg, cells, props, swat_init_, press); const int np = props.numPhases(); for (int p = 0; p < np; ++p) { diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index ba2df1ed4..31215de58 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -579,7 +579,8 @@ namespace Opm std::vector< std::vector > phaseSaturations(const Region& reg, const CellRange& cells, - const BlackoilPropertiesInterface& props, + BlackoilPropertiesInterface& props, + const std::vector swat_init, std::vector< std::vector >& phase_pressures) { const double z0 = reg.datum(); @@ -610,8 +611,14 @@ namespace Opm double sw = 0.0; if (water) { const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index]; - sw = satFromPc(props, waterpos, cell, pcov); - phase_saturations[waterpos][local_index] = sw; + if (swat_init.empty()) { // Invert Pc to find sw + sw = satFromPc(props, 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); + phase_saturations[waterpos][local_index] = sw; + } } double sg = 0.0; if (gas) { @@ -736,7 +743,7 @@ namespace Opm * \param[in] gravity Acceleration of gravity, assumed to be in Z direction. */ void initStateEquil(const UnstructuredGrid& grid, - const BlackoilPropertiesInterface& props, + BlackoilPropertiesInterface& props, const Opm::DeckConstPtr deck, const double gravity, BlackoilState& state)