OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1049 from akva2/import_adaptive_time_stepping

Browse Source 
Import adaptive time stepping and associated classes
This commit is contained in:
Atgeirr Flø Rasmussen 2017-02-10 15:17:02 +01:00 committed by GitHub
commit 12adcdf218
18 changed files with 1906 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
CMakeLists_files.cmake
View File

@ -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
)

144
examples/wells_example.cpp Normal file
View File

@ -0,0 +1,144 @@
#include "config.h"
#include <iostream>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/wells.h>
#include <opm/core/grid.h>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
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<std::string > ("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<double>("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<double> 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<int> 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<double> 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<double> 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<double> 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;
}

2
opm/autodiff/FlowMain.hpp
View File

@ -48,7 +48,7 @@
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/simulators/thresholdPressures.hpp> // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix.

6
opm/autodiff/SimulatorBase.hpp
View File

@ -40,8 +40,8 @@
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
@ -49,7 +49,7 @@
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>

2
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -32,7 +32,7 @@
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/SimFIBODetails.hpp>
#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
#include <opm/core/utility/initHydroCarbonState.hpp>
#include <opm/core/utility/StopWatch.hpp>

4
opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp
View File

@ -42,14 +42,14 @@
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
//#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
//#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
//#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
//#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
//#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>

166
opm/simulators/timestepping/AdaptiveSimulatorTimer.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <cassert>
#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
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<steps_.size(); ++i )
{
os << " step[ " << i << " ] = " << unit::convert::to( steps_[ i ], unit::day ) << " (days)" << std::endl;
}
os << "sub steps end time = " << unit::convert::to( simulationTimeElapsed(), unit::day ) << " (days)" << std::endl;
}
boost::posix_time::ptime AdaptiveSimulatorTimer::startDateTime() const
{
return start_date_time_;
}
/// return copy of object
std::unique_ptr< SimulatorTimerInterface >
AdaptiveSimulatorTimer::clone() const
{
return std::unique_ptr< SimulatorTimerInterface > (new AdaptiveSimulatorTimer( *this ));
}
} // namespace Opm

123
opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ADAPTIVESIMULATORTIMER_HEADER_INCLUDED
#define OPM_ADAPTIVESIMULATORTIMER_HEADER_INCLUDED
#include <cassert>
#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric>
#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
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<double>::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

112
opm/simulators/timestepping/AdaptiveTimeStepping.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SUBSTEPPING_HEADER_INCLUDED
#define OPM_SUBSTEPPING_HEADER_INCLUDED
#include <iostream>
#include <utility>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
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 <class Solver, class State, class WellState>
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 <class Solver, class State, class WellState, class Output>
SimulatorReport step( const SimulatorTimer& timer,
Solver& solver, State& state, WellState& well_state,
Output& outputWriter,
const std::vector<int>* fipnum = nullptr);
protected:
template <class Solver, class State, class WellState, class Output>
SimulatorReport stepImpl( const SimulatorTimer& timer,
Solver& solver, State& state, WellState& well_state,
Output* outputWriter,
const std::vector<int>* 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 <opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp>
#endif

356
opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
#define OPM_ADAPTIVETIMESTEPPING_IMPL_HEADER_INCLUDED
#include <iostream>
#include <string>
#include <utility>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
#include <opm/simulators/timestepping/TimeStepControl.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <dune/istl/istlexception.hh>
#include <dune/istl/ilu.hh> // For MatrixBlockException
#include <opm/parser/eclipse/EclipseState/Schedule/Tuning.hpp>
namespace Opm {
namespace detail
{
template <class Solver, class State>
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<class E>
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 <class Solver, class State, class WellState>
SimulatorReport AdaptiveTimeStepping::
step( const SimulatorTimer& simulatorTimer, Solver& solver, State& state, WellState& well_state )
{
return stepImpl( simulatorTimer, solver, state, well_state, nullptr, nullptr );
}
template <class Solver, class State, class WellState, class Output>
SimulatorReport AdaptiveTimeStepping::
step( const SimulatorTimer& simulatorTimer,
Solver& solver, State& state, WellState& well_state,
Output& outputWriter,
const std::vector<int>* fipnum)
{
return stepImpl( simulatorTimer, solver, state, well_state, &outputWriter, fipnum );
}
// implementation of the step method
template <class Solver, class State, class WState, class Output >
SimulatorReport AdaptiveTimeStepping::
stepImpl( const SimulatorTimer& simulatorTimer,
Solver& solver, State& state, WState& well_state,
Output* outputWriter,
const std::vector<int>* 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

171
opm/simulators/timestepping/SimulatorTimer.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <ostream>
#include <numeric>
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

131
opm/simulators/timestepping/SimulatorTimer.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SIMULATORTIMER_HEADER_INCLUDED
#define OPM_SIMULATORTIMER_HEADER_INCLUDED
#include <opm/parser/eclipse/EclipseState/Schedule/TimeMap.hpp>
#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
#include <iosfwd>
#include <vector>
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<double> timesteps_;
int current_step_;
double current_time_;
double total_time_;
boost::gregorian::date start_date_;
};
} // namespace Opm
#endif // OPM_SIMULATORTIMER_HEADER_INCLUDED

114
opm/simulators/timestepping/SimulatorTimerInterface.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SIMULATORTIMERINTERFACE_HEADER_INCLUDED
#define OPM_SIMULATORTIMERINTERFACE_HEADER_INCLUDED
#include <memory>
#include <boost/date_time/gregorian/gregorian.hpp>
#include <boost/date_time/posix_time/posix_time_types.hpp>
#include <boost/date_time/posix_time/conversion.hpp>
namespace Opm
{
namespace parameter { class ParameterGroup; }
/// Interface class for SimulatorTimer objects, to be improved.
class SimulatorTimerInterface
{
protected:
/// Default constructor, protected to not allow explicit instances of this class.
SimulatorTimerInterface() {}
public:
/// destructor
virtual ~SimulatorTimerInterface() {}
/// Current step number. This is the number of timesteps that
/// has been completed from the start of the run. The time
/// after initialization but before the simulation has started
/// is timestep number zero.
virtual int currentStepNum() const = 0;
/// Current report step number. This might differ from currentStepNum in case of sub stepping
virtual int reportStepNum() const { return currentStepNum(); }
/// Current step length. This is the length of the step
/// the simulator will take in the next iteration.
///
/// @note if done(), it is an error to call currentStepLength().
virtual double currentStepLength() const = 0;
/// Previous step length. This is the length of the step that
/// was taken to arrive at this time.
///
/// @note if no increments have been done (i.e. the timer is
/// still in its constructed state and currentStepNum() == 0),
/// it is an error to call stepLengthTaken().
virtual double stepLengthTaken () const = 0;
/// Previous report step length. This is the length of the step that
/// was taken to arrive at this report step time.
///
/// @note if no increments have been done (i.e. the timer is
/// still in its constructed state and reportStepNum() == 0),
/// it is an error to call stepLengthTaken().
virtual double reportStepLengthTaken () const { return stepLengthTaken(); }
/// Time elapsed since the start of the simulation until the
/// beginning of the current time step [s].
virtual double simulationTimeElapsed() const = 0;
/// 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

192
opm/simulators/timestepping/TimeStepControl.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <cassert>
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <string>
#include <fstream>
#include <iostream>
#include <opm/common/ErrorMacros.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/simulators/timestepping/TimeStepControl.hpp>
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

147
opm/simulators/timestepping/TimeStepControl.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TIMESTEPCONTROL_HEADER_INCLUDED
#define OPM_TIMESTEPCONTROL_HEADER_INCLUDED
#include <vector>
#include <boost/any.hpp>
#include <boost/range/iterator_range.hpp>
#include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
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<double> subStepTime_;
};
} // end namespace Opm
#endif

66
opm/simulators/timestepping/TimeStepControlInterface.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#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

49
tests/TESTTIMER.DATA Normal file
View File

@ -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 ==

115
tests/test_timer.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#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 <boost/test/unit_test.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <string>
#include <iostream>
#include <vector>
#include <memory>
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) );
}