first attempt to implement a time step control.

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -25,6 +25,7 @@
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
+#include <opm/autodiff/TimeStepControl.hpp>
 
 struct UnstructuredGrid;
 struct Wells;
@@ -114,7 +115,8 @@ namespace Opm {
         /// \param[in] dt        time step size
         /// \param[in] state     reservoir state
         /// \param[in] wstate    well state
-        void
+        /// \return suggested time step for next step call
+        double
         step(const double   dt    ,
              BlackoilState& state ,
              WellStateFullyImplicitBlackoil&     wstate);
@@ -189,6 +191,8 @@ namespace Opm {
 
         std::vector<int>         primalVariable_;
 
+        IterationCountTimeStepControl timeStepControl_;
+
         // Private methods.
         SolutionState
         constantState(const BlackoilState& x,

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -235,6 +235,7 @@ namespace {
         , residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
                         ADB::null(),
                         ADB::null() } )
+        , timeStepControl_()
     {
     }
 
@@ -262,7 +263,7 @@ namespace {
 
 
     template<class T>
-    void
+    double
     FullyImplicitBlackoilSolver<T>::
     step(const double   dt,
          BlackoilState& x ,
@@ -304,10 +305,14 @@ namespace {
         bool isOscillate = false;
         bool isStagnate = false;
         const enum RelaxType relaxtype = relaxType();
+        int linearIterations = 0;
 
         while ((!converged) && (it < maxIter())) {
             V dx = solveJacobianSystem();
 
+            // store number of linear iterations used
+            linearIterations += linsolver_.iterations();
+
             detectNewtonOscillations(residual_history, it, relaxRelTol(), isOscillate, isStagnate);
 
             if (isOscillate) {
@@ -337,6 +342,9 @@ namespace {
             std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
             // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
         }
+
+        std::cout << "Linear iterations: " << linearIterations << std::endl;
+        return timeStepControl_.computeTimeStepSize( dt, linearIterations );
     }
 
 

opm/autodiff/NewtonIterationBlackoilCPR.cpp

@@ -108,6 +108,7 @@ namespace Opm
 
     /// Construct a system solver.
     NewtonIterationBlackoilCPR::NewtonIterationBlackoilCPR(const parameter::ParameterGroup& param)
+        : iterations_( 0 )
     {
         use_amg_ = param.getDefault("cpr_use_amg", false);
         use_bicgstab_ = param.getDefault("cpr_use_bicgstab", true);
@@ -199,6 +200,9 @@ namespace Opm
         Dune::InverseOperatorResult result;
         linsolve.apply(x, istlb, result);
 
+        // store number of iterations
+        iterations_ = result.iterations;
+
         // Check for failure of linear solver.
         if (!result.converged) {
             OPM_THROW(std::runtime_error, "Convergence failure for linear solver.");

opm/autodiff/NewtonIterationBlackoilCPR.hpp

@@ -54,7 +54,10 @@ namespace Opm
         /// \return               the solution x
         virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const;
 
+        /// \copydoc NewtonIterationBlackoilInterface::iterations
+        virtual int iterations () const { return iterations_; }
     private:
+        mutable int iterations_;
         bool use_amg_;
         bool use_bicgstab_;
     };

opm/autodiff/NewtonIterationBlackoilInterface.hpp

@@ -1,5 +1,6 @@
 /*
   Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 IRIS AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -38,6 +39,9 @@ namespace Opm
         /// \param[in] residual   residual object containing A and b.
         /// \return               the solution x
         virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const = 0;
+
+        /// \return number of iterations used during last call of computeNewtonIncrement
+        virtual int iterations () const = 0;
     };
 
 } // namespace Opm

opm/autodiff/NewtonIterationBlackoilSimple.cpp

@@ -30,6 +30,7 @@ namespace Opm
     /// Construct a system solver.
     /// \param[in] linsolver   linear solver to use
     NewtonIterationBlackoilSimple::NewtonIterationBlackoilSimple(const parameter::ParameterGroup& param)
+        : iterations_( 0 )
     {
         linsolver_.reset(new LinearSolverFactory(param));
     }
@@ -58,6 +59,10 @@ namespace Opm
             = linsolver_->solve(matr.rows(), matr.nonZeros(),
                                 matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
                                 total_residual.value().data(), dx.data());
+
+        // store iterations 
+        iterations_ = rep.iterations;
+
         if (!rep.converged) {
             OPM_THROW(std::runtime_error,
                       "FullyImplicitBlackoilSolver::solveJacobianSystem(): "

opm/autodiff/NewtonIterationBlackoilSimple.hpp

@@ -1,5 +1,6 @@
 /*
   Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 IRIS AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -48,8 +49,12 @@ namespace Opm
         /// \return               the solution x
         virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const;
 
+        /// \copydoc NewtonIterationBlackoilInterface::iterations
+        virtual int iterations () const { return iterations_; }
+
     private:
         std::unique_ptr<LinearSolverInterface> linsolver_;
+        mutable int iterations_;
     };
 
 } // namespace Opm

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -1,5 +1,6 @@
 /*
   Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 IRIS AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -291,6 +292,8 @@ namespace Opm
         std::string tstep_filename = output_dir_ + "/step_timing.txt";
         std::ofstream tstep_os(tstep_filename.c_str());
 
+        double lastSubStep = timer.currentStepLength();
+
         typename FullyImplicitBlackoilSolver<T>::SolverParameter solverParam( param_ );
 
         // Main simulation loop.
@@ -340,13 +343,40 @@ namespace Opm
             // Compute reservoir volumes for RESV controls.
             computeRESV(timer.currentStepNum(), wells, state, well_state);
 
-            // Run a single step of the solver.
+            // Run a multiple steps of the solver depending on the time step control.
             solver_timer.start();
+
             FullyImplicitBlackoilSolver<T> solver(solverParam, grid_, props_, geo_, rock_comp_props_, *wells, solver_, has_disgas_, has_vapoil_);
             if (!threshold_pressures_by_face_.empty()) {
                 solver.setThresholdPressures(threshold_pressures_by_face_);
             }
-            solver.step(timer.currentStepLength(), state, well_state);
+
+            const bool subStepping = false;
+            if( subStepping )
+            {
+
+                // create sub step simulator timer with previously used sub step size
+                SubStepSimulatorTimer subStepper( timer, lastSubStep );
+
+                while( ! subStepper.done() )
+                {
+                    const double dt_new = solver.step(subStepper.currentStepLength(), state, well_state);
+                    subStepper.next( dt_new );
+                }
+
+                subStepper.report( std::cout );
+
+                // store last small time step for next reportStep
+                lastSubStep = subStepper.currentStepLength();
+
+                std::cout << "Last suggested step size = " << lastSubStep << std::endl;
+                if( lastSubStep != lastSubStep )
+                    lastSubStep = timer.currentStepLength();
+            }
+            else
+                solver.step(timer.currentStepLength(), state, well_state);
+
+            // take time that was used to solve system for this reportStep
             solver_timer.stop();
 
             // Report timing.

opm/autodiff/TimeStepControl.hpp

@@ -0,0 +1,103 @@
+/*
+  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_TIMESTEPCONTROL_HEADER_INCLUDED
+#define OPM_TIMESTEPCONTROL_HEADER_INCLUDED
+
+namespace Opm
+{
+
+    class TimeStepControlInterface 
+    {
+    protected:    
+        TimeStepControlInterface() {}
+    public:
+        virtual double computeTimeStepSize( const double dt, const int iterations ) const = 0;
+        virtual ~TimeStepControlInterface () {}
+    };
+
+    class IterationCountTimeStepControl : public TimeStepControlInterface
+    {
+    protected:    
+        mutable double prevDt_;
+        mutable int prevIterations_;
+        const int targetIterationCount_;
+        const double adjustmentFactor_;
+
+        const int upperTargetIterationCount_;
+        const int lowerTargetIterationCount_;
+
+    public:
+        IterationCountTimeStepControl() 
+            : prevDt_( 0.0 ), prevIterations_( 0 ), 
+              targetIterationCount_( 100 ), adjustmentFactor_( 1.25 ),
+              upperTargetIterationCount_( 200 ), lowerTargetIterationCount_( 30 )
+        {}
+
+        double computeTimeStepSize( const double dt, const int iterations ) const
+        {
+            // make sure dt is somewhat reliable
+            assert( dt > 0 && dt == dt );
+            double newDt = dt;
+            double derivation = double(std::abs( iterations - targetIterationCount_ )) / double(targetIterationCount_);
+
+            if( derivation > 0.1 )
+            {
+                if( iterations < targetIterationCount_ ) 
+                {
+                    newDt = dt * adjustmentFactor_;
+                }
+                else 
+                    newDt = dt / adjustmentFactor_;
+            }
+                           
+            /*
+            if( prevDt_ > 0 && std::abs( dt - prevDt_ ) > 1e-12 ) {
+                const double dFdt  = double(iterations - prevIterations_) / ( dt - prevDt_ );
+                if( std::abs( dFdt ) > 1e-12 )
+                    newDt = dt + (targetIterationCount_ - iterations) / dFdt;
+                else
+                    // if iterations was the same or dts were the same, do some magic
+                    newDt = dt * double( targetIterationCount_ ) / double(targetIterationCount_ - iterations);
+            }
+
+            if( newDt < 0 || ! (prevDt_ > 0) || ( iterations == prevIterations_) ) 
+            {
+                if( iterations > upperTargetIterationCount_ )
+                    newDt = dt / adjustmentFactor_;
+                else if( iterations < lowerTargetIterationCount_ )
+                    newDt = dt * adjustmentFactor_;
+                else
+                    newDt = dt;
+            }
+            */
+
+            assert( newDt == newDt );
+
+            //std::cout << "dt = " << dt << "  " << prevDt_ << std::endl;
+
+            prevDt_ = dt;
+            prevIterations_ = iterations;
+
+            return newDt;
+        }
+    };
+
+} // end namespace OPM
+
+#endif