From a723a01f7258d8e5bec611618dcf610dca88c19b Mon Sep 17 00:00:00 2001 From: Robert K Date: Fri, 3 Oct 2014 14:18:31 +0200 Subject: [PATCH] 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