From d03f9411b673e1c84d7a4fc5e7bd0457dbd5c53a Mon Sep 17 00:00:00 2001 From: Robert K Date: Wed, 1 Oct 2014 13:50:08 +0200 Subject: [PATCH 01/19] first attempt to implement a time step control. --- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 6 +- .../FullyImplicitBlackoilSolver_impl.hpp | 10 +- opm/autodiff/NewtonIterationBlackoilCPR.cpp | 4 + opm/autodiff/NewtonIterationBlackoilCPR.hpp | 3 + .../NewtonIterationBlackoilInterface.hpp | 4 + .../NewtonIterationBlackoilSimple.cpp | 5 + .../NewtonIterationBlackoilSimple.hpp | 5 + .../SimulatorFullyImplicitBlackoil_impl.hpp | 34 +++++- opm/autodiff/TimeStepControl.hpp | 103 ++++++++++++++++++ 9 files changed, 170 insertions(+), 4 deletions(-) create mode 100644 opm/autodiff/TimeStepControl.hpp diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 0e3ec3cef..8858f91ea 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -25,6 +25,7 @@ #include #include #include +#include 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 primalVariable_; + IterationCountTimeStepControl timeStepControl_; + // Private methods. SolutionState constantState(const BlackoilState& x, diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 1f7492cab..1b2a621b2 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -235,6 +235,7 @@ namespace { , residual_ ( { std::vector(fluid.numPhases(), ADB::null()), ADB::null(), ADB::null() } ) + , timeStepControl_() { } @@ -262,7 +263,7 @@ namespace { template - void + double FullyImplicitBlackoilSolver:: 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 ); } diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index f7d3b2281..218dfa306 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/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."); diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.hpp b/opm/autodiff/NewtonIterationBlackoilCPR.hpp index fb2bd2b73..c06310393 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.hpp +++ b/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_; }; diff --git a/opm/autodiff/NewtonIterationBlackoilInterface.hpp b/opm/autodiff/NewtonIterationBlackoilInterface.hpp index 277b36d9b..c105cc783 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterface.hpp +++ b/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 diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.cpp b/opm/autodiff/NewtonIterationBlackoilSimple.cpp index 4ced9bcb5..c8d2472d9 100644 --- a/opm/autodiff/NewtonIterationBlackoilSimple.cpp +++ b/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(): " diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.hpp b/opm/autodiff/NewtonIterationBlackoilSimple.hpp index 7fceaa28a..6bdaf7b07 100644 --- a/opm/autodiff/NewtonIterationBlackoilSimple.hpp +++ b/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 linsolver_; + mutable int iterations_; }; } // namespace Opm diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 4cad4ed41..99a07053d 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/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::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 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. diff --git a/opm/autodiff/TimeStepControl.hpp b/opm/autodiff/TimeStepControl.hpp new file mode 100644 index 000000000..912e42e33 --- /dev/null +++ b/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 . +*/ +#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 From f535761a17ce10ed751ae73cb82b93d706ec6d9b Mon Sep 17 00:00:00 2001 From: Robert K Date: Wed, 1 Oct 2014 15:45:11 +0200 Subject: [PATCH 02/19] only warn when non-diagonal block is found. --- opm/autodiff/NewtonIterationBlackoilCPR.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index 218dfa306..aeb026cb0 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -192,7 +192,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 linsolve(opA, sp, precond, tolerance, restart, maxit, verbosity); @@ -245,7 +245,8 @@ namespace Opm // std::cout << "++++++++++++++++++++++++++++++++++++++++++++\n" // << D // << "++++++++++++++++++++++++++++++++++++++++++++\n" << std::endl; - OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); + std::cerr << "WARNING (ignored): Cannot do Schur complement with respect to non-diagonal block." << std::endl; + //OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); } V diag = D.diagonal(); Eigen::DiagonalMatrix invD = (1.0 / diag).matrix().asDiagonal(); @@ -300,7 +301,8 @@ namespace Opm // Find inv(D). const M& D = equation.derivative()[n]; if (!isDiagonal(D)) { - OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); + std::cerr << "WARNING (ignored): Cannot do Schur complement with respect to non-diagonal block." << std::endl; + //OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); } V diag = D.diagonal(); Eigen::DiagonalMatrix invD = (1.0 / diag).matrix().asDiagonal(); From 2602ae7baf05088e06a03e32ee38fa12c3b1f696 Mon Sep 17 00:00:00 2001 From: Robert K Date: Wed, 1 Oct 2014 15:45:27 +0200 Subject: [PATCH 03/19] enable substepping. --- opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 99a07053d..74a707473 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -351,7 +351,7 @@ namespace Opm solver.setThresholdPressures(threshold_pressures_by_face_); } - const bool subStepping = false; + const bool subStepping = true; if( subStepping ) { From d3bc836536e776456e8415a5d58858bfdf05c7ad Mon Sep 17 00:00:00 2001 From: Robert K Date: Wed, 1 Oct 2014 16:36:38 +0200 Subject: [PATCH 04/19] make sub stepping a parameter. --- opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 74a707473..f76b83472 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -296,6 +296,8 @@ namespace Opm typename FullyImplicitBlackoilSolver::SolverParameter solverParam( param_ ); + const bool subStepping = param_.getDefault("substepping", bool(false) ); + // Main simulation loop. while (!timer.done()) { // Report timestep. @@ -351,7 +353,6 @@ namespace Opm solver.setThresholdPressures(threshold_pressures_by_face_); } - const bool subStepping = true; if( subStepping ) { From fcf6cd5f90178c5b35561e75d861760d19fd9676 Mon Sep 17 00:00:00 2001 From: Robert K Date: Thu, 2 Oct 2014 14:04:59 +0200 Subject: [PATCH 05/19] implemented the PID controler, seems to work fine. More testing needed. --- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 3 +- .../FullyImplicitBlackoilSolver_impl.hpp | 4 +- .../SimulatorFullyImplicitBlackoil_impl.hpp | 6 +- opm/autodiff/TimeStepControl.hpp | 97 ++++++++++++++++++- 4 files changed, 104 insertions(+), 6 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 8858f91ea..9b00f53d6 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -191,7 +191,8 @@ namespace Opm { std::vector primalVariable_; - IterationCountTimeStepControl timeStepControl_; + //IterationCountTimeStepControl timeStepControl_; + PIDTimeStepControl timeStepControl_; // Private methods. SolutionState diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 1b2a621b2..5f6154cba 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -279,6 +279,8 @@ namespace { computeWellConnectionPressures(state, xw); } + // initialize time step control + timeStepControl_.initialize( x ); std::vector> residual_history; @@ -344,7 +346,7 @@ namespace { } std::cout << "Linear iterations: " << linearIterations << std::endl; - return timeStepControl_.computeTimeStepSize( dt, linearIterations ); + return timeStepControl_.computeTimeStepSize( dt, linearIterations, x ); } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index f76b83472..0d7d62fe3 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -355,7 +355,6 @@ namespace Opm if( subStepping ) { - // create sub step simulator timer with previously used sub step size SubStepSimulatorTimer subStepper( timer, lastSubStep ); @@ -368,7 +367,10 @@ namespace Opm subStepper.report( std::cout ); // store last small time step for next reportStep - lastSubStep = subStepper.currentStepLength(); + //lastSubStep = subStepper.maxStepLength(); + //lastSubStep = subStepper.averageStepLength(); + //lastSubStep = subStepper.suggestedMax(); + lastSubStep = subStepper.suggestedAverage(); std::cout << "Last suggested step size = " << lastSubStep << std::endl; if( lastSubStep != lastSubStep ) diff --git a/opm/autodiff/TimeStepControl.hpp b/opm/autodiff/TimeStepControl.hpp index 912e42e33..e18d2946b 100644 --- a/opm/autodiff/TimeStepControl.hpp +++ b/opm/autodiff/TimeStepControl.hpp @@ -27,7 +27,8 @@ namespace Opm protected: TimeStepControlInterface() {} public: - virtual double computeTimeStepSize( const double dt, const int iterations ) const = 0; + virtual void initialize( const SimulatorState& state ) {} + virtual double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& ) const = 0; virtual ~TimeStepControlInterface () {} }; @@ -49,7 +50,7 @@ namespace Opm upperTargetIterationCount_( 200 ), lowerTargetIterationCount_( 30 ) {} - double computeTimeStepSize( const double dt, const int iterations ) const + double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& ) const { // make sure dt is somewhat reliable assert( dt > 0 && dt == dt ); @@ -98,6 +99,98 @@ namespace Opm } }; + class PIDTimeStepControl : public TimeStepControlInterface + { + protected: + mutable std::vector p0_; + mutable std::vector sat0_; + + mutable double prevDt_; + mutable int prevIterations_; + const int targetIterationCount_; + const double adjustmentFactor_; + + const int upperTargetIterationCount_; + const int lowerTargetIterationCount_; + + const double tol_; + + mutable std::vector< double > errors_; + + public: + PIDTimeStepControl( const double tol = 8e-4 ) + : p0_(), sat0_(), prevDt_( 0.0 ), prevIterations_( 0 ), + targetIterationCount_( 100 ), adjustmentFactor_( 1.25 ), + upperTargetIterationCount_( 200 ), lowerTargetIterationCount_( 30 ), + tol_( tol ), + errors_( 3, tol_ ) + {} + + void initialize( const SimulatorState& state ) + { + // store current state + p0_ = state.pressure(); + sat0_ = state.saturation(); + } + + template + double inner_product( Iterator it, const Iterator end ) const + { + double product = 0.0 ; + for( ; it != end; ++it ) + product += ( *it * *it ); + return product; + } + + double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const + { + const size_t size = p0_.size(); + assert( state.pressure().size() == size ); + // compute u^n - u^n+1 + for( size_t i=0; i tol_ ) + { + // adjust dt by given tolerance + std::cout << "Computed step size (tol): " << (dt * tol_ / error ) << std::endl; + return (dt * tol_ / error ); + } + else + { + const double kP = 0.075 ; + const double kI = 0.175 ; + const double kD = 0.01 ; + double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) * + std::pow( tol_ / errors_[ 2 ], kI ) * + std::pow( (errors_[0]*errors_[0]/errors_[ 1 ]*errors_[ 2 ]), kD )); + std::cout << "Computed step size (pow): " << newDt << std::endl; + return newDt; + } + } + }; + } // end namespace OPM #endif From a723a01f7258d8e5bec611618dcf610dca88c19b Mon Sep 17 00:00:00 2001 From: Robert K Date: Fri, 3 Oct 2014 14:18:31 +0200 Subject: [PATCH 06/19] some revision, time step control is now completly in the Simulator run method. The solver simply returns a number of iterations. --- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 8 +- .../FullyImplicitBlackoilSolver_impl.hpp | 16 +- .../SimulatorFullyImplicitBlackoil_impl.hpp | 61 ++++- opm/autodiff/TimeStepControl.hpp | 245 +++++++++--------- 4 files changed, 180 insertions(+), 150 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 9b00f53d6..f0eb90111 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -25,7 +25,6 @@ #include #include #include -#include struct UnstructuredGrid; struct Wells; @@ -115,8 +114,8 @@ namespace Opm { /// \param[in] dt time step size /// \param[in] state reservoir state /// \param[in] wstate well state - /// \return suggested time step for next step call - double + /// \return number of iterations used + int step(const double dt , BlackoilState& state , WellStateFullyImplicitBlackoil& wstate); @@ -191,9 +190,6 @@ namespace Opm { std::vector primalVariable_; - //IterationCountTimeStepControl timeStepControl_; - PIDTimeStepControl timeStepControl_; - // Private methods. SolutionState constantState(const BlackoilState& x, diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 5f6154cba..2a1079072 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -235,7 +235,6 @@ namespace { , residual_ ( { std::vector(fluid.numPhases(), ADB::null()), ADB::null(), ADB::null() } ) - , timeStepControl_() { } @@ -263,7 +262,7 @@ namespace { template - double + int FullyImplicitBlackoilSolver:: step(const double dt, BlackoilState& x , @@ -279,9 +278,6 @@ namespace { computeWellConnectionPressures(state, xw); } - // initialize time step control - timeStepControl_.initialize( x ); - std::vector> residual_history; assemble(pvdt, x, xw); @@ -335,18 +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"; + std::cerr << "ERROR: Failed to compute converged solution in " << it << " iterations." << std::endl; // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); + return -1; } - std::cout << "Linear iterations: " << linearIterations << std::endl; - return timeStepControl_.computeTimeStepSize( dt, linearIterations, x ); + return linearIterations; } @@ -1704,6 +1701,7 @@ namespace { for (; quantityIt != endQuantityIt; ++quantityIt) { const double quantityResid = (*quantityIt).value().matrix().norm(); if (!std::isfinite(quantityResid)) { + //std::cout << quantityResid << " quantity" << std::endl; const int trouble_phase = quantityIt - residual_.material_balance_eq.begin(); OPM_THROW(Opm::NumericalProblem, "Encountered a non-finite residual in material balance equation " diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 0d7d62fe3..ccf186ab4 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -18,7 +18,7 @@ along with OPM. If not, see . */ -#include +#include #include #include #include @@ -28,6 +28,7 @@ #include #include #include +#include #include #include @@ -37,6 +38,7 @@ #include #include #include +#include #include #include #include @@ -298,6 +300,9 @@ namespace Opm const bool subStepping = param_.getDefault("substepping", bool(false) ); + std::unique_ptr< TimeStepControlInterface > + timeStepControl( new PIDAndIterationCountTimeStepControl( 1e-3, 50 ) ); + // Main simulation loop. while (!timer.done()) { // Report timestep. @@ -353,31 +358,67 @@ namespace Opm solver.setThresholdPressures(threshold_pressures_by_face_); } + // 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 if( subStepping ) { // create sub step simulator timer with previously used sub step size - SubStepSimulatorTimer subStepper( timer, lastSubStep ); + const double start_time = timer.simulationTimeElapsed(); + const double end_time = start_time + timer.currentStepLength(); + AdaptiveSimulatorTimer subStepper( start_time, end_time, lastSubStep ); + // copy states in case solver has to be restarted + //BlackoilState last_state( state ); + //WellStateFullyImplicitBlackoil last_well_state( well_state ); + + // sub step time loop while( ! subStepper.done() ) { - const double dt_new = solver.step(subStepper.currentStepLength(), state, well_state); - subStepper.next( dt_new ); + // initialize time step control in case current state is needed later + timeStepControl->initialize( state ); + + // keep last state for solver restart and time step control + // (linearIterations < 0 means on convergence in solver) + const int linearIterations = solver.step(subStepper.currentStepLength(), state, well_state); + + // (linearIterations < 0 means on convergence in solver) + if( linearIterations >= 0 ) + { + // advance by current dt + subStepper.advance(); + + // compute new time step estimate + const double dtEstimate = + timeStepControl->computeTimeStepSize( subStepper.currentStepLength(), linearIterations, state ); + // set new time step length + subStepper.provideTimeStepEstimate( dtEstimate ); + } + else // in case of no convergence + { + // we need to revise this + abort(); + subStepper.halfTimeStep(); + std::cerr << "Solver convergence failed, restarting solver with half time step ("<< subStepper.currentStepLength()<<" days)." << std::endl; + // reset states + // state = last_state; + // well_state = last_well_state; + } } subStepper.report( std::cout ); // store last small time step for next reportStep - //lastSubStep = subStepper.maxStepLength(); - //lastSubStep = subStepper.averageStepLength(); - //lastSubStep = subStepper.suggestedMax(); lastSubStep = subStepper.suggestedAverage(); + std::cout << "Last suggested step size = " << lastSubStep/86400.0 << " (days)" << std::endl; - std::cout << "Last suggested step size = " << lastSubStep << std::endl; - if( lastSubStep != lastSubStep ) + if( ! std::isfinite( lastSubStep ) ) // check for NaN lastSubStep = timer.currentStepLength(); } - else + 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(); diff --git a/opm/autodiff/TimeStepControl.hpp b/opm/autodiff/TimeStepControl.hpp index e18d2946b..9e8ee1511 100644 --- a/opm/autodiff/TimeStepControl.hpp +++ b/opm/autodiff/TimeStepControl.hpp @@ -22,117 +22,130 @@ namespace Opm { + /////////////////////////////////////////////////////////////////// + /// + /// 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 state new solution state + /// + /// \return suggested time step size for the next step virtual double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& ) const = 0; + + /// virtual destructor (empty) 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 SimulatorState& ) 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; - } - }; - + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /// + /// PID controller based adaptive time step control as suggested in: + /// Turek and Kuzmin. Algebraic Flux Correction III. Incompressible Flow Problems. Uni Dortmund. + /// + /// See also: + /// D. Kuzmin and S.Turek. Numerical simulation of turbulent bubbly flows. Techreport Uni Dortmund. 2004 + /// + /// and the original article: + /// Valli, Coutinho, and Carey. Adaptive Control for Time Step Selection in Finite Element + /// Simulation of Coupled Viscous Flow and Heat Transfer. Proc of the 10th + /// International Conference on Numerical Methods in Fluids. 1998. + /// + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// class PIDTimeStepControl : public TimeStepControlInterface { protected: mutable std::vector p0_; mutable std::vector sat0_; - mutable double prevDt_; - mutable int prevIterations_; - const int targetIterationCount_; - const double adjustmentFactor_; - - const int upperTargetIterationCount_; - const int lowerTargetIterationCount_; - const double tol_; - mutable std::vector< double > errors_; + const bool verbose_; public: - PIDTimeStepControl( const double tol = 8e-4 ) - : p0_(), sat0_(), prevDt_( 0.0 ), prevIterations_( 0 ), - targetIterationCount_( 100 ), adjustmentFactor_( 1.25 ), - upperTargetIterationCount_( 200 ), lowerTargetIterationCount_( 30 ), - tol_( tol ), - errors_( 3, tol_ ) + /// \brief constructor + /// \param tol tolerance for the relative changes of the numerical solution to be accepted + /// in one time step (default is 1e-3) + PIDTimeStepControl( const double tol = 1e-3, const bool verbose = false ) + : p0_() + , sat0_() + , tol_( tol ) + , errors_( 3, tol_ ) + , verbose_( verbose ) {} + /// \brief \copydoc TimeStepControlInterface::initialize void initialize( const SimulatorState& state ) { - // store current state + // store current state for later time step computation p0_ = state.pressure(); sat0_ = state.saturation(); } + /// \brief \copydoc TimeStepControlInterface::computeTimeStepSize + double computeTimeStepSize( const double dt, const int /* iterations */, const SimulatorState& state ) const + { + const size_t size = p0_.size(); + assert( state.pressure().size() == size ); + assert( state.saturation().size() == size ); + assert( sat0_.size() == size ); + + // compute u^n - u^n+1 + for( size_t i=0; i tol_ ) + { + // adjust dt by given tolerance + if( verbose_ ) + std::cout << "Computed step size (tol): " << (dt * tol_ / error )/86400.0 << " (days)" << std::endl; + return (dt * tol_ / error ); + } + else + { + // values taking from turek time stepping paper + const double kP = 0.075 ; + const double kI = 0.175 ; + const double kD = 0.01 ; + double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) * + std::pow( tol_ / errors_[ 2 ], kI ) * + std::pow( (errors_[0]*errors_[0]/errors_[ 1 ]*errors_[ 2 ]), kD )); + if( verbose_ ) + std::cout << "Computed step size (pow): " << newDt/86400.0 << " (days)" << std::endl; + return newDt; + } + } + + protected: + // return inner product for given container, here std::vector template double inner_product( Iterator it, const Iterator end ) const { @@ -142,55 +155,37 @@ namespace Opm return product; } + }; + + class PIDAndIterationCountTimeStepControl : public PIDTimeStepControl + { + typedef PIDTimeStepControl BaseType; + protected: + const int targetIterationCount_; + + public: + PIDAndIterationCountTimeStepControl( const int target_iterations = 20, + const double tol = 1e-3, + const bool verbose = false) + : BaseType( tol, verbose ) + , targetIterationCount_( target_iterations ) + {} + double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const { - const size_t size = p0_.size(); - assert( state.pressure().size() == size ); - // compute u^n - u^n+1 - for( size_t i=0; i targetIterationCount_ ) { - p0_[ i ] -= state.pressure()[ i ]; - sat0_[ i ] -= state.saturation()[ i ]; + // if iterations was the same or dts were the same, do some magic + dtEstimate *= double( targetIterationCount_ ) / double(iterations); } - // compute || u^n - u^n+1 || - double stateN0 = inner_product( p0_.begin(), p0_.end() ) + - inner_product( sat0_.begin(), sat0_.end() ); - - // compute || u^n+1 || - double stateN = inner_product( state.pressure().begin(), state.pressure().end() ) + - inner_product( state.saturation().begin(), state.saturation().end() ); - - - for( int i=0; i<2; ++i ) - errors_[ i ] = errors_[i+1]; - - double error = stateN0 / stateN ; - errors_[ 2 ] = error ; - - prevDt_ = dt; - prevIterations_ = iterations; - - if( error > tol_ ) - { - // adjust dt by given tolerance - std::cout << "Computed step size (tol): " << (dt * tol_ / error ) << std::endl; - return (dt * tol_ / error ); - } - else - { - const double kP = 0.075 ; - const double kI = 0.175 ; - const double kD = 0.01 ; - double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) * - std::pow( tol_ / errors_[ 2 ], kI ) * - std::pow( (errors_[0]*errors_[0]/errors_[ 1 ]*errors_[ 2 ]), kD )); - std::cout << "Computed step size (pow): " << newDt << std::endl; - return newDt; - } + return dtEstimate; } }; -} // end namespace OPM +} // end namespace OPM #endif From 9e9ef0155ca90c591ff5ce127566a68df9399ac1 Mon Sep 17 00:00:00 2001 From: Robert K Date: Fri, 3 Oct 2014 14:31:57 +0200 Subject: [PATCH 07/19] moved TimeStepControl to Simulator::run and make it work again. --- opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp | 5 ++++- opm/autodiff/TimeStepControl.hpp | 2 +- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index ccf186ab4..a9a83e316 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -300,8 +300,9 @@ namespace Opm const bool subStepping = param_.getDefault("substepping", bool(false) ); + // create time step control object, TODO introduce parameter std::unique_ptr< TimeStepControlInterface > - timeStepControl( new PIDAndIterationCountTimeStepControl( 1e-3, 50 ) ); + timeStepControl( new PIDAndIterationCountTimeStepControl( 100, 1e-3 ) ); // Main simulation loop. while (!timer.done()) { @@ -391,6 +392,8 @@ namespace Opm // compute new time step estimate const double dtEstimate = timeStepControl->computeTimeStepSize( subStepper.currentStepLength(), linearIterations, state ); + std::cout << "Suggested time step size = " << dtEstimate/86400.0 << " (days)" << std::endl; + // set new time step length subStepper.provideTimeStepEstimate( dtEstimate ); } diff --git a/opm/autodiff/TimeStepControl.hpp b/opm/autodiff/TimeStepControl.hpp index 9e8ee1511..75f94fe1f 100644 --- a/opm/autodiff/TimeStepControl.hpp +++ b/opm/autodiff/TimeStepControl.hpp @@ -164,7 +164,7 @@ namespace Opm const int targetIterationCount_; public: - PIDAndIterationCountTimeStepControl( const int target_iterations = 20, + explicit PIDAndIterationCountTimeStepControl( const int target_iterations = 20, const double tol = 1e-3, const bool verbose = false) : BaseType( tol, verbose ) From 2295718f59beaf4a835f890897f9390a9adf1007 Mon Sep 17 00:00:00 2001 From: Robert K Date: Fri, 3 Oct 2014 16:01:59 +0200 Subject: [PATCH 08/19] enabled solver restart again. --- .../SimulatorFullyImplicitBlackoil_impl.hpp | 35 +++++++++++++------ opm/autodiff/TimeStepControl.hpp | 2 +- 2 files changed, 25 insertions(+), 12 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index a9a83e316..77028cc66 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -298,11 +298,12 @@ namespace Opm typename FullyImplicitBlackoilSolver::SolverParameter solverParam( param_ ); + // sub stepping flag, if false just the normal time steps will be used const bool subStepping = param_.getDefault("substepping", bool(false) ); // create time step control object, TODO introduce parameter std::unique_ptr< TimeStepControlInterface > - timeStepControl( new PIDAndIterationCountTimeStepControl( 100, 1e-3 ) ); + timeStepControl( new PIDAndIterationCountTimeStepControl( 50, 8e-4 ) ); // Main simulation loop. while (!timer.done()) { @@ -369,9 +370,9 @@ namespace Opm const double end_time = start_time + timer.currentStepLength(); AdaptiveSimulatorTimer subStepper( start_time, end_time, lastSubStep ); - // copy states in case solver has to be restarted - //BlackoilState last_state( state ); - //WellStateFullyImplicitBlackoil last_well_state( well_state ); + // copy states in case solver has to be restarted (to be revised) + BlackoilState last_state( state ); + WellStateFullyImplicitBlackoil last_well_state( well_state ); // sub step time loop while( ! subStepper.done() ) @@ -379,9 +380,18 @@ namespace Opm // initialize time step control in case current state is needed later timeStepControl->initialize( state ); - // keep last state for solver restart and time step control - // (linearIterations < 0 means on convergence in solver) - const int linearIterations = solver.step(subStepper.currentStepLength(), state, well_state); + int linearIterations = -1; + try { + // (linearIterations < 0 means on convergence in solver) + linearIterations = solver.step(subStepper.currentStepLength(), state, well_state); + + // report number of linear iterations + std::cout << "Overall linear iterations used: " << linearIterations << std::endl; + } + catch (Opm::NumericalProblem) + { + // since linearIterations is < 0 this will restart the solver + } // (linearIterations < 0 means on convergence in solver) if( linearIterations >= 0 ) @@ -396,16 +406,19 @@ namespace Opm // set new time step length subStepper.provideTimeStepEstimate( dtEstimate ); + + // update states + last_state = state ; + last_well_state = well_state; } else // in case of no convergence { // we need to revise this - abort(); - subStepper.halfTimeStep(); + subStepper.provideTimeStepEstimate( 0.1 * subStepper.currentStepLength() ); std::cerr << "Solver convergence failed, restarting solver with half time step ("<< subStepper.currentStepLength()<<" days)." << std::endl; // reset states - // state = last_state; - // well_state = last_well_state; + state = last_state; + well_state = last_well_state; } } diff --git a/opm/autodiff/TimeStepControl.hpp b/opm/autodiff/TimeStepControl.hpp index 75f94fe1f..cca302614 100644 --- a/opm/autodiff/TimeStepControl.hpp +++ b/opm/autodiff/TimeStepControl.hpp @@ -137,7 +137,7 @@ namespace Opm const double kD = 0.01 ; double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) * std::pow( tol_ / errors_[ 2 ], kI ) * - std::pow( (errors_[0]*errors_[0]/errors_[ 1 ]*errors_[ 2 ]), kD )); + std::pow( errors_[0]*errors_[0]/errors_[ 1 ]/errors_[ 2 ], kD )); if( verbose_ ) std::cout << "Computed step size (pow): " << newDt/86400.0 << " (days)" << std::endl; return newDt; From d4802121d3ddd88dcd395da11bedee1725b83f09 Mon Sep 17 00:00:00 2001 From: Robert K Date: Mon, 6 Oct 2014 13:59:21 +0200 Subject: [PATCH 09/19] moved the utility classes to opm-core. --- .../SimulatorFullyImplicitBlackoil_impl.hpp | 87 ++------ opm/autodiff/TimeStepControl.hpp | 191 ------------------ 2 files changed, 14 insertions(+), 264 deletions(-) delete mode 100644 opm/autodiff/TimeStepControl.hpp diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 77028cc66..b6e779bac 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -28,7 +28,6 @@ #include #include #include -#include #include #include @@ -47,6 +46,7 @@ #include #include +#include #include #include @@ -54,6 +54,7 @@ #include #include + #include #include @@ -294,12 +295,14 @@ 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::SolverParameter solverParam( param_ ); - // sub stepping flag, if false just the normal time steps will be used - const bool subStepping = param_.getDefault("substepping", bool(false) ); + // sub stepping + std::unique_ptr< AdaptiveTimeStepping > subStepping; + if( param_.getDefault("timestep.adaptive", bool(false) ) ) + { + subStepping = std::unique_ptr< AdaptiveTimeStepping > (new AdaptiveTimeStepping( param_ )); + } // create time step control object, TODO introduce parameter std::unique_ptr< TimeStepControlInterface > @@ -362,74 +365,12 @@ namespace Opm // 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 - if( subStepping ) - { - // create sub step simulator timer with previously used sub step size - const double start_time = timer.simulationTimeElapsed(); - const double end_time = start_time + timer.currentStepLength(); - AdaptiveSimulatorTimer subStepper( start_time, end_time, lastSubStep ); - - // copy states in case solver has to be restarted (to be revised) - BlackoilState last_state( state ); - WellStateFullyImplicitBlackoil last_well_state( well_state ); - - // sub step time loop - while( ! subStepper.done() ) - { - // initialize time step control in case current state is needed later - timeStepControl->initialize( state ); - - int linearIterations = -1; - try { - // (linearIterations < 0 means on convergence in solver) - linearIterations = solver.step(subStepper.currentStepLength(), state, well_state); - - // report number of linear iterations - std::cout << "Overall linear iterations used: " << linearIterations << std::endl; - } - catch (Opm::NumericalProblem) - { - // since linearIterations is < 0 this will restart the solver - } - - // (linearIterations < 0 means on convergence in solver) - if( linearIterations >= 0 ) - { - // advance by current dt - subStepper.advance(); - - // compute new time step estimate - const double dtEstimate = - timeStepControl->computeTimeStepSize( subStepper.currentStepLength(), linearIterations, state ); - std::cout << "Suggested time step size = " << dtEstimate/86400.0 << " (days)" << std::endl; - - // set new time step length - subStepper.provideTimeStepEstimate( dtEstimate ); - - // update states - last_state = state ; - last_well_state = well_state; - } - else // in case of no convergence - { - // we need to revise this - subStepper.provideTimeStepEstimate( 0.1 * subStepper.currentStepLength() ); - std::cerr << "Solver convergence failed, restarting solver with half time step ("<< subStepper.currentStepLength()<<" days)." << std::endl; - // reset states - state = last_state; - well_state = last_well_state; - } - } - - subStepper.report( std::cout ); - - // store last small time step for next reportStep - lastSubStep = subStepper.suggestedAverage(); - std::cout << "Last suggested step size = " << lastSubStep/86400.0 << " (days)" << std::endl; - - if( ! std::isfinite( lastSubStep ) ) // check for NaN - lastSubStep = timer.currentStepLength(); + // + // \Note: The report steps are met in any case + // \Note: The sub stepping will require a copy of the state variables + if( subStepping ) { + subStepping->step( solver, state, well_state, + timer.simulationTimeElapsed(), timer.currentStepLength() ); } else { // solve for complete report step diff --git a/opm/autodiff/TimeStepControl.hpp b/opm/autodiff/TimeStepControl.hpp deleted file mode 100644 index cca302614..000000000 --- a/opm/autodiff/TimeStepControl.hpp +++ /dev/null @@ -1,191 +0,0 @@ -/* - 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 . -*/ -#ifndef OPM_TIMESTEPCONTROL_HEADER_INCLUDED -#define OPM_TIMESTEPCONTROL_HEADER_INCLUDED - -namespace Opm -{ - - /////////////////////////////////////////////////////////////////// - /// - /// 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 state new solution state - /// - /// \return suggested time step size for the next step - virtual double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& ) const = 0; - - /// virtual destructor (empty) - virtual ~TimeStepControlInterface () {} - }; - - /////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /// - /// PID controller based adaptive time step control as suggested in: - /// Turek and Kuzmin. Algebraic Flux Correction III. Incompressible Flow Problems. Uni Dortmund. - /// - /// See also: - /// D. Kuzmin and S.Turek. Numerical simulation of turbulent bubbly flows. Techreport Uni Dortmund. 2004 - /// - /// and the original article: - /// Valli, Coutinho, and Carey. Adaptive Control for Time Step Selection in Finite Element - /// Simulation of Coupled Viscous Flow and Heat Transfer. Proc of the 10th - /// International Conference on Numerical Methods in Fluids. 1998. - /// - /////////////////////////////////////////////////////////////////////////////////////////////////////////////// - class PIDTimeStepControl : public TimeStepControlInterface - { - protected: - mutable std::vector p0_; - mutable std::vector sat0_; - - const double tol_; - mutable std::vector< double > errors_; - - const bool verbose_; - public: - /// \brief constructor - /// \param tol tolerance for the relative changes of the numerical solution to be accepted - /// in one time step (default is 1e-3) - PIDTimeStepControl( const double tol = 1e-3, const bool verbose = false ) - : p0_() - , sat0_() - , tol_( tol ) - , errors_( 3, tol_ ) - , verbose_( verbose ) - {} - - /// \brief \copydoc TimeStepControlInterface::initialize - void initialize( const SimulatorState& state ) - { - // store current state for later time step computation - p0_ = state.pressure(); - sat0_ = state.saturation(); - } - - /// \brief \copydoc TimeStepControlInterface::computeTimeStepSize - double computeTimeStepSize( const double dt, const int /* iterations */, const SimulatorState& state ) const - { - const size_t size = p0_.size(); - assert( state.pressure().size() == size ); - assert( state.saturation().size() == size ); - assert( sat0_.size() == size ); - - // compute u^n - u^n+1 - for( size_t i=0; i tol_ ) - { - // adjust dt by given tolerance - if( verbose_ ) - std::cout << "Computed step size (tol): " << (dt * tol_ / error )/86400.0 << " (days)" << std::endl; - return (dt * tol_ / error ); - } - else - { - // values taking from turek time stepping paper - const double kP = 0.075 ; - const double kI = 0.175 ; - const double kD = 0.01 ; - double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) * - std::pow( tol_ / errors_[ 2 ], kI ) * - std::pow( errors_[0]*errors_[0]/errors_[ 1 ]/errors_[ 2 ], kD )); - if( verbose_ ) - std::cout << "Computed step size (pow): " << newDt/86400.0 << " (days)" << std::endl; - return newDt; - } - } - - protected: - // return inner product for given container, here std::vector - template - double inner_product( Iterator it, const Iterator end ) const - { - double product = 0.0 ; - for( ; it != end; ++it ) - product += ( *it * *it ); - return product; - } - - }; - - class PIDAndIterationCountTimeStepControl : public PIDTimeStepControl - { - typedef PIDTimeStepControl BaseType; - protected: - const int targetIterationCount_; - - public: - explicit PIDAndIterationCountTimeStepControl( const int target_iterations = 20, - const double tol = 1e-3, - const bool verbose = false) - : BaseType( tol, verbose ) - , targetIterationCount_( target_iterations ) - {} - - double computeTimeStepSize( const double dt, const int iterations, const SimulatorState& state ) const - { - double dtEstimate = BaseType :: computeTimeStepSize( dt, iterations, state ); - - // further reduce step size if to many iterations were used - if( iterations > targetIterationCount_ ) - { - // if iterations was the same or dts were the same, do some magic - dtEstimate *= double( targetIterationCount_ ) / double(iterations); - } - - return dtEstimate; - } - }; - - -} // end namespace OPM -#endif From c42eeffdeb3c7fce6087a4cd3509ce33869f860f Mon Sep 17 00:00:00 2001 From: Robert K Date: Mon, 6 Oct 2014 14:27:55 +0200 Subject: [PATCH 10/19] remove unused output --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index cb6138dac..725c9b003 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -1702,7 +1702,6 @@ namespace { for (; quantityIt != endQuantityIt; ++quantityIt) { const double quantityResid = (*quantityIt).value().matrix().norm(); if (!std::isfinite(quantityResid)) { - //std::cout << quantityResid << " quantity" << std::endl; const int trouble_phase = quantityIt - residual_.material_balance_eq.begin(); OPM_THROW(Opm::NumericalProblem, "Encountered a non-finite residual in material balance equation " From a8c0f7df9214d2c30ce2347eccecd134b0507043 Mon Sep 17 00:00:00 2001 From: Robert K Date: Mon, 6 Oct 2014 15:53:17 +0200 Subject: [PATCH 11/19] apply changes made in opm-core. --- .../SimulatorFullyImplicitBlackoil_impl.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index b6e779bac..76b76f2e4 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -297,15 +297,15 @@ namespace Opm typename FullyImplicitBlackoilSolver::SolverParameter solverParam( param_ ); - // sub stepping - std::unique_ptr< AdaptiveTimeStepping > subStepping; + // adaptive time stepping + std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; if( param_.getDefault("timestep.adaptive", bool(false) ) ) { - subStepping = std::unique_ptr< AdaptiveTimeStepping > (new AdaptiveTimeStepping( param_ )); + adaptiveTimeStepping = std::unique_ptr< AdaptiveTimeStepping > (new AdaptiveTimeStepping( param_ )); } // create time step control object, TODO introduce parameter - std::unique_ptr< TimeStepControlInterface > + std::unique_ptr< TimeStepControlInterface > timeStepControl( new PIDAndIterationCountTimeStepControl( 50, 8e-4 ) ); // Main simulation loop. @@ -366,10 +366,10 @@ namespace Opm // 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 report steps are met in any case // \Note: The sub stepping will require a copy of the state variables - if( subStepping ) { - subStepping->step( solver, state, well_state, + if( adaptiveTimeStepping ) { + adaptiveTimeStepping->step( solver, state, well_state, timer.simulationTimeElapsed(), timer.currentStepLength() ); } else { From c2e6b368ae7781e8d3ca5fada9598964d4e5d5bd Mon Sep 17 00:00:00 2001 From: Robert K Date: Mon, 6 Oct 2014 15:59:01 +0200 Subject: [PATCH 12/19] revert Schur fix. --- opm/autodiff/NewtonIterationBlackoilCPR.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp index cd5f59098..4b211af8f 100644 --- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp +++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -247,8 +247,7 @@ namespace Opm // std::cout << "++++++++++++++++++++++++++++++++++++++++++++\n" // << D // << "++++++++++++++++++++++++++++++++++++++++++++\n" << std::endl; - std::cerr << "WARNING (ignored): Cannot do Schur complement with respect to non-diagonal block." << std::endl; - //OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); + OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); } V diag = D.diagonal(); Eigen::DiagonalMatrix invD = (1.0 / diag).matrix().asDiagonal(); @@ -303,8 +302,7 @@ namespace Opm // Find inv(D). const M& D = equation.derivative()[n]; if (!isDiagonal(D)) { - std::cerr << "WARNING (ignored): Cannot do Schur complement with respect to non-diagonal block." << std::endl; - //OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); + OPM_THROW(std::logic_error, "Cannot do Schur complement with respect to non-diagonal block."); } V diag = D.diagonal(); Eigen::DiagonalMatrix invD = (1.0 / diag).matrix().asDiagonal(); From cd7b6ce7f05f6aacd0544d92e684451b120789da Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Fri, 17 Oct 2014 12:25:13 +0200 Subject: [PATCH 13/19] remove blank at end of line. --- opm/autodiff/NewtonIterationBlackoilSimple.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.cpp b/opm/autodiff/NewtonIterationBlackoilSimple.cpp index c8d2472d9..228225b4c 100644 --- a/opm/autodiff/NewtonIterationBlackoilSimple.cpp +++ b/opm/autodiff/NewtonIterationBlackoilSimple.cpp @@ -60,7 +60,7 @@ namespace Opm matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(), total_residual.value().data(), dx.data()); - // store iterations + // store iterations iterations_ = rep.iterations; if (!rep.converged) { From fb32376d8f8421a3a133869d43cbbaa185c05b19 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Fri, 17 Oct 2014 12:40:25 +0200 Subject: [PATCH 14/19] throw exception when convergence failed, also in NewtonSolver. --- opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 4190de1fb..a13289b67 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -338,8 +338,8 @@ namespace { } if (!converged) { - std::cerr << "ERROR: Failed to compute converged solution in " << it << " iterations." << std::endl; - // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); + // the runtime_error is caught by the AdaptiveTimeStepping + OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); return -1; } From 42e4043c18679e0f32a83d627aeacedb18857616 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Mon, 20 Oct 2014 14:47:33 +0200 Subject: [PATCH 15/19] remove unused variable. --- opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 76b76f2e4..cb3b7028e 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -301,13 +301,9 @@ namespace Opm std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; if( param_.getDefault("timestep.adaptive", bool(false) ) ) { - adaptiveTimeStepping = std::unique_ptr< AdaptiveTimeStepping > (new AdaptiveTimeStepping( param_ )); + adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_ ) ); } - // create time step control object, TODO introduce parameter - std::unique_ptr< TimeStepControlInterface > - timeStepControl( new PIDAndIterationCountTimeStepControl( 50, 8e-4 ) ); - // Main simulation loop. while (!timer.done()) { // Report timestep. From e68c58fb5971ea833064e02adbeb51f4c0f62e14 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Mon, 20 Oct 2014 14:47:45 +0200 Subject: [PATCH 16/19] added linear to docu. --- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 2 +- opm/autodiff/NewtonIterationBlackoilInterface.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index f0eb90111..dfc85b715 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -114,7 +114,7 @@ namespace Opm { /// \param[in] dt time step size /// \param[in] state reservoir state /// \param[in] wstate well state - /// \return number of iterations used + /// \return number of linear iterations used int step(const double dt , BlackoilState& state , diff --git a/opm/autodiff/NewtonIterationBlackoilInterface.hpp b/opm/autodiff/NewtonIterationBlackoilInterface.hpp index c105cc783..2e26db07e 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterface.hpp +++ b/opm/autodiff/NewtonIterationBlackoilInterface.hpp @@ -40,7 +40,7 @@ namespace Opm /// \return the solution x virtual SolutionVector computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const = 0; - /// \return number of iterations used during last call of computeNewtonIncrement + /// \return number of linear iterations used during last call of computeNewtonIncrement virtual int iterations () const = 0; }; From 937555bb3be199cf409fd2ab93b13bda3e1c1731 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Mon, 20 Oct 2014 15:18:14 +0200 Subject: [PATCH 17/19] use group for timestep parameters. --- opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index cb3b7028e..3a11cdea0 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -299,9 +299,10 @@ namespace Opm // adaptive time stepping std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; - if( param_.getDefault("timestep.adaptive", bool(false) ) ) + const parameter::ParameterGroup& timeStepParam = param_.getGroup( "timestep" ); + if( timeStepParam.getDefault("adaptive", bool(false) ) ) { - adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_ ) ); + adaptiveTimeStepping.reset( new AdaptiveTimeStepping( timeStepParam ) ); } // Main simulation loop. From 2a7f951e259854a801141ed95687a2a357ee959b Mon Sep 17 00:00:00 2001 From: Robert K Date: Wed, 22 Oct 2014 15:00:10 +0200 Subject: [PATCH 18/19] provide default for parameter group. --- opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 3a11cdea0..6b7da7672 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -299,7 +299,7 @@ namespace Opm // adaptive time stepping std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; - const parameter::ParameterGroup& timeStepParam = param_.getGroup( "timestep" ); + const parameter::ParameterGroup& timeStepParam = param_.getDefault( "timestep", param_ ); if( timeStepParam.getDefault("adaptive", bool(false) ) ) { adaptiveTimeStepping.reset( new AdaptiveTimeStepping( timeStepParam ) ); From 05608a61455baa40f212def30043201112da65b0 Mon Sep 17 00:00:00 2001 From: Robert K Date: Wed, 22 Oct 2014 15:29:20 +0200 Subject: [PATCH 19/19] stick to previous method of parameter extraction. grouping is not what we wanted here. --- opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp index 6b7da7672..cb3b7028e 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp @@ -299,10 +299,9 @@ namespace Opm // adaptive time stepping std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; - const parameter::ParameterGroup& timeStepParam = param_.getDefault( "timestep", param_ ); - if( timeStepParam.getDefault("adaptive", bool(false) ) ) + if( param_.getDefault("timestep.adaptive", bool(false) ) ) { - adaptiveTimeStepping.reset( new AdaptiveTimeStepping( timeStepParam ) ); + adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_ ) ); } // Main simulation loop.