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

remove files

Browse Source
This commit is contained in:
Arne Morten Kvarving
2018-11-14 14:42:52 +01:00
parent 78158bf44a
commit f027262ec4
259 changed files with 1 additions and 72573 deletions
Download Patch File Download Diff File Expand all files Collapse all files

139
opm/simulators/timestepping/AdaptiveTimeStepping.hpp
View File

@@ -1,139 +0,0 @@
/*
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/common/utility/parameters/ParameterGroup.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/simulators/timestepping/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 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 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
\param event event status for possible tuning
*/
template <class Solver, class State, class WellState>
SimulatorReport step( const SimulatorTimer& timer,
Solver& solver, State& state, WellState& well_state,
const bool event);
/** \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 event event status for possible tuning
\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,
const bool event,
Output& outputWriter,
const std::vector<int>* fipnum = nullptr);
/** \brief Returns the simulator report for the failed substeps of the last
* report step.
*/
const SimulatorReport& failureReport() const { return failureReport_; };
double suggestedNextStep() const { return suggested_next_timestep_; }
void setSuggestedNextStep(const double x) { suggested_next_timestep_ = x; }
void updateTUNING(const Tuning& tuning, size_t time_step) {
restart_factor_ = tuning.getTSFCNV(time_step);
growth_factor_ = tuning.getTFDIFF(time_step);
max_growth_ = tuning.getTSFMAX(time_step);
max_time_step_ = tuning.getTSMAXZ(time_step);
suggested_next_timestep_ = tuning.getTSINIT(time_step);
timestep_after_event_ = tuning.getTMAXWC(time_step);
}
protected:
template <class Solver, class State, class WellState, class Output>
SimulatorReport stepImpl( const SimulatorTimer& timer,
Solver& solver, State& state, WellState& well_state,
const bool event,
Output* outputWriter,
const std::vector<int>* fipnum);
void init(const ParameterGroup& param);
typedef std::unique_ptr< TimeStepControlInterface > TimeStepControlType;
SimulatorReport failureReport_; //!< statistics for the failed substeps of the last timestep
TimeStepControlType timeStepControl_; //!< time step control object
double restart_factor_; //!< factor to multiply time step with when solver fails to converge
double growth_factor_; //!< factor to multiply time step when solver recovered from failed convergence
double max_growth_; //!< factor that limits the maximum growth of a time step
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
double timestep_after_event_; //!< suggested size of timestep after an event
bool use_newton_iteration_; //!< use newton iteration count for adaptive time step control
};
}
#include <opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp>
#endif

394
opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp
View File

@@ -1,394 +0,0 @@
/*
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/simulators/timestepping/SimulatorTimer.hpp>
#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
#include <opm/simulators/timestepping/TimeStepControl.hpp>
#include <opm/grid/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::string message;
message = "Caught Exception: ";
message += exception.what();
OpmLog::debug(message);
}
}
}
// AdaptiveTimeStepping
//---------------------
inline AdaptiveTimeStepping::AdaptiveTimeStepping( const Tuning& tuning,
size_t time_step,
const 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) ) )
, timestep_after_event_( tuning.getTMAXWC(time_step))
, use_newton_iteration_(false)
{
init(param);
}
inline AdaptiveTimeStepping::AdaptiveTimeStepping( const 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) ) )
, timestep_after_event_( unit::convert::from(param.getDefault("timestep.timestep_in_days_after_event", -1.0 ), unit::day))
, use_newton_iteration_(false)
{
init(param);
}
inline void AdaptiveTimeStepping::
init(const 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 int defaultTargetNewtonIterations = 8;
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 == "pid+newtoniteration" )
{
const int iterations = param.getDefault("timestep.control.targetiteration", defaultTargetNewtonIterations );
timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
use_newton_iteration_ = true;
}
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, const bool event )
{
return stepImpl( simulatorTimer, solver, state, well_state, event, nullptr, nullptr );
}
template <class Solver, class State, class WellState, class Output>
SimulatorReport AdaptiveTimeStepping::
step( const SimulatorTimer& simulatorTimer,
Solver& solver, State& state, WellState& well_state,
const bool event,
Output& outputWriter,
const std::vector<int>* fipnum)
{
return stepImpl( simulatorTimer, solver, state, well_state, event, &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,
const bool event,
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;
}
// use seperate time step after event
if (event && timestep_after_event_ > 0) {
suggested_next_timestep_ = timestep_after_event_;
}
// 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 );
// reset the statistics for the failed substeps
failureReport_ = SimulatorReport();
// 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 <<"\nTime step " << substepTimer.currentStepNum() << ", stepsize "
<< unit::convert::to(substepTimer.currentStepLength(), unit::day) << " days.";
OpmLog::info(ss.str());
}
SimulatorReport substepReport;
std::string cause_of_failure = "";
try {
substepReport = solver.step( substepTimer, state, well_state);
report += substepReport;
if( solver_verbose_ ) {
// report number of linear iterations
OpmLog::debug("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
}
}
catch (const Opm::TooManyIterations& e) {
substepReport += solver.failureReport();
cause_of_failure = "Solver convergence failure - Iteration limit reached";
detail::logException(e, solver_verbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const Opm::LinearSolverProblem& e) {
substepReport += solver.failureReport();
cause_of_failure = "Linear solver convergence failure";
detail::logException(e, solver_verbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const Opm::NumericalIssue& e) {
substepReport += solver.failureReport();
cause_of_failure = "Solver convergence failure - Numerical problem encountered";
detail::logException(e, solver_verbose_);
// since linearIterations is < 0 this will restart the solver
}
catch (const std::runtime_error& e) {
substepReport += solver.failureReport();
detail::logException(e, solver_verbose_);
// also catch linear solver not converged
}
catch (const Dune::ISTLError& e) {
substepReport += solver.failureReport();
detail::logException(e, solver_verbose_);
// also catch errors in ISTL AMG that occur when time step is too large
}
catch (const Dune::MatrixBlockError& e) {
substepReport += solver.failureReport();
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
const int iterations = use_newton_iteration_ ? substepReport.total_newton_iterations
: substepReport.total_linear_iterations;
double dtEstimate = timeStepControl_->computeTimeStepSize( dt, 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;
substepReport.reportStep(ss);
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, -1.0, substepReport);
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)
{
substepTimer.setLastStepFailed(true);
failureReport_ += substepReport;
// increase restart counter
if( restarts >= solver_restart_max_ ) {
const auto msg = std::string("Solver failed to converge after cutting timestep ")
+ std::to_string(restarts) + " times.";
if (solver_verbose_) {
OpmLog::error(msg);
}
OPM_THROW_NOLOG(Opm::NumericalIssue, msg);
}
const double newTimeStep = restart_factor_ * dt;
// we need to revise this
substepTimer.provideTimeStepEstimate( newTimeStep );
if( solver_verbose_ ) {
std::string msg;
msg = cause_of_failure + "\nTimestep chopped to "
+ std::to_string(unit::convert::to( substepTimer.currentStepLength(), 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::debug(ss.str());
}
if( ! std::isfinite( suggested_next_timestep_ ) ) { // check for NaN
suggested_next_timestep_ = timestep;
}
return report;
}
}
#endif