mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-09 10:15:34 -06:00
AdaptiveTimeStepping: pass object to compute time error to time step control. This
allows us to shift the computation of the error to the physical model.
This commit is contained in:
Robert Kloefkorn
parent
2399629ebc
commit
005a0c2930
@ -43,7 +43,6 @@ namespace Opm {
|
||||
//! \param pinfo The information about the data distribution
|
||||
//! and communication for a parallel run.
|
||||
AdaptiveTimeStepping( const parameter::ParameterGroup& param,
|
||||
const boost::any& pinfo=boost::any(),
|
||||
const bool terminal_output = true );
|
||||
|
||||
/** \brief step method that acts like the solver::step method
|
||||
|
||||
@ -31,11 +31,35 @@
|
||||
|
||||
namespace Opm {
|
||||
|
||||
namespace detail
|
||||
{
|
||||
template <class Solver, class State>
|
||||
class SolutionTimeErrorSolverWrapper : public SolutionTimeErrorInterface
|
||||
{
|
||||
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 timeError() const
|
||||
{
|
||||
return solver_.model().computeTimeError( previous_, current_ );
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
// AdaptiveTimeStepping
|
||||
//---------------------
|
||||
|
||||
AdaptiveTimeStepping::AdaptiveTimeStepping( const parameter::ParameterGroup& param,
|
||||
const boost::any& parallel_information,
|
||||
const bool terminal_output )
|
||||
: timeStepControl_()
|
||||
, restart_factor_( param.getDefault("solver.restartfactor", double(0.33) ) )
|
||||
@ -56,12 +80,12 @@ namespace Opm {
|
||||
|
||||
const double tol = param.getDefault("timestep.control.tol", double(1e-1) );
|
||||
if( control == "pid" ) {
|
||||
timeStepControl_ = TimeStepControlType( new PIDTimeStepControl( tol, parallel_information ) );
|
||||
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, parallel_information ) );
|
||||
timeStepControl_ = TimeStepControlType( new PIDAndIterationCountTimeStepControl( iterations, tol ) );
|
||||
}
|
||||
else if ( control == "iterationcount" )
|
||||
{
|
||||
@ -130,9 +154,6 @@ namespace Opm {
|
||||
// get current delta t
|
||||
const double dt = substepTimer.currentStepLength() ;
|
||||
|
||||
// initialize time step control in case current state is needed later
|
||||
timeStepControl_->initialize( state );
|
||||
|
||||
if( timestep_verbose_ )
|
||||
{
|
||||
std::cout <<"Substep( " << substepTimer.currentStepNum() << " ), try with stepsize "
|
||||
@ -172,9 +193,13 @@ namespace Opm {
|
||||
// advance by current dt
|
||||
++substepTimer;
|
||||
|
||||
// create object to compute the time error, simply forwards the call to the model
|
||||
detail::SolutionTimeErrorSolverWrapper< Solver, State >
|
||||
timeError( solver, last_state, state );
|
||||
|
||||
// compute new time step estimate
|
||||
double dtEstimate =
|
||||
timeStepControl_->computeTimeStepSize( dt, linearIterations, state );
|
||||
timeStepControl_->computeTimeStepSize( dt, linearIterations, timeError );
|
||||
|
||||
// limit the growth of the timestep size by the growth factor
|
||||
dtEstimate = std::min( dtEstimate, double(max_growth_ * dt) );
|
||||
|
||||
@ -55,7 +55,7 @@ namespace Opm
|
||||
}
|
||||
|
||||
double SimpleIterationCountTimeStepControl::
|
||||
computeTimeStepSize( const double dt, const int iterations, const SimulatorState& /* state */ ) const
|
||||
computeTimeStepSize( const double dt, const int iterations, const SolutionTimeErrorInterface& /* timeError */ ) const
|
||||
{
|
||||
double dtEstimate = dt ;
|
||||
|
||||
@ -82,58 +82,24 @@ namespace Opm
|
||||
//
|
||||
////////////////////////////////////////////////////////
|
||||
|
||||
PIDTimeStepControl::PIDTimeStepControl( const double tol, const boost::any& pinfo,
|
||||
PIDTimeStepControl::PIDTimeStepControl( const double tol,
|
||||
const bool verbose )
|
||||
: p0_()
|
||||
, sat0_()
|
||||
, tol_( tol )
|
||||
: tol_( tol )
|
||||
, errors_( 3, tol_ )
|
||||
, verbose_( verbose )
|
||||
, parallel_information_(pinfo)
|
||||
{}
|
||||
|
||||
void PIDTimeStepControl::initialize( const SimulatorState& state )
|
||||
{
|
||||
// store current state for later time step computation
|
||||
p0_ = state.pressure();
|
||||
sat0_ = state.saturation();
|
||||
}
|
||||
|
||||
double PIDTimeStepControl::
|
||||
computeTimeStepSize( const double dt, const int /* iterations */, const SimulatorState& state ) const
|
||||
computeTimeStepSize( const double dt, const int /* iterations */, const SolutionTimeErrorInterface& errorObj ) const
|
||||
{
|
||||
const std::size_t pSize = p0_.size();
|
||||
assert( state.pressure().size() == pSize );
|
||||
const std::size_t satSize = sat0_.size();
|
||||
assert( state.saturation().size() == satSize );
|
||||
|
||||
// compute u^n - u^n+1
|
||||
for( std::size_t i=0; i<pSize; ++i ) {
|
||||
p0_[ i ] -= state.pressure()[ i ];
|
||||
}
|
||||
|
||||
for( std::size_t i=0; i<satSize; ++i ) {
|
||||
sat0_[ i ] -= state.saturation()[ i ];
|
||||
}
|
||||
|
||||
// compute || u^n - u^n+1 ||
|
||||
const double stateOld = euclidianNormSquared( p0_.begin(), p0_.end() ) +
|
||||
euclidianNormSquared( sat0_.begin(), sat0_.end(),
|
||||
state.numPhases() );
|
||||
|
||||
// compute || u^n+1 ||
|
||||
const double stateNew = euclidianNormSquared( state.pressure().begin(), state.pressure().end() ) +
|
||||
euclidianNormSquared( state.saturation().begin(), state.saturation().end(),
|
||||
state.numPhases() );
|
||||
|
||||
// shift errors
|
||||
for( int i=0; i<2; ++i ) {
|
||||
errors_[ i ] = errors_[i+1];
|
||||
}
|
||||
|
||||
// store new error
|
||||
const double error = stateOld / stateNew;
|
||||
errors_[ 2 ] = error ;
|
||||
const double error = errorObj.timeError();
|
||||
errors_[ 2 ] = error;
|
||||
|
||||
if( error > tol_ )
|
||||
{
|
||||
@ -169,16 +135,15 @@ namespace Opm
|
||||
PIDAndIterationCountTimeStepControl::
|
||||
PIDAndIterationCountTimeStepControl( const int target_iterations,
|
||||
const double tol,
|
||||
const boost::any& pinfo,
|
||||
const bool verbose)
|
||||
: BaseType( tol, pinfo, verbose )
|
||||
: BaseType( tol, verbose )
|
||||
, target_iterations_( target_iterations )
|
||||
{}
|
||||
|
||||
double PIDAndIterationCountTimeStepControl::
|
||||
computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const
|
||||
computeTimeStepSize( const double dt, const int iterations, const SolutionTimeErrorInterface& errorObj ) const
|
||||
{
|
||||
double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, state );
|
||||
double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, errorObj );
|
||||
|
||||
// further reduce step size if to many iterations were used
|
||||
if( iterations > target_iterations_ )
|
||||
|
||||
@ -49,7 +49,7 @@ namespace Opm
|
||||
const bool verbose = false);
|
||||
|
||||
/// \brief \copydoc TimeStepControlInterface::computeTimeStepSize
|
||||
double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const;
|
||||
double computeTimeStepSize( const double dt, const int iterations, const SolutionTimeErrorInterface& /* timeError */ ) const;
|
||||
|
||||
protected:
|
||||
const int target_iterations_;
|
||||
@ -82,60 +82,18 @@ namespace Opm
|
||||
/// compute parallel scalarproducts.
|
||||
/// \param verbose if true get some output (default = false)
|
||||
PIDTimeStepControl( const double tol = 1e-3,
|
||||
const boost::any& pinfo = boost::any(),
|
||||
const bool verbose = false );
|
||||
|
||||
/// \brief \copydoc TimeStepControlInterface::initialize
|
||||
void initialize( const SimulatorState& state );
|
||||
|
||||
/// \brief \copydoc TimeStepControlInterface::computeTimeStepSize
|
||||
double computeTimeStepSize( const double dt, const int /* iterations */, const SimulatorState& state ) const;
|
||||
double computeTimeStepSize( const double dt, const int /* iterations */, const SolutionTimeErrorInterface& timeError ) const;
|
||||
|
||||
protected:
|
||||
template <class Iterator>
|
||||
double euclidianNormSquared( Iterator it, const Iterator end, int num_components = 1 ) const
|
||||
{
|
||||
static_cast<void>(num_components); // Suppress warning in the serial case.
|
||||
#if HAVE_MPI
|
||||
if ( parallel_information_.type() == typeid(ParallelISTLInformation) )
|
||||
{
|
||||
const ParallelISTLInformation& info =
|
||||
boost::any_cast<const ParallelISTLInformation&>(parallel_information_);
|
||||
int size_per_component = (end - it) / num_components;
|
||||
assert((end - it) == num_components * size_per_component);
|
||||
double component_product = 0.0;
|
||||
for( int i = 0; i < num_components; ++i )
|
||||
{
|
||||
auto component_container =
|
||||
boost::make_iterator_range(it + i * size_per_component,
|
||||
it + (i + 1) * size_per_component);
|
||||
info.computeReduction(component_container,
|
||||
Opm::Reduction::makeInnerProductFunctor<double>(),
|
||||
component_product);
|
||||
}
|
||||
return component_product;
|
||||
}
|
||||
else
|
||||
#endif
|
||||
{
|
||||
double product = 0.0 ;
|
||||
for( ; it != end; ++it ) {
|
||||
product += ( *it * *it );
|
||||
}
|
||||
return product;
|
||||
}
|
||||
}
|
||||
|
||||
protected:
|
||||
mutable std::vector<double> p0_;
|
||||
mutable std::vector<double> sat0_;
|
||||
|
||||
const double tol_;
|
||||
mutable std::vector< double > errors_;
|
||||
|
||||
const bool verbose_;
|
||||
private:
|
||||
const boost::any parallel_information_;
|
||||
// const boost::any parallel_information_;
|
||||
};
|
||||
|
||||
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
@ -158,11 +116,10 @@ namespace Opm
|
||||
/// \param verbose if true get some output (default = false)
|
||||
PIDAndIterationCountTimeStepControl( const int target_iterations = 20,
|
||||
const double tol = 1e-3,
|
||||
const boost::any& = boost::any(),
|
||||
const bool verbose = false);
|
||||
|
||||
/// \brief \copydoc TimeStepControlInterface::computeTimeStepSize
|
||||
double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const;
|
||||
double computeTimeStepSize( const double dt, const int iterations, const SolutionTimeErrorInterface& timeError ) const;
|
||||
|
||||
protected:
|
||||
const int target_iterations_;
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
Copyright 2014 IRIS AS
|
||||
Copyright 2014 IRIS AS
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
@ -19,31 +19,45 @@
|
||||
#ifndef OPM_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED
|
||||
#define OPM_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/simulator/SimulatorState.hpp>
|
||||
#include <opm/core/simulator/SimulatorState.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
///////////////////////////////////////////////////////////////////
|
||||
///
|
||||
/// TimeStepControlInterface
|
||||
///
|
||||
/// TimeStepControlInterface
|
||||
///
|
||||
///////////////////////////////////////////////////////////////////
|
||||
class TimeStepControlInterface
|
||||
class SolutionTimeErrorInterface
|
||||
{
|
||||
protected:
|
||||
protected:
|
||||
SolutionTimeErrorInterface() {}
|
||||
public:
|
||||
/// \return || u^n+1 - u^n || / || u^n+1 ||
|
||||
virtual double timeError() const = 0;
|
||||
|
||||
/// virtual destructor (empty)
|
||||
virtual ~SolutionTimeErrorInterface () {}
|
||||
};
|
||||
|
||||
///////////////////////////////////////////////////////////////////
|
||||
///
|
||||
/// TimeStepControlInterface
|
||||
///
|
||||
///////////////////////////////////////////////////////////////////
|
||||
class TimeStepControlInterface
|
||||
{
|
||||
protected:
|
||||
TimeStepControlInterface() {}
|
||||
public:
|
||||
/// \param state simulation state before computing update in the solver (default is empty)
|
||||
virtual void initialize( const SimulatorState& /*state*/ ) {}
|
||||
|
||||
/// 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 state new solution state
|
||||
/// \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 SimulatorState& ) const = 0;
|
||||
virtual double computeTimeStepSize( const double dt, const int iterations, const SolutionTimeErrorInterface& timeError ) const = 0;
|
||||
|
||||
/// virtual destructor (empty)
|
||||
virtual ~TimeStepControlInterface () {}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user