OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-10-12 08:54:53 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1775 from andlaus/update_ewoms_simulator_in_flow

Browse Source 
flow: make sure that the ebosSimulator's episode machinery is up to date
This commit is contained in:
Tor Harald Sandve authored and GitHub committed 2019-05-07 10:51:52 +02:00
2 parents e19508e599 dd74fe1027
commit 71cd83bb2b
18 files changed
+288 -228
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

ebos/eclproblem.hh
+191 -178
View File
@@ -582,15 +582,17 @@ public:
const auto& vanguard = simulator.vanguard();
SolventModule::initFromDeck(vanguard.deck(), vanguard.eclState());
PolymerModule::initFromDeck(vanguard.deck(), vanguard.eclState());
if (EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput))
// create the ECL writer
eclWriter_.reset(new EclWriterType(simulator));
// create the ECL writer
eclWriter_.reset(new EclWriterType(simulator));
if (enableExperiments)
enableDriftCompensation_ = EWOMS_GET_PARAM(TypeTag, bool, EclEnableDriftCompensation);
else
enableDriftCompensation_ = false;
enableEclOutput_ = EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput);
enableTuning_ = EWOMS_GET_PARAM(TypeTag, bool, EclEnableTuning);
initialTimeStepSize_ = EWOMS_GET_PARAM(TypeTag, Scalar, InitialTimeStepSize);
minTimeStepSize_ = EWOMS_GET_PARAM(TypeTag, Scalar, MinTimeStepSize);
@@ -608,6 +610,19 @@ public:
ParentType::finishInit();
auto& simulator = this->simulator();
const auto& eclState = simulator.vanguard().eclState();
const auto& schedule = simulator.vanguard().schedule();
const auto& timeMap = schedule.getTimeMap();
// Set the start time of the simulation
simulator.setStartTime(timeMap.getStartTime(/*reportStepIdx=*/0));
// We want the episode index to be the same as the report step index to make
// things simpler, so we have to set the episode index to -1 because it is
// incremented by endEpisode(). The size of the initial time step and
// length of the initial episode is set to zero for the same reason.
simulator.setEpisodeIndex(-1);
simulator.setEpisodeLength(0.0);
// the "NOGRAV" keyword from Frontsim or setting the EnableGravity to false
// disables gravity, else the standard value of the gravity constant at sea level
@@ -632,7 +647,6 @@ public:
}
// deal with DRSDT
const auto& eclState = simulator.vanguard().eclState();
unsigned ntpvt = eclState.runspec().tabdims().getNumPVTTables();
maxDRs_.resize(ntpvt, 1e30);
dRsDtOnlyFreeGas_.resize(ntpvt, false);
@@ -649,34 +663,17 @@ public:
readThermalParameters_();
transmissibilities_.finishInit();
const auto& timeMap = simulator.vanguard().schedule().getTimeMap();
readInitialCondition_();
// Set the start time of the simulation
simulator.setStartTime(timeMap.getStartTime(/*timeStepIdx=*/0));
// We want the episode index to be the same as the report step index to make
// things simpler, so we have to set the episode index to -1 because it is
// incremented inside beginEpisode(). The size of the initial time step and
// length of the initial episode is set to zero for the same reason.
simulator.setEpisodeIndex(-1);
simulator.setEpisodeLength(0.0);
simulator.setTimeStepSize(0.0);
updatePffDofData_();
if (GET_PROP_VALUE(TypeTag, EnablePolymer)) {
const auto& vanguard = this->simulator().vanguard();
const auto& vanguard = simulator.vanguard();
const auto& gridView = vanguard.gridView();
int numElements = gridView.size(/*codim=*/0);
maxPolymerAdsorption_.resize(numElements, 0.0);
}
if (eclWriter_) {
eclWriter_->writeInit();
this->simulator().vanguard().releaseGlobalTransmissibilities();
}
tracerModel_.init();
readBoundaryConditions_();
@@ -689,6 +686,16 @@ public:
if (enableExperiments)
checkDeckCompatibility_();
// write the static output files (EGRID, INIT, SMSPEC, etc.)
if (enableEclOutput_)
eclWriter_->writeInit();
simulator.vanguard().releaseGlobalTransmissibilities();
// after finishing the initialization and writing the initial solution, we move
// to the first "real" episode/report step
simulator.startNextEpisode(timeMap.getTimeStepLength(0));
simulator.setEpisodeIndex(0);
}
void prefetch(const Element& elem) const
@@ -737,37 +744,21 @@ public:
void beginEpisode(bool isOnRestart = false)
{
// Proceed to the next report step
Simulator& simulator = this->simulator();
auto& eclState = this->simulator().vanguard().eclState();
const auto& schedule = this->simulator().vanguard().schedule();
auto& simulator = this->simulator();
auto& eclState = simulator.vanguard().eclState();
const auto& schedule = simulator.vanguard().schedule();
const auto& events = schedule.getEvents();
const auto& timeMap = schedule.getTimeMap();
int episodeIdx = simulator.episodeIndex();
// The first thing to do in the morning of an episode is update update the
// eclState and the deck if they need to be changed.
int nextEpisodeIdx = simulator.episodeIndex();
if (enableExperiments && this->gridView().comm().rank() == 0) {
boost::posix_time::ptime curDateTime =
boost::posix_time::from_time_t(timeMap.getStartTime(nextEpisodeIdx+1));
std::cout << "Report step " << nextEpisodeIdx + 2
<< "/" << timeMap.numTimesteps()
<< " at day " << timeMap.getTimePassedUntil(nextEpisodeIdx+1)/(24*3600)
<< "/" << timeMap.getTotalTime()/(24*3600)
<< ", date = " << curDateTime.date()
<< "\n ";
}
if (nextEpisodeIdx > 0 &&
events.hasEvent(Opm::ScheduleEvents::GEO_MODIFIER, nextEpisodeIdx))
{
if (episodeIdx >= 0 && events.hasEvent(Opm::ScheduleEvents::GEO_MODIFIER, episodeIdx)) {
// bring the contents of the keywords to the current state of the SCHEDULE
// section
// section.
//
// TODO (?): make grid topology changes possible (depending on what exactly
// has changed, the grid may need be re-created which has some serious
// implications on e.g., the solution of the simulation.)
const auto& miniDeck = schedule.getModifierDeck(nextEpisodeIdx);
const auto& miniDeck = schedule.getModifierDeck(episodeIdx);
eclState.applyModifierDeck(miniDeck);
// re-compute all quantities which may possibly be affected.
@@ -776,44 +767,32 @@ public:
updatePffDofData_();
}
// Opm::TimeMap deals with points in time, so the number of time intervals (i.e.,
// report steps) is one less!
int numReportSteps = timeMap.size() - 1;
// start the next episode if there are additional report steps, else finish the
// simulation
while (nextEpisodeIdx < numReportSteps &&
simulator.time() >= timeMap.getTimePassedUntil(nextEpisodeIdx + 1)*(1 - 1e-10))
{
++ nextEpisodeIdx;
if (enableExperiments && this->gridView().comm().rank() == 0 && episodeIdx >= 0) {
// print some useful information in experimental mode. (the production
// simulator does this externally.)
boost::posix_time::ptime curDateTime =
boost::posix_time::from_time_t(timeMap.getStartTime(episodeIdx));
std::cout << "Report step " << episodeIdx + 1
<< "/" << timeMap.numTimesteps()
<< " at day " << timeMap.getTimePassedUntil(episodeIdx)/(24*3600)
<< "/" << timeMap.getTotalTime()/(24*3600)
<< ", date = " << curDateTime.date()
<< "\n ";
}
// react to TUNING changes
bool tuningEvent = false;
if (nextEpisodeIdx > 0
&& enableTuning_
&& events.hasEvent(Opm::ScheduleEvents::TUNING_CHANGE, nextEpisodeIdx))
if (episodeIdx > 0 && enableTuning_ && events.hasEvent(Opm::ScheduleEvents::TUNING_CHANGE, episodeIdx))
{
const auto& tuning = schedule.getTuning();
initialTimeStepSize_ = tuning.getTSINIT(nextEpisodeIdx);
maxTimeStepAfterWellEvent_ = tuning.getTMAXWC(nextEpisodeIdx);
maxTimeStepSize_ = tuning.getTSMAXZ(nextEpisodeIdx);
restartShrinkFactor_ = 1./tuning.getTSFCNV(nextEpisodeIdx);
minTimeStepSize_ = tuning.getTSMINZ(nextEpisodeIdx);
initialTimeStepSize_ = tuning.getTSINIT(episodeIdx);
maxTimeStepAfterWellEvent_ = tuning.getTMAXWC(episodeIdx);
maxTimeStepSize_ = tuning.getTSMAXZ(episodeIdx);
restartShrinkFactor_ = 1./tuning.getTSFCNV(episodeIdx);
minTimeStepSize_ = tuning.getTSMINZ(episodeIdx);
tuningEvent = true;
}
Scalar episodeLength = timeMap.getTimeStepLength(nextEpisodeIdx);
simulator.startNextEpisode(episodeLength);
Scalar dt = limitNextTimeStepSize_(episodeLength);
if (nextEpisodeIdx == 0 || tuningEvent)
// allow the size of the initial time step to be set via an external parameter
// if TUNING is enabled, also limit the time step size after a tuning event to TSINIT
dt = std::min(dt, initialTimeStepSize_);
simulator.setTimeStepSize(dt);
const bool invalidateFromHyst = updateHysteresis_();
const bool invalidateFromMaxOilSat = updateMaxOilSaturation_();
const bool doInvalidate = invalidateFromHyst || invalidateFromMaxOilSat;
@@ -822,14 +801,19 @@ public:
updateMaxPolymerAdsorption_();
// set up the wells for the next episode.
//
// TODO: the first two arguments seem to be unnecessary
wellModel_.beginEpisode(this->simulator().vanguard().eclState(),
this->simulator().vanguard().schedule(),
isOnRestart);
wellModel_.beginEpisode();
// set up the aquifers for the next episode.
aquiferModel_.beginEpisode();
// set the size of the initial time step of the episode
Scalar dt = limitNextTimeStepSize_(simulator.episodeLength());
if (episodeIdx == 0 || tuningEvent)
// allow the size of the initial time step to be set via an external parameter
// if TUNING is enabled, also limit the time step size after a tuning event to TSINIT
dt = std::min(dt, initialTimeStepSize_);
simulator.setTimeStepSize(dt);
if (doInvalidate)
this->model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
}
@@ -839,17 +823,18 @@ public:
*/
void beginTimeStep()
{
int epsiodeIdx = this->simulator().episodeIndex();
const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
const auto& simulator = this->simulator();
int epsiodeIdx = simulator.episodeIndex();
const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
if (drsdtActive_())
// DRSDT is enabled
for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx)
maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*this->simulator().timeStepSize();
maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*simulator.timeStepSize();
if (drvdtActive_())
// DRVDT is enabled
for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx)
maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*this->simulator().timeStepSize();
maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*simulator.timeStepSize();
wellModel_.beginTimeStep();
aquiferModel_.beginTimeStep();
@@ -896,10 +881,13 @@ public:
// the right thing (i.e., the mass change inside the whole reservoir must be
// equivalent to the fluxes over the grid's boundaries plus the source rates
// specified by the problem)
std::cout << "checking conservativeness of solution\n";
this->model().checkConservativeness(/*tolerance=*/-1, /*verbose=*/true);
std::cout << "solution is sufficiently conservative\n";
}
#endif // NDEBUG
const auto& simulator = this->simulator();
wellModel_.endTimeStep();
aquiferModel_.endTimeStep();
tracerModel_.endTimeStep();
@@ -911,7 +899,7 @@ public:
const auto& residual = this->model().linearizer().residual();
for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
drift_[globalDofIdx] = residual[globalDofIdx];
drift_[globalDofIdx] *= this->simulator().timeStepSize();
drift_[globalDofIdx] *= simulator.timeStepSize();
if (GET_PROP_VALUE(TypeTag, UseVolumetricResidual))
drift_[globalDofIdx] *= this->model().dofTotalVolume(globalDofIdx);
}
@@ -925,15 +913,18 @@ public:
{
auto& simulator = this->simulator();
const auto& schedule = simulator.vanguard().schedule();
const auto& timeMap = schedule.getTimeMap();
int episodeIdx = simulator.episodeIndex();
const auto& timeMap = schedule.getTimeMap();
int numReportSteps = timeMap.size() - 1;
if (episodeIdx + 1 >= numReportSteps) {
// check if we're finished ...
if (episodeIdx + 1 >= static_cast<int>(timeMap.numTimesteps())) {
simulator.setFinished(true);
return;
}
// .. if we're not yet done, start the next episode (report step)
simulator.startNextEpisode(timeMap.getTimeStepLength(episodeIdx + 1));
}
/*!
@@ -945,27 +936,25 @@ public:
*/
bool shouldWriteOutput() const
{
if (this->simulator().timeStepIndex() < 0)
const auto& simulator = this->simulator();
if (simulator.timeStepIndex() < 0)
// always write the initial solution
return true;
if (EWOMS_GET_PARAM(TypeTag, bool, EnableWriteAllSolutions))
return true;
return this->simulator().episodeWillBeOver();
return simulator.episodeWillBeOver();
}
/*!
* \brief Returns true if an eWoms restart file should be written to disk.
*
* The EclProblem does not write any restart files using the ad-hoc format, only ones
* using the ECL format.
*/
bool shouldWriteRestartFile() const
{
unsigned n = EWOMS_GET_PARAM(TypeTag, unsigned, RestartWritingInterval);
unsigned i = this->simulator().timeStepIndex();
if (i > 0 && (i%n) == 0)
return true; // we don't write a restart file for the initial condition
return false;
}
{ return false; }
/*!
* \brief Write the requested quantities of the current solution into the output
@@ -977,10 +966,8 @@ public:
// write the desired VTK files.
ParentType::writeOutput(isSubStep, verbose);
if (!eclWriter_)
return;
eclWriter_->writeOutput(isSubStep);
if (enableEclOutput_)
eclWriter_->writeOutput(isSubStep);
}
/*!
@@ -1443,21 +1430,20 @@ public:
*/
void initialSolutionApplied()
{
const auto& eclState = this->simulator().vanguard().eclState();
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& eclState = vanguard.eclState();
// initialize the wells. Note that this needs to be done after initializing the
// intrinsic permeabilities and the after applying the initial solution because
// the well model uses these...
wellModel_.init(eclState, this->simulator().vanguard().schedule());
wellModel_.init();
// let the object for threshold pressures initialize itself. this is done only at
// this point, because determining the threshold pressures may require to access
// the initial solution.
thresholdPressures_.finishInit();
// release the memory of the EQUIL grid since it's no longer needed after this point
this->simulator().vanguard().releaseEquilGrid();
updateCompositionChangeLimits_();
aquiferModel_.initialSolutionApplied();
@@ -1613,8 +1599,9 @@ public:
bool vapparsActive() const
{
int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0);
const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
const auto& simulator = this->simulator();
int epsiodeIdx = std::max(simulator.episodeIndex(), 0);
const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
return (oilVaporizationControl.getType() == Opm::OilVaporizationEnum::VAPPARS);
}
@@ -1633,7 +1620,8 @@ public:
if (this->nextTimeStepSize_ > 0.0)
return this->nextTimeStepSize_;
int episodeIdx = this->simulator().episodeIndex();
const auto& simulator = this->simulator();
int episodeIdx = simulator.episodeIndex();
// for the initial episode, we use a fixed time step size
if (episodeIdx < 0)
@@ -1643,7 +1631,6 @@ public:
// well the previous time step converged. After that, apply the runtime time
// stepping constraints.
const auto& newtonMethod = this->model().newtonMethod();
const auto& simulator = this->simulator();
return limitNextTimeStepSize_(newtonMethod.suggestTimeStepSize(simulator.timeStepSize()));
}
@@ -1745,15 +1732,17 @@ private:
bool drsdtActive_() const
{
int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0);
const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
const auto& simulator = this->simulator();
int epsiodeIdx = std::max(simulator.episodeIndex(), 0);
const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
return (oilVaporizationControl.drsdtActive());
}
bool drvdtActive_() const
{
int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0);
const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
const auto& simulator = this->simulator();
int epsiodeIdx = std::max(simulator.episodeIndex(), 0);
const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
return (oilVaporizationControl.drvdtActive());
}
@@ -1774,7 +1763,8 @@ private:
void updateElementDepths_()
{
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& gridView = vanguard.gridView();
const auto& elemMapper = this->elementMapper();;
@@ -1796,12 +1786,13 @@ private:
{
// update the "last Rs" values for all elements, including the ones in the ghost
// and overlap regions
int epsiodeIdx = std::max(this->simulator().episodeIndex(), 0);
const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
const auto& simulator = this->simulator();
int epsiodeIdx = std::max(simulator.episodeIndex(), 0);
const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
if (oilVaporizationControl.drsdtActive()) {
ElementContext elemCtx(this->simulator());
const auto& vanguard = this->simulator().vanguard();
ElementContext elemCtx(simulator);
const auto& vanguard = simulator.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
@@ -1830,8 +1821,8 @@ private:
// update the "last Rv" values for all elements, including the ones in the ghost
// and overlap regions
if (drvdtActive_()) {
ElementContext elemCtx(this->simulator());
const auto& vanguard = this->simulator().vanguard();
ElementContext elemCtx(simulator);
const auto& vanguard = simulator.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
@@ -1856,10 +1847,12 @@ private:
bool updateMaxOilSaturation_()
{
const auto& simulator = this->simulator();
// we use VAPPARS
if (vapparsActive()) {
ElementContext elemCtx(this->simulator());
const auto& vanguard = this->simulator().vanguard();
ElementContext elemCtx(simulator);
const auto& vanguard = simulator.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
@@ -1887,9 +1880,10 @@ private:
void readRockParameters_()
{
const auto& deck = this->simulator().vanguard().deck();
const auto& eclState = this->simulator().vanguard().eclState();
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
const auto& deck = simulator.vanguard().deck();
const auto& eclState = simulator.vanguard().eclState();
const auto& vanguard = simulator.vanguard();
// the ROCK keyword has not been specified, so we don't need
// to read rock parameters
@@ -1945,7 +1939,8 @@ private:
void readMaterialParameters_()
{
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& deck = vanguard.deck();
const auto& eclState = vanguard.eclState();
@@ -1980,7 +1975,8 @@ private:
if (!enableEnergy)
return;
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& deck = vanguard.deck();
const auto& eclState = vanguard.eclState();
@@ -1996,7 +1992,8 @@ private:
void updatePorosity_()
{
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& eclState = vanguard.eclState();
const auto& eclGrid = eclState.getInputGrid();
const auto& props = eclState.get3DProperties();
@@ -2044,7 +2041,7 @@ private:
// we define the porosity as the accumulated pore volume divided by the
// geometric volume of the element. Note that -- in pathetic cases -- it can
// be larger than 1.0!
Scalar dofVolume = this->simulator().model().dofTotalVolume(dofIdx);
Scalar dofVolume = simulator.model().dofTotalVolume(dofIdx);
assert(dofVolume > 0.0);
porosity_[dofIdx] = poreVolume/dofVolume;
}
@@ -2052,15 +2049,17 @@ private:
void initFluidSystem_()
{
const auto& deck = this->simulator().vanguard().deck();
const auto& eclState = this->simulator().vanguard().eclState();
const auto& simulator = this->simulator();
const auto& deck = simulator.vanguard().deck();
const auto& eclState = simulator.vanguard().eclState();
FluidSystem::initFromDeck(deck, eclState);
}
void readInitialCondition_()
{
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& deck = vanguard.deck();
if (!deck.hasKeyword("EQUIL"))
@@ -2074,9 +2073,11 @@ private:
void readEquilInitialCondition_()
{
const auto& simulator = this->simulator();
// initial condition corresponds to hydrostatic conditions.
typedef Ewoms::EclEquilInitializer<TypeTag> EquilInitializer;
EquilInitializer equilInitializer(this->simulator(), *materialLawManager_);
EquilInitializer equilInitializer(simulator, *materialLawManager_);
size_t numElems = this->model().numGridDof();
initialFluidStates_.resize(numElems);
@@ -2089,17 +2090,22 @@ private:
void readEclRestartSolution_()
{
// Set the start time of the simulation
const auto& schedule = this->simulator().vanguard().schedule();
const auto& eclState = this->simulator().vanguard().eclState();
auto& simulator = this->simulator();
const auto& schedule = simulator.vanguard().schedule();
const auto& eclState = simulator.vanguard().eclState();
const auto& timeMap = schedule.getTimeMap();
const auto& initconfig = eclState.getInitConfig();
int episodeIdx = initconfig.getRestartStep() - 1;
this->simulator().setStartTime(timeMap.getStartTime(/*timeStepIdx=*/0));
this->simulator().setTime(timeMap.getTimePassedUntil(episodeIdx));
this->simulator().setEpisodeIndex(episodeIdx);
this->simulator().setEpisodeLength(timeMap.getTimeStepLength(episodeIdx));
this->simulator().setTimeStepSize(eclWriter_->restartTimeStepSize());
simulator.setStartTime(timeMap.getStartTime(/*timeStepIdx=*/0));
simulator.setTime(timeMap.getTimePassedUntil(episodeIdx));
simulator.startNextEpisode(simulator.startTime() + simulator.time(),
timeMap.getTimeStepLength(episodeIdx));
simulator.setEpisodeIndex(episodeIdx);
Scalar dt = std::min(eclWriter_->restartTimeStepSize(), simulator.episodeLength());
simulator.setTimeStepSize(dt);
eclWriter_->beginRestart();
@@ -2125,7 +2131,7 @@ private:
for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
auto& elemFluidState = initialFluidStates_[elemIdx];
elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
eclWriter_->eclOutputModule().initHysteresisParams(this->simulator(), elemIdx);
eclWriter_->eclOutputModule().initHysteresisParams(simulator, elemIdx);
eclWriter_->eclOutputModule().assignToFluidState(elemFluidState, elemIdx);
processRestartSaturations_(elemFluidState);
@@ -2139,17 +2145,17 @@ private:
// if we need to restart for polymer molecular weight simulation, we need to add related here
}
const int epsiodeIdx = this->simulator().episodeIndex();
const auto& oilVaporizationControl = this->simulator().vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
const int epsiodeIdx = simulator.episodeIndex();
const auto& oilVaporizationControl = simulator.vanguard().schedule().getOilVaporizationProperties(epsiodeIdx);
if (drsdtActive_())
// DRSDT is enabled
for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRs_.size(); ++pvtRegionIdx)
maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*this->simulator().timeStepSize();
maxDRs_[pvtRegionIdx] = oilVaporizationControl.getMaxDRSDT(pvtRegionIdx)*simulator.timeStepSize();
if (drvdtActive_())
// DRVDT is enabled
for (size_t pvtRegionIdx = 0; pvtRegionIdx < maxDRv_.size(); ++pvtRegionIdx)
maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*this->simulator().timeStepSize();
maxDRv_[pvtRegionIdx] = oilVaporizationControl.getMaxDRVDT(pvtRegionIdx)*simulator.timeStepSize();
if (tracerModel().numTracers() > 0)
std::cout << "Warning: Restart is not implemented for the tracer model, it will initialize with initial tracer concentration" << std::endl;
@@ -2159,7 +2165,7 @@ private:
// need to be correct for stuff like boundary conditions.
auto& sol = this->model().solution(/*timeIdx=*/0);
const auto& gridView = this->gridView();
ElementContext elemCtx(this->simulator());
ElementContext elemCtx(simulator);
auto elemIt = gridView.template begin</*codim=*/0>();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
@@ -2183,7 +2189,7 @@ private:
initialFluidStates_ = tmpInitialFs;
// make sure that the stuff which needs to be done at the beginning of an episode
// is run.
this->beginEpisode(/*isOnRestart=*/true);
this->beginEpisode();
eclWriter_->endRestart();
}
@@ -2215,7 +2221,8 @@ private:
void readExplicitInitialCondition_()
{
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& eclState = vanguard.eclState();
const auto& eclProps = eclState.get3DProperties();
@@ -2241,7 +2248,7 @@ private:
initialFluidStates_.resize(numDof);
const auto& cartSize = this->simulator().vanguard().cartesianDimensions();
const auto& cartSize = simulator.vanguard().cartesianDimensions();
size_t numCartesianCells = cartSize[0] * cartSize[1] * cartSize[2];
std::vector<double> waterSaturationData;
@@ -2359,7 +2366,8 @@ private:
void readBlackoilExtentionsInitialConditions_()
{
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
const auto& eclState = vanguard.eclState();
size_t numDof = this->model().numGridDof();
@@ -2405,8 +2413,9 @@ private:
// we need to update the hysteresis data for _all_ elements (i.e., not just the
// interior ones) to avoid desynchronization of the processes in the parallel case!
ElementContext elemCtx(this->simulator());
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
ElementContext elemCtx(simulator);
const auto& vanguard = simulator.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
@@ -2425,8 +2434,9 @@ private:
void updateMaxPolymerAdsorption_()
{
// we need to update the max polymer adsoption data for all elements
ElementContext elemCtx(this->simulator());
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
ElementContext elemCtx(simulator);
const auto& vanguard = simulator.vanguard();
auto elemIt = vanguard.gridView().template begin</*codim=*/0>();
const auto& elemEndIt = vanguard.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++elemIt) {
@@ -2444,14 +2454,15 @@ private:
void updatePvtnum_()
{
const auto& eclState = this->simulator().vanguard().eclState();
const auto& simulator = this->simulator();
const auto& eclState = simulator.vanguard().eclState();
const auto& eclProps = eclState.get3DProperties();
if (!eclProps.hasDeckIntGridProperty("PVTNUM"))
return;
const auto& pvtnumData = eclProps.getIntGridProperty("PVTNUM").getData();
const auto& vanguard = this->simulator().vanguard();
const auto& vanguard = simulator.vanguard();
unsigned numElems = vanguard.gridView().size(/*codim=*/0);
pvtnum_.resize(numElems);
@@ -2463,14 +2474,15 @@ private:
void updateSatnum_()
{
const auto& eclState = this->simulator().vanguard().eclState();
const auto& simulator = this->simulator();
const auto& eclState = simulator.vanguard().eclState();
const auto& eclProps = eclState.get3DProperties();
if (!eclProps.hasDeckIntGridProperty("SATNUM"))
return;
const auto& satnumData = eclProps.getIntGridProperty("SATNUM").getData();
const auto& vanguard = this->simulator().vanguard();
const auto& vanguard = simulator.vanguard();
unsigned numElems = vanguard.gridView().size(/*codim=*/0);
satnum_.resize(numElems);
@@ -2482,14 +2494,15 @@ private:
void updateMiscnum_()
{
const auto& eclState = this->simulator().vanguard().eclState();
const auto& simulator = this->simulator();
const auto& eclState = simulator.vanguard().eclState();
const auto& eclProps = eclState.get3DProperties();
if (!eclProps.hasDeckIntGridProperty("MISCNUM"))
return;
const auto& miscnumData = eclProps.getIntGridProperty("MISCNUM").getData();
const auto& vanguard = this->simulator().vanguard();
const auto& vanguard = simulator.vanguard();
unsigned numElems = vanguard.gridView().size(/*codim=*/0);
miscnum_.resize(numElems);
@@ -2501,14 +2514,15 @@ private:
void updatePlmixnum_()
{
const auto& eclState = this->simulator().vanguard().eclState();
const auto& simulator = this->simulator();
const auto& eclState = simulator.vanguard().eclState();
const auto& eclProps = eclState.get3DProperties();
if (!eclProps.hasDeckIntGridProperty("PLMIXNUM"))
return;
const auto& plmixnumData = eclProps.getIntGridProperty("PLMIXNUM").getData();
const auto& vanguard = this->simulator().vanguard();
const auto& vanguard = simulator.vanguard();
unsigned numElems = vanguard.gridView().size(/*codim=*/0);
plmixnum_.resize(numElems);
@@ -2551,7 +2565,8 @@ private:
void readBoundaryConditions_()
{
nonTrivialBoundaryConditions_ = false;
const auto& vanguard = this->simulator().vanguard();
const auto& simulator = this->simulator();
const auto& vanguard = simulator.vanguard();
if (vanguard.deck().hasKeyword("BC")) {
nonTrivialBoundaryConditions_ = true;
@@ -2684,31 +2699,27 @@ private:
}
// this method applies the runtime constraints specified via the deck and/or command
// line parameters for the size of the next time step size.
// line parameters for the size of the next time step.
Scalar limitNextTimeStepSize_(Scalar dtNext) const
{
const auto& events = this->simulator().vanguard().schedule().getEvents();
int episodeIdx = this->simulator().episodeIndex();
if (!enableExperiments)
return dtNext;
const auto& simulator = this->simulator();
const auto& events = simulator.vanguard().schedule().getEvents();
int episodeIdx = simulator.episodeIndex();
// first thing in the morning, limit the time step size to the maximum size
dtNext = std::min(dtNext, maxTimeStepSize_);
// use at least slightly more than half of the maximum time step size by default
if (dtNext < maxTimeStepSize_ && maxTimeStepSize_ < dtNext*2)
dtNext = 1.01 * maxTimeStepSize_/2.0;
Scalar remainingEpisodeTime =
this->simulator().episodeStartTime()
+ this->simulator().episodeLength()
- this->simulator().time()
- this->simulator().timeStepSize();
// limit the time step size to the remaining time of the episode
dtNext = std::min(dtNext, remainingEpisodeTime);
simulator.episodeStartTime() + simulator.episodeLength()
- (simulator.startTime() + simulator.time());
assert(remainingEpisodeTime >= 0.0);
// if we would have a small amount of time left over in the current episode, make
// two equal time steps instead of a big and a small one
if (dtNext < remainingEpisodeTime*(1.0 - 1e-5) && dtNext*1.25 > remainingEpisodeTime)
if (remainingEpisodeTime/2.0 < dtNext && dtNext < remainingEpisodeTime*(1.0 - 1e-5))
// note: limiting to the maximum time step size here is probably not strictly
// necessary, but it should not hurt and is more fool-proof
dtNext = std::min(maxTimeStepSize_, remainingEpisodeTime/2.0);
@@ -2768,6 +2779,8 @@ private:
EclWellModel wellModel_;
EclAquiferModel aquiferModel_;
bool enableEclOutput_;
std::unique_ptr<EclWriterType> eclWriter_;
PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
ebos/eclthresholdpressure.hh
+14 -7
View File
@@ -91,10 +91,6 @@ public:
enableThresholdPressure_ = false;
}
void setFromRestart(const std::vector<Scalar>& values)
{ thpres_ = values; }
/*!
* \brief Actually compute the threshold pressures over a face as a pre-compute step.
*/
@@ -201,10 +197,21 @@ public:
return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx];
}
const std::vector<Scalar>& data() const {
return thpres_;
}
/*!
* \brief Return the raw array with the threshold pressures
*
* This is used for the restart capability.
*/
const std::vector<Scalar>& data() const
{ return thpres_; }
/*!
* \brief Set the threshold pressures from a raw array
*
* This is used for the restart capability.
*/
void setFromRestart(const std::vector<Scalar>& values)
{ thpres_ = values; }
private:
// compute the defaults of the threshold pressures using the initial condition
ebos/eclwellmanager.hh
+7 -4
View File
@@ -101,9 +101,10 @@ public:
*
* I.e., well positions, names etc...
*/
void init(const Opm::EclipseState& eclState OPM_UNUSED,
const Opm::Schedule& deckSchedule)
void init()
{
const Opm::Schedule& deckSchedule = simulator_.vanguard().schedule();
// create the wells which intersect with the current process' grid
for (size_t deckWellIdx = 0; deckWellIdx < deckSchedule.numWells(); ++deckWellIdx)
{
@@ -130,8 +131,10 @@ public:
* \brief This should be called the problem before each simulation
* episode to adapt the well controls.
*/
void beginEpisode(const Opm::EclipseState& eclState, const Opm::Schedule& deckSchedule, bool wasRestarted=false)
void beginEpisode(bool wasRestarted=false)
{
const Opm::EclipseState& eclState = simulator_.vanguard().eclState();
const Opm::Schedule& deckSchedule = simulator_.vanguard().schedule();
unsigned episodeIdx = simulator_.episodeIndex();
WellConnectionsMap wellCompMap;
@@ -589,7 +592,7 @@ public:
void deserialize(Restarter& res OPM_UNUSED)
{
// initialize the wells for the current episode
beginEpisode(simulator_.vanguard().eclState(), simulator_.vanguard().schedule(), /*wasRestarted=*/true);
beginEpisode(/*wasRestarted=*/true);
}
/*!
ebos/eclwriter.hh
+8 -8
View File
@@ -160,11 +160,11 @@ public:
// output using eclWriter if enabled
Opm::data::Wells localWellData = simulator_.problem().wellModel().wellData();
int episodeIdx = simulator_.episodeIndex() + 1;
int reportStepNum = simulator_.episodeIndex() + 1;
const auto& gridView = simulator_.vanguard().gridView();
int numElements = gridView.size(/*codim=*/0);
bool log = collectToIORank_.isIORank();
eclOutputModule_.allocBuffers(numElements, episodeIdx, isSubStep, log);
eclOutputModule_.allocBuffers(numElements, reportStepNum, isSubStep, log);
ElementContext elemCtx(simulator_);
ElementIterator elemIt = gridView.template begin</*codim=*/0>();
@@ -184,7 +184,7 @@ public:
// add cell data to perforations for Rft output
if (!isSubStep)
eclOutputModule_.addRftDataToWells(localWellData, episodeIdx);
eclOutputModule_.addRftDataToWells(localWellData, reportStepNum);
if (collectToIORank_.isParallel())
collectToIORank_.collect(localCellData, eclOutputModule_.getBlockData(), localWellData);
@@ -221,7 +221,7 @@ public:
// first, create a tasklet to write the data for the current time step to disk
auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(*eclIO_,
episodeIdx,
reportStepNum,
isSubStep,
curTime,
restartValue,
@@ -437,7 +437,7 @@ private:
: public TaskletInterface
{
Opm::EclipseIO& eclIO_;
int episodeIdx_;
int reportStepNum_;
bool isSubStep_;
double secondsElapsed_;
Opm::RestartValue restartValue_;
@@ -447,7 +447,7 @@ private:
bool writeDoublePrecision_;
explicit EclWriteTasklet(Opm::EclipseIO& eclIO,
int episodeIdx,
int reportStepNum,
bool isSubStep,
double secondsElapsed,
Opm::RestartValue restartValue,
@@ -456,7 +456,7 @@ private:
const std::map<std::pair<std::string, int>, double>& blockSummaryValues,
bool writeDoublePrecision)
: eclIO_(eclIO)
, episodeIdx_(episodeIdx)
, reportStepNum_(reportStepNum)
, isSubStep_(isSubStep)
, secondsElapsed_(secondsElapsed)
, restartValue_(restartValue)
@@ -469,7 +469,7 @@ private:
// callback to eclIO serial writeTimeStep method
void run()
{
eclIO_.writeTimeStep(episodeIdx_,
eclIO_.writeTimeStep(reportStepNum_,
isSubStep_,
secondsElapsed_,
restartValue_,
opm/autodiff/BlackoilModelEbos.hpp
+1 -7
View File
@@ -182,7 +182,6 @@ namespace Opm {
/// \param[in] timer simulation timer
void prepareStep(const SimulatorTimerInterface& timer)
{
// update the solution variables in ebos
if ( timer.lastStepFailed() ) {
ebosSimulator_.model().updateFailed();
@@ -194,13 +193,8 @@ namespace Opm {
// know the report step/episode index because of timing dependend data
// despide the fact that flow uses its own time stepper. (The length of the
// episode does not matter, though.)
Scalar t = timer.simulationTimeElapsed();
ebosSimulator_.startNextEpisode(/*episodeStartTime=*/t, /*episodeLength=*/1e30);
ebosSimulator_.setEpisodeIndex(timer.reportStepNum());
ebosSimulator_.setTime(t);
ebosSimulator_.setTime(timer.simulationTimeElapsed());
ebosSimulator_.setTimeStepSize(timer.currentStepLength());
ebosSimulator_.setTimeStepIndex(ebosSimulator_.timeStepIndex() + 1);
ebosSimulator_.problem().beginTimeStep();
unsigned numDof = ebosSimulator_.model().numGridDof();
opm/autodiff/BlackoilWellModel.hpp
+3 -12
View File
@@ -112,7 +112,7 @@ namespace Opm {
BlackoilWellModel(Simulator& ebosSimulator);
void init(const Opm::EclipseState& eclState, const Opm::Schedule& schedule);
void init();
/////////////
// <eWoms auxiliary module stuff>
@@ -153,18 +153,9 @@ namespace Opm {
// TODO (?)
}
void beginEpisode(const Opm::EclipseState& /* eclState */,
const Opm::Schedule& /* schedule */,
bool isRestart)
void beginEpisode()
{
size_t episodeIdx = ebosSimulator_.episodeIndex();
// beginEpisode in eclProblem advances the episode index
// we don't want this when we are at the beginning of an
// restart.
if (isRestart)
episodeIdx -= 1;
beginReportStep(episodeIdx);
beginReportStep(ebosSimulator_.episodeIndex());
}
void beginTimeStep();
opm/autodiff/BlackoilWellModel_impl.hpp
+3 -1
View File
@@ -38,8 +38,10 @@ namespace Opm {
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
init(const Opm::EclipseState& eclState, const Opm::Schedule& /* schedule */)
init()
{
const Opm::EclipseState& eclState = ebosSimulator_.vanguard().eclState();
gravity_ = ebosSimulator_.problem().gravity()[2];
extractLegacyCellPvtRegionIndex_();
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
+31 -8
View File
@@ -133,6 +133,8 @@ public:
{
failureReport_ = SimulatorReport();
ebosSimulator_.setEpisodeIndex(-1);
// handle restarts
std::unique_ptr<RestartValue> restartValues;
if (isRestart()) {
@@ -164,6 +166,21 @@ public:
double suggestedStepSize = -1.0;
if (isRestart()) {
// Set the start time of the simulation
const auto& schedule = ebosSimulator_.vanguard().schedule();
const auto& eclState = ebosSimulator_.vanguard().eclState();
const auto& timeMap = schedule.getTimeMap();
const auto& initconfig = eclState.getInitConfig();
int episodeIdx = initconfig.getRestartStep() - 1;
ebosSimulator_.setStartTime(timeMap.getStartTime(/*timeStepIdx=*/0));
ebosSimulator_.setTime(timeMap.getTimePassedUntil(episodeIdx));
ebosSimulator_.startNextEpisode(ebosSimulator_.startTime() + ebosSimulator_.time(),
timeMap.getTimeStepLength(episodeIdx));
ebosSimulator_.setEpisodeIndex(episodeIdx);
// This is a restart, determine the time step size from the restart data
if (restartValues->hasExtra("OPMEXTRA")) {
std::vector<double> opmextra = restartValues->getExtra("OPMEXTRA");
@@ -206,12 +223,18 @@ public:
Dune::Timer perfTimer;
perfTimer.start();
ebosSimulator_.setEpisodeIndex(-1);
ebosSimulator_.setEpisodeLength(0.0);
ebosSimulator_.setTimeStepSize(0.0);
wellModel_().beginReportStep(timer.currentStepNum());
ebosSimulator_.problem().writeOutput(false);
report.output_write_time += perfTimer.stop();
}
ebosSimulator_.setEpisodeIndex(timer.currentStepNum());
// Run a multiple steps of the solver depending on the time step control.
solverTimer.start();
@@ -265,6 +288,14 @@ public:
}
}
// write simulation state at the report stage
Dune::Timer perfTimer;
perfTimer.start();
const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
ebosSimulator_.problem().setNextTimeStepSize(nextstep);
ebosSimulator_.problem().writeOutput(false);
report.output_write_time += perfTimer.stop();
solver->model().endReportStep();
// take time that was used to solve system for this reportStep
@@ -284,14 +315,6 @@ public:
}
}
// write simulation state at the report stage
Dune::Timer perfTimer;
perfTimer.start();
const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
ebosSimulator_.problem().setNextTimeStepSize(nextstep);
ebosSimulator_.problem().writeOutput(false);
report.output_write_time += perfTimer.stop();
if (terminalOutput_) {
std::string msg =
"Time step took " + std::to_string(solverTimer.secsSinceStart()) + " seconds; "
tests/equil_base.DATA
+3
View File
@@ -113,3 +113,6 @@ PRESSURE
SUMMARY
SCHEDULE
TSTEP
1 /
tests/equil_capillary.DATA
+3 -1
View File
@@ -84,4 +84,6 @@ EQUIL
-------------------------------------
SCHEDULE
-- empty section
TSTEP
1 /
tests/equil_capillary_overlap.DATA
+3
View File
@@ -121,6 +121,9 @@ SEPARATE
SCHEDULE ======
TSTEP
1 /
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RESTART=3' 'NEWTON=2' /
tests/equil_capillary_swatinit.DATA
+4 -1
View File
@@ -99,4 +99,7 @@ RPTSOL
/
-------------------------------------
SCHEDULE
-- empty section
TSTEP
1 /
tests/equil_deadfluids.DATA
+4 -1
View File
@@ -85,4 +85,7 @@ EQUIL
-------------------------------------
SCHEDULE
-- empty section
TSTEP
1 /
tests/equil_livegas.DATA
+3
View File
@@ -124,6 +124,9 @@ SEPARATE
SCHEDULE ======
TSTEP
1 /
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
tests/equil_liveoil.DATA
+3
View File
@@ -133,6 +133,9 @@ SEPARATE
SCHEDULE ======
TSTEP
1 /
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
tests/equil_pbvd_and_pdvd.DATA
+3
View File
@@ -145,6 +145,9 @@ SEPARATE
SCHEDULE ======
TSTEP
1 /
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
tests/equil_rsvd_and_rvvd.DATA
+3
View File
@@ -144,6 +144,9 @@ SEPARATE
SCHEDULE ======
TSTEP
1 /
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
tests/test_ecl_output.cc
+1
View File
@@ -131,6 +131,7 @@ void test_summary()
simulator->model().applyInitialSolution();
Opm::data::Wells dw;
bool substep = false;
simulator->startNextEpisode(0.0, 1e30);
simulator->setEpisodeIndex(0);
eclWriter->writeOutput(substep);
simulator->setEpisodeIndex(1);