SWATINIT: Initialisation and capillary pressure scaling.

This commit is contained in:
osae 2014-06-10 14:02:22 +02:00 committed by Andreas Lauser
parent 897f64c21a
commit 10f5b07915
2 changed files with 36 additions and 12 deletions

View File

@ -150,8 +150,9 @@ namespace Opm
std::vector< std::vector<double> > std::vector< std::vector<double> >
phaseSaturations(const Region& reg, phaseSaturations(const Region& reg,
const CellRange& cells, const CellRange& cells,
const BlackoilPropertiesInterface& props, BlackoilPropertiesInterface& props,
const std::vector< std::vector<double> >& phase_pressures); const std::vector<double> swat_init,
std::vector< std::vector<double> >& phase_pressures);
@ -254,7 +255,7 @@ namespace Opm
class InitialStateComputer { class InitialStateComputer {
public: public:
InitialStateComputer(const BlackoilPropertiesInterface& props, InitialStateComputer(BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr deck, const Opm::DeckConstPtr deck,
const UnstructuredGrid& G , const UnstructuredGrid& G ,
const double grav = unit::gravity) const double grav = unit::gravity)
@ -333,6 +334,19 @@ namespace Opm
rv_func_.push_back(std::make_shared<Miscibility::NoMixing>()); rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
} }
} }
// Check for presence of kw SWATINIT
if (deck->hasKeyword("SWATINIT")) {
const std::vector<double>& 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. // Compute pressures, saturations, rs and rv factors.
calcPressSatRsRv(eqlmap, rec, props, G, grav); calcPressSatRsRv(eqlmap, rec, props, G, grav);
@ -360,15 +374,18 @@ namespace Opm
PVec sat_; PVec sat_;
Vec rs_; Vec rs_;
Vec rv_; Vec rv_;
Vec swat_init_;
template <class RMap> template <class RMap>
void void
calcPressSatRsRv(const RMap& reg , calcPressSatRsRv(const RMap& reg ,
const std::vector< EquilRecord >& rec , const std::vector< EquilRecord >& rec ,
const Opm::BlackoilPropertiesInterface& props, Opm::BlackoilPropertiesInterface& props,
const UnstructuredGrid& G , const UnstructuredGrid& G ,
const double grav) const double grav)
{ {
typedef Miscibility::NoMixing NoMix;
for (typename RMap::RegionId for (typename RMap::RegionId
r = 0, nr = reg.numRegions(); r = 0, nr = reg.numRegions();
r < nr; ++r) r < nr; ++r)
@ -383,7 +400,7 @@ namespace Opm
PVec press = phasePressures(G, eqreg, cells, grav); 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(); const int np = props.numPhases();
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {

View File

@ -579,7 +579,8 @@ namespace Opm
std::vector< std::vector<double> > std::vector< std::vector<double> >
phaseSaturations(const Region& reg, phaseSaturations(const Region& reg,
const CellRange& cells, const CellRange& cells,
const BlackoilPropertiesInterface& props, BlackoilPropertiesInterface& props,
const std::vector<double> swat_init,
std::vector< std::vector<double> >& phase_pressures) std::vector< std::vector<double> >& phase_pressures)
{ {
const double z0 = reg.datum(); const double z0 = reg.datum();
@ -610,8 +611,14 @@ namespace Opm
double sw = 0.0; double sw = 0.0;
if (water) { if (water) {
const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index]; const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
sw = satFromPc(props, waterpos, cell, pcov); if (swat_init.empty()) { // Invert Pc to find sw
phase_saturations[waterpos][local_index] = 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; double sg = 0.0;
if (gas) { if (gas) {
@ -736,7 +743,7 @@ namespace Opm
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction. * \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
*/ */
void initStateEquil(const UnstructuredGrid& grid, void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props, BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr deck, const Opm::DeckConstPtr deck,
const double gravity, const double gravity,
BlackoilState& state) BlackoilState& state)