diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp index 732e4006e..3725ea7c5 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping.hpp @@ -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 diff --git a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp index 0cfe655b3..354fe9108 100644 --- a/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp +++ b/opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp @@ -31,11 +31,35 @@ namespace Opm { + namespace detail + { + template + 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) ); diff --git a/opm/simulators/timestepping/TimeStepControl.cpp b/opm/simulators/timestepping/TimeStepControl.cpp index cd0fb73f2..7ac85d2cd 100644 --- a/opm/simulators/timestepping/TimeStepControl.cpp +++ b/opm/simulators/timestepping/TimeStepControl.cpp @@ -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 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_ ) diff --git a/opm/simulators/timestepping/TimeStepControl.hpp b/opm/simulators/timestepping/TimeStepControl.hpp index eee92d5f2..b3ebc3b5c 100644 --- a/opm/simulators/timestepping/TimeStepControl.hpp +++ b/opm/simulators/timestepping/TimeStepControl.hpp @@ -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 - double euclidianNormSquared( Iterator it, const Iterator end, int num_components = 1 ) const - { - static_cast(num_components); // Suppress warning in the serial case. -#if HAVE_MPI - if ( parallel_information_.type() == typeid(ParallelISTLInformation) ) - { - const ParallelISTLInformation& info = - boost::any_cast(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(), - component_product); - } - return component_product; - } - else -#endif - { - double product = 0.0 ; - for( ; it != end; ++it ) { - product += ( *it * *it ); - } - return product; - } - } - - protected: - mutable std::vector p0_; - mutable std::vector 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_; diff --git a/opm/simulators/timestepping/TimeStepControlInterface.hpp b/opm/simulators/timestepping/TimeStepControlInterface.hpp index a32e1285d..0da1637a0 100644 --- a/opm/simulators/timestepping/TimeStepControlInterface.hpp +++ b/opm/simulators/timestepping/TimeStepControlInterface.hpp @@ -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 +#include 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 () {}