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

BlackOilModelBase: revision of adaptive time stepping.

opm/autodiff/BlackoilModelBase.hpp

@@ -178,6 +178,10 @@ namespace Opm {
         /// and afterwards the norm of the residual of the well flux and the well equation.
         std::vector<double> computeResidualNorms() const;
 
+        /// \brief compute the relative change between to simulation states
+        //  \return || u^n+1 - u^n || / || u^n+1 ||
+        double relativeChange( const SimulatorState& previous, const SimulatorState& current ) const;
+
         /// The size (number of unknowns) of the nonlinear system of equations.
         int sizeNonLinear() const;
 

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -1800,6 +1800,47 @@ namespace detail {
             }
         }
 
+        /// \brief Compute the Euclidian norm of a vector
+        /// \param it              begin iterator for the given vector
+        /// \param end             end iterator for the given vector
+        /// \param num_components  number of components (i.e. phases) in the vector
+        /// \param pinfo           In a parallel this holds the information about the data distribution.
+        template <class Iterator>
+        inline
+        double euclidianNormSquared( Iterator it, const Iterator end, int num_components, const boost::any& pinfo = boost::any() )
+        {
+            static_cast<void>(num_components); // Suppress warning in the serial case.
+            static_cast<void>(pinfo); // Suppress warning in non-MPI case.
+#if HAVE_MPI
+            if ( pinfo.type() == typeid(ParallelISTLInformation) )
+            {
+                const ParallelISTLInformation& info =
+                    boost::any_cast<const ParallelISTLInformation&>(pinfo);
+                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;
+            }
+        }
+
         /// \brief Compute the L-infinity norm of a vector representing a well equation.
         /// \param a The container to compute the infinity norm on.
         /// \param info In a parallel this holds the information about the data distribution.
@@ -2343,8 +2384,46 @@ namespace detail {
     }
 
 
+    template <class Grid, class Implementation>
+    double
+    BlackoilModelBase<Grid, Implementation>::
+    relativeChange(const SimulatorState& previous,
+                   const SimulatorState& current ) const
+    {
+        std::vector< double > p0  ( previous.pressure() );
+        std::vector< double > sat0( previous.saturation() );
 
+        const std::size_t pSize = p0.size();
+        const std::size_t satSize = sat0.size();
 
+        // compute u^n - u^n+1
+        for( std::size_t i=0; i<pSize; ++i ) {
+            p0[ i ] -= current.pressure()[ i ];
+        }
+
+        for( std::size_t i=0; i<satSize; ++i ) {
+            sat0[ i ] -= current.saturation()[ i ];
+        }
+
+        // compute || u^n - u^n+1 ||
+        const double stateOld  = detail::euclidianNormSquared( p0.begin(),   p0.end(), 1, linsolver_.parallelInformation() ) +
+                                 detail::euclidianNormSquared( sat0.begin(), sat0.end(),
+                                                               current.numPhases(),
+                                                               linsolver_.parallelInformation() );
+
+        // compute || u^n+1 ||
+        const double stateNew  = detail::euclidianNormSquared( current.pressure().begin(),   current.pressure().end(), 1, linsolver_.parallelInformation() ) +
+                                 detail::euclidianNormSquared( current.saturation().begin(), current.saturation().end(),
+                                                               current.numPhases(),
+                                                               linsolver_.parallelInformation() );
+
+        if( stateNew > 0.0 ) {
+            return stateOld / stateNew ;
+        }
+        else {
+            return 0.0;
+        }
+    }
 
     template <class Grid, class Implementation>
     double

opm/autodiff/NewtonSolver.hpp

@@ -96,6 +96,9 @@ namespace Opm {
         /// Number of linear solver iterations used in the last call to step().
         unsigned int linearIterationsLastStep() const;
 
+        /// return reference to physical model
+        const PhysicalModel& model() const;
+
     private:
         // ---------  Data members  ---------
         SolverParameters param_;

opm/autodiff/NewtonSolver_impl.hpp

@@ -49,6 +49,12 @@ namespace Opm
         return linearIterations_;
     }
 
+    template <class PhysicalModel>
+    const PhysicalModel& NewtonSolver<PhysicalModel>::model() const
+    {
+        assert( model_ );
+        return *model_;
+    }
 
     template <class PhysicalModel>
     int

opm/autodiff/SimulatorBase_impl.hpp

@@ -92,7 +92,7 @@ namespace Opm
         std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
         if( param_.getDefault("timestep.adaptive", true ) )
         {
-            adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, solver_.parallelInformation(), terminal_output_ ) );
+            adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
         }
 
         // init output writer