From f3400f75cd15eae68edfae043f0e4c411836a80f Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Thu, 20 Apr 2017 18:09:29 +0200
Subject: [PATCH] [dt] Added a member of NormType
---
NumLib/ODESolver/TimeDiscretization.cpp | 20 ++++++---
.../Algorithms/EvolutionaryPIDcontroller.cpp | 24 +++++++++--
.../Algorithms/EvolutionaryPIDcontroller.h | 13 +-----
.../Algorithms/FixedTimeStepping.cpp | 7 +++-
.../Algorithms/FixedTimeStepping.h | 1 -
.../Algorithms/ITimeStepAlgorithm.h | 41 +++++++++++++++----
...erationNumberBasedAdaptiveTimeStepping.cpp | 6 +--
7 files changed, 76 insertions(+), 36 deletions(-)
diff --git a/NumLib/ODESolver/TimeDiscretization.cpp b/NumLib/ODESolver/TimeDiscretization.cpp
index 4e84ac8ccf9..a3236ca1c5c 100644
--- a/NumLib/ODESolver/TimeDiscretization.cpp
+++ b/NumLib/ODESolver/TimeDiscretization.cpp
@@ -25,6 +25,9 @@ namespace NumLib
static double computeRelativeError(GlobalVector const& x, GlobalVector& x_old,
MathLib::VecNormType norm_type)
{
+ if (norm_type == MathLib::VecNormType::INVALID)
+ return 0.;
+
const double norm2_x = MathLib::LinAlg::norm(x, norm_type);
assert(norm2_x > std::numeric_limits<double>::epsilon());
@@ -33,25 +36,30 @@ static double computeRelativeError(GlobalVector const& x, GlobalVector& x_old,
return MathLib::LinAlg::norm(x_old, norm_type) / norm2_x;
}
-double BackwardEuler::getRelativeError(GlobalVector const& x, MathLib::VecNormType norm_type)
+double BackwardEuler::getRelativeError(GlobalVector const& x,
+ MathLib::VecNormType norm_type)
{
return computeRelativeError(x, _x_old, norm_type);
}
-double ForwardEuler::getRelativeError(GlobalVector const& x, MathLib::VecNormType norm_type)
+double ForwardEuler::getRelativeError(GlobalVector const& x,
+ MathLib::VecNormType norm_type)
{
return computeRelativeError(x, _x_old, norm_type);
}
-double CrankNicolson::getRelativeError(GlobalVector const& x, MathLib::VecNormType norm_type)
+double CrankNicolson::getRelativeError(GlobalVector const& x,
+ MathLib::VecNormType norm_type)
{
return computeRelativeError(x, _x_old, norm_type);
}
-double BackwardDifferentiationFormula::getRelativeError(GlobalVector const& x, MathLib::VecNormType norm_type)
+double BackwardDifferentiationFormula::getRelativeError(
+ GlobalVector const& x, MathLib::VecNormType norm_type)
{
- return (_xs_old.size() < _num_steps) ? 0. : computeRelativeError(
- x, *_xs_old[_offset], norm_type);
+ return (_xs_old.size() < _num_steps)
+ ? 0.
+ : computeRelativeError(x, *_xs_old[_offset], norm_type);
}
void BackwardDifferentiationFormula::pushState(const double,
diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
index fda6f69bb33..df1af35e162 100644
--- a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
+++ b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
@@ -42,8 +42,6 @@ bool EvolutionaryPIDcontroller::next(const double solution_error)
// step accepted.
_is_accepted = true;
- if (_specific_times.size() > 0)
- _specific_times.pop_back();
if (_ts_current.steps() == 0)
{
@@ -81,10 +79,11 @@ bool EvolutionaryPIDcontroller::next(const double solution_error)
h_new = limitStepSize(h_new, h_n);
}
- _dt_vector.push_back(h_new);
+ const double checked_h_new = checkSpecificTimeReached(h_new);
+ _dt_vector.push_back(checked_h_new);
_ts_prev = _ts_current;
- _ts_current += h_new;
+ _ts_current += checked_h_new;
_e_n_minus2 = _e_n_minus1;
_e_n_minus1 = e_n;
@@ -93,6 +92,23 @@ bool EvolutionaryPIDcontroller::next(const double solution_error)
return true;
}
+double EvolutionaryPIDcontroller::checkSpecificTimeReached(const double h_new)
+{
+ if (_specific_times.empty())
+ return h_new;
+
+ const double specific_time = _specific_times.back();
+ const double zero_threshlod = std::numeric_limits<double>::epsilon();
+ if ((specific_time > _ts_current.current()) &&
+ (_ts_current.current() + h_new - specific_time > zero_threshlod))
+ {
+ _specific_times.pop_back();
+ return specific_time - _ts_current.current();
+ }
+
+ return h_new;
+}
+
/// Create an EvolutionaryPIDcontroller time stepper from the given
/// configuration
std::unique_ptr<ITimeStepAlgorithm> createEvolutionaryPIDcontroller(
diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h
index 0f167dbce79..9cdc115542b 100644
--- a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h
+++ b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h
@@ -16,8 +16,6 @@
#include "ITimeStepAlgorithm.h"
-#include "MathLib/LinAlg/LinAlg.h"
-
namespace BaseLib
{
class ConfigTree;
@@ -60,7 +58,7 @@ public:
const std::vector<double>&& specific_times,
const double tol,
const MathLib::VecNormType norm_type)
- : ITimeStepAlgorithm(t0, t_end),
+ : ITimeStepAlgorithm(t0, t_end, norm_type),
_h0(h0),
_h_min(h_min),
_h_max(h_max),
@@ -68,7 +66,6 @@ public:
_rel_h_max(rel_h_max),
_specific_times(std::move(specific_times)),
_tol(tol),
- _norm_type(norm_type),
_e_n_minus1(0.),
_e_n_minus2(0.),
_is_accepted(true)
@@ -103,7 +100,6 @@ private:
std::vector<double> _specific_times;
const double _tol;
- const MathLib::VecNormType _norm_type;
double _e_n_minus1; ///< \f$e_{n-1}\f$.
double _e_n_minus2; ///< \f$e_{n-2}\f$.
@@ -117,12 +113,7 @@ private:
std::min(h_in_range, _rel_h_max * h_n));
}
- double checkSpecificTimeReached(const double h_new) const
- {
- if (_ts_current.current() + h_new;
- return std::max(_rel_h_min * h_n,
- std::min(h_in_range, _rel_h_max * h_n));
- }
+ double checkSpecificTimeReached(const double h_new);
};
/// Create an EvolutionaryPIDcontroller time stepper from the given
diff --git a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
index fa74b69f7ad..cc6024d1ecb 100644
--- a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
@@ -26,12 +26,15 @@ namespace NumLib
FixedTimeStepping::FixedTimeStepping(double t0,
double tn,
const std::vector<double>& vec_all_dt)
- : ITimeStepAlgorithm(t0, computeEnd(t0, tn, vec_all_dt), vec_all_dt)
+ : ITimeStepAlgorithm(t0,
+ computeEnd(t0, tn, vec_all_dt),
+ vec_all_dt,
+ MathLib::VecNormType::INVALID)
{
}
FixedTimeStepping::FixedTimeStepping(double t0, double tn, double dt)
- : ITimeStepAlgorithm(t0, tn, dt)
+ : ITimeStepAlgorithm(t0, tn, dt, MathLib::VecNormType::INVALID)
{
}
diff --git a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h
index 40e0e3cf179..db8d8608a59 100644
--- a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h
+++ b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h
@@ -71,7 +71,6 @@ public:
/// return if current time step is accepted
bool accepted() const override { return true; }
-
private:
/// determine true end time
static double computeEnd(double t_initial, double t_end,
diff --git a/NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h b/NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h
index e43b2cee496..c5c3e56acf2 100644
--- a/NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h
+++ b/NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h
@@ -16,6 +16,8 @@
#include "NumLib/TimeStepping/TimeStep.h"
+#include "MathLib/LinAlg/LinAlg.h" // For MathLib::VecNormType
+
namespace NumLib
{
/**
@@ -25,24 +27,42 @@ namespace NumLib
class ITimeStepAlgorithm
{
public:
- ITimeStepAlgorithm(const double t0, const double t_end)
- : _t_initial(t0), _t_end(t_end), _ts_prev(t0), _ts_current(t0)
+ ITimeStepAlgorithm(const double t0, const double t_end,
+ const MathLib::VecNormType norm_type)
+ : _t_initial(t0),
+ _t_end(t_end),
+ _ts_prev(t0),
+ _ts_current(t0),
+ _norm_type(norm_type)
{
}
- ITimeStepAlgorithm(const double t0, const double t_end, const double dt)
- : _t_initial(t0), _t_end(t_end), _ts_prev(t0), _ts_current(t0),
- _dt_vector(static_cast<std::size_t>(std::ceil((t_end - t0) / dt)), dt)
+ ITimeStepAlgorithm(const double t0, const double t_end, const double dt,
+ const MathLib::VecNormType norm_type)
+ : _t_initial(t0),
+ _t_end(t_end),
+ _ts_prev(t0),
+ _ts_current(t0),
+ _dt_vector(static_cast<std::size_t>(std::ceil((t_end - t0) / dt)),
+ dt),
+ _norm_type(norm_type)
{
}
ITimeStepAlgorithm(const double t0, const double t_end,
- const std::vector<double>& all_step_sizes)
- : _t_initial(t0), _t_end(t_end), _ts_prev(t0), _ts_current(t0),
- _dt_vector(all_step_sizes)
+ const std::vector<double>& all_step_sizes,
+ const MathLib::VecNormType norm_type)
+ : _t_initial(t0),
+ _t_end(t_end),
+ _ts_prev(t0),
+ _ts_current(t0),
+ _dt_vector(all_step_sizes),
+ _norm_type(norm_type)
{
}
+ virtual ~ITimeStepAlgorithm() = default;
+
/// return the beginning of time steps
double begin() const { return _t_initial; }
/// return the end of time steps
@@ -70,7 +90,7 @@ public:
return _dt_vector;
}
- virtual ~ITimeStepAlgorithm() = default;
+ MathLib::VecNormType getSolutionNormType() const { return _norm_type; }
protected:
/// initial time
const double _t_initial;
@@ -84,6 +104,9 @@ protected:
/// a vector of time step sizes
std::vector<double> _dt_vector;
+
+ /// Type of the norm of the solution vector.
+ const MathLib::VecNormType _norm_type;
};
} // NumLib
diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp
index 001f569af1f..a136671ca6c 100644
--- a/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp
@@ -24,9 +24,9 @@ IterationNumberBasedAdaptiveTimeStepping::
double t0, double tn, double min_ts, double max_ts, double initial_ts,
std::vector<std::size_t>& iter_times_vector,
std::vector<double>& multiplier_vector)
- : ITimeStepAlgorithm(t0, tn),
- _iter_times_vector(std::move(iter_times_vector)),
- _multiplier_vector(std::move(multiplier_vector)),
+ : ITimeStepAlgorithm(t0, tn, MathLib::VecNormType::INVALID),
+ _iter_times_vector(iter_times_vector),
+ _multiplier_vector(multiplier_vector),
_min_ts(min_ts),
_max_ts(max_ts),
_initial_ts(initial_ts),
--
GitLab