diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index c65af9dc4..3abd6613e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -67,6 +67,9 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/SimulatorIncompTwophase.cpp opm/simulators/WellSwitchingLogger.cpp opm/simulators/vtk/writeVtkData.cpp + opm/simulators/timestepping/TimeStepControl.cpp + opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp + opm/simulators/timestepping/SimulatorTimer.cpp ) @@ -89,6 +92,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_multisegmentwells.cpp # tests/test_thresholdpressure.cpp tests/test_wellswitchlogger.cpp + tests/test_timer.cpp ) list (APPEND TEST_DATA_FILES @@ -96,6 +100,7 @@ list (APPEND TEST_DATA_FILES tests/VFPPROD1 tests/VFPPROD2 tests/msw.data + tests/TESTTIMER.DATA ) @@ -117,6 +122,7 @@ list (APPEND EXAMPLE_SOURCE_FILES examples/sim_poly2p_incomp_reorder.cpp examples/sim_poly_fi2p_comp_ad.cpp examples/flow_polymer.cpp + examples/wells_example.cpp ) # programs listed here will not only be compiled, but also marked for @@ -259,5 +265,12 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/thresholdPressures.hpp opm/simulators/WellSwitchingLogger.hpp opm/simulators/vtk/writeVtkData.hpp + opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp + opm/simulators/timestepping/AdaptiveTimeStepping.hpp + opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp + opm/simulators/timestepping/TimeStepControl.hpp + opm/simulators/timestepping/TimeStepControlInterface.hpp + opm/simulators/timestepping/SimulatorTimer.hpp + opm/simulators/timestepping/SimulatorTimerInterface.hpp ) diff --git a/examples/wells_example.cpp b/examples/wells_example.cpp new file mode 100644 index 000000000..53c6f548e --- /dev/null +++ b/examples/wells_example.cpp @@ -0,0 +1,144 @@ +#include "config.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +int main(int argc, char** argv) +try +{ + using namespace Opm::parameter; + using namespace Opm; + ParameterGroup parameters(argc, argv, false); + std::string file_name = parameters.getDefault ("inputdeck", "data.data"); + + SimulatorTimer simtimer; + simtimer.init(parameters); + + // Read input file + ParseContext parseContext; + Opm::Parser parser; + Opm::Deck deck = parser.parseFile(file_name , parseContext); + Opm::EclipseState eclipseState(deck , parseContext); + std::cout << "Done!" << std::endl; + + // Setup grid + GridManager grid(eclipseState.getInputGrid()); + + // Define rock and fluid properties + IncompPropertiesFromDeck incomp_properties(deck, eclipseState, *grid.c_grid()); + RockCompressibility rock_comp(eclipseState); + + // Finally handle the wells + WellsManager wells(eclipseState , 0 , *grid.c_grid()); + + double gravity[3] = {0.0, 0.0, parameters.getDefault("gravity", 0.0)}; + Opm::LinearSolverFactory linsolver(parameters); + double nl_pressure_residual_tolerance = 1e-8; + double nl_pressure_change_tolerance = 0.0; + int nl_pressure_maxiter = 100; + if (rock_comp.isActive()) { + nl_pressure_residual_tolerance = parameters.getDefault("nl_pressure_residual_tolerance", 1e-8); + nl_pressure_change_tolerance = parameters.getDefault("nl_pressure_change_tolerance", 1.0); // in Pascal + nl_pressure_maxiter = parameters.getDefault("nl_pressure_maxiter", 10); + } + + std::vector src; + Opm::FlowBCManager bcs; + + // EXPERIMENT_ISTL + IncompTpfa pressure_solver(*grid.c_grid(), incomp_properties, &rock_comp, linsolver, + nl_pressure_residual_tolerance, nl_pressure_change_tolerance, nl_pressure_maxiter, + gravity, wells.c_wells(), src, bcs.c_bcs()); + + + std::vector all_cells; + for (int i = 0; i < grid.c_grid()->number_of_cells; i++) { + all_cells.push_back(i); + } + + Opm::TwophaseState state( grid.c_grid()->number_of_cells , grid.c_grid()->number_of_faces ); + + initStateFromDeck(*grid.c_grid(), incomp_properties, deck, gravity[2], state); + + Opm::WellState well_state; + well_state.init(wells.c_wells(), state); + + pressure_solver.solve(simtimer.currentStepLength(), state, well_state); + + const int np = incomp_properties.numPhases(); + std::vector fractional_flows(grid.c_grid()->number_of_cells*np, 0.0); + computeFractionalFlow(incomp_properties, all_cells, state.saturation(), fractional_flows); + + // This will be refactored into a separate function once done + std::vector well_resflows(wells.c_wells()->number_of_wells*np, 0.0); + computePhaseFlowRatesPerWell(*wells.c_wells(), well_state.perfRates(), fractional_flows, well_resflows); + // We approximate (for _testing_ that resflows = surfaceflows) + for (int wc_iter = 0; wc_iter < 10 && !wells.conditionsMet(well_state.bhp(), well_resflows, well_resflows); ++wc_iter) { + std::cout << "Conditions not met for well, trying again" << std::endl; + pressure_solver.solve(simtimer.currentStepLength(), state, well_state); + std::cout << "Solved" << std::endl; + + computePhaseFlowRatesPerWell(*wells.c_wells(), well_state.perfRates(), fractional_flows, well_resflows); + } + +#if 0 + std::vector porevol; + computePorevolume(*grid->c_grid(), incomp_properties, porevol); + + + + TwophaseFluid fluid(incomp_properties); + TransportContextl model(fluid, *grid->c_grid(), porevol, gravity[2], true); + + TransportSolver tsolver(model); + + TransportSource* tsrc = create_transport_source(2, 2); + double ssrc[] = {1.0, 0.0}; + double ssink[] = {0.0, 1.0}; + double zdummy[] = {0.0, 0.0}; + + { + int well_cell_index = 0; + for (int well = 0; well < wells.c_wells()->number_of_wells; ++well) { + for (int cell = wells.c_wells()->well_connpos[well]; cell < wells.c_wells()->well_connpos[well + 1]; ++cell) { + if (well_rate_per_cell[well_cell_index] > 0.0) { + append_transport_source(well_cell_index, 2, 0, + well_rate_per_cell[well_cell_index], ssrc, zdummy, tsrc); + } else if (well_rate_per_cell[well_cell_index] < 0.0) { + append_transport_source(well_cell_index, 2, 0, + well_rate_per_cell[well_cell_index], ssink, zdummy, tsrc); + } + } + } + } + + tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt); + + Opm::computeInjectedProduced(*props, state.saturation(), src, stepsize, injected, produced); +#endif + return 0; +} +catch (const std::exception &e) { + std::cerr << "Program threw an exception: " << e.what() << "\n"; + throw; +} diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index d9b5dc8ff..62877d03f 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -48,7 +48,7 @@ #include #include #include -#include +#include #include #include #include // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix. diff --git a/opm/autodiff/SimulatorBase.hpp b/opm/autodiff/SimulatorBase.hpp index dc5591b22..2a09f3506 100644 --- a/opm/autodiff/SimulatorBase.hpp +++ b/opm/autodiff/SimulatorBase.hpp @@ -40,8 +40,8 @@ #include #include -#include -#include +#include +#include #include #include #include @@ -49,7 +49,7 @@ #include #include -#include +#include #include #include diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 871062fa0..cb83aad67 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -32,7 +32,7 @@ #include #include -#include +#include #include #include diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp index 1cba6695a..81d095c1e 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp @@ -42,14 +42,14 @@ #include #include -//#include +//#include #include #include #include #include -//#include +//#include //#include #include diff --git a/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp b/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp new file mode 100644 index 000000000..577cbbe42 --- /dev/null +++ b/opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp @@ -0,0 +1,166 @@ +/* + 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 . +*/ + +#include +#include +#include + +#include +#include + +#include +#include + +namespace Opm +{ + AdaptiveSimulatorTimer:: + AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer, + const double lastStepTaken, + const double maxTimeStep ) + : start_date_time_( timer.startDateTime() ) + , start_time_( timer.simulationTimeElapsed() ) + , total_time_( start_time_ + timer.currentStepLength() ) + , report_step_( timer.reportStepNum() ) + , max_time_step_( maxTimeStep ) + , current_time_( start_time_ ) + , dt_( 0.0 ) + , current_step_( 0 ) + , steps_() + , lastStepFailed_( false ) + { + // reserve memory for sub steps + steps_.reserve( 10 ); + + // set appropriate value for dt_ + provideTimeStepEstimate( lastStepTaken ); + } + + AdaptiveSimulatorTimer& AdaptiveSimulatorTimer::operator++ () + { + ++current_step_; + current_time_ += dt_; + // store used time step sizes + steps_.push_back( dt_ ); + return *this; + } + + void AdaptiveSimulatorTimer:: + provideTimeStepEstimate( const double dt_estimate ) + { + double remaining = (total_time_ - current_time_); + // apply max time step if it was set + dt_ = std::min( dt_estimate, max_time_step_ ); + + if( remaining > 0 ) { + + // set new time step (depending on remaining time) + if( 1.05 * dt_ > remaining ) { + dt_ = remaining; + // check max time step again and use half remaining if to large + if( dt_ > max_time_step_ ) { + dt_ = 0.5 * remaining; + } + return; + } + + // check for half interval step to avoid very small step at the end + // remaining *= 0.5; + + if( 1.5 * dt_ > remaining ) { + dt_ = 0.5 * remaining; + return; + } + } + } + + int AdaptiveSimulatorTimer:: + currentStepNum () const { return current_step_; } + + int AdaptiveSimulatorTimer:: + reportStepNum () const { return report_step_; } + + double AdaptiveSimulatorTimer::currentStepLength () const + { + 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_; } + + bool AdaptiveSimulatorTimer::done () const { return (current_time_ >= total_time_) ; } + + double AdaptiveSimulatorTimer::averageStepLength() const + { + const int size = steps_.size(); + if( size == 0 ) return 0.0; + + const double sum = std::accumulate(steps_.begin(), steps_.end(), 0.0); + return sum / double(size); + } + + /// \brief return max step length used so far + double AdaptiveSimulatorTimer::maxStepLength () const + { + if( steps_.size() == 0 ) return 0.0; + return *(std::max_element( steps_.begin(), steps_.end() )); + } + + /// \brief return min step length used so far + double AdaptiveSimulatorTimer::minStepLength () const + { + if( steps_.size() == 0 ) return 0.0; + return *(std::min_element( steps_.begin(), steps_.end() )); + } + + /// \brief report start and end time as well as used steps so far + void AdaptiveSimulatorTimer:: + report(std::ostream& os) const + { + os << "Sub steps started at time = " << unit::convert::to( start_time_, unit::day ) << " (days)" << std::endl; + for( size_t i=0; i + AdaptiveSimulatorTimer::clone() const + { + return std::unique_ptr< SimulatorTimerInterface > (new AdaptiveSimulatorTimer( *this )); + } + + + +} // namespace Opm diff --git a/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp b/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp new file mode 100644 index 000000000..fd34af8c1 --- /dev/null +++ b/opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp @@ -0,0 +1,123 @@ +/* + 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_ADAPTIVESIMULATORTIMER_HEADER_INCLUDED +#define OPM_ADAPTIVESIMULATORTIMER_HEADER_INCLUDED + +#include +#include +#include + +#include +#include + +#include + +namespace Opm +{ + + ///////////////////////////////////////////////////////// + /// + /// \brief Simulation timer for adaptive time stepping + /// + ///////////////////////////////////////////////////////// + class AdaptiveSimulatorTimer : public SimulatorTimerInterface + { + public: + /// \brief constructor taking a simulator timer to determine start and end time + /// \param timer in case of sub stepping this is the outer timer + /// \param lastStepTaken last suggested time step + /// \param maxTimeStep maximum time step allowed + AdaptiveSimulatorTimer( const SimulatorTimerInterface& timer, + const double lastStepTaken, + const double maxTimeStep = std::numeric_limits::max() ); + + /// \brief advance time by currentStepLength + AdaptiveSimulatorTimer& operator++ (); + + /// \brief advance time by currentStepLength + void advance() { this->operator++ (); } + + /// \brief provide and estimate for new time step size + void provideTimeStepEstimate( const double dt_estimate ); + + /// \brief \copydoc SimulationTimer::currentStepNum + int currentStepNum () const; + + /// \brief return current report step + int reportStepNum() const; + + /// \brief \copydoc SimulationTimer::currentStepLength + double currentStepLength () const; + + /// \brief \copydoc SimulationTimer::totalTime + double totalTime() const; + + /// \brief \copydoc SimulationTimer::simulationTimeElapsed + double simulationTimeElapsed() const; + + /// \brief \copydoc SimulationTimer::done + bool done () const; + + /// \brief return average step length used so far + double averageStepLength() const; + + /// \brief return max step length used so far + double maxStepLength () const; + + /// \brief return min step length used so far + double minStepLength () 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; + + /// \brief Return true if last time step failed + bool lastStepFailed() const {return lastStepFailed_;} + + /// \brief tell the timestepper whether timestep failed or not + void setLastStepFailed(bool lastStepFailed) {lastStepFailed_ = lastStepFailed;} + + /// return copy of object + virtual std::unique_ptr< SimulatorTimerInterface > clone() const; + + protected: + const boost::posix_time::ptime start_date_time_; + const double start_time_; + const double total_time_; + const int report_step_; + const double max_time_step_; + + double current_time_; + double dt_; + int current_step_; + + std::vector< double > steps_; + bool lastStepFailed_; + + }; + +} // namespace Opm + +#endif // OPM_SIMULATORTIMER_HEADER_INCLUDED diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp new file mode 100644 index 000000000..74cfb99f4 --- /dev/null +++ b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp @@ -0,0 +1,112 @@ +/* + Copyright 2014 IRIS AS + Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015 Statoil 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 + +#include +#include + +#include +#include +#include +#include + +namespace Opm { + + + // AdaptiveTimeStepping + //--------------------- + + class AdaptiveTimeStepping + { + public: + //! \brief contructor taking parameter object + //! \param param The parameter object + //! \param pinfo The information about the data distribution + //! and communication for a parallel run. + AdaptiveTimeStepping( const parameter::ParameterGroup& param, + const bool terminal_output = true ); + + //! \brief contructor taking parameter object + //! \param tuning Pointer to ecl TUNING keyword + //! \param time_step current report step + //! \param param The parameter object + //! \param pinfo The information about the data distribution + //! and communication for a parallel run. + AdaptiveTimeStepping( const Tuning& tuning, + size_t time_step, + const parameter::ParameterGroup& param, + const bool terminal_output = true ); + + /** \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 + SimulatorReport 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 fipnum Fluid-in-place numbering array + \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 + SimulatorReport step( const SimulatorTimer& timer, + Solver& solver, State& state, WellState& well_state, + Output& outputWriter, + const std::vector* fipnum = nullptr); + + protected: + template + SimulatorReport stepImpl( const SimulatorTimer& timer, + Solver& solver, State& state, WellState& well_state, + Output* outputWriter, + const std::vector* fipnum); + + void init(const parameter::ParameterGroup& param); + + typedef std::unique_ptr< TimeStepControlInterface > TimeStepControlType; + + TimeStepControlType timeStepControl_; //!< time step control object + 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 double max_growth_; //!< factor that limits the maximum growth of a time step + const double max_time_step_; //!< maximal allowed time step size + const int solver_restart_max_; //!< how many restart of solver are allowed + const bool solver_verbose_; //!< solver verbosity + const bool timestep_verbose_; //!< timestep verbosity + double suggested_next_timestep_; //!< suggested size of next timestep + bool full_timestep_initially_; //!< beginning with the size of the time step from data file + }; +} + +#include +#endif diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp new file mode 100644 index 000000000..440555c5f --- /dev/null +++ b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp @@ -0,0 +1,356 @@ +/* + 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 + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include // For MatrixBlockException +#include + +namespace Opm { + + namespace detail + { + template + class SolutionTimeErrorSolverWrapper : public RelativeChangeInterface + { + const Solver& solver_; + const State& previous_; + const State& current_; + public: + SolutionTimeErrorSolverWrapper( const Solver& solver, + const State& previous, + const State& current ) + : solver_( solver ), + previous_( previous ), + current_( current ) + {} + + /// return || u^n+1 - u^n || / || u^n+1 || + double relativeChange() const + { + return solver_.model().relativeChange( previous_, current_ ); + } + }; + + template + void logException(const E& exception, bool verbose) + { + if( verbose ) + { + std::ostringstream message; + message << "Caught Exception: " << exception.what(); + OpmLog::debug(message.str()); + } + } + } + + // AdaptiveTimeStepping + //--------------------- + + AdaptiveTimeStepping::AdaptiveTimeStepping( const Tuning& tuning, + size_t time_step, + const parameter::ParameterGroup& param, + const bool terminal_output ) + : timeStepControl_() + , restart_factor_( tuning.getTSFCNV(time_step) ) + , growth_factor_(tuning.getTFDIFF(time_step) ) + , max_growth_( tuning.getTSFMAX(time_step) ) + // default is 1 year, convert to seconds + , max_time_step_( tuning.getTSMAXZ(time_step) ) + , solver_restart_max_( param.getDefault("solver.restart", int(10) ) ) + , solver_verbose_( param.getDefault("solver.verbose", bool(true) ) && terminal_output ) + , timestep_verbose_( param.getDefault("timestep.verbose", bool(true) ) && terminal_output ) + , suggested_next_timestep_( tuning.getTSINIT(time_step) ) + , full_timestep_initially_( param.getDefault("full_timestep_initially", bool(false) ) ) + { + init(param); + + } + + AdaptiveTimeStepping::AdaptiveTimeStepping( const parameter::ParameterGroup& param, + const bool terminal_output ) + : timeStepControl_() + , restart_factor_( param.getDefault("solver.restartfactor", double(0.33) ) ) + , growth_factor_( param.getDefault("solver.growthfactor", double(2) ) ) + , max_growth_( param.getDefault("timestep.control.maxgrowth", double(3.0) ) ) + // default is 1 year, convert to seconds + , max_time_step_( unit::convert::from(param.getDefault("timestep.max_timestep_in_days", 365.0 ), unit::day) ) + , solver_restart_max_( param.getDefault("solver.restart", int(10) ) ) + , solver_verbose_( param.getDefault("solver.verbose", bool(true) ) && terminal_output ) + , timestep_verbose_( param.getDefault("timestep.verbose", bool(true) ) && terminal_output ) + , suggested_next_timestep_( unit::convert::from(param.getDefault("timestep.initial_timestep_in_days", -1.0 ), unit::day) ) + , full_timestep_initially_( param.getDefault("full_timestep_initially", bool(false) ) ) + { + init(param); + } + + void AdaptiveTimeStepping:: + init(const parameter::ParameterGroup& param) + { + // valid are "pid" and "pid+iteration" + std::string control = param.getDefault("timestep.control", std::string("pid") ); + // iterations is the accumulation of all linear iterations over all newton steops per time step + const int defaultTargetIterations = 30; + + const double tol = param.getDefault("timestep.control.tol", double(1e-1) ); + if( control == "pid" ) { + timeStepControl_ = TimeStepControlType( new PIDTimeStepControl( tol ) ); + } + else if ( control == "pid+iteration" ) + { + const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations ); + timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) ); + } + else if ( control == "iterationcount" ) + { + const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetIterations ); + const double decayrate = param.getDefault("timestep.control.decayrate", double(0.75) ); + const double growthrate = param.getDefault("timestep.control.growthrate", double(1.25) ); + timeStepControl_ = TimeStepControlType( new SimpleIterationCountTimeStepControl( iterations, decayrate, growthrate ) ); + } else if ( control == "hardcoded") { + const std::string filename = param.getDefault("timestep.control.filename", std::string("timesteps")); + timeStepControl_ = TimeStepControlType( new HardcodedTimeStepControl( filename ) ); + + } + else + OPM_THROW(std::runtime_error,"Unsupported time step control selected "<< control ); + + // make sure growth factor is something reasonable + assert( growth_factor_ >= 1.0 ); + } + + + + template + SimulatorReport AdaptiveTimeStepping:: + step( const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state ) + { + return stepImpl( simulatorTimer, solver, state, well_state, nullptr, nullptr ); + } + + template + SimulatorReport AdaptiveTimeStepping:: + step( const SimulatorTimer& simulatorTimer, + Solver& solver, State& state, WellState& well_state, + Output& outputWriter, + const std::vector* fipnum) + { + return stepImpl( simulatorTimer, solver, state, well_state, &outputWriter, fipnum ); + } + + + // implementation of the step method + template + SimulatorReport AdaptiveTimeStepping:: + stepImpl( const SimulatorTimer& simulatorTimer, + Solver& solver, State& state, WState& well_state, + Output* outputWriter, + const std::vector* fipnum) + { + SimulatorReport report; + const double timestep = simulatorTimer.currentStepLength(); + + // init last time step as a fraction of the given time step + if( suggested_next_timestep_ < 0 ) { + suggested_next_timestep_ = restart_factor_ * timestep; + } + + if (full_timestep_initially_) { + suggested_next_timestep_ = timestep; + } + + // TODO + // take change in well state into account + + // create adaptive step timer with previously used sub step size + AdaptiveSimulatorTimer substepTimer( simulatorTimer, suggested_next_timestep_, max_time_step_ ); + + // copy states in case solver has to be restarted (to be revised) + State last_state( state ); + WState last_well_state( well_state ); + + // counter for solver restarts + int restarts = 0; + + // sub step time loop + while( ! substepTimer.done() ) + { + // get current delta t + const double dt = substepTimer.currentStepLength() ; + if( timestep_verbose_ ) + { + std::ostringstream ss; + ss <<" Substep " << substepTimer.currentStepNum() << ", stepsize " + << unit::convert::to(substepTimer.currentStepLength(), unit::day) << " days."; + OpmLog::info(ss.str()); + } + + SimulatorReport substepReport; + try { + substepReport = solver.step( substepTimer, state, well_state); + report += substepReport; + + if( solver_verbose_ ) { + // report number of linear iterations + OpmLog::note("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations)); + } + } + catch (const Opm::NumericalProblem& e) { + detail::logException(e, solver_verbose_); + // since linearIterations is < 0 this will restart the solver + } + catch (const std::runtime_error& e) { + detail::logException(e, solver_verbose_); + // also catch linear solver not converged + } + catch (const Dune::ISTLError& e) { + detail::logException(e, solver_verbose_); + // also catch errors in ISTL AMG that occur when time step is too large + } + catch (const Dune::MatrixBlockError& e) { + detail::logException(e, solver_verbose_); + // this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError + } + + if( substepReport.converged ) + { + // advance by current dt + ++substepTimer; + + // create object to compute the time error, simply forwards the call to the model + detail::SolutionTimeErrorSolverWrapper< Solver, State > + relativeChange( solver, last_state, state ); + + // compute new time step estimate + double dtEstimate = + timeStepControl_->computeTimeStepSize( dt, substepReport.total_linear_iterations, relativeChange, substepTimer.simulationTimeElapsed()); + + // limit the growth of the timestep size by the growth factor + dtEstimate = std::min( dtEstimate, double(max_growth_ * dt) ); + + // further restrict time step size growth after convergence problems + if( restarts > 0 ) { + dtEstimate = std::min( growth_factor_ * dt, dtEstimate ); + // solver converged, reset restarts counter + restarts = 0; + } + + if( timestep_verbose_ ) + { + std::ostringstream ss; + ss << " Substep summary: "; + if (report.total_well_iterations != 0) { + ss << "well iterations = " << report.total_well_iterations << ", "; + } + ss << "newton iterations = " << report.total_newton_iterations << ", " + << "linearizations = " << report.total_linearizations + << " (" << report.assemble_time << " sec), " + << "linear iterations = " << report.total_linear_iterations + << " (" << report.linear_solve_time << " sec)"; + OpmLog::info(ss.str()); + } + + // write data if outputWriter was provided + // if the time step is done we do not need + // to write it as this will be done by the simulator + // anyway. + if( outputWriter && !substepTimer.done() ) { + if (fipnum) { + solver.computeFluidInPlace(state, *fipnum); + } + Opm::time::StopWatch perfTimer; + perfTimer.start(); + bool substep = true; + const auto& physicalModel = solver.model(); + outputWriter->writeTimeStep( substepTimer, state, well_state, physicalModel, substep); + report.output_write_time += perfTimer.secsSinceStart(); + } + + // set new time step length + substepTimer.provideTimeStepEstimate( dtEstimate ); + + // update states + last_state = state ; + last_well_state = well_state; + + report.converged = substepTimer.done(); + substepTimer.setLastStepFailed(false); + + } + else // in case of no convergence (linearIterations < 0) + { + report.converged = false; + substepTimer.setLastStepFailed(true); + // increase restart counter + if( restarts >= solver_restart_max_ ) { + const auto msg = std::string("Solver failed to converge after ") + + std::to_string(restarts) + " restarts."; + if (solver_verbose_) { + OpmLog::error(msg); + } + OPM_THROW_NOLOG(Opm::NumericalProblem, msg); + } + + const double newTimeStep = restart_factor_ * dt; + // we need to revise this + substepTimer.provideTimeStepEstimate( newTimeStep ); + if( solver_verbose_ ) { + std::string msg; + msg = "Solver convergence failed, restarting solver with new time step (" + + std::to_string(unit::convert::to( newTimeStep, unit::day )) + " days).\n"; + OpmLog::problem(msg); + } + // reset states + state = last_state; + well_state = last_well_state; + + ++restarts; + } + } + + + // store estimated time step for next reportStep + suggested_next_timestep_ = substepTimer.currentStepLength(); + if( timestep_verbose_ ) + { + std::ostringstream ss; + substepTimer.report(ss); + ss << "Suggested next step size = " << unit::convert::to( suggested_next_timestep_, unit::day ) << " (days)" << std::endl; + OpmLog::note(ss.str()); + } + + if( ! std::isfinite( suggested_next_timestep_ ) ) { // check for NaN + suggested_next_timestep_ = timestep; + } + return report; + } +} + +#endif diff --git a/opm/simulators/timestepping/SimulatorTimer.cpp b/opm/simulators/timestepping/SimulatorTimer.cpp new file mode 100644 index 000000000..cedcaeefe --- /dev/null +++ b/opm/simulators/timestepping/SimulatorTimer.cpp @@ -0,0 +1,171 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + 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 . +*/ + +#include "config.h" +#include +#include +#include +#include +#include + +namespace Opm +{ + + /// Default constructor. + SimulatorTimer::SimulatorTimer() + : current_step_(0), + current_time_(0.0), + start_date_(2012,1,1) // A really arbitrary default starting value?! + { + } + + /// Initialize from parameters. Accepts the following: + /// num_psteps (default 1) + /// stepsize_days (default 1) + void SimulatorTimer::init(const parameter::ParameterGroup& param) + { + const int num_psteps = param.getDefault("num_psteps", 1); + const double stepsize_days = param.getDefault("stepsize_days", 1.0); + const double stepsize = Opm::unit::convert::from(stepsize_days, Opm::unit::day); + timesteps_.clear(); + timesteps_.resize(num_psteps, stepsize); + total_time_ = num_psteps*stepsize; + } + + /// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap + void SimulatorTimer::init(const TimeMap& timeMap, size_t report_step) + { + total_time_ = timeMap.getTotalTime(); + timesteps_.resize(timeMap.numTimesteps()); + for ( size_t i = 0; i < timeMap.numTimesteps(); ++i ) { + timesteps_[i] = timeMap.getTimeStepLength(i); + } + + setCurrentStepNum(report_step); + + boost::posix_time::ptime start_time = timeMap.getStartTime(0); + start_date_ = start_time.date(); + } + + /// Whether the current step is the first step. + bool SimulatorTimer::initialStep() const + { + return (current_step_ == 0); + } + + /// Total number of steps. + int SimulatorTimer::numSteps() const + { + return timesteps_.size(); + } + + /// Current step number. + int SimulatorTimer::currentStepNum() const + { + return current_step_; + } + + /// Set current step number. + void SimulatorTimer::setCurrentStepNum(int step) + { + current_step_ = step; + current_time_ = std::accumulate(timesteps_.begin(), timesteps_.begin() + step, 0.0); + } + + + /// Current step length. + double SimulatorTimer::currentStepLength() const + { + assert(!done()); + return timesteps_[current_step_]; + } + + double SimulatorTimer::stepLengthTaken() const + { + assert(current_step_ > 0); + return timesteps_[current_step_ - 1]; + } + + /// time elapsed since the start of the simulation [s]. + double SimulatorTimer::simulationTimeElapsed() const + { + return current_time_; + } + + boost::posix_time::ptime SimulatorTimer::startDateTime() const + { + return boost::posix_time::ptime(start_date_); + } + + + boost::posix_time::ptime SimulatorTimer::currentDateTime() const + { + return startDateTime() + boost::posix_time::seconds( (int) simulationTimeElapsed()); + } + + /// Total time. + double SimulatorTimer::totalTime() const + { + return total_time_; + } + + /// Set total time. + /// This is primarily intended for multi-epoch schedules, + /// where a timer for a given epoch does not have + /// access to later timesteps. + void SimulatorTimer::setTotalTime(double time) + { + total_time_ = time; + } + + /// Print a report with current and total time etc. + void SimulatorTimer::report(std::ostream& os) const + { + os << "\n\n--------------- Simulation step number " << currentStepNum() << " ---------------" + << "\n Current time (days) " << Opm::unit::convert::to(simulationTimeElapsed(), Opm::unit::day) + << "\n Current stepsize (days) " << Opm::unit::convert::to(currentStepLength(), Opm::unit::day) + << "\n Total time (days) " << Opm::unit::convert::to(totalTime(), Opm::unit::day) + << "\n" << std::endl; + } + + /// Next step. + SimulatorTimer& SimulatorTimer::operator++() + { + assert(!done()); + current_time_ += timesteps_[current_step_]; + ++current_step_; + return *this; + } + + /// Return true if op++() has been called numSteps() times. + bool SimulatorTimer::done() const + { + return int(timesteps_.size()) == current_step_; + } + + /// return copy of object + std::unique_ptr< SimulatorTimerInterface > + SimulatorTimer::clone() const + { + return std::unique_ptr< SimulatorTimerInterface > (new SimulatorTimer( *this )); + } + + + +} // namespace Opm diff --git a/opm/simulators/timestepping/SimulatorTimer.hpp b/opm/simulators/timestepping/SimulatorTimer.hpp new file mode 100644 index 000000000..0caf225a5 --- /dev/null +++ b/opm/simulators/timestepping/SimulatorTimer.hpp @@ -0,0 +1,131 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + 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_SIMULATORTIMER_HEADER_INCLUDED +#define OPM_SIMULATORTIMER_HEADER_INCLUDED + +#include +#include + +#include +#include + +namespace Opm +{ + + namespace parameter { class ParameterGroup; } + + class SimulatorTimer : public SimulatorTimerInterface + { + public: + // use default implementation of these methods + using SimulatorTimerInterface::currentDateTime; + using SimulatorTimerInterface::currentPosixTime; + + /// Default constructor. + SimulatorTimer(); + + /// Initialize from parameters. Accepts the following: + /// num_psteps (default 1) + /// stepsize_days (default 1) + void init(const parameter::ParameterGroup& param); + + /// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap + void init(const TimeMap& timeMap, size_t report_step = 0); + + /// Whether the current step is the first step. + bool initialStep() const; + + /// Total number of steps. + int numSteps() const; + + /// 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. + int currentStepNum() const; + + /// Set current step number. + void setCurrentStepNum(int step); + + /// 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(). + double currentStepLength() const; + + /// 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(). + double stepLengthTaken () const; + + /// Time elapsed since the start of the simulation until the + /// beginning of the current time step [s]. + double simulationTimeElapsed() const; + + /// Total time. + double totalTime() const; + + /// Return start date of simulation + boost::posix_time::ptime startDateTime() const; + + /// Return current date. + boost::posix_time::ptime currentDateTime() const; + + /// Set total time. + /// This is primarily intended for multi-epoch schedules, + /// where a timer for a given epoch does not have + /// access to later timesteps. + void setTotalTime(double time); + + /// Print a report with current and total time etc. + /// Note: if done(), it is an error to call report(). + void report(std::ostream& os) const; + + /// advance time by currentStepLength + SimulatorTimer& operator++(); + + /// advance time by currentStepLength + void advance() { this->operator++(); } + + /// Return true if op++() has been called numSteps() times. + bool done() const; + + /// Always return false. Timestep failures is handled in the + /// substepTimer + bool lastStepFailed() const {return false;} + + /// return copy of object + virtual std::unique_ptr< SimulatorTimerInterface > clone() const; + + private: + std::vector timesteps_; + int current_step_; + double current_time_; + double total_time_; + boost::gregorian::date start_date_; + }; + + +} // namespace Opm + +#endif // OPM_SIMULATORTIMER_HEADER_INCLUDED diff --git a/opm/simulators/timestepping/SimulatorTimerInterface.hpp b/opm/simulators/timestepping/SimulatorTimerInterface.hpp new file mode 100644 index 000000000..ae30d381c --- /dev/null +++ b/opm/simulators/timestepping/SimulatorTimerInterface.hpp @@ -0,0 +1,114 @@ +/* + 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 +#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; + + /// advance time by currentStepLength + virtual void advance() = 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); + } + + /// Return true if last time step failed + virtual bool lastStepFailed() const = 0; + + /// return copy of current timer instance + virtual std::unique_ptr< SimulatorTimerInterface > clone () const = 0; + }; + + +} // namespace Opm + +#endif // OPM_SIMULATORTIMER_HEADER_INCLUDED diff --git a/opm/simulators/timestepping/TimeStepControl.cpp b/opm/simulators/timestepping/TimeStepControl.cpp new file mode 100644 index 000000000..ebbfaa428 --- /dev/null +++ b/opm/simulators/timestepping/TimeStepControl.cpp @@ -0,0 +1,192 @@ +/* + Copyright 2014 IRIS AS + Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015 Statoil 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 . +*/ +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace Opm +{ + //////////////////////////////////////////////////////// + // + // InterationCountTimeStepControl Implementation + // + //////////////////////////////////////////////////////// + + SimpleIterationCountTimeStepControl:: + SimpleIterationCountTimeStepControl( const int target_iterations, + const double decayrate, + const double growthrate, + const bool verbose) + : target_iterations_( target_iterations ) + , decayrate_( decayrate ) + , growthrate_( growthrate ) + , verbose_( verbose ) + { + if( decayrate_ > 1.0 ) { + OPM_THROW(std::runtime_error,"SimpleIterationCountTimeStepControl: decay should be <= 1 " << decayrate_ ); + } + if( growthrate_ < 1.0 ) { + OPM_THROW(std::runtime_error,"SimpleIterationCountTimeStepControl: growth should be >= 1 " << growthrate_ ); + } + } + + double SimpleIterationCountTimeStepControl:: + computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */, const double /*simulationTimeElapsed */) const + { + double dtEstimate = dt ; + + // reduce the time step size if we exceed the number of target iterations + if( iterations > target_iterations_ ) + { + // scale dtEstimate down with a given rate + dtEstimate *= decayrate_; + } + // increase the time step size if we are below the number of target iterations + else if ( iterations < target_iterations_-1 ) + { + dtEstimate *= growthrate_; + } + + return dtEstimate; + } + + //////////////////////////////////////////////////////// + // + // HardcodedTimeStepControl Implementation + // + //////////////////////////////////////////////////////// + + HardcodedTimeStepControl:: + HardcodedTimeStepControl( const std::string& filename) + { + std::ifstream infile (filename); + if (!infile.is_open()) { + OPM_THROW(std::runtime_error,"Incorrect or no filename is provided to the hardcodedTimeStep. Use timestep.control.filename=your_file_name"); + } + std::string::size_type sz; + std::string line; + while ( std::getline(infile, line)) { + if( line[0] != '-') { // ignore lines starting with '-' + const double time = std::stod(line,&sz); // read the first number i.e. the actual substep time + subStepTime_.push_back( time * unit::day ); + } + + } + } + + double HardcodedTimeStepControl:: + computeTimeStepSize( const double /*dt */, const int /*iterations */, const RelativeChangeInterface& /* relativeChange */ , const double simulationTimeElapsed) const + { + auto nextTime = std::upper_bound(subStepTime_.begin(), subStepTime_.end(), simulationTimeElapsed); + return (*nextTime - simulationTimeElapsed); + } + + + + //////////////////////////////////////////////////////// + // + // PIDTimeStepControl Implementation + // + //////////////////////////////////////////////////////// + + PIDTimeStepControl::PIDTimeStepControl( const double tol, + const bool verbose ) + : tol_( tol ) + , errors_( 3, tol_ ) + , verbose_( verbose ) + {} + + double PIDTimeStepControl:: + computeTimeStepSize( const double dt, const int /* iterations */, const RelativeChangeInterface& relChange, const double /*simulationTimeElapsed */) const + { + // shift errors + for( int i=0; i<2; ++i ) { + errors_[ i ] = errors_[i+1]; + } + + // store new error + const double error = relChange.relativeChange(); + errors_[ 2 ] = error; + + if( error > tol_ ) + { + // adjust dt by given tolerance + const double newDt = dt * tol_ / error; + if( verbose_ ) + std::cout << "Computed step size (tol): " << unit::convert::to( newDt, unit::day ) << " (days)" << std::endl; + return newDt; + } + else + { + // values taking from turek time stepping paper + const double kP = 0.075 ; + const double kI = 0.175 ; + const double kD = 0.01 ; + const double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) * + std::pow( tol_ / errors_[ 2 ], kI ) * + std::pow( errors_[0]*errors_[0]/errors_[ 1 ]/errors_[ 2 ], kD )); + if( verbose_ ) + std::cout << "Computed step size (pow): " << unit::convert::to( newDt, unit::day ) << " (days)" << std::endl; + return newDt; + } + } + + + + //////////////////////////////////////////////////////////// + // + // PIDAndIterationCountTimeStepControl Implementation + // + //////////////////////////////////////////////////////////// + + PIDAndIterationCountTimeStepControl:: + PIDAndIterationCountTimeStepControl( const int target_iterations, + const double tol, + const bool verbose) + : BaseType( tol, verbose ) + , target_iterations_( target_iterations ) + {} + + double PIDAndIterationCountTimeStepControl:: + computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relChange, const double simulationTimeElapsed ) const + { + double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, relChange, simulationTimeElapsed); + + // further reduce step size if to many iterations were used + if( iterations > target_iterations_ ) + { + // if iterations was the same or dts were the same, do some magic + dtEstimate *= double( target_iterations_ ) / double(iterations); + } + + return dtEstimate; + } + +} // end namespace Opm diff --git a/opm/simulators/timestepping/TimeStepControl.hpp b/opm/simulators/timestepping/TimeStepControl.hpp new file mode 100644 index 000000000..e0e202dbd --- /dev/null +++ b/opm/simulators/timestepping/TimeStepControl.hpp @@ -0,0 +1,147 @@ +/* + Copyright 2014 IRIS AS + Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015 Statoil 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_TIMESTEPCONTROL_HEADER_INCLUDED +#define OPM_TIMESTEPCONTROL_HEADER_INCLUDED + +#include + +#include +#include +#include + +namespace Opm +{ + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// + /// A simple iteration count based adaptive time step control. + // + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + class SimpleIterationCountTimeStepControl : public TimeStepControlInterface + { + public: + /// \brief constructor + /// \param target_iterations number of desired iterations (e.g. Newton iterations) per time step in one time step + // \param decayrate decayrate of time step when target iterations are not met (should be <= 1) + // \param growthrate growthrate of time step when target iterations are not met (should be >= 1) + /// \param verbose if true get some output (default = false) + SimpleIterationCountTimeStepControl( const int target_iterations, + const double decayrate, + const double growthrate, + const bool verbose = false); + + /// \brief \copydoc TimeStepControlInterface::computeTimeStepSize + double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */, const double /*simulationTimeElapsed */ ) const; + + protected: + const int target_iterations_; + const double decayrate_; + const double growthrate_; + const bool verbose_; + }; + + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// + /// PID controller based adaptive time step control as suggested in: + /// Turek and Kuzmin. Algebraic Flux Correction III. Incompressible Flow Problems. Uni Dortmund. + /// + /// See also: + /// D. Kuzmin and S.Turek. Numerical simulation of turbulent bubbly flows. Techreport Uni Dortmund. 2004 + /// + /// and the original article: + /// Valli, Coutinho, and Carey. Adaptive Control for Time Step Selection in Finite Element + /// Simulation of Coupled Viscous Flow and Heat Transfer. Proc of the 10th + /// International Conference on Numerical Methods in Fluids. 1998. + /// + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + class PIDTimeStepControl : public TimeStepControlInterface + { + public: + /// \brief constructor + /// \param tol tolerance for the relative changes of the numerical solution to be accepted + /// in one time step (default is 1e-3) + /// \param verbose if true get some output (default = false) + PIDTimeStepControl( const double tol = 1e-3, + const bool verbose = false ); + + /// \brief \copydoc TimeStepControlInterface::computeTimeStepSize + double computeTimeStepSize( const double dt, const int /* iterations */, const RelativeChangeInterface& relativeChange, const double /*simulationTimeElapsed */ ) const; + + protected: + const double tol_; + mutable std::vector< double > errors_; + + const bool verbose_; + }; + + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// + /// PID controller based adaptive time step control as above that also takes + /// an target iteration into account. + // + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + class PIDAndIterationCountTimeStepControl : public PIDTimeStepControl + { + typedef PIDTimeStepControl BaseType; + public: + /// \brief constructor + /// \param target_iterations number of desired iterations per time step + /// \param tol tolerance for the relative changes of the numerical solution to be accepted + /// in one time step (default is 1e-3) + /// \param verbose if true get some output (default = false) + PIDAndIterationCountTimeStepControl( const int target_iterations = 20, + const double tol = 1e-3, + const bool verbose = false); + + /// \brief \copydoc TimeStepControlInterface::computeTimeStepSize + double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relativeChange, const double /*simulationTimeElapsed */ ) const; + + protected: + const int target_iterations_; + }; + + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// + /// HardcodedTimeStepControl + /// Input generated from summary file using the ert application: + /// + /// ecl_summary DECK TIME > filename + /// + /// Assumes time is given in days + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + class HardcodedTimeStepControl : public TimeStepControlInterface + { + public: + /// \brief constructor + /// \param filename filename contaning the timesteps + explicit HardcodedTimeStepControl( const std::string& filename); + + /// \brief \copydoc TimeStepControlInterface::computeTimeStepSize + double computeTimeStepSize( const double dt, const int /* iterations */, const RelativeChangeInterface& /*relativeChange */, const double simulationTimeElapsed) const; + + protected: + // store the time (in days) of the substeps the simulator should use + std::vector subStepTime_; + }; + + +} // end namespace Opm +#endif + diff --git a/opm/simulators/timestepping/TimeStepControlInterface.hpp b/opm/simulators/timestepping/TimeStepControlInterface.hpp new file mode 100644 index 000000000..bdab5e3a4 --- /dev/null +++ b/opm/simulators/timestepping/TimeStepControlInterface.hpp @@ -0,0 +1,66 @@ +/* + 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_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED +#define OPM_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED + + +namespace Opm +{ + + /////////////////////////////////////////////////////////////////// + /// + /// RelativeChangeInterface + /// + /////////////////////////////////////////////////////////////////// + class RelativeChangeInterface + { + protected: + RelativeChangeInterface() {} + public: + /// \return || u^n+1 - u^n || / || u^n+1 || + virtual double relativeChange() const = 0; + + /// virtual destructor (empty) + virtual ~RelativeChangeInterface() {} + }; + + /////////////////////////////////////////////////////////////////// + /// + /// TimeStepControlInterface + /// + /////////////////////////////////////////////////////////////////// + class TimeStepControlInterface + { + protected: + TimeStepControlInterface() {} + public: + /// compute new time step size suggestions based on the PID controller + /// \param dt time step size used in the current step + /// \param iterations number of iterations used (linear/nonlinear) + /// \param timeError object to compute || u^n+1 - u^n || / || u^n+1 || + /// + /// \return suggested time step size for the next step + virtual double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relativeChange , const double simulationTimeElapsed) const = 0; + + /// virtual destructor (empty) + virtual ~TimeStepControlInterface () {} + }; + +} +#endif diff --git a/tests/TESTTIMER.DATA b/tests/TESTTIMER.DATA new file mode 100644 index 000000000..efe90c9e6 --- /dev/null +++ b/tests/TESTTIMER.DATA @@ -0,0 +1,49 @@ +RUNSPEC + +START + 26 'MAR' 2014 / + +TSTEP +1.0 2*5.0 + / + +TSTEP + 7*10.0 14*25.0 + / + +TSTEP + 19.0 + / + +TSTEP + 18*13.0 + / + +TSTEP + 17*10.0 + / + +TSTEP + 13.0 + / + +TSTEP + 18.0 + / + +TSTEP + 11.0 + / + +TSTEP + 17*5.0 + / + +TSTEP + 19*6.0 + / + +TSTEP + 21*5.0 / + +END == diff --git a/tests/test_timer.cpp b/tests/test_timer.cpp new file mode 100644 index 000000000..68511f444 --- /dev/null +++ b/tests/test_timer.cpp @@ -0,0 +1,115 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + 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 . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE OPM-TimerTest +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +BOOST_AUTO_TEST_CASE(CreateTimer) +{ + const std::string filename1 = "TESTTIMER.DATA"; + Opm::ParseContext parseContext; + Opm::Parser parser; + Opm::Deck parserDeck = parser.parseFile( filename1 , parseContext); + + Opm::TimeMap timeMap( parserDeck ); + Opm::SimulatorTimer simtimer; + + boost::gregorian::date defaultStartDate( 2012, 1, 1); + BOOST_CHECK_EQUAL( boost::posix_time::ptime(defaultStartDate), simtimer.currentDateTime() ); + + simtimer.init(timeMap); + boost::gregorian::date startDate( 2014, 3, 26); + BOOST_CHECK_EQUAL( boost::posix_time::ptime(startDate), simtimer.currentDateTime() ); + + BOOST_CHECK_EQUAL( 0, simtimer.currentStepNum() ); + BOOST_CHECK_EQUAL( 0., simtimer.simulationTimeElapsed() ); + BOOST_CHECK_EQUAL( 120, simtimer.numSteps() ); + BOOST_CHECK_EQUAL( 1200., Opm::unit::convert::to(simtimer.totalTime(), Opm::unit::day) ); + BOOST_CHECK_EQUAL( 0., Opm::unit::convert::to(simtimer.simulationTimeElapsed(), Opm::unit::day) ); + + double testCurrentTime = 0.; + BOOST_CHECK_EQUAL( Opm::unit::convert::to(testCurrentTime, Opm::unit::day), + Opm::unit::convert::to(simtimer.simulationTimeElapsed(), Opm::unit::day) ); + + for ( int i = 0; i < simtimer.numSteps(); ++i ) { + BOOST_CHECK_EQUAL( i, simtimer.currentStepNum() ); + BOOST_CHECK_EQUAL( Opm::unit::convert::to(testCurrentTime, Opm::unit::minute), + Opm::unit::convert::to(simtimer.simulationTimeElapsed(), Opm::unit::minute) ); + testCurrentTime += simtimer.currentStepLength(); + ++simtimer; + } + + for ( int i = 0; i <= simtimer.numSteps(); ++i ) { + simtimer.setCurrentStepNum(i); + BOOST_CHECK_EQUAL( i, simtimer.currentStepNum() ); + } + + BOOST_CHECK_EQUAL( true, simtimer.done() ); + simtimer.setCurrentStepNum(0); + BOOST_CHECK_EQUAL( false, simtimer.done() ); + BOOST_CHECK_EQUAL( 0., Opm::unit::convert::to(simtimer.simulationTimeElapsed(), Opm::unit::day) ); + + simtimer.setCurrentStepNum(120); + BOOST_CHECK_EQUAL( Opm::unit::convert::to(simtimer.simulationTimeElapsed(), Opm::unit::day), + Opm::unit::convert::to(simtimer.totalTime(), Opm::unit::day) ); + + int i = 0; + double testCurrentTime1 = 0.; + double testCurrentTime2 = 0.; + simtimer.setCurrentStepNum(0); + + while (!simtimer.done()) { + testCurrentTime1 += simtimer.currentStepLength(); + BOOST_CHECK_EQUAL( i, simtimer.currentStepNum() ); + ++i; + ++simtimer; + testCurrentTime2 += simtimer.stepLengthTaken(); + BOOST_CHECK_EQUAL( Opm::unit::convert::to(testCurrentTime1, Opm::unit::minute), + Opm::unit::convert::to(simtimer.simulationTimeElapsed(), Opm::unit::minute) ); + BOOST_CHECK_EQUAL( Opm::unit::convert::to(testCurrentTime2, Opm::unit::minute), + Opm::unit::convert::to(simtimer.simulationTimeElapsed(), Opm::unit::minute) ); + } + + BOOST_CHECK_EQUAL( true, simtimer.done() ); + BOOST_CHECK_EQUAL( Opm::unit::convert::to(testCurrentTime1, Opm::unit::minute), + Opm::unit::convert::to(simtimer.totalTime(), Opm::unit::minute) ); + BOOST_CHECK_EQUAL( Opm::unit::convert::to(testCurrentTime2, Opm::unit::minute), + Opm::unit::convert::to(simtimer.totalTime(), Opm::unit::minute) ); + +}