diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 3a1ce398..149f73f9 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -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 diff --git a/opm/core/io/OutputWriter.cpp b/opm/core/io/OutputWriter.cpp index fcdc2b5d..1431ad46 100644 --- a/opm/core/io/OutputWriter.cpp +++ b/opm/core/io/OutputWriter.cpp @@ -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); } diff --git a/opm/core/io/OutputWriter.hpp b/opm/core/io/OutputWriter.hpp index 84d2591b..27f4dd3b 100644 --- a/opm/core/io/OutputWriter.hpp +++ b/opm/core/io/OutputWriter.hpp @@ -21,6 +21,7 @@ #define OPM_OUTPUT_WRITER_HPP #include // unique_ptr, shared_ptr +#include 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; diff --git a/opm/core/io/eclipse/EclipseWriter.cpp b/opm/core/io/eclipse/EclipseWriter.cpp index a83bf3a0..9eb48e28 100644 --- a/opm/core/io/eclipse/EclipseWriter.cpp +++ b/opm/core/io/eclipse/EclipseWriter.cpp @@ -27,7 +27,7 @@ #include #include #include -#include +#include #include #include #include @@ -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& 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& 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& 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& 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 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 wells_ptr = eclipseState_->getSchedule()->getWells(timer.currentStepNum()); + std::vector wells_ptr = eclipseState_->getSchedule()->getWells(timer.reportStepNum()); std::vector iwell_data; std::vector zwell_data; std::vector 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::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("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 diff --git a/opm/core/io/eclipse/EclipseWriter.hpp b/opm/core/io/eclipse/EclipseWriter.hpp index 5e715274..5ff4e7f3 100644 --- a/opm/core/io/eclipse/EclipseWriter.hpp +++ b/opm/core/io/eclipse/EclipseWriter.hpp @@ -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 #include #include // WellType +#include #include @@ -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 diff --git a/opm/core/simulator/AdaptiveSimulatorTimer.cpp b/opm/core/simulator/AdaptiveSimulatorTimer.cpp index 576e4fae..83716dfe 100644 --- a/opm/core/simulator/AdaptiveSimulatorTimer.cpp +++ b/opm/core/simulator/AdaptiveSimulatorTimer.cpp @@ -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 { diff --git a/opm/core/simulator/AdaptiveSimulatorTimer.hpp b/opm/core/simulator/AdaptiveSimulatorTimer.hpp index 3336a69c..c2082619 100644 --- a/opm/core/simulator/AdaptiveSimulatorTimer.hpp +++ b/opm/core/simulator/AdaptiveSimulatorTimer.hpp @@ -16,7 +16,6 @@ You should have received a copy of the GNU General Public License along with OPM. If not, see . */ - #ifndef OPM_ADAPTIVESIMULATORTIMER_HEADER_INCLUDED #define OPM_ADAPTIVESIMULATORTIMER_HEADER_INCLUDED @@ -27,6 +26,8 @@ #include #include +#include + 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_; diff --git a/opm/core/simulator/AdaptiveTimeStepping.hpp b/opm/core/simulator/AdaptiveTimeStepping.hpp index e73790db..fe5f860a 100644 --- a/opm/core/simulator/AdaptiveTimeStepping.hpp +++ b/opm/core/simulator/AdaptiveTimeStepping.hpp @@ -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 . +*/ #ifndef OPM_SUBSTEPPING_HEADER_INCLUDED #define OPM_SUBSTEPPING_HEADER_INCLUDED @@ -6,33 +24,53 @@ #include #include +#include #include 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 - 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 + void step( const SimulatorTimer& timer, + Solver& solver, State& state, WellState& well_state, + OutputWriter& outputWriter ); protected: + template + 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 }; } diff --git a/opm/core/simulator/AdaptiveTimeStepping_impl.hpp b/opm/core/simulator/AdaptiveTimeStepping_impl.hpp index 8e1bdda4..b1d5ed71 100644 --- a/opm/core/simulator/AdaptiveTimeStepping_impl.hpp +++ b/opm/core/simulator/AdaptiveTimeStepping_impl.hpp @@ -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 . +*/ #ifndef OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED #define OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED @@ -5,15 +23,16 @@ #include #include +#include #include #include 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 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 + 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 + 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; } diff --git a/opm/core/simulator/SimulatorTimer.cpp b/opm/core/simulator/SimulatorTimer.cpp index add8b738..75785751 100644 --- a/opm/core/simulator/SimulatorTimer.cpp +++ b/opm/core/simulator/SimulatorTimer.cpp @@ -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 diff --git a/opm/core/simulator/SimulatorTimer.hpp b/opm/core/simulator/SimulatorTimer.hpp index 61faef09..8f9e7b0b 100644 --- a/opm/core/simulator/SimulatorTimer.hpp +++ b/opm/core/simulator/SimulatorTimer.hpp @@ -21,20 +21,23 @@ #define OPM_SIMULATORTIMER_HEADER_INCLUDED #include +#include #include #include -#include -#include 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 diff --git a/opm/core/simulator/SimulatorTimerInterface.hpp b/opm/core/simulator/SimulatorTimerInterface.hpp new file mode 100644 index 00000000..25118fd2 --- /dev/null +++ b/opm/core/simulator/SimulatorTimerInterface.hpp @@ -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 . +*/ + +#ifndef OPM_SIMULATORTIMERINTERFACE_HEADER_INCLUDED +#define OPM_SIMULATORTIMERINTERFACE_HEADER_INCLUDED + +#include +#include +#include + +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