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