diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
index 71b17e2304683c26c90a4d62182d36465a622ab1..54e4ae16e69b69524994a97885af469998f03f9c 100644
--- a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
+++ b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
@@ -22,6 +22,8 @@ namespace NumLib
{
bool EvolutionaryPIDcontroller::next(const double solution_error)
{
+ const bool is_previous_step_accepted = _is_accepted;
+
const double e_n = solution_error;
const double zero_threshlod = std::numeric_limits<double>::epsilon();
// step rejected.
@@ -32,19 +34,22 @@ bool EvolutionaryPIDcontroller::next(const double solution_error)
double h_new = (e_n > zero_threshlod) ? _ts_current.dt() * _tol / e_n
: 0.5 * _ts_current.dt();
- h_new = limitStepSize(h_new, _ts_current.dt());
+ h_new =
+ limitStepSize(h_new, _ts_current.dt(), is_previous_step_accepted);
h_new = checkSpecificTimeReached(h_new);
_ts_current = _ts_prev;
_ts_current += h_new;
- 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);
+ 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 false;
}
@@ -90,7 +95,7 @@ bool EvolutionaryPIDcontroller::next(const double solution_error)
}
}
- h_new = limitStepSize(h_new, h_n);
+ h_new = limitStepSize(h_new, h_n, is_previous_step_accepted);
h_new = checkSpecificTimeReached(h_new);
_dt_vector.push_back(h_new);
@@ -104,6 +109,35 @@ bool EvolutionaryPIDcontroller::next(const double solution_error)
return true;
}
+double EvolutionaryPIDcontroller::limitStepSize(
+ const double h_new, const double h_n,
+ const bool previous_step_accepted) const
+{
+ const double h_in_range = std::max(_h_min, std::min(h_new, _h_max));
+ 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::fabs(limited_h - _ts_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 > _ts_current.dt())
+ limited_h = std::max(_h_min, 0.5 * _ts_current.dt());
+ }
+ return limited_h;
+}
+
double EvolutionaryPIDcontroller::checkSpecificTimeReached(const double h_new)
{
if (_specific_times.empty())
diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h
index bd22e0080f0dfa6a46ea334b12cb4dfc78c2e619..5f0581b96710051a5048ed639e5d18c7e1130acd 100644
--- a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h
+++ b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h
@@ -81,10 +81,16 @@ public:
/// return if current time step is accepted
bool accepted() const override { return _is_accepted; }
+ /// Set the status of the step
+ /// \param accepted A flag for whether the step is accepted or not
+ void setAcceptedOrNot(const bool accepted) override
+ {
+ _is_accepted = accepted;
+ }
+
/// Get a flag to indicate that this algorithm need to compute
/// solution error.
bool isSolutionErrorComputationNeeded() override { return true; }
-
private:
const double _kP = 0.075; ///< Parameter. \see EvolutionaryPIDcontroller
const double _kI = 0.175; ///< Parameter. \see EvolutionaryPIDcontroller
@@ -109,12 +115,8 @@ private:
bool _is_accepted;
- double limitStepSize(const double h_new, const double h_n) const
- {
- const double h_in_range = std::max(_h_min, std::min(h_new, _h_max));
- return std::max(_rel_h_min * h_n,
- std::min(h_in_range, _rel_h_max * h_n));
- }
+ double limitStepSize(const double h_new, const double h_n,
+ const bool previous_step_accepted) const;
double checkSpecificTimeReached(const double h_new);
};
diff --git a/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h b/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h
index a74cd0f6a9fea415f7eefab749f3c88650a8e49a..659dd07c6145cbeb70208d1d3131599f31e0f6db 100644
--- a/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h
+++ b/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h
@@ -72,6 +72,10 @@ public:
/// return if current time step is accepted or not
virtual bool accepted() const = 0;
+ /// Set the status of the step. A boolean parameter is needed to indicated
+ /// whether the step is accepted or not.
+ virtual void setAcceptedOrNot(const bool /*accepted*/){};
+
/// return a history of time step sizes
const std::vector<double>& getTimeStepSizeHistory() const
{
@@ -81,7 +85,6 @@ public:
/// Get a flag to indicate whether this algorithm needs to compute
/// solution error. The default return value is false.
virtual bool isSolutionErrorComputationNeeded() { return false; }
-
protected:
/// initial time
const double _t_initial;
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index 594525cfbdc00eedf8692fb0f98350b44e0ecad6..9f73bb1e06a319b6beabacf217b73ca904b82c39 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -576,6 +576,9 @@ double UncoupledProcessesTimeLoop::computeTimeStepping(
x, norm_type))
: 0.;
+ if (!ppd.nonlinear_solver_converged)
+ timestepper->setAcceptedOrNot(false);
+
if (!timestepper->next(solution_error) &&
// In case of FixedTimeStepping, which makes timestepper->next(...)
// return false when the ending time is reached.