diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index c46205c61..5b7c66a3c 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -105,6 +105,12 @@ namespace Opm { std::vector soMax; // Maximum oil saturation + //Hysteresis parameters + std::vector krnswdc_ow; + std::vector krnswdc_go; + std::vector pcswmdc_ow; + std::vector pcswmdc_go; + std::array fip; }; diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 1d425f9c7..1e5d82792 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -464,6 +464,10 @@ typedef Eigen::ArrayoilWaterHysteresisParams( - pcSwMdc_ow[cellIdx], - krnSwMdc_ow[cellIdx], - cellIdx); - matLawManager->gasOilHysteresisParams( - pcSwMdc_go[cellIdx], - krnSwMdc_go[cellIdx], - 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(); diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.cpp b/opm/autodiff/BlackoilPropsAdFromDeck.cpp index 17492e098..9c3c481dd 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.cpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.cpp @@ -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& pcswmdc, + const std::vector& krnswdc, + const std::vector& 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& pcswmdc, + std::vector& krnswdc, + const std::vector& 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& pcswmdc, + const std::vector& krnswdc, + const std::vector& 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& pcswmdc, + std::vector& krnswdc, + const std::vector& 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& saturation) { @@ -895,6 +947,11 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& return satOilMax_; } + void BlackoilPropsAdFromDeck::setSatOilMax(const std::vector& 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. diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.hpp b/opm/autodiff/BlackoilPropsAdFromDeck.hpp index 41277b2bd..9293081ca 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.hpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.hpp @@ -359,12 +359,45 @@ namespace Opm void updateSatHyst(const std::vector& saturation, const std::vector& 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& pcswmdc, + const std::vector& krnswdc, + const std::vector& 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& pcswmdc, + std::vector& krnswdc, + const std::vector& 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& pcswmdc, + const std::vector& krnswdc, + const std::vector& 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& pcswmdc, + std::vector& krnswdc, + const std::vector& cells) const; + /// Update for max oil saturation. + /// \param[in] saturation Saturations for all phases void updateSatOilMax(const std::vector& saturation); /// Returns the max oil saturation const std::vector& 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& 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. diff --git a/opm/autodiff/SimulatorBase.hpp b/opm/autodiff/SimulatorBase.hpp index 2a09f3506..eaffbd599 100644 --- a/opm/autodiff/SimulatorBase.hpp +++ b/opm/autodiff/SimulatorBase.hpp @@ -186,6 +186,7 @@ namespace Opm const WellState& well_state, DynamicListEconLimited& list_econ_limited) const; + void initHysteresisParams(ReservoirState& state); // Data. typedef RateConverter:: diff --git a/opm/autodiff/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index e3fc5a508..e7d35a121 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -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 + void + SimulatorBase:: + initHysteresisParams(ReservoirState& state) + { + typedef std::vector 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 diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 63c677243..de2a1089c 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -797,22 +797,25 @@ protected: ebosSimulator_.model().setMaxOilSaturation(somax[cellIdx], cellIdx); } - 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) { + if (ebosSimulator_.problem().materialLawManager()->enableHysteresis()) { auto matLawManager = ebosSimulator_.problem().materialLawManager(); - matLawManager->setOilWaterHysteresisParams( - pcSwMdc_ow[cellIdx], - krnSwMdc_ow[cellIdx], - cellIdx); - matLawManager->setGasOilHysteresisParams( - pcSwMdc_go[cellIdx], - krnSwMdc_go[cellIdx], - cellIdx); + + 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); + } } } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index 2c0156a28..3704e5de2 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -556,6 +556,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; } diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp index 5ef511292..0661e4723 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp @@ -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 soMax; + std::vector krnswdc_ow; + std::vector krnswdc_go; + std::vector pcswmdc_ow; + std::vector pcswmdc_go; std::array fip; };