Detailed Description

This class gives an adaptive algorithm whose time step control is evolutionary PID controller. With an definition of relative solution change \(e_n=\frac{\|u^{n+1}-u^n\|}{\|u^{n+1}\|}\), the algorithm gives a time step size estimation as

\[ h_{n+1} = \left(\frac{e_{n-1}}{e_n}\right)^{k_P} \left(\frac{TOL}{e_n}\right)^{k_I} \left(\frac{e^2_{n-1}}{e_n e_{n-2}}\right)^{k_D} \]

where \(k_P=0.075\), \(k_I=0.175\), \(k_D=0.01\) are empirical PID parameters.

In the computation, \( e_n\) is calculated firstly. If \(e_n>TOL\), the current time step is rejected and repeated with a new time step size of \(h=\frac{TOL}{e_n} h_n\).

Limits of the time step size are given as

\[ h_{\mbox{min}} \leq h_{n+1} \leq h_{\mbox{max}}, l \leq \frac{h_{n+1}}{h_n} \leq L \]

Similar algorithm can be found in [1] .

Definition at line 46 of file EvolutionaryPIDcontroller.h.

#include <EvolutionaryPIDcontroller.h>

Inheritance diagram for NumLib::EvolutionaryPIDcontroller:

[legend]

Collaboration diagram for NumLib::EvolutionaryPIDcontroller:

[legend]

Public Member Functions
	EvolutionaryPIDcontroller (const double t0, const double t_end, const double h0, const double h_min, const double h_max, const double rel_h_min, const double rel_h_max, const double tol, std::vector< double > const &fixed_times_for_output)

std::tuple< bool, double >	next (double solution_error, int number_iterations, NumLib::TimeStep &timestep_previous, NumLib::TimeStep &timestep_current) override

bool	isSolutionErrorComputationNeeded () const override

virtual bool	canReduceTimestepSize (NumLib::TimeStep const &timestep_previous, NumLib::TimeStep const &timestep_current) const override
	Query the timestepper if further time step size reduction is possible.

Public Member Functions inherited from NumLib::TimeStepAlgorithm
	TimeStepAlgorithm (const double t0, const double t_end)

virtual	~TimeStepAlgorithm ()=default

double	begin () const
	return the beginning of time steps

double	end () const
	return the end of time steps

virtual void	resetCurrentTimeStep (const double, TimeStep &, TimeStep &)
	reset the current step size from the previous time

Private Member Functions
double	limitStepSize (const double h_new, const bool previous_step_accepted, NumLib::TimeStep const &timestep_current) const

Private Attributes
const double	_kP = 0.075
	Parameter.

const double	_kI = 0.175
	Parameter.

const double	_kD = 0.01
	Parameter.

const double	_h0
	initial time step size.

const double	_h_min
	minimum step size.

const double	_h_max
	maximum step size.

const double	_rel_h_min
	\(l\) in \( h_{\mbox{min}} \leq h_{n+1} \leq h_{\mbox{max}},\)

const double	_rel_h_max
	\(L\) in \( h_{\mbox{min}} \leq h_{n+1} \leq h_{\mbox{max}},\)

const double	_tol

double	_e_n_minus1
	\(e_{n-1}\).

double	_e_n_minus2
	\(e_{n-2}\).

std::vector< double > const	_fixed_times_for_output

Additional Inherited Members
Protected Attributes inherited from NumLib::TimeStepAlgorithm
const double	_t_initial
	initial time

const double	_t_end
	end time

Constructor & Destructor Documentation

◆ EvolutionaryPIDcontroller()

NumLib::EvolutionaryPIDcontroller::EvolutionaryPIDcontroller	(	const double	t0,
		const double	t_end,
		const double	h0,
		const double	h_min,
		const double	h_max,
		const double	rel_h_min,
		const double	rel_h_max,
		const double	tol,
		std::vector< double > const &	fixed_times_for_output )

inline

Definition at line 49 of file EvolutionaryPIDcontroller.h.

        : TimeStepAlgorithm(t0, t_end),
          _h0(h0),
          _h_min(h_min),
          _h_max(h_max),
          _rel_h_min(rel_h_min),
          _rel_h_max(rel_h_max),
          _tol(tol),
          _e_n_minus1(0.),
          _e_n_minus2(0.),
          _fixed_times_for_output(fixed_times_for_output)
    {
    }

Member Function Documentation

◆ canReduceTimestepSize()

bool NumLib::EvolutionaryPIDcontroller::canReduceTimestepSize	(	NumLib::TimeStep const &	,
		NumLib::TimeStep const &	) const

overridevirtual

Query the timestepper if further time step size reduction is possible.

Reimplemented from NumLib::TimeStepAlgorithm.

Definition at line 170 of file EvolutionaryPIDcontroller.cpp.

{
    return NumLib::canReduceTimestepSize(timestep_previous, timestep_current,
                                         _h_min);
}

References _h_min, and NumLib::canReduceTimestepSize().

◆ isSolutionErrorComputationNeeded()

bool NumLib::EvolutionaryPIDcontroller::isSolutionErrorComputationNeeded ( ) const

inlineoverridevirtual

Get a flag to indicate whether this algorithm needs to compute solution error. The default return value is false.

Reimplemented from NumLib::TimeStepAlgorithm.

Definition at line 72 of file EvolutionaryPIDcontroller.h.

72{ return true; }

◆ limitStepSize()

double NumLib::EvolutionaryPIDcontroller::limitStepSize	(	const double	h_new,
		const bool	previous_step_accepted,
		NumLib::TimeStep const &	timestep_current ) const

private

Force the computed time step size in the given range (see the formulas in the documentation of the class) or use the half of the previous time step size under some other constrains.

Parameters

h_new	The computed time step size.
previous_step_accepted	An indicator for whether the previous time step is rejected.
timestep_current	the current time step

Returns: The new time step after apply the constrains.

Definition at line 105 of file EvolutionaryPIDcontroller.cpp.

{
    const double h_n = timestep_current.dt();
    // Forced the computed time step size in the given range
    // (see the formulas in the documentation of the class)
    const double h_in_range = std::max(_h_min, std::min(h_new, _h_max));
    // Limit the step size change ratio.
    double limited_h =
        std::max(_rel_h_min * h_n, std::min(h_in_range, _rel_h_max * h_n));
 
    if (!previous_step_accepted)
    {
        // If the last time step was rejected and the new predicted time step
        // size is identical to that of the previous rejected step, the new
        // step size is then reduced by half.
        if (std::abs(limited_h - timestep_current.dt()) <
            std::numeric_limits<double>::min())
        {
            limited_h = std::max(_h_min, 0.5 * limited_h);
        }
 
        // If the last time step was rejected and the new predicted time step
        // size is larger than the step size of the rejected step, the new step
        // size takes the half of the size of the rejected step. This could
        // happen when a time step is rejected due to a diverged non-linear
        // solver. In such case, this algorithm may give a large time step size
        // by using the diverged solution.
        if (limited_h > timestep_current.dt())
        {
            limited_h = std::max(_h_min, 0.5 * timestep_current.dt());
        }
    }
 
    if (_fixed_times_for_output.empty())
    {
        return limited_h;
    }
 
    // find first fixed timestep for output larger than the current time, i.e.,
    // current time < fixed output time
    auto fixed_output_time_it = std::find_if(
        std::begin(_fixed_times_for_output), std::end(_fixed_times_for_output),
        [&timestep_current](auto const fixed_output_time)
        { return timestep_current.current() < fixed_output_time; });
 
    if (fixed_output_time_it != _fixed_times_for_output.end())
    {
        // check if the fixed output time is in the interval
        // (current time, current time + limited_h)
        if (*fixed_output_time_it < timestep_current.current() + limited_h)
        {
            // check if the potential adjusted time step is larger than zero
            if (std::abs(*fixed_output_time_it - timestep_current.current()) >
                std::numeric_limits<double>::epsilon() *
                    timestep_current.current())
            {
                return *fixed_output_time_it - timestep_current.current();
            }
        }
    }
    return limited_h;
}

References _fixed_times_for_output, _h_max, _h_min, _rel_h_max, _rel_h_min, NumLib::TimeStep::current(), and NumLib::TimeStep::dt().

Referenced by next().

◆ next()

std::tuple< bool, double > NumLib::EvolutionaryPIDcontroller::next	(	double	solution_error,
		int	number_iterations,
		NumLib::TimeStep &	ts_previous,
		NumLib::TimeStep &	ts_current )

overridevirtual

Move to the next time step

Parameters

solution_error	Solution error \(e_n\) between two successive time steps.
number_iterations	Number of non-linear iterations used.
ts_previous	the previous time step used to compute the size of the next step
ts_current	the current time step used to compute the size of the next step

Returns: A step acceptance flag and the computed step size.

Implements NumLib::TimeStepAlgorithm.

Definition at line 23 of file EvolutionaryPIDcontroller.cpp.

{
    const bool is_previous_step_accepted = timestep_current.isAccepted();
 
    const double e_n = solution_error;
    const double zero_threshold = std::numeric_limits<double>::epsilon();
    // step rejected.
    if (e_n > _tol)
    {
        timestep_current.setAccepted(false);
 
        double h_new = (e_n > zero_threshold)
                           ? timestep_current.dt() * _tol / e_n
                           : 0.5 * timestep_current.dt();
 
        h_new =
            limitStepSize(h_new, is_previous_step_accepted, timestep_current);
 
        WARN(
            "This step is rejected due to the relative change from the"
            " solution of the previous\n"
            "\t time step to the current solution exceeds the given tolerance"
            " of {:g}.\n"
            "\t This time step will be repeated with a new time step size of"
            " {:g}\n"
            "\t or the simulation will be halted.",
            _tol, h_new);
 
        return std::make_tuple(timestep_current.isAccepted(), h_new);
    }
 
    // step accepted.
    timestep_current.setAccepted(true);
 
    if (timestep_current.timeStepNumber() == 0)
    {
        _e_n_minus1 = e_n;
 
        return std::make_tuple(timestep_current.isAccepted(), _h0);
    }
    else
    {
        const double h_n = timestep_current.dt();
        double h_new = h_n;
 
        if (e_n > zero_threshold)
        {
            if (_e_n_minus1 > zero_threshold)
            {
                if (_e_n_minus2 > zero_threshold)
                {
                    h_new = std::pow(_e_n_minus1 / e_n, _kP) *
                            std::pow(_tol / e_n, _kI) *
                            std::pow(
                                _e_n_minus1 * _e_n_minus1 / (e_n * _e_n_minus2),
                                _kD) *
                            h_n;
                }
                else
                {
                    h_new = std::pow(_e_n_minus1 / e_n, _kP) *
                            std::pow(_tol / e_n, _kI) * h_n;
                }
            }
            else
            {
                h_new = std::pow(_tol / e_n, _kI) * h_n;
            }
        }
 
        h_new =
            limitStepSize(h_new, is_previous_step_accepted, timestep_current);
 
        _e_n_minus2 = _e_n_minus1;
        _e_n_minus1 = e_n;
 
        return std::make_tuple(timestep_current.isAccepted(), h_new);
    }
}