Merge pull request #1084 from babrodtk/hysteresis_output

Adds hysteresis output and input (for restarting)
This commit is contained in:
Atgeirr Flø Rasmussen 2017-04-11 07:44:59 +02:00 committed by GitHub
commit 88e4646b71
13 changed files with 263 additions and 20 deletions

View File

@ -105,6 +105,12 @@ namespace Opm {
std::vector<double> soMax; // Maximum oil saturation
//Hysteresis parameters
std::vector<double> krnswdc_ow;
std::vector<double> krnswdc_go;
std::vector<double> pcswmdc_ow;
std::vector<double> pcswmdc_go;
std::array<V, fipValues> fip;
};

View File

@ -464,6 +464,10 @@ typedef Eigen::Array<double,
, rsSat(ADB::null())
, rvSat(ADB::null())
, soMax()
, krnswdc_ow()
, krnswdc_go()
, pcswmdc_ow()
, pcswmdc_go()
, fip()
{
}
@ -664,6 +668,8 @@ typedef Eigen::Array<double,
state.rv = sd_.rvSat;
}
sd_.soMax = fluid_.satOilMax();
fluid_.getGasOilHystParams(sd_.krnswdc_go, sd_.pcswmdc_go, cells_);
fluid_.getOilWaterHystParams(sd_.krnswdc_ow, sd_.pcswmdc_ow, cells_);
}
}
else {

View File

@ -1250,6 +1250,20 @@ namespace Opm {
simData.registerCellData( "SOMAX", 1 );
VectorType& somax = simData.getCellData( "SOMAX" );
// Two components for hysteresis parameters
// pcSwMdc/krnSwMdc, one for oil-water and one for gas-oil
simData.registerCellData( "PCSWMDC_GO", 1 );
simData.registerCellData( "KRNSWMDC_GO", 1 );
simData.registerCellData( "PCSWMDC_OW", 1 );
simData.registerCellData( "KRNSWMDC_OW", 1 );
VectorType& pcSwMdc_go = simData.getCellData( "PCSWMDC_GO" );
VectorType& krnSwMdc_go = simData.getCellData( "KRNSWMDC_GO" );
VectorType& pcSwMdc_ow = simData.getCellData( "PCSWMDC_OW" );
VectorType& krnSwMdc_ow = simData.getCellData( "KRNSWMDC_OW" );
std::vector<int> failed_cells_pb;
std::vector<int> failed_cells_pd;
const auto& gridView = ebosSimulator().gridView();
@ -1276,6 +1290,18 @@ namespace Opm {
somax[cellIdx] = ebosSimulator().model().maxOilSaturation(cellIdx);
const auto& matLawManager = ebosSimulator().problem().materialLawManager();
if (matLawManager->enableHysteresis()) {
matLawManager->oilWaterHysteresisParams(
pcSwMdc_ow[cellIdx],
krnSwMdc_ow[cellIdx],
cellIdx);
matLawManager->gasOilHysteresisParams(
pcSwMdc_go[cellIdx],
krnSwMdc_go[cellIdx],
cellIdx);
}
if (aqua_active) {
saturation[ satIdx + aqua_pos ] = fs.saturation(FluidSystem::waterPhaseIdx).value();
bWater[cellIdx] = fs.invB(FluidSystem::waterPhaseIdx).value();

View File

@ -876,6 +876,58 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
satprops_->updateSatHyst(n, cells.data(), saturation.data());
}
/// Set gas-oil hysteresis parameters
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
void BlackoilPropsAdFromDeck::setGasOilHystParams(const std::vector<double>& pcswmdc,
const std::vector<double>& krnswdc,
const std::vector<int>& cells)
{
const int n = cells.size();
assert(pcswmdc.size() == n);
assert(krnswdc.size() == n);
satprops_->setGasOilHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
}
/// Get gas-oil hysteresis parameters
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
void BlackoilPropsAdFromDeck::getGasOilHystParams(std::vector<double>& pcswmdc,
std::vector<double>& krnswdc,
const std::vector<int>& cells) const
{
const int n = cells.size();
pcswmdc.resize(n);
krnswdc.resize(n);
satprops_->getGasOilHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
}
/// Set oil-water hysteresis parameters
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
void BlackoilPropsAdFromDeck::setOilWaterHystParams(const std::vector<double>& pcswmdc,
const std::vector<double>& krnswdc,
const std::vector<int>& cells)
{
const int n = cells.size();
assert(pcswmdc.size() == n);
assert(krnswdc.size() == n);
satprops_->setOilWaterHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
}
/// Get oil-water hysteresis parameters
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
void BlackoilPropsAdFromDeck::getOilWaterHystParams(std::vector<double>& pcswmdc,
std::vector<double>& krnswdc,
const std::vector<int>& cells) const
{
const int n = cells.size();
pcswmdc.resize(n);
krnswdc.resize(n);
satprops_->getOilWaterHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
}
/// Update for max oil saturation.
void BlackoilPropsAdFromDeck::updateSatOilMax(const std::vector<double>& saturation)
{
@ -895,6 +947,11 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
return satOilMax_;
}
void BlackoilPropsAdFromDeck::setSatOilMax(const std::vector<double>& max_sat) {
assert(satOilMax_.size() == max_sat.size());
satOilMax_ = max_sat;
}
/// Set capillary pressure scaling according to pressure diff. and initial water saturation.
/// \param[in] saturation Array of n*numPhases saturation values.
/// \param[in] pc Array of n*numPhases capillary pressure values.

View File

@ -359,12 +359,45 @@ namespace Opm
void updateSatHyst(const std::vector<double>& saturation,
const std::vector<int>& cells);
/// Set gas-oil hysteresis parameters
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
void setGasOilHystParams(const std::vector<double>& pcswmdc,
const std::vector<double>& krnswdc,
const std::vector<int>& cells);
/// Get gas-oil hysteresis parameters
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
void getGasOilHystParams(std::vector<double>& pcswmdc,
std::vector<double>& krnswdc,
const std::vector<int>& cells) const;
/// Set oil-water hysteresis parameters
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
void setOilWaterHystParams(const std::vector<double>& pcswmdc,
const std::vector<double>& krnswdc,
const std::vector<int>& cells);
/// Set oil-water hysteresis parameters
/// \param[in] pcswmdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
/// \param[in] krnswdc Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
void getOilWaterHystParams(std::vector<double>& pcswmdc,
std::vector<double>& krnswdc,
const std::vector<int>& cells) const;
/// Update for max oil saturation.
/// \param[in] saturation Saturations for all phases
void updateSatOilMax(const std::vector<double>& saturation);
/// Returns the max oil saturation
const std::vector<double>& satOilMax() const;
/// Set max oil saturation (for restarting)
/// \param[in] max_sat Max oil saturations. Note that this is *only* oil saturations
void setSatOilMax(const std::vector<double>& max_sat);
/// Set capillary pressure scaling according to pressure diff. and initial water saturation.
/// \param[in] saturation Array of n*numPhases saturation values.
/// \param[in] pc Array of n*numPhases capillary pressure values.

View File

@ -123,6 +123,7 @@ data::Solution simToSolution( const SimulationDataContainer& reservoir,
void solutionToSim( const data::Solution& sol,
const std::map<std::string,std::vector<double> >& extra,
PhaseUsage phases,
SimulationDataContainer& state ) {
@ -170,13 +171,33 @@ void solutionToSim( const data::Solution& sol,
state.getCellData( "RV" ) = sol.data( "RV" );
}
if (sol.has( "SSOL" ) ) {
if ( sol.has( "SSOL" ) ) {
state.getCellData("SSOL") = sol.data("SSOL");
}
if ( sol.has( "SOMAX" ) ) {
if ( extra.find("SOMAX") != extra.end() ) {
state.registerCellData("SOMAX", 1);
state.getCellData("SOMAX") = sol.data("SOMAX");
state.getCellData("SOMAX") = extra.find("SOMAX")->second;
}
if ( extra.find("PCSWM_OW") != extra.end() ) {
state.registerCellData("PCSWMDC_OW", 1);
state.getCellData("PCSWMDC_OW") = extra.find("PCSWM_OW")->second;
}
if ( extra.find("KRNSW_OW") != extra.end() ) {
state.registerCellData("KRNSWMDC_OW", 1);
state.getCellData("KRNSWMDC_OW") = extra.find("KRNSW_OW")->second;
}
if ( extra.find("PCSWM_GO") != extra.end() ) {
state.registerCellData("PCSWMDC_GO", 1);
state.getCellData("PCSWMDC_GO") = extra.find("PCSWM_GO")->second;
}
if ( extra.find("KRNSW_GO") != extra.end() ) {
state.registerCellData("KRNSWMDC_GO", 1);
state.getCellData("KRNSWMDC_GO") = extra.find("KRNSW_GO")->second;
}
}

View File

@ -59,7 +59,9 @@ namespace Opm {
/// Copies the following fields from sol into state (all conditionally):
/// PRESSURE, TEMP, SWAT, SGAS, RS, RV, SSOL
/// Also handles extra data such as hysteresis parameters, SOMAX, etc.
void solutionToSim( const data::Solution& sol,
const std::map<std::string,std::vector<double> >& extra,
PhaseUsage phases,
SimulationDataContainer& state );

View File

@ -665,7 +665,8 @@ namespace Opm
if( isIORank() )
{
// Update values in the globalReservoirState
solutionToSim(*globalCellData_, phaseUsage_, *globalReservoirState_);
const std::map<std::string, std::vector<double> > no_extra_data;
solutionToSim(*globalCellData_, no_extra_data, phaseUsage_, *globalReservoirState_);
}
return isIORank();
}

View File

@ -186,6 +186,7 @@ namespace Opm
const WellState& well_state,
DynamicListEconLimited& list_econ_limited) const;
void initHysteresisParams(ReservoirState& state);
// Data.
typedef RateConverter::

View File

@ -94,6 +94,7 @@ namespace Opm
// This is a restart, populate WellState and ReservoirState state objects from restart file
output_writer_.initFromRestartFile(props_.phaseUsage(), grid_, state, prev_well_state, extra);
initHydroCarbonState(state, props_.phaseUsage(), Opm::UgGridHelpers::numCells(grid_), has_disgas_, has_vapoil_);
initHysteresisParams(state);
}
// Create timers and file for writing timing info.
@ -844,4 +845,25 @@ namespace Opm
well_state, list_econ_limited);
}
template <class Implementation>
void
SimulatorBase<Implementation>::
initHysteresisParams(ReservoirState& state)
{
typedef std::vector<double> VectorType;
const VectorType& somax = state.getCellData( "SOMAX" );
VectorType& pcSwMdc_ow = state.getCellData( "PCSWMDC_OW" );
VectorType& krnSwMdc_ow = state.getCellData( "KRNSWMDC_OW" );
VectorType& pcSwMdc_go = state.getCellData( "PCSWMDC_GO" );
VectorType& krnSwMdc_go = state.getCellData( "KRNSWMDC_GO" );
props_.setSatOilMax(somax);
props_.setOilWaterHystParams(pcSwMdc_ow, krnSwMdc_ow, allcells_);
props_.setGasOilHystParams(pcSwMdc_go, krnSwMdc_go, allcells_);
}
} // namespace Opm

View File

@ -151,17 +151,7 @@ public:
// This is a restart, populate WellState and ReservoirState state objects from restart file
output_writer_.initFromRestartFile(props_.phaseUsage(), grid(), state, prev_well_state, extra);
initHydroCarbonState(state, props_.phaseUsage(), Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
{
const int num_cells = Opm::UgGridHelpers::numCells(grid());
typedef std::vector<double> VectorType;
const VectorType& somax = state.getCellData( "SOMAX" );
for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
ebosSimulator_.model().setMaxOilSaturation(somax[cellIdx], cellIdx);
}
}
initHysteresisParams(state);
}
// Create timers and file for writing timing info.
@ -799,6 +789,40 @@ protected:
}
}
void initHysteresisParams(ReservoirState& state) {
const int num_cells = Opm::UgGridHelpers::numCells(grid());
typedef std::vector<double> VectorType;
const VectorType& somax = state.getCellData( "SOMAX" );
for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
ebosSimulator_.model().setMaxOilSaturation(somax[cellIdx], cellIdx);
}
if (ebosSimulator_.problem().materialLawManager()->enableHysteresis()) {
auto matLawManager = ebosSimulator_.problem().materialLawManager();
VectorType& pcSwMdc_ow = state.getCellData( "PCSWMDC_OW" );
VectorType& krnSwMdc_ow = state.getCellData( "KRNSWMDC_OW" );
VectorType& pcSwMdc_go = state.getCellData( "PCSWMDC_GO" );
VectorType& krnSwMdc_go = state.getCellData( "KRNSWMDC_GO" );
for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
matLawManager->setOilWaterHysteresisParams(
pcSwMdc_ow[cellIdx],
krnSwMdc_ow[cellIdx],
cellIdx);
matLawManager->setGasOilHysteresisParams(
pcSwMdc_go[cellIdx],
krnSwMdc_go[cellIdx],
cellIdx);
}
}
}
// Data.
Simulator& ebosSimulator_;

View File

@ -412,9 +412,15 @@ namespace Opm
{"SGAS" , UnitSystem::measure::identity},
{"TEMP" , UnitSystem::measure::temperature},
{"RS" , UnitSystem::measure::gas_oil_ratio},
{"RV" , UnitSystem::measure::oil_gas_ratio},
{"SOMAX", UnitSystem::measure::identity}};
std::map<std::string, bool> extra_keys {{"OPMEXTRA" , false}};
{"RV" , UnitSystem::measure::oil_gas_ratio}};
std::map<std::string, bool> extra_keys {
{"OPMEXTRA" , false},
{"SOMAX", false},
{"PCSWM_OW", false},
{"KRNSW_OW", false},
{"PCSWM_GO", false},
{"KRNSW_GO", false}
};
if (restart_double_si_) {
// Avoid any unit conversions, treat restart input as SI units.
@ -446,7 +452,7 @@ namespace Opm
wellstate.resize(wells, simulatorstate, phaseUsage ); //Resize for restart step
auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys);
solutionToSim( restart_values.solution, phaseUsage, simulatorstate );
solutionToSim( restart_values.solution, restart_values.extra, phaseUsage, simulatorstate );
wellsToState( restart_values.wells, phaseUsage, wellstate );
const auto opmextra_iter = restart_values.extra.find("OPMEXTRA");
@ -540,6 +546,10 @@ namespace Opm
addToSimData( simData, "RVSAT", sd.rvSat );
addToSimData( simData, "SOMAX", sd.soMax );
addToSimData( simData, "PCSWMDC_OW", sd.pcswmdc_ow);
addToSimData( simData, "KRNSWMDC_OW", sd.krnswdc_ow);
addToSimData( simData, "PCSWMDC_GO", sd.pcswmdc_go);
addToSimData( simData, "KRNSWMDC_GO", sd.krnswdc_go);
return simData;
}
@ -772,6 +782,32 @@ namespace Opm
data::TargetType::RESTART_AUXILIARY);
}
if (sd.hasCellData("PCSWMDC_OW")) {
output.insert("PCSWM_OW", //FIXME: Eight-long variable name
Opm::UnitSystem::measure::identity,
std::move( sd.getCellData("PCSWMDC_OW") ),
data::TargetType::RESTART_AUXILIARY);
}
if (sd.hasCellData("KRNSWMDC_OW")) {
output.insert("KRNSW_OW",
Opm::UnitSystem::measure::identity,
std::move( sd.getCellData("KRNSWMDC_OW") ),
data::TargetType::RESTART_AUXILIARY);
}
if (sd.hasCellData("PCSWMDC_GO")) {
output.insert("PCSWM_GO", //FIXME: Eight-long variable name
Opm::UnitSystem::measure::identity,
std::move( sd.getCellData("PCSWMDC_GO") ),
data::TargetType::RESTART_AUXILIARY);
}
if (sd.hasCellData("KRNSWMDC_GO")) {
output.insert("KRNSW_GO",
Opm::UnitSystem::measure::identity,
std::move( sd.getCellData("KRNSWMDC_GO") ),
data::TargetType::RESTART_AUXILIARY);
}
//Warn for any unhandled keyword
if (log) {
for (auto& keyValue : rstKeywords) {

View File

@ -82,7 +82,11 @@ namespace Opm {
: rq(num_phases)
, rsSat(ADB::null())
, rvSat(ADB::null())
, soMax()
, soMax() // FIXME: Not handled properly
, krnswdc_ow() // FIXME: Not handled properly
, krnswdc_go() // FIXME: Not handled properly
, pcswmdc_ow() // FIXME: Not handled properly
, pcswmdc_go() // FIXME: Not handled properly
, fip()
{
}
@ -94,6 +98,10 @@ namespace Opm {
ADB rsSat;
ADB rvSat;
std::vector<double> soMax;
std::vector<double> krnswdc_ow;
std::vector<double> krnswdc_go;
std::vector<double> pcswmdc_ow;
std::vector<double> pcswmdc_go;
std::array<V, fipValues> fip;
};