Merge pull request #216 from dr-robertk/timestepcontrol

Timestepcontrol

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -114,7 +114,8 @@ namespace Opm {
         /// \param[in] dt        time step size
         /// \param[in] state     reservoir state
         /// \param[in] wstate    well state
-        void
+        /// \return              number of linear iterations used
+        int
         step(const double   dt    ,
              BlackoilState& state ,
              WellStateFullyImplicitBlackoil&     wstate);

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -262,7 +262,7 @@ namespace {
 
 
     template<class T>
-    void
+    int
     FullyImplicitBlackoilSolver<T>::
     step(const double   dt,
          BlackoilState& x ,
@@ -278,7 +278,6 @@ namespace {
             computeWellConnectionPressures(state, xw);
         }
 
-
         std::vector<std::vector<double>> residual_history;
 
         assemble(pvdt, x, xw);
@@ -304,10 +303,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) {
@@ -328,15 +331,19 @@ namespace {
 
             converged = getConvergence(dt);
 
-            it += 1;
+            // increase iteration counter
+            ++it;
             std::cout << std::setw(9) << it << std::setprecision(9)
                       << std::setw(18) << r << std::endl;
         }
 
         if (!converged) {
-            std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
+            // the runtime_error is caught by the AdaptiveTimeStepping
-            // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
+            OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
+            return -1;
         }
+
+        return linearIterations;
     }
 
 

opm/autodiff/NewtonIterationBlackoilCPR.cpp

@@ -110,6 +110,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);
@@ -193,7 +194,7 @@ namespace Opm
         // Construct linear solver.
         const double tolerance = 1e-3;
         const int maxit = 5000;
-        const int verbosity = 1;
+        const int verbosity = 0;
         const int restart = 40;
         Dune::RestartedGMResSolver<Vector> linsolve(opA, sp, precond, tolerance, restart, maxit, verbosity);
 
@@ -201,6 +202,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 linear 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).
 
@@ -17,7 +18,7 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include  <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
+#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
 #include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
@@ -36,6 +37,7 @@
 #include <opm/core/io/eclipse/EclipseWriter.hpp>
 #include <opm/core/simulator/SimulatorReport.hpp>
 #include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
 #include <opm/core/utility/StopWatch.hpp>
 #include <opm/core/io/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
@@ -44,6 +46,7 @@
 #include <opm/core/props/rock/RockCompressibility.hpp>
 
 #include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/simulator/AdaptiveTimeStepping.hpp>
 #include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
@@ -51,6 +54,7 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
 
+
 #include <boost/filesystem.hpp>
 #include <boost/lexical_cast.hpp>
 
@@ -293,6 +297,13 @@ namespace Opm
 
         typename FullyImplicitBlackoilSolver<T>::SolverParameter solverParam( param_ );
 
+        // adaptive time stepping
+        std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
+        if( param_.getDefault("timestep.adaptive", bool(false) ) )
+        {
+            adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_ ) );
+        }
+
         // Main simulation loop.
         while (!timer.done()) {
             // Report timestep.
@@ -340,13 +351,29 @@ 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);
+
+            // If sub stepping is enabled allow the solver to sub cycle
+            // in case the report steps are to large for the solver to converge
+            //
+            // \Note: The report steps are met in any case
+            // \Note: The sub stepping will require a copy of the state variables
+            if( adaptiveTimeStepping ) {
+                adaptiveTimeStepping->step( solver, state, well_state,
+                        timer.simulationTimeElapsed(), timer.currentStepLength() );
+            }
+            else {
+                // solve for complete report step
+                solver.step(timer.currentStepLength(), state, well_state);
+            }
+
+            // take time that was used to solve system for this reportStep
             solver_timer.stop();
 
             // Report timing.