Merge pull request #921 from dr-robertk/PR/adaptive-timestepping-revision

AdaptiveTimeStepping: allow for more flexible error computation in the future.

opm/core/simulator/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

opm/core/simulator/AdaptiveTimeStepping_impl.hpp

@@ -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 relativeChange() const
+            {
+                return solver_.model().relativeChange( 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 >
+                    relativeChange( solver, last_state, state );
+
                 // compute new time step estimate
                 double dtEstimate =
-                    timeStepControl_->computeTimeStepSize( dt, linearIterations, state );
+                    timeStepControl_->computeTimeStepSize( dt, linearIterations, relativeChange );
 
                 // limit the growth of the timestep size by the growth factor
                 dtEstimate = std::min( dtEstimate, double(max_growth_ * dt) );

opm/core/simulator/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 RelativeChangeInterface& /* relativeChange */ ) 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_()
+        : tol_( tol )
-        , sat0_()
-        , 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 RelativeChangeInterface& relChange ) 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;
+        const double error = relChange.relativeChange();
-        errors_[ 2 ] =  error ;
+        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 RelativeChangeInterface& relChange ) const
     {
-        double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, state );
+        double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, relChange );
 
         // further reduce step size if to many iterations were used
         if( iterations > target_iterations_ )

opm/core/simulator/TimeStepControl.hpp

@@ -26,7 +26,6 @@
 #include <boost/any.hpp>
 #include <boost/range/iterator_range.hpp>
 #include <opm/core/simulator/TimeStepControlInterface.hpp>
-#include <opm/core/linalg/ParallelIstlInformation.hpp>
 
 namespace Opm
 {
@@ -49,7 +48,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 RelativeChangeInterface& /* relativeChange */ ) const;
 
     protected:
         const int     target_iterations_;
@@ -78,64 +77,18 @@ namespace Opm
         /// \brief constructor
         /// \param tol      tolerance for the relative changes of the numerical solution to be accepted
         ///                 in one time step (default is 1e-3)
-        /// \paramm pinfo   The information about the parallel information. Needed to
-        ///                 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 RelativeChangeInterface& relativeChange ) 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_;
     };
 
     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -152,17 +105,13 @@ namespace Opm
         /// \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 maxgrowth  max growth factor for new time step in relation of old time step (default = 3.0)
-        /// \paramm pinfo     The information about the parallel information. Needed to
-        ///                   compute parallel scalarproducts.
         /// \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 RelativeChangeInterface& relativeChange ) const;
 
     protected:
         const int     target_iterations_;

opm/core/simulator/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 <opm/core/simulator/SimulatorState.hpp> 
+#include <opm/core/simulator/SimulatorState.hpp>
 
 namespace Opm
 {
 
     ///////////////////////////////////////////////////////////////////
     ///
-    ///  TimeStepControlInterface 
+    ///  RelativeChangeInterface
-    /// 
+    ///
     ///////////////////////////////////////////////////////////////////
-    class TimeStepControlInterface 
+    class RelativeChangeInterface
     {
-    protected:    
+    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:
-        /// \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 iterations  number of iterations used (linear/nonlinear)
-        /// \param state       new solution state
+        /// \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 RelativeChangeInterface& relativeChange ) const = 0;
 
         /// virtual destructor (empty)
         virtual ~TimeStepControlInterface () {}