EclipseWriter: allow for writing of substeps in addition to report steps.

This commit is contained in:
Robert K 2014-12-10 11:20:29 +01:00
parent 1a430f76db
commit 10cffa770b
6 changed files with 254 additions and 55 deletions

View File

@ -40,6 +40,15 @@ struct MultiWriter : public OutputWriter {
}
}
virtual void writeTimeStep(const SimulatorTimer& timer,
const AdaptiveSimulatorTimer& substepTimer,
const SimulatorState& reservoirState,
const WellState& wellState) {
for (it_t it = writers_->begin (); it != writers_->end(); ++it) {
(*it)->writeTimeStep (timer, substepTimer, reservoirState, wellState);
}
}
private:
ptr_t writers_;
};

View File

@ -31,6 +31,7 @@ class EclipseState;
namespace parameter { class ParameterGroup; }
class SimulatorState;
class SimulatorTimer;
class AdaptiveSimulatorTimer;
class WellState;
struct PhaseUsage;
@ -79,8 +80,9 @@ public:
* \brief Write a blackoil reservoir state to disk for later inspection with
* visualization tools like ResInsight
*
* \param[in] timer The timer providing time, time step, etc. information
* \param[in] reservoirState The thermodynamic state of the reservoir
* \param[in] wellState The production/injection data for all wells
* \param[in] wellState The production/injection data for all wells
*
* This routine should be called after the timestep has been advanced,
* i.e. timer.currentStepNum () > 0.
@ -89,6 +91,23 @@ public:
const SimulatorState& reservoirState,
const WellState& wellState) = 0;
/*!
* \brief Write a blackoil reservoir state to disk for later inspection with
* visualization tools like ResInsight
*
* \param[in] timer The timer providing time, time step, etc. information
* \param[in] subStepTimer The timer providing sub step time information
* \param[in] reservoirState The thermodynamic state of the reservoir
* \param[in] wellState The production/injection data for all wells
*
* This routine should be called after the timestep has been advanced,
* i.e. timer.currentStepNum () > 0.
*/
virtual void writeTimeStep(const SimulatorTimer& timer,
const AdaptiveSimulatorTimer& subStepTimer,
const SimulatorState& reservoirState,
const WellState& wellState) = 0;
/*!
* Create a suitable set of output formats based on configuration.
*

View File

@ -28,6 +28,7 @@
#include <opm/core/grid/cpgpreprocess/preprocess.h>
#include <opm/core/simulator/SimulatorState.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/parameters/Parameter.hpp>
@ -125,6 +126,38 @@ int ertPhaseMask(const PhaseUsage uses)
}
// wrapper class to make EclipseWriter work with either the SimulatorTimer or a
// combination of SimulatorTimer and AdaptiveSimulatorTimer.
struct WriterTimer
{
const double time_;
const double stepLength_;
const time_t posixTime_;
const int deckReportStep_;
// copy values from SimulatorTimer
WriterTimer( const SimulatorTimer& timer )
: time_( timer.simulationTimeElapsed() ),
stepLength_( timer.stepLengthTaken() ),
posixTime_( timer.currentPosixTime() ),
deckReportStep_( timer.currentStepNum() )
{}
// copy values from SimulatorTimer and add values from
// AdaptiveSimulatorTimer, except the deck's reportStep
WriterTimer( const SimulatorTimer& timer, const AdaptiveSimulatorTimer& subStepTimer )
: time_( subStepTimer.simulationTimeElapsed() ),
stepLength_( subStepTimer.currentStepLength() ),
posixTime_( timer.currentPosixTime() + time_t(time_ - timer.simulationTimeElapsed()) ),
deckReportStep_( timer.currentStepNum() )
{}
int deckReportStep () const { return deckReportStep_; }
double simulationTimeElapsed() const { return time_; }
/// time elapsed since the start of the POSIX epoch (Jan 1st, 1970) [s].
time_t currentPosixTime() const { return posixTime_; }
double stepLengthTaken () const { return stepLength_; }
};
/**
@ -372,7 +405,7 @@ public:
ecl_rst_file_close(restartFileHandle_);
}
void writeHeader(const SimulatorTimer& timer,
void writeHeader(const WriterTimer& timer,
int reportStepIdx,
ecl_rsthead_type * rsthead_data)
{
@ -427,7 +460,7 @@ class Summary : private boost::noncopyable
public:
Summary(const std::string& outputDir,
const std::string& baseName,
const SimulatorTimer& timer,
const WriterTimer& timer,
int nx,
int ny,
int nz)
@ -464,7 +497,7 @@ public:
void addAllWells(Opm::EclipseStateConstPtr eclipseState,
const PhaseUsage& uses);
void writeTimeStep(int reportStepIdx,
const SimulatorTimer& timer,
const WriterTimer& timer,
const WellState& wellState);
ecl_sum_type *ertHandle() const
@ -482,7 +515,7 @@ class SummaryTimeStep : private boost::noncopyable
public:
SummaryTimeStep(Summary& summaryHandle,
int reportStepIdx,
const SimulatorTimer &timer)
const WriterTimer &timer)
{
ertHandle_ = ecl_sum_add_tstep(summaryHandle.ertHandle(),
reportStepIdx,
@ -537,7 +570,7 @@ public:
void writeHeader(int numCells,
const int* compressedToCartesianCellIdx,
const SimulatorTimer& timer,
const WriterTimer& timer,
Opm::EclipseStateConstPtr eclipseState,
const PhaseUsage uses)
{
@ -621,7 +654,8 @@ protected:
public:
/// Retrieve the value which the monitor is supposed to write to the summary file
/// according to the state of the well.
virtual double retrieveValue(const SimulatorTimer& timer,
virtual double retrieveValue(const int reportStepIdx,
const WriterTimer& timer,
const WellState& wellState,
const std::map<std::string, int>& nameToIdxMap) = 0;
@ -752,7 +786,8 @@ public:
"SM3/DAY" /* surf. cub. m. per day */)
{ }
virtual double retrieveValue(const SimulatorTimer& timer,
virtual double retrieveValue(const int reportStepIdx,
const WriterTimer& timer,
const WellState& wellState,
const std::map<std::string, int>& wellNameToIdxMap)
{
@ -763,7 +798,7 @@ public:
return 0.0;
}
if (well_->getStatus(timer.currentStepNum()) == WellCommon::SHUT) {
if (well_->getStatus(timer.deckReportStep()) == WellCommon::SHUT) {
// well is shut in the current time step
return 0.0;
}
@ -797,17 +832,18 @@ public:
, total_(0.)
{ }
virtual double retrieveValue(const SimulatorTimer& timer,
virtual double retrieveValue(const int reportStepIdx,
const WriterTimer& timer,
const WellState& wellState,
const std::map<std::string, int>& wellNameToIdxMap)
{
if (timer.currentStepNum() == 0) {
if (reportStepIdx == 0) {
// We are at the initial state.
// No step has been taken yet.
return 0.0;
}
if (well_->getStatus(timer.currentStepNum()) == WellCommon::SHUT) {
if (well_->getStatus(timer.deckReportStep()) == WellCommon::SHUT) {
// well is shut in the current time step
return 0.0;
}
@ -855,7 +891,8 @@ public:
"Pascal")
{ }
virtual double retrieveValue(const SimulatorTimer& timer,
virtual double retrieveValue(const int reportStepIdx,
const WriterTimer& timer,
const WellState& wellState,
const std::map<std::string, int>& wellNameToIdxMap)
{
@ -865,7 +902,7 @@ public:
// well not active in current time step
return 0.0;
}
if (well_->getStatus(timer.currentStepNum()) == WellCommon::SHUT) {
if (well_->getStatus(timer.deckReportStep()) == WellCommon::SHUT) {
// well is shut in the current time step
return 0.0;
}
@ -877,16 +914,16 @@ public:
// no inline implementation of this since it depends on the
// WellReport type being completed first
void Summary::writeTimeStep(int reportStepIdx,
const SimulatorTimer& timer,
const WriterTimer& timer,
const WellState& wellState)
{
// create a name -> well index map
const Opm::ScheduleConstPtr schedule = eclipseState_->getSchedule();
const auto& timeStepWells = schedule->getWells(reportStepIdx);
const auto& timeStepWells = schedule->getWells(timer.deckReportStep());
std::map<std::string, int> wellNameToIdxMap;
int openWellIdx = 0;
for (size_t tsWellIdx = 0; tsWellIdx < timeStepWells.size(); ++tsWellIdx) {
if (timeStepWells[tsWellIdx]->getStatus(timer.currentStepNum()) != WellCommon::SHUT ) {
if (timeStepWells[tsWellIdx]->getStatus(timer.deckReportStep()) != WellCommon::SHUT ) {
wellNameToIdxMap[timeStepWells[tsWellIdx]->name()] = openWellIdx;
openWellIdx++;
}
@ -898,7 +935,7 @@ void Summary::writeTimeStep(int reportStepIdx,
for (auto varIt = summaryReportVars_.begin(); varIt != summaryReportVars_.end(); ++varIt) {
ecl_sum_tstep_iset(tstep.ertHandle(),
smspec_node_get_params_index((*varIt)->ertHandle()),
(*varIt)->retrieveValue(timer, wellState, wellNameToIdxMap));
(*varIt)->retrieveValue(reportStepIdx, timer, wellState, wellNameToIdxMap));
}
// write the summary file to disk
@ -1014,6 +1051,12 @@ int EclipseWriter::eclipseWellStatusMask(WellCommon::StatusEnum wellStatus)
void EclipseWriter::writeInit(const SimulatorTimer &timer)
{
EclipseWriterDetails::WriterTimer writerTimer( timer );
writeInit( writerTimer );
}
void EclipseWriter::writeInit(const EclipseWriterDetails::WriterTimer &timer)
{
// if we don't want to write anything, this method becomes a
// no-op...
@ -1062,6 +1105,24 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer)
void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
const SimulatorState& reservoirState,
const WellState& wellState)
{
EclipseWriterDetails::WriterTimer writerTimer ( timer );
writeTimeStep( writerTimer, reservoirState, wellState );
}
void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
const AdaptiveSimulatorTimer& subStepTimer,
const SimulatorState& reservoirState,
const WellState& wellState)
{
EclipseWriterDetails::WriterTimer writerTimer ( timer, subStepTimer );
writeTimeStep( writerTimer, reservoirState, wellState );
}
// implementation of the writeTimeStep method
void EclipseWriter::writeTimeStep(const EclipseWriterDetails::WriterTimer& timer,
const SimulatorState& reservoirState,
const WellState& wellState)
{
// if we don't want to write anything, this method becomes a
// no-op...
@ -1075,7 +1136,7 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
}
std::vector<WellConstPtr> wells_ptr = eclipseState_->getSchedule()->getWells(timer.currentStepNum());
std::vector<WellConstPtr> wells_ptr = eclipseState_->getSchedule()->getWells(timer.deckReportStep());
std::vector<int> iwell_data;
std::vector<const char*> zwell_data;
std::vector<int> icon_data;
@ -1086,7 +1147,7 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
rsthead_data.nx = cartesianSize_[0];
rsthead_data.ny = cartesianSize_[1];
rsthead_data.nz = cartesianSize_[2];
rsthead_data.nwells = eclipseState_->getSchedule()->numWells(timer.currentStepNum());
rsthead_data.nwells = eclipseState_->getSchedule()->numWells(timer.deckReportStep());
rsthead_data.niwelz = 0;
rsthead_data.nzwelz = 0;
rsthead_data.niconz = 0;
@ -1098,10 +1159,10 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
for (std::vector<WellConstPtr>::const_iterator c_iter = wells_ptr.begin(); c_iter != wells_ptr.end(); ++c_iter) {
WellConstPtr well_ptr = *c_iter;
rsthead_data.ncwmax = eclipseState_->getSchedule()->getMaxNumCompletionsForWells(timer.currentStepNum());
restartHandle.addRestartFileIwelData(iwell_data, timer.currentStepNum(), well_ptr);
restartHandle.addRestartFileZwelData(zwell_data, timer.currentStepNum(), well_ptr);
restartHandle.addRestartFileIconData(icon_data, timer.currentStepNum(), rsthead_data.ncwmax, well_ptr);
rsthead_data.ncwmax = eclipseState_->getSchedule()->getMaxNumCompletionsForWells(timer.deckReportStep());
restartHandle.addRestartFileIwelData(iwell_data, timer.deckReportStep(), well_ptr);
restartHandle.addRestartFileZwelData(zwell_data, timer.deckReportStep(), well_ptr);
restartHandle.addRestartFileIconData(icon_data, timer.deckReportStep(), rsthead_data.ncwmax, well_ptr);
rsthead_data.niwelz = EclipseWriterDetails::Restart::NIWELZ;
rsthead_data.nzwelz = EclipseWriterDetails::Restart::NZWELZ;

View File

@ -1,6 +1,7 @@
/*
Copyright (c) 2013 Andreas Lauser
Copyright (c) 2013 Uni Research AS
Copyright (c) 2014 IRIS AS
This file is part of the Open Porous Media project (OPM).
@ -39,10 +40,12 @@ namespace Opm {
// forward declarations
namespace EclipseWriterDetails {
class Summary;
struct WriterTimer;
}
class SimulatorState;
class SimulatorTimer;
class AdaptiveSimulatorTimer;
class WellState;
namespace parameter { class ParameterGroup; }
@ -91,13 +94,36 @@ public:
* ERT or ECLIPSE. Note that calling this method is only
* meaningful after the first time step has been completed.
*
* \param[in] timer The timer providing time step and time information
* \param[in] reservoirState The thermodynamic state of the reservoir
* \param[in] wellState The production/injection data for all wells
* \param[in] wellState The production/injection data for all wells
*/
virtual void writeTimeStep(const SimulatorTimer& timer,
const SimulatorState& reservoirState,
const WellState& wellState);
/*!
* \brief Write a reservoir state and summary information to disk.
*
*
* The reservoir state can be inspected with visualization tools like
* ResInsight.
*
* The summary information can then be visualized using tools from
* ERT or ECLIPSE. Note that calling this method is only
* meaningful after the first time step has been completed.
*
* \param[in] timer The timer providing time step and time information
* \param[in] subStepTimer The timer providing sub step information
* \param[in] reservoirState The thermodynamic state of the reservoir
* \param[in] wellState The production/injection data for all wells
*/
virtual void writeTimeStep(const SimulatorTimer& timer,
const AdaptiveSimulatorTimer& subStepTimer,
const SimulatorState& reservoirState,
const WellState& wellState);
static int eclipseWellTypeMask(WellType wellType, WellInjector::TypeEnum injectorType);
static int eclipseWellStatusMask(WellCommon::StatusEnum wellStatus);
@ -117,6 +143,13 @@ private:
std::shared_ptr<EclipseWriterDetails::Summary> summary_;
void init(const parameter::ParameterGroup& params);
// implementation of writeInit
void writeInit(const EclipseWriterDetails::WriterTimer &timer);
// implementation of writeTimeStep
void writeTimeStep(const EclipseWriterDetails::WriterTimer& timer,
const SimulatorState& reservoirState,
const WellState& wellState);
};
typedef std::shared_ptr<EclipseWriter> EclipseWriterPtr;

View File

@ -6,33 +6,66 @@
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/simulator/TimeStepControlInterface.hpp>
namespace Opm {
// AdaptiveTimeStepping
// AdaptiveTimeStepping
//---------------------
class AdaptiveTimeStepping
{
public:
public:
//! \brief contructor taking parameter object
AdaptiveTimeStepping( const parameter::ParameterGroup& param );
/** \brief step method that acts like the solver::step method
in a sub cycle of time steps
in a sub cycle of time steps
\param solver solver object that must implement a method step( dt, state, well_state )
\param state current state of the solution variables
\param state current state of the solution variables
\param well_state additional well state object
\param time current simulation time
\param timestep current time step length that is to be sub cycled
*/
\param timestep current time step length that is to be sub cycled
*/
template <class Solver, class State, class WellState>
void step( Solver& solver, State& state, WellState& well_state,
const double time, const double timestep );
/** \brief step method that acts like the solver::step method
in a sub cycle of time steps
\param timer simulator timer providing time and timestep
\param solver solver object that must implement a method step( dt, state, well_state )
\param state current state of the solution variables
\param well_state additional well state object
*/
template <class Solver, class State, class WellState>
void step( const SimulatorTimer& timer,
Solver& solver, State& state, WellState& well_state );
/** \brief step method that acts like the solver::step method
in a sub cycle of time steps
\param timer simulator timer providing time and timestep
\param solver solver object that must implement a method step( dt, state, well_state )
\param state current state of the solution variables
\param well_state additional well state object
\param outputWriter writer object to write sub steps
*/
template <class Solver, class State, class WellState>
void step( const SimulatorTimer& timer,
Solver& solver, State& state, WellState& well_state,
OutputWriter& outputWriter );
protected:
template <class Solver, class State, class WellState>
void stepImpl( Solver& solver, State& state, WellState& well_state,
const double time, const double timestep,
const SimulatorTimer* timer, OutputWriter* outputWriter);
typedef std::unique_ptr< TimeStepControlInterface > TimeStepControlType;
TimeStepControlType timeStepControl_; //!< time step control object
@ -40,8 +73,8 @@ namespace Opm {
const double restart_factor_; //!< factor to multiply time step with when solver fails to converge
const double growth_factor_; //!< factor to multiply time step when solver recovered from failed convergence
const int solver_restart_max_; //!< how many restart of solver are allowed
const bool solver_verbose_; //!< solver verbosity
const bool timestep_verbose_; //!< timestep verbosity
const bool solver_verbose_; //!< solver verbosity
const bool timestep_verbose_; //!< timestep verbosity
double last_timestep_; //!< size of last timestep
};
}

View File

@ -5,15 +5,16 @@
#include <string>
#include <utility>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
#include <opm/core/simulator/PIDTimeStepControl.hpp>
namespace Opm {
// AdaptiveTimeStepping
// AdaptiveTimeStepping
//---------------------
AdaptiveTimeStepping::AdaptiveTimeStepping( const parameter::ParameterGroup& param )
AdaptiveTimeStepping::AdaptiveTimeStepping( const parameter::ParameterGroup& param )
: timeStepControl_()
, initial_fraction_( param.getDefault("solver.initialfraction", double(0.25) ) )
, restart_factor_( param.getDefault("solver.restartfactor", double(0.1) ) )
@ -25,17 +26,17 @@ namespace Opm {
{
// valid are "pid" and "pid+iteration"
std::string control = param.getDefault("timestep.control", std::string("pid") );
const double tol = param.getDefault("timestep.control.tol", double(1e-3) );
if( control == "pid" ) {
timeStepControl_ = TimeStepControlType( new PIDTimeStepControl( tol ) );
}
else if ( control == "pid+iteration" )
else if ( control == "pid+iteration" )
{
const int iterations = param.getDefault("timestep.control.targetiteration", int(25) );
timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
}
else
else
OPM_THROW(std::runtime_error,"Unsupported time step control selected "<< control );
// make sure growth factor is something reasonable
@ -43,10 +44,44 @@ namespace Opm {
}
template <class Solver, class State, class WellState>
void AdaptiveTimeStepping::
step( const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state )
{
const double time = simulatorTimer.simulationTimeElapsed();
const double timestep = simulatorTimer.currentStepLength();
step( solver, state, well_state, time, timestep );
}
template <class Solver, class State, class WellState>
void AdaptiveTimeStepping::
step( const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state,
OutputWriter& outputWriter )
{
const double time = simulatorTimer.simulationTimeElapsed();
const double timestep = simulatorTimer.currentStepLength();
stepImpl( solver, state, well_state, time, timestep, &simulatorTimer, &outputWriter );
}
// implementation of the step method
template <class Solver, class State, class WellState>
void AdaptiveTimeStepping::
step( Solver& solver, State& state, WellState& well_state,
const double time, const double timestep )
{
stepImpl( solver, state, well_state, time, timestep,
(SimulatorTimer *) 0, (OutputWriter *) 0 );
}
// implementation of the step method
template <class Solver, class State, class WState>
void AdaptiveTimeStepping::
stepImpl( Solver& solver, State& state, WState& well_state,
const double time, const double timestep,
const SimulatorTimer* simulatorTimer,
OutputWriter* outputWriter )
{
// init last time step as a fraction of the given time step
if( last_timestep_ < 0 ) {
@ -54,26 +89,26 @@ namespace Opm {
}
// create adaptive step timer with previously used sub step size
AdaptiveSimulatorTimer timer( time, time+timestep, last_timestep_ );
AdaptiveSimulatorTimer substepTimer( time, time+timestep, last_timestep_ );
// copy states in case solver has to be restarted (to be revised)
State last_state( state );
WellState last_well_state( well_state );
State last_state( state );
WState last_well_state( well_state );
// counter for solver restarts
int restarts = 0;
// sub step time loop
while( ! timer.done() )
while( ! substepTimer.done() )
{
// get current delta t
const double dt = timer.currentStepLength() ;
const double dt = substepTimer.currentStepLength() ;
// initialize time step control in case current state is needed later
timeStepControl_->initialize( state );
int linearIterations = -1;
try {
try {
// (linearIterations < 0 means on convergence in solver)
linearIterations = solver.step( dt, state, well_state);
@ -95,7 +130,7 @@ namespace Opm {
if( linearIterations >= 0 )
{
// advance by current dt
++timer;
++substepTimer;
// compute new time step estimate
double dtEstimate =
@ -109,14 +144,23 @@ namespace Opm {
}
if( timestep_verbose_ )
{
std::cout << "Substep[ " << substepTimer.currentStepNum() << " ] " << unit::convert::to(substepTimer.simulationTimeElapsed(),unit::day) << std::endl;
std::cout << "Suggested time step size = " << unit::convert::to(dtEstimate, unit::day) << " (days)" << std::endl;
}
// set new time step length
timer.provideTimeStepEstimate( dtEstimate );
substepTimer.provideTimeStepEstimate( dtEstimate );
// update states
// update states
last_state = state ;
last_well_state = well_state;
// write data if outputWriter was provided
if( outputWriter ) {
assert( simulatorTimer );
outputWriter->writeTimeStep( *simulatorTimer, substepTimer, state, well_state );
}
}
else // in case of no convergence (linearIterations < 0)
{
@ -127,12 +171,12 @@ namespace Opm {
const double newTimeStep = restart_factor_ * dt;
// we need to revise this
timer.provideTimeStepEstimate( newTimeStep );
if( solver_verbose_ )
substepTimer.provideTimeStepEstimate( newTimeStep );
if( solver_verbose_ )
std::cerr << "Solver convergence failed, restarting solver with new time step ("
<< unit::convert::to( newTimeStep, unit::day ) <<" days)." << std::endl;
// reset states
// reset states
state = last_state;
well_state = last_well_state;
@ -142,10 +186,10 @@ namespace Opm {
// store last small time step for next reportStep
last_timestep_ = timer.suggestedAverage();
last_timestep_ = substepTimer.suggestedAverage();
if( timestep_verbose_ )
{
timer.report( std::cout );
substepTimer.report( std::cout );
std::cout << "Last suggested step size = " << unit::convert::to( last_timestep_, unit::day ) << " (days)" << std::endl;
}