OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #713 from dr-robertk/PR/EclipseWriter-revision-to-write-substeps

Browse Source 
Introduce interface for SimulatorTimers and support writing of substeps with EclipseWriter.
This commit is contained in:
Atgeirr Flø Rasmussen 2015-01-15 13:52:22 +01:00
commit 75ae3c33a2
12 changed files with 350 additions and 132 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -367,6 +367,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/simulator/SimulatorOutput.hpp
opm/core/simulator/SimulatorReport.hpp
opm/core/simulator/SimulatorState.hpp
opm/core/simulator/SimulatorTimerInterface.hpp
opm/core/simulator/SimulatorTimer.hpp
opm/core/simulator/TimeStepControlInterface.hpp
opm/core/simulator/TwophaseState.hpp

8
opm/core/io/OutputWriter.cpp
View File

@ -26,15 +26,15 @@ struct MultiWriter : public OutputWriter {
MultiWriter (ptr_t writers) : writers_ (std::move (writers)) { }
/// Forward the call to all writers
virtual void writeInit(const SimulatorTimer &timer) {
virtual void writeInit(const SimulatorTimerInterface &timer) {
for (it_t it = writers_->begin (); it != writers_->end (); ++it) {
(*it)->writeInit (timer);
}
}
virtual void writeTimeStep(const SimulatorTimer& timer,
const SimulatorState& reservoirState,
const WellState& wellState) {
virtual void writeTimeStep(const SimulatorTimerInterface& timer,
const SimulatorState& reservoirState,
const WellState& wellState) {
for (it_t it = writers_->begin (); it != writers_->end(); ++it) {
(*it)->writeTimeStep (timer, reservoirState, wellState);
}

9
opm/core/io/OutputWriter.hpp
View File

@ -21,6 +21,7 @@
#define OPM_OUTPUT_WRITER_HPP
#include <memory> // unique_ptr, shared_ptr
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
struct UnstructuredGrid;
@ -30,7 +31,6 @@ namespace Opm {
class EclipseState;
namespace parameter { class ParameterGroup; }
class SimulatorState;
class SimulatorTimer;
class WellState;
struct PhaseUsage;
@ -73,19 +73,20 @@ public:
* This routine should be called before the first timestep (i.e. when
* timer.currentStepNum () == 0)
*/
virtual void writeInit(const SimulatorTimer &timer) = 0;
virtual void writeInit(const SimulatorTimerInterface &timer) = 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] 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.
*/
virtual void writeTimeStep(const SimulatorTimer& timer,
virtual void writeTimeStep(const SimulatorTimerInterface& timer,
const SimulatorState& reservoirState,
const WellState& wellState) = 0;

104
opm/core/io/eclipse/EclipseWriter.cpp
View File

@ -27,7 +27,7 @@
#include <opm/core/grid.h>
#include <opm/core/grid/cpgpreprocess/preprocess.h>
#include <opm/core/simulator/SimulatorState.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/parameters/Parameter.hpp>
@ -125,8 +125,6 @@ int ertPhaseMask(const PhaseUsage uses)
}
/**
* Eclipse "keyword" (i.e. named data) for a vector.
*/
@ -237,13 +235,13 @@ public:
FileName(const std::string& outputDir,
const std::string& baseName,
ecl_file_enum type,
int reportStepIdx)
int writeStepIdx)
{
ertHandle_ = ecl_util_alloc_filename(outputDir.c_str(),
baseName.c_str(),
type,
false, // formatted?
reportStepIdx);
writeStepIdx);
}
~FileName()
@ -278,15 +276,15 @@ public:
Restart(const std::string& outputDir,
const std::string& baseName,
int reportStepIdx)
int writeStepIdx)
{
restartFileName_ = ecl_util_alloc_filename(outputDir.c_str(),
baseName.c_str(),
/*type=*/ECL_UNIFIED_RESTART_FILE,
false, // use formatted instead of binary output?
reportStepIdx);
writeStepIdx);
if (reportStepIdx == 0) {
if (writeStepIdx == 0) {
restartFileHandle_ = ecl_rst_file_open_write(restartFileName_);
}
else {
@ -372,13 +370,13 @@ public:
ecl_rst_file_close(restartFileHandle_);
}
void writeHeader(const SimulatorTimer& timer,
int reportStepIdx,
void writeHeader(const SimulatorTimerInterface& timer,
int writeStepIdx,
ecl_rsthead_type * rsthead_data)
{
ecl_rst_file_fwrite_header(restartFileHandle_,
reportStepIdx,
writeStepIdx,
rsthead_data);
}
@ -427,7 +425,7 @@ class Summary : private boost::noncopyable
public:
Summary(const std::string& outputDir,
const std::string& baseName,
const SimulatorTimer& timer,
const SimulatorTimerInterface& timer,
int nx,
int ny,
int nz)
@ -463,8 +461,8 @@ public:
// add rate variables for each of the well in the input file
void addAllWells(Opm::EclipseStateConstPtr eclipseState,
const PhaseUsage& uses);
void writeTimeStep(int reportStepIdx,
const SimulatorTimer& timer,
void writeTimeStep(int writeStepIdx,
const SimulatorTimerInterface& timer,
const WellState& wellState);
ecl_sum_type *ertHandle() const
@ -481,11 +479,11 @@ class SummaryTimeStep : private boost::noncopyable
{
public:
SummaryTimeStep(Summary& summaryHandle,
int reportStepIdx,
const SimulatorTimer &timer)
int writeStepIdx,
const SimulatorTimerInterface &timer)
{
ertHandle_ = ecl_sum_add_tstep(summaryHandle.ertHandle(),
reportStepIdx,
writeStepIdx,
Opm::unit::convert::to(timer.simulationTimeElapsed(),
Opm::unit::day));
}
@ -510,16 +508,16 @@ class Init : private boost::noncopyable
public:
Init(const std::string& outputDir,
const std::string& baseName,
int reportStepIdx)
int writeStepIdx)
: egridFileName_(outputDir,
baseName,
ECL_EGRID_FILE,
reportStepIdx)
writeStepIdx)
{
FileName initFileName(outputDir,
baseName,
ECL_INIT_FILE,
reportStepIdx);
writeStepIdx);
bool isFormatted;
if (!ecl_util_fmt_file(initFileName.ertHandle(), &isFormatted)) {
@ -537,7 +535,7 @@ public:
void writeHeader(int numCells,
const int* compressedToCartesianCellIdx,
const SimulatorTimer& timer,
const SimulatorTimerInterface& timer,
Opm::EclipseStateConstPtr eclipseState,
const PhaseUsage uses)
{
@ -621,7 +619,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 writeStepIdx,
const SimulatorTimerInterface& timer,
const WellState& wellState,
const std::map<std::string, int>& nameToIdxMap) = 0;
@ -752,7 +751,8 @@ public:
"SM3/DAY" /* surf. cub. m. per day */)
{ }
virtual double retrieveValue(const SimulatorTimer& timer,
virtual double retrieveValue(const int writeStepIdx,
const SimulatorTimerInterface& timer,
const WellState& wellState,
const std::map<std::string, int>& wellNameToIdxMap)
{
@ -763,7 +763,7 @@ public:
return 0.0;
}
if (well_->getStatus(timer.currentStepNum()) == WellCommon::SHUT) {
if (well_->getStatus(timer.reportStepNum()) == WellCommon::SHUT) {
// well is shut in the current time step
return 0.0;
}
@ -797,17 +797,18 @@ public:
, total_(0.)
{ }
virtual double retrieveValue(const SimulatorTimer& timer,
virtual double retrieveValue(const int writeStepIdx,
const SimulatorTimerInterface& timer,
const WellState& wellState,
const std::map<std::string, int>& wellNameToIdxMap)
{
if (timer.currentStepNum() == 0) {
if (writeStepIdx == 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.reportStepNum()) == WellCommon::SHUT) {
// well is shut in the current time step
return 0.0;
}
@ -855,7 +856,8 @@ public:
"Pascal")
{ }
virtual double retrieveValue(const SimulatorTimer& timer,
virtual double retrieveValue(const int writeStepIdx,
const SimulatorTimerInterface& timer,
const WellState& wellState,
const std::map<std::string, int>& wellNameToIdxMap)
{
@ -865,7 +867,7 @@ public:
// well not active in current time step
return 0.0;
}
if (well_->getStatus(timer.currentStepNum()) == WellCommon::SHUT) {
if (well_->getStatus(timer.reportStepNum()) == WellCommon::SHUT) {
// well is shut in the current time step
return 0.0;
}
@ -876,29 +878,29 @@ public:
// no inline implementation of this since it depends on the
// WellReport type being completed first
void Summary::writeTimeStep(int reportStepIdx,
const SimulatorTimer& timer,
void Summary::writeTimeStep(int writeStepIdx,
const SimulatorTimerInterface& 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.reportStepNum());
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.reportStepNum()) != WellCommon::SHUT ) {
wellNameToIdxMap[timeStepWells[tsWellIdx]->name()] = openWellIdx;
openWellIdx++;
}
}
// internal view; do not move this code out of Summary!
SummaryTimeStep tstep(*this, reportStepIdx, timer);
SummaryTimeStep tstep(*this, writeStepIdx, timer);
// write all the variables
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(writeStepIdx, timer, wellState, wellNameToIdxMap));
}
// write the summary file to disk
@ -1013,7 +1015,7 @@ int EclipseWriter::eclipseWellStatusMask(WellCommon::StatusEnum wellStatus)
}
void EclipseWriter::writeInit(const SimulatorTimer &timer)
void EclipseWriter::writeInit(const SimulatorTimerInterface &timer)
{
// if we don't want to write anything, this method becomes a
// no-op...
@ -1021,7 +1023,7 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer)
return;
}
reportStepIdx_ = 0;
writeStepIdx_ = 0;
EclipseWriterDetails::Init fortio(outputDir_, baseName_, /*stepIdx=*/0);
fortio.writeHeader(numCells_,
@ -1058,8 +1060,8 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer)
summary_->addAllWells(eclipseState_, phaseUsage_);
}
void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
// implementation of the writeTimeStep method
void EclipseWriter::writeTimeStep(const SimulatorTimerInterface& timer,
const SimulatorState& reservoirState,
const WellState& wellState)
{
@ -1070,12 +1072,12 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
}
// respected the output_interval parameter
if (reportStepIdx_ % outputInterval_ != 0) {
if (writeStepIdx_ % outputInterval_ != 0) {
return;
}
std::vector<WellConstPtr> wells_ptr = eclipseState_->getSchedule()->getWells(timer.currentStepNum());
std::vector<WellConstPtr> wells_ptr = eclipseState_->getSchedule()->getWells(timer.reportStepNum());
std::vector<int> iwell_data;
std::vector<const char*> zwell_data;
std::vector<int> icon_data;
@ -1086,22 +1088,22 @@ 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.reportStepNum());
rsthead_data.niwelz = 0;
rsthead_data.nzwelz = 0;
rsthead_data.niconz = 0;
rsthead_data.ncwmax = 0;
rsthead_data.phase_sum = Opm::EclipseWriterDetails::ertPhaseMask(phaseUsage_);
EclipseWriterDetails::Restart restartHandle(outputDir_, baseName_, reportStepIdx_);
EclipseWriterDetails::Restart restartHandle(outputDir_, baseName_, writeStepIdx_);
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.reportStepNum());
restartHandle.addRestartFileIwelData(iwell_data, timer.reportStepNum(), well_ptr);
restartHandle.addRestartFileZwelData(zwell_data, timer.reportStepNum(), well_ptr);
restartHandle.addRestartFileIconData(icon_data, timer.reportStepNum(), rsthead_data.ncwmax, well_ptr);
rsthead_data.niwelz = EclipseWriterDetails::Restart::NIWELZ;
rsthead_data.nzwelz = EclipseWriterDetails::Restart::NZWELZ;
@ -1111,7 +1113,7 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
rsthead_data.sim_days = Opm::unit::convert::to(timer.simulationTimeElapsed(), Opm::unit::day); //data for doubhead
restartHandle.writeHeader(timer,
reportStepIdx_,
writeStepIdx_,
&rsthead_data);
@ -1163,9 +1165,9 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
// instead of creating a temporary EclipseWriterDetails::Summary in this function
// every time it is called. This has been changed so that the final summary file
// will contain data from the whole simulation, instead of just the last step.
summary_->writeTimeStep(reportStepIdx_, timer, wellState);
summary_->writeTimeStep(writeStepIdx_, timer, wellState);
++reportStepIdx_;
++writeStepIdx_;
}
@ -1241,7 +1243,7 @@ void EclipseWriter::init(const parameter::ParameterGroup& params)
outputDir_ = params.getDefault<std::string>("output_dir", ".");
// set the index of the first time step written to 0...
reportStepIdx_ = 0;
writeStepIdx_ = 0;
if (enableOutput_) {
// make sure that the output directory exists, if not try to create it

13
opm/core/io/eclipse/EclipseWriter.hpp
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).
@ -24,6 +25,7 @@
#include <opm/core/io/OutputWriter.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/wells.h> // WellType
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@ -42,7 +44,6 @@ class Summary;
}
class SimulatorState;
class SimulatorTimer;
class WellState;
namespace parameter { class ParameterGroup; }
@ -78,7 +79,7 @@ public:
/**
* Write the static eclipse data (grid, PVT curves, etc) to disk.
*/
virtual void writeInit(const SimulatorTimer &timer);
virtual void writeInit(const SimulatorTimerInterface &timer);
/*!
* \brief Write a reservoir state and summary information to disk.
@ -91,13 +92,15 @@ 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,
virtual void writeTimeStep(const SimulatorTimerInterface& timer,
const SimulatorState& reservoirState,
const WellState& wellState);
static int eclipseWellTypeMask(WellType wellType, WellInjector::TypeEnum injectorType);
static int eclipseWellStatusMask(WellCommon::StatusEnum wellStatus);
@ -110,7 +113,7 @@ private:
double deckToSiPressure_;
bool enableOutput_;
int outputInterval_;
int reportStepIdx_;
int writeStepIdx_;
std::string outputDir_;
std::string baseName_;
PhaseUsage phaseUsage_; // active phases in the input deck

28
opm/core/simulator/AdaptiveSimulatorTimer.cpp
View File

@ -1,5 +1,5 @@
/*
Copyright 2014 IRIS AS
Copyright (c) 2014 IRIS AS
This file is part of the Open Porous Media project (OPM).
@ -30,11 +30,13 @@
namespace Opm
{
AdaptiveSimulatorTimer::
AdaptiveSimulatorTimer( const double start_time, const double total_time, const double lastDt )
: start_time_( start_time )
, total_time_( total_time )
AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer, const double lastStepTaken )
: start_date_time_( timer.startDateTime() )
, start_time_( timer.simulationTimeElapsed() )
, total_time_( start_time_ + timer.currentStepLength() )
, report_step_( timer.reportStepNum() )
, current_time_( start_time_ )
, dt_( computeInitialTimeStep( lastDt ) )
, dt_( computeInitialTimeStep( lastStepTaken ) )
, current_step_( 0 )
, steps_()
, suggestedMax_( 0.0 )
@ -86,12 +88,23 @@ namespace Opm
int AdaptiveSimulatorTimer::
currentStepNum () const { return current_step_; }
int AdaptiveSimulatorTimer::
reportStepNum () const { return report_step_; }
double AdaptiveSimulatorTimer::currentStepLength () const
{
assert( ! done () );
return dt_;
}
double AdaptiveSimulatorTimer::stepLengthTaken() const
{
assert( ! steps_.empty() );
return steps_.back();
}
double AdaptiveSimulatorTimer::totalTime() const { return total_time_; }
double AdaptiveSimulatorTimer::simulationTimeElapsed() const { return current_time_; }
@ -143,6 +156,11 @@ namespace Opm
std::cout << "sub steps end time = " << unit::convert::to( simulationTimeElapsed(), unit::day ) << " (days)" << std::endl;
}
boost::posix_time::ptime AdaptiveSimulatorTimer::startDateTime() const
{
return start_date_time_;
}
double AdaptiveSimulatorTimer::
computeInitialTimeStep( const double lastDt ) const
{

24
opm/core/simulator/AdaptiveSimulatorTimer.hpp
View File

@ -16,7 +16,6 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ADAPTIVESIMULATORTIMER_HEADER_INCLUDED
#define OPM_ADAPTIVESIMULATORTIMER_HEADER_INCLUDED
@ -27,6 +26,8 @@
#include <algorithm>
#include <numeric>
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
namespace Opm
{
@ -35,14 +36,12 @@ namespace Opm
/// \brief Simulation timer for adaptive time stepping
///
/////////////////////////////////////////////////////////
class AdaptiveSimulatorTimer
class AdaptiveSimulatorTimer : public SimulatorTimerInterface
{
public:
/// \brief constructor taking a simulator timer to determine start and end time
/// \param start_time start time of timer
/// \param total_time total time of timer
/// \param lastDt last suggested length of time step interval
AdaptiveSimulatorTimer( const double start_time, const double total_time, const double lastDt );
/// \param timer in case of sub stepping this is the outer timer
explicit AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer, const double lastStepTaken );
/// \brief advance time by currentStepLength
AdaptiveSimulatorTimer& operator++ ();
@ -53,6 +52,9 @@ namespace Opm
/// \brief \copydoc SimulationTimer::currentStepNum
int currentStepNum () const;
/// \brief return current report step
int reportStepNum() const;
/// \brief \copydoc SimulationTimer::currentStepLength
double currentStepLength () const;
@ -80,12 +82,22 @@ namespace Opm
/// \brief return average suggested step length
double suggestedAverage () const;
/// \brief Previous step length. This is the length of the step that
/// was taken to arrive at this time.
double stepLengthTaken () const;
/// \brief report start and end time as well as used steps so far
void report(std::ostream& os) const;
/// \brief start date time of simulation
boost::posix_time::ptime startDateTime() const;
protected:
const boost::posix_time::ptime start_date_time_;
const double start_time_;
const double total_time_;
const int report_step_;
double current_time_;
double dt_;
int current_step_;

64
opm/core/simulator/AdaptiveTimeStepping.hpp
View File

@ -1,3 +1,21 @@
/*
Copyright 2014 IRIS AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SUBSTEPPING_HEADER_INCLUDED
#define OPM_SUBSTEPPING_HEADER_INCLUDED
@ -6,33 +24,53 @@
#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 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 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
*/
*/
template <class Solver, class State, class WellState>
void step( Solver& solver, State& state, WellState& well_state,
const double time, const double timestep );
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( const SimulatorTimer& timer,
Solver& solver, State& state, WellState& well_state,
OutputWriter* outputWriter);
typedef std::unique_ptr< TimeStepControlInterface > TimeStepControlType;
TimeStepControlType timeStepControl_; //!< time step control object
@ -40,8 +78,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
};
}

97
opm/core/simulator/AdaptiveTimeStepping_impl.hpp
View File

@ -1,3 +1,21 @@
/*
Copyright 2014 IRIS AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
#define OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
@ -5,15 +23,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 +44,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
@ -45,35 +64,54 @@ namespace Opm {
template <class Solver, class State, class WellState>
void AdaptiveTimeStepping::
step( Solver& solver, State& state, WellState& well_state,
const double time, const double timestep )
step( const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state )
{
stepImpl( simulatorTimer, solver, state, well_state );
}
template <class Solver, class State, class WellState>
void AdaptiveTimeStepping::
step( const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state,
OutputWriter& outputWriter )
{
stepImpl( simulatorTimer, solver, state, well_state, &outputWriter );
}
// implementation of the step method
template <class Solver, class State, class WState>
void AdaptiveTimeStepping::
stepImpl( const SimulatorTimer& simulatorTimer,
Solver& solver, State& state, WState& well_state,
OutputWriter* outputWriter )
{
const double timestep = simulatorTimer.currentStepLength();
// init last time step as a fraction of the given time step
if( last_timestep_ < 0 ) {
last_timestep_ = initial_fraction_ * timestep ;
last_timestep_ = initial_fraction_ * timestep;
}
// create adaptive step timer with previously used sub step size
AdaptiveSimulatorTimer timer( time, time+timestep, last_timestep_ );
AdaptiveSimulatorTimer substepTimer( simulatorTimer, 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 +133,7 @@ namespace Opm {
if( linearIterations >= 0 )
{
// advance by current dt
++timer;
++substepTimer;
// compute new time step estimate
double dtEstimate =
@ -109,14 +147,25 @@ namespace Opm {
}
if( timestep_verbose_ )
std::cout << "Suggested time step size = " << unit::convert::to(dtEstimate, unit::day) << " (days)" << std::endl;
{
std::cout << std::endl
<<"Substep( " << substepTimer.currentStepNum()
<< " ): Current time (days) " << unit::convert::to(substepTimer.simulationTimeElapsed(),unit::day) << std::endl
<< " Current stepsize est (days) " << unit::convert::to(dtEstimate, unit::day) << std::endl;
}
// write data if outputWriter was provided
if( outputWriter ) {
outputWriter->writeTimeStep( substepTimer, state, well_state );
}
// set new time step length
timer.provideTimeStepEstimate( dtEstimate );
substepTimer.provideTimeStepEstimate( dtEstimate );
// update states
// update states
last_state = state ;
last_well_state = well_state;
}
else // in case of no convergence (linearIterations < 0)
{
@ -127,12 +176,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 +191,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;
}

12
opm/core/simulator/SimulatorTimer.cpp
View File

@ -100,19 +100,11 @@ namespace Opm
return current_time_;
}
/// time elapsed since the start of the POSIX epoch (Jan 1st, 1970) [s].
time_t SimulatorTimer::currentPosixTime() const
boost::posix_time::ptime SimulatorTimer::startDateTime() const
{
tm t = boost::posix_time::to_tm(currentDateTime());
return std::mktime(&t);
return boost::posix_time::ptime(start_date_);
}
boost::posix_time::ptime SimulatorTimer::currentDateTime() const
{
return boost::posix_time::ptime(start_date_) + boost::posix_time::seconds( (int) current_time_ );
}
/// Total time.
double SimulatorTimer::totalTime() const

19
opm/core/simulator/SimulatorTimer.hpp
View File

@ -21,20 +21,23 @@
#define OPM_SIMULATORTIMER_HEADER_INCLUDED
#include <opm/parser/eclipse/EclipseState/Schedule/TimeMap.hpp>
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
#include <iosfwd>
#include <vector>
#include <boost/date_time/gregorian/gregorian.hpp>
#include <boost/date_time/posix_time/posix_time_types.hpp>
namespace Opm
{
namespace parameter { class ParameterGroup; }
class SimulatorTimer
class SimulatorTimer : public SimulatorTimerInterface
{
public:
// use default implementation of these methods
using SimulatorTimerInterface::currentDateTime;
using SimulatorTimerInterface::currentPosixTime;
/// Default constructor.
SimulatorTimer();
@ -72,20 +75,16 @@ namespace Opm
/// it is an error to call stepLengthTaken().
double stepLengthTaken () const;
/// Time elapsed since the start of the POSIX epoch (Jan 1st,
/// 1970) until the current time step begins [s].
time_t currentPosixTime() const;
/// Time elapsed since the start of the simulation until the
/// beginning of the current time step [s].
double simulationTimeElapsed() const;
/// Return the current time as a posix time object.
boost::posix_time::ptime currentDateTime() const;
/// Total time.
double totalTime() const;
/// Return start date of simulation
boost::posix_time::ptime startDateTime() const;
/// Set total time.
/// This is primarily intended for multi-epoch schedules,
/// where a timer for a given epoch does not have

103
opm/core/simulator/SimulatorTimerInterface.hpp Normal file
View File

@ -0,0 +1,103 @@
/*
Copyright (c) 2014 IRIS AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SIMULATORTIMERINTERFACE_HEADER_INCLUDED
#define OPM_SIMULATORTIMERINTERFACE_HEADER_INCLUDED
#include <boost/date_time/gregorian/gregorian.hpp>
#include <boost/date_time/posix_time/posix_time_types.hpp>
#include <boost/date_time/posix_time/conversion.hpp>
namespace Opm
{
namespace parameter { class ParameterGroup; }
/// Interface class for SimulatorTimer objects, to be improved.
class SimulatorTimerInterface
{
protected:
/// Default constructor, protected to not allow explicit instances of this class.
SimulatorTimerInterface() {}
public:
/// destructor
virtual ~SimulatorTimerInterface() {}
/// Current step number. This is the number of timesteps that
/// has been completed from the start of the run. The time
/// after initialization but before the simulation has started
/// is timestep number zero.
virtual int currentStepNum() const = 0;
/// Current report step number. This might differ from currentStepNum in case of sub stepping
virtual int reportStepNum() const { return currentStepNum(); }
/// Current step length. This is the length of the step
/// the simulator will take in the next iteration.
///
/// @note if done(), it is an error to call currentStepLength().
virtual double currentStepLength() const = 0;
/// Previous step length. This is the length of the step that
/// was taken to arrive at this time.
///
/// @note if no increments have been done (i.e. the timer is
/// still in its constructed state and currentStepNum() == 0),
/// it is an error to call stepLengthTaken().
virtual double stepLengthTaken () const = 0;
/// Previous report step length. This is the length of the step that
/// was taken to arrive at this report step time.
///
/// @note if no increments have been done (i.e. the timer is
/// still in its constructed state and reportStepNum() == 0),
/// it is an error to call stepLengthTaken().
virtual double reportStepLengthTaken () const { return stepLengthTaken(); }
/// Time elapsed since the start of the simulation until the
/// beginning of the current time step [s].
virtual double simulationTimeElapsed() const = 0;
/// Return true if timer indicates that simulation of timer interval is finished
virtual bool done() const = 0;
/// Return start date of simulation
virtual boost::posix_time::ptime startDateTime() const = 0;
/// Return the current time as a posix time object.
virtual boost::posix_time::ptime currentDateTime() const
{
return startDateTime() + boost::posix_time::seconds( (int) simulationTimeElapsed());
//boost::posix_time::ptime(startDate()) + boost::posix_time::seconds( (int) simulationTimeElapsed());
}
/// Time elapsed since the start of the POSIX epoch (Jan 1st,
/// 1970) until the current time step begins [s].
virtual time_t currentPosixTime() const
{
tm t = boost::posix_time::to_tm(currentDateTime());
return std::mktime(&t);
}
};
} // namespace Opm
#endif // OPM_SIMULATORTIMER_HEADER_INCLUDED