From d03f9411b673e1c84d7a4fc5e7bd0457dbd5c53a Mon Sep 17 00:00:00 2001 From: Robert K Date: Wed, 1 Oct 2014 13:50:08 +0200 Subject: [PATCH] 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