Merge pull request #328 from dr-robertk/PR/tolerance-wells

Output Newton and Lienar solver iterations

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -72,6 +72,7 @@ namespace Opm {
             double                          max_residual_allowed_;
             double                          tolerance_mb_;
             double                          tolerance_cnv_;
+            double                          tolerance_wells_;
             int                             max_iter_;
 
             SolverParameter( const parameter::ParameterGroup& param );
@@ -100,7 +101,8 @@ namespace Opm {
                                     const Wells*                    wells,
                                     const NewtonIterationBlackoilInterface& linsolver,
                                     const bool has_disgas,
-                                    const bool has_vapoil );
+                                    const bool has_vapoil,
+                                    const bool terminal_output);
 
         /// \brief Set threshold pressures that prevent or reduce flow.
         /// This prevents flow across faces if the potential
@@ -127,6 +129,9 @@ namespace Opm {
              BlackoilState& state ,
              WellStateFullyImplicitBlackoil&     wstate);
 
+        unsigned int newtonIterations () const { return newtonIterations_; }
+        unsigned int linearIterations () const { return linearIterations_; }
+
     private:
         // Types and enums
         typedef AutoDiffBlock<double> ADB;
@@ -200,6 +205,8 @@ namespace Opm {
 
         /// \brief Whether we print something to std::cout
         bool terminal_output_;
+        unsigned int newtonIterations_;
+        unsigned int linearIterations_;
 
         std::vector<int>         primalVariable_;
 

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -147,6 +147,7 @@ namespace detail {
         max_residual_allowed_  = std::numeric_limits< double >::max();
         tolerance_mb_    = 1.0e-7;
         tolerance_cnv_   = 1.0e-3;
+        tolerance_wells_ = 1./Opm::unit::day;
     }
 
     template<class T>
@@ -174,6 +175,7 @@ namespace detail {
 
         tolerance_mb_    = param.getDefault("tolerance_mb", tolerance_mb_);
         tolerance_cnv_   = param.getDefault("tolerance_cnv", tolerance_cnv_);
+        tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ );
 
         std::string relaxation_type = param.getDefault("relax_type", std::string("dampen"));
         if (relaxation_type == "dampen") {
@@ -196,7 +198,8 @@ namespace detail {
                                 const Wells*                    wells,
                                 const NewtonIterationBlackoilInterface&    linsolver,
                                 const bool has_disgas,
-                                const bool has_vapoil)
+                                const bool has_vapoil,
+                                const bool terminal_output)
         : grid_  (grid)
         , fluid_ (fluid)
         , geo_   (geo)
@@ -217,16 +220,20 @@ namespace detail {
         , residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
                         ADB::null(),
                         ADB::null() } )
-        , terminal_output_ (true)
+        , terminal_output_ (terminal_output)
+        , newtonIterations_( 0 )
+        , linearIterations_( 0 )
     {
 #if HAVE_MPI
+        if ( terminal_output_ ) {
             if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
             {
                 const ParallelISTLInformation& info =
                     boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
-            // Only rank 0 does print to std::cout
+                // Only rank 0 does print to std::cout if terminal_output is enabled
                 terminal_output_ = (info.communicator().rank()==0);
             }
+        }
 #endif
     }
 
@@ -330,6 +337,9 @@ namespace detail {
             return -1;
         }
 
+        linearIterations_ += linearIterations;
+        newtonIterations_ += it;
+
         return linearIterations;
     }
 
@@ -1972,6 +1982,7 @@ namespace detail {
     {
         const double tol_mb    = param_.tolerance_mb_;
         const double tol_cnv   = param_.tolerance_cnv_;
+        const double tol_wells = param_.tolerance_wells_;
 
         const int nc = Opm::AutoDiffGrid::numCells(grid_);
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
@@ -2016,7 +2027,7 @@ namespace detail {
 
         const double residualWellFlux = detail::infinityNorm(residual_.well_flux_eq);
         const double residualWell     = detail::infinityNorm(residual_.well_eq);
-        const bool   converged_Well   = (residualWellFlux < 1./Opm::unit::day) && (residualWell < Opm::unit::barsa);
+        const bool   converged_Well   = (residualWellFlux < tol_wells) && (residualWell < Opm::unit::barsa);
         const bool   converged        = converged_MB && converged_CNV && converged_Well;
 
         // if one of the residuals is NaN, throw exception, so that the solver can be restarted

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -185,7 +185,7 @@ namespace Opm
           solver_(linsolver),
           has_disgas_(has_disgas),
           has_vapoil_(has_vapoil),
-          terminal_output_(true),
+          terminal_output_(param.getDefault("output_terminal", true)),
           eclipse_state_(eclipse_state),
           output_writer_(output_writer),
           rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0)),
@@ -198,6 +198,7 @@ namespace Opm
             allcells_[cell] = cell;
         }
 #if HAVE_MPI
+        if ( terminal_output_ ) {
             if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
             {
                 const ParallelISTLInformation& info =
@@ -205,6 +206,7 @@ namespace Opm
                 // Only rank 0 does print to std::cout
                 terminal_output_= (info.communicator().rank()==0);
             }
+        }
 #endif
     }
 
@@ -246,6 +248,9 @@ namespace Opm
             output_writer_.restore( timer, state, prev_well_state, restorefilename, desiredRestoreStep );
         }
 
+        unsigned int totalNewtonIterations = 0;
+        unsigned int totalLinearIterations = 0;
+
         // Main simulation loop.
         while (!timer.done()) {
             // Report timestep.
@@ -283,7 +288,7 @@ namespace Opm
             // 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_);
+            FullyImplicitBlackoilSolver<T> solver(solverParam, grid_, props_, geo_, rock_comp_props_, wells, solver_, has_disgas_, has_vapoil_, terminal_output_);
             if (!threshold_pressures_by_face_.empty()) {
                 solver.setThresholdPressures(threshold_pressures_by_face_);
             }
@@ -304,6 +309,10 @@ namespace Opm
             // take time that was used to solve system for this reportStep
             solver_timer.stop();
 
+            // accumulate the number of Newton and Linear Iterations
+            totalNewtonIterations += solver.newtonIterations();
+            totalLinearIterations += solver.linearIterations();
+
             // Report timing.
             const double st = solver_timer.secsSinceStart();
 
@@ -336,6 +345,8 @@ namespace Opm
         report.pressure_time = stime;
         report.transport_time = 0.0;
         report.total_time = total_timer.secsSinceStart();
+        report.total_newton_iterations = totalNewtonIterations;
+        report.total_linear_iterations = totalLinearIterations;
         return report;
     }