diff --git a/Documentation/ProjectFile/prj/time_loop/i_time_loop.md b/Documentation/ProjectFile/prj/time_loop/i_time_loop.md
index 576add64edccf5d980a0ff86143dc43a7109a788..aa3cece571f10e2567a06028b73a3626306799d2 100644
--- a/Documentation/ProjectFile/prj/time_loop/i_time_loop.md
+++ b/Documentation/ProjectFile/prj/time_loop/i_time_loop.md
@@ -1 +1 @@
-\ogs_missing_documentation
+Time discretization loop.
diff --git a/Documentation/ProjectFile/prj/time_loop/processes/process/t_time_stepping.md b/Documentation/ProjectFile/prj/time_loop/processes/process/t_time_stepping.md
new file mode 100644
index 0000000000000000000000000000000000000000..2c33d8289fdb551170cc9460b7246af6acdeecf4
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/processes/process/t_time_stepping.md
@@ -0,0 +1 @@
+Define a time stepping method for the modelling of a process.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/c_EvolutionaryPIDcontroller.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/c_EvolutionaryPIDcontroller.md
new file mode 100644
index 0000000000000000000000000000000000000000..730f3fe09a1c9d2661ad05c7eefea7edf7e2d72d
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/c_EvolutionaryPIDcontroller.md
@@ -0,0 +1 @@
+\copydoc NumLib::EvolutionaryPIDcontroller
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_dt_guess.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_dt_guess.md
new file mode 100644
index 0000000000000000000000000000000000000000..96fc64a7031a9a2b48ad4b8860842cdf3b216418
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_dt_guess.md
@@ -0,0 +1 @@
+The initial guess of time step size.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_dt_max.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_dt_max.md
new file mode 100644
index 0000000000000000000000000000000000000000..da4f14e93324fe19019153f66e3cf7926175e4f3
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_dt_max.md
@@ -0,0 +1 @@
+The maximum restriction of time step size.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_dt_min.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_dt_min.md
new file mode 100644
index 0000000000000000000000000000000000000000..c975b5939b85682055c0d50049ff03b2a74ade5a
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_dt_min.md
@@ -0,0 +1 @@
+The minimum restriction of time step size.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_rel_dt_max.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_rel_dt_max.md
new file mode 100644
index 0000000000000000000000000000000000000000..b83f827b1b9a03f6a3fe474d03e904c3202ac72a
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_rel_dt_max.md
@@ -0,0 +1 @@
+The maximum restriction of time step size ratio against the previous time step.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_rel_dt_min.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_rel_dt_min.md
new file mode 100644
index 0000000000000000000000000000000000000000..5d846abd2b589b661dfa5d19d20fd3b3ce40767b
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_rel_dt_min.md
@@ -0,0 +1 @@
+The minimum restriction of time step size ratio against the previous time step.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_specific_times.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_specific_times.md
new file mode 100644
index 0000000000000000000000000000000000000000..92fff5fadbdc54f4f17ebfe9d85beb72be5ee700
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_specific_times.md
@@ -0,0 +1,2 @@
+The specified times that must be reached in the time stepping.
+
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_t_end.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_t_end.md
new file mode 100644
index 0000000000000000000000000000000000000000..ce96f332203a067242130456016be3ff2455f3d4
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_t_end.md
@@ -0,0 +1 @@
+The end time.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_t_initial.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_t_initial.md
new file mode 100644
index 0000000000000000000000000000000000000000..7ca53cb25af157b6f0b483ed4128adac11267fc8
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_t_initial.md
@@ -0,0 +1 @@
+The begin time.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_tol.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_tol.md
new file mode 100644
index 0000000000000000000000000000000000000000..804855455c24028e0b0f842e004cd0b71583e33d
--- /dev/null
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/EvolutionaryPIDcontroller/t_tol.md
@@ -0,0 +1 @@
+The tolerance of the relative error that controls the time stepping.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/c_FixedTimeStepping.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/c_FixedTimeStepping.md
index 576add64edccf5d980a0ff86143dc43a7109a788..6c73dd28077ea87e2a77268299c7d816776654b0 100644
--- a/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/c_FixedTimeStepping.md
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/c_FixedTimeStepping.md
@@ -1 +1 @@
-\ogs_missing_documentation
+Defines a fixed time step method.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/t_t_end.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/t_t_end.md
index 576add64edccf5d980a0ff86143dc43a7109a788..ce96f332203a067242130456016be3ff2455f3d4 100644
--- a/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/t_t_end.md
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/t_t_end.md
@@ -1 +1 @@
-\ogs_missing_documentation
+The end time.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/t_t_initial.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/t_t_initial.md
index 576add64edccf5d980a0ff86143dc43a7109a788..7ca53cb25af157b6f0b483ed4128adac11267fc8 100644
--- a/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/t_t_initial.md
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/FixedTimeStepping/t_t_initial.md
@@ -1 +1 @@
-\ogs_missing_documentation
+The begin time.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/SingleStep/c_SingleStep.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/SingleStep/c_SingleStep.md
index 576add64edccf5d980a0ff86143dc43a7109a788..116071e10a23f299e7fa086986f3eac260d1c677 100644
--- a/Documentation/ProjectFile/prj/time_loop/time_stepping/SingleStep/c_SingleStep.md
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/SingleStep/c_SingleStep.md
@@ -1 +1 @@
-\ogs_missing_documentation
+Defines a SingleStep time stepper.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/i_time_stepping.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/i_time_stepping.md
index 576add64edccf5d980a0ff86143dc43a7109a788..3e5fafed0952b5cc42ed90619f990818ab02883a 100644
--- a/Documentation/ProjectFile/prj/time_loop/time_stepping/i_time_stepping.md
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/i_time_stepping.md
@@ -1 +1 @@
-\ogs_missing_documentation
+Defines a time stepper.
diff --git a/Documentation/ProjectFile/prj/time_loop/time_stepping/t_type.md b/Documentation/ProjectFile/prj/time_loop/time_stepping/t_type.md
index 576add64edccf5d980a0ff86143dc43a7109a788..2cac2c7c0deff1b4f42f6fa915daca19e8c51d31 100644
--- a/Documentation/ProjectFile/prj/time_loop/time_stepping/t_type.md
+++ b/Documentation/ProjectFile/prj/time_loop/time_stepping/t_type.md
@@ -1 +1 @@
-\ogs_missing_documentation
+The type of a time stepping method.
diff --git a/Documentation/bibliography.bib b/Documentation/bibliography.bib
index c101c7cb1f04389646fda183d0e4026aae003ad9..696f7b6de36ad5615b4aff1f1c1349a901dc0166 100644
--- a/Documentation/bibliography.bib
+++ b/Documentation/bibliography.bib
@@ -110,3 +110,13 @@
keywords = "Mont Terri ",
abstract = "Abstract The constitutive modeling of shales is an important topic in the geomechanics community as it is often encountered in advanced applications such as nuclear waste storage, \{CO2\} sequestration, and unconventional oil and gas exploitation. The goal of this work is to describe, within a unique plastic-damage framework, the full mechanical behavior of shales such as pre-peak hardening plasticity, non-linearity in the onset of inelastic strains, dilatancy, post-peak softening, and degradation of the elastic parameters. The model validation against three sets of triaxial experimental results on shale demonstrates its capability to reproduce the main mechanical characteristics. "
}
+
+@article{ahmed2015adaptive,
+ title={Adaptive time step control for higher order variational time discretizations applied to convection--diffusion--reaction equations},
+ author={Ahmed, Naveed and John, Volker},
+ journal={Computer Methods in Applied Mechanics and Engineering},
+ volume={285},
+ pages={83--101},
+ year={2015},
+ publisher={Elsevier}
+}
diff --git a/NumLib/ODESolver/ConvergenceCriterion.h b/NumLib/ODESolver/ConvergenceCriterion.h
index 154b23b4de4c4d7514a4507d1ad2670d45495528..59babdd8450e4a88623e1a2c5d183b7372f56b4d 100644
--- a/NumLib/ODESolver/ConvergenceCriterion.h
+++ b/NumLib/ODESolver/ConvergenceCriterion.h
@@ -10,6 +10,7 @@
#pragma once
#include <memory>
+#include "MathLib/LinAlg/LinAlg.h" // For MathLib::VecNormType
#include "NumLib/NumericsConfig.h"
namespace BaseLib {
@@ -27,6 +28,11 @@ namespace NumLib
class ConvergenceCriterion
{
public:
+ ConvergenceCriterion(const MathLib::VecNormType norm_type)
+ : _norm_type(norm_type)
+ {
+ }
+
//! Tells if the change of the solution between iterations is checked.
//!
//! \remark
@@ -76,11 +82,14 @@ public:
//! Tell if the convergence criterion is satisfied.
virtual bool isSatisfied() const { return _satisfied; };
+ MathLib::VecNormType getVectorNormType () const { return _norm_type; }
+
virtual ~ConvergenceCriterion() = default;
protected:
bool _satisfied = true;
bool _is_first_iteration = true;
+ const MathLib::VecNormType _norm_type;
};
//! Creates a convergence criterion from the given configuration.
diff --git a/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp b/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
index 88de56434264079d10a16d5fb528cb703bf0a6ee..1bbece408622b75a3d9f16963e343b707dfc350c 100644
--- a/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
@@ -18,10 +18,9 @@ namespace NumLib
ConvergenceCriterionDeltaX::ConvergenceCriterionDeltaX(
boost::optional<double>&& absolute_tolerance,
boost::optional<double>&& relative_tolerance,
- MathLib::VecNormType norm_type)
- : _abstol(std::move(absolute_tolerance)),
- _reltol(std::move(relative_tolerance)),
- _norm_type(norm_type)
+ const MathLib::VecNormType norm_type)
+ : ConvergenceCriterion(norm_type), _abstol(std::move(absolute_tolerance)),
+ _reltol(std::move(relative_tolerance))
{
if ((!_abstol) && (!_reltol))
OGS_FATAL(
diff --git a/NumLib/ODESolver/ConvergenceCriterionDeltaX.h b/NumLib/ODESolver/ConvergenceCriterionDeltaX.h
index 75ede7ec830da063c8cb22d4b466c8e17e12a87f..71715ca1950c2e25469e5ddace7d37d9a6ca1efd 100644
--- a/NumLib/ODESolver/ConvergenceCriterionDeltaX.h
+++ b/NumLib/ODESolver/ConvergenceCriterionDeltaX.h
@@ -27,7 +27,7 @@ public:
ConvergenceCriterionDeltaX(
boost::optional<double>&& absolute_tolerance,
boost::optional<double>&& relative_tolerance,
- MathLib::VecNormType norm_type);
+ const MathLib::VecNormType norm_type);
bool hasDeltaXCheck() const override { return true; }
bool hasResidualCheck() const override { return false; }
@@ -40,7 +40,6 @@ public:
private:
const boost::optional<double> _abstol;
const boost::optional<double> _reltol;
- const MathLib::VecNormType _norm_type;
};
std::unique_ptr<ConvergenceCriterionDeltaX> createConvergenceCriterionDeltaX(
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponent.h b/NumLib/ODESolver/ConvergenceCriterionPerComponent.h
index 0ff42fdebf1c5c6f302110ce0df9158242e64f90..8da246ad60d1856e59dec0f516a5061665b79194 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponent.h
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponent.h
@@ -11,6 +11,8 @@
#include "ConvergenceCriterion.h"
+#include "MathLib/LinAlg/LinAlg.h" // For MathLib::VecNormType
+
namespace MeshLib
{
class Mesh;
@@ -26,6 +28,9 @@ class LocalToGlobalIndexMap;
class ConvergenceCriterionPerComponent : public ConvergenceCriterion
{
public:
+ ConvergenceCriterionPerComponent(const MathLib::VecNormType norm_type)
+ :ConvergenceCriterion(norm_type) {}
+
//! Sets the d.o.f. table used to extract data for a specific component.
virtual void setDOFTable(NumLib::LocalToGlobalIndexMap const& dof_table,
MeshLib::Mesh const& mesh) = 0;
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
index b9faaeed41ed0f4ae6e37633a7185bee5e30b930..f5ef433565e803e7d0f9a1ee6141cdbe1ee3f51e 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
@@ -19,10 +19,10 @@ namespace NumLib
ConvergenceCriterionPerComponentDeltaX::ConvergenceCriterionPerComponentDeltaX(
std::vector<double>&& absolute_tolerances,
std::vector<double>&& relative_tolerances,
- MathLib::VecNormType norm_type)
- : _abstols(std::move(absolute_tolerances)),
- _reltols(std::move(relative_tolerances)),
- _norm_type(norm_type)
+ const MathLib::VecNormType norm_type)
+ : ConvergenceCriterionPerComponent(norm_type),
+ _abstols(std::move(absolute_tolerances)),
+ _reltols(std::move(relative_tolerances))
{
if (_abstols.size() != _reltols.size())
OGS_FATAL(
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.h b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.h
index a35ea708527e8c83dc54775637b05e33692b2180..3622a32e286d9b979749181d94a014e4b2eebdd6 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.h
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.h
@@ -29,7 +29,7 @@ public:
ConvergenceCriterionPerComponentDeltaX(
std::vector<double>&& absolute_tolerances,
std::vector<double>&& relative_tolerances,
- MathLib::VecNormType norm_type);
+ const MathLib::VecNormType norm_type);
bool hasDeltaXCheck() const override { return true; }
bool hasResidualCheck() const override { return false; }
@@ -46,7 +46,6 @@ public:
private:
const std::vector<double> _abstols;
const std::vector<double> _reltols;
- const MathLib::VecNormType _norm_type;
LocalToGlobalIndexMap const* _dof_table = nullptr;
MeshLib::Mesh const* _mesh = nullptr;
};
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
index 29600a819952605a9669b730bef056398eaa357e..acd515f7ccab2dc91f921428ac29f554d13ff568 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
@@ -20,10 +20,10 @@ ConvergenceCriterionPerComponentResidual::
ConvergenceCriterionPerComponentResidual(
std::vector<double>&& absolute_tolerances,
std::vector<double>&& relative_tolerances,
- MathLib::VecNormType norm_type)
- : _abstols(std::move(absolute_tolerances)),
+ const MathLib::VecNormType norm_type)
+ : ConvergenceCriterionPerComponent(norm_type),
+ _abstols(std::move(absolute_tolerances)),
_reltols(std::move(relative_tolerances)),
- _norm_type(norm_type),
_residual_norms_0(_abstols.size())
{
if (_abstols.size() != _reltols.size())
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.h b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.h
index f54c1a6a932e04919a0d4f4dcc8f5d5f7dd755a2..743a6fc5da5b63e9428de7d94cd63d312e721257 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.h
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.h
@@ -30,7 +30,7 @@ public:
ConvergenceCriterionPerComponentResidual(
std::vector<double>&& absolute_tolerances,
std::vector<double>&& relative_tolerances,
- MathLib::VecNormType norm_type);
+ const MathLib::VecNormType norm_type);
bool hasDeltaXCheck() const override { return true; }
bool hasResidualCheck() const override { return true; }
@@ -47,7 +47,6 @@ public:
private:
const std::vector<double> _abstols;
const std::vector<double> _reltols;
- const MathLib::VecNormType _norm_type;
LocalToGlobalIndexMap const* _dof_table = nullptr;
MeshLib::Mesh const* _mesh = nullptr;
std::vector<double> _residual_norms_0;
diff --git a/NumLib/ODESolver/ConvergenceCriterionResidual.cpp b/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
index dc437fc66fe3ecf6dbe14e6e31dfd42f1332973f..e160ac6a4bf7a4385d711f0056fca15d41ddf2bb 100644
--- a/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
@@ -18,10 +18,9 @@ namespace NumLib
ConvergenceCriterionResidual::ConvergenceCriterionResidual(
boost::optional<double>&& absolute_tolerance,
boost::optional<double>&& relative_tolerance,
- MathLib::VecNormType norm_type)
- : _abstol(std::move(absolute_tolerance)),
- _reltol(std::move(relative_tolerance)),
- _norm_type(norm_type)
+ const MathLib::VecNormType norm_type)
+ : ConvergenceCriterion(norm_type), _abstol(std::move(absolute_tolerance)),
+ _reltol(std::move(relative_tolerance))
{
if ((!_abstol) && (!_reltol))
OGS_FATAL(
@@ -29,8 +28,8 @@ ConvergenceCriterionResidual::ConvergenceCriterionResidual(
"specified.");
}
-void ConvergenceCriterionResidual::checkDeltaX(const GlobalVector& minus_delta_x,
- GlobalVector const& x)
+void ConvergenceCriterionResidual::checkDeltaX(
+ const GlobalVector& minus_delta_x, GlobalVector const& x)
{
auto error_dx = MathLib::LinAlg::norm(minus_delta_x, _norm_type);
auto norm_x = MathLib::LinAlg::norm(x, _norm_type);
@@ -43,21 +42,38 @@ void ConvergenceCriterionResidual::checkResidual(const GlobalVector& residual)
{
auto norm_res = MathLib::LinAlg::norm(residual, _norm_type);
- if (_is_first_iteration) {
+ if (_is_first_iteration)
+ {
INFO("Convergence criterion: |r0|=%.4e", norm_res);
_residual_norm_0 = norm_res;
- } else {
- INFO("Convergence criterion: |r|=%.4e |r0|=%.4e |r|/|r0|=%.4e",
- norm_res, _residual_norm_0, norm_res / _residual_norm_0);
+ }
+ else
+ {
+ _residual_norm_0 =
+ (_residual_norm_0 < std::numeric_limits<double>::epsilon())
+ ? norm_res
+ : _residual_norm_0;
+ if (_residual_norm_0 < std::numeric_limits<double>::epsilon())
+ {
+ INFO("Convergence criterion: |r|=%.4e |r0|=%.4e",
+ norm_res, _residual_norm_0);
+ }
+ else
+ {
+ INFO("Convergence criterion: |r|=%.4e |r0|=%.4e |r|/|r0|=%.4e",
+ norm_res, _residual_norm_0, norm_res / _residual_norm_0);
+ }
}
bool satisfied_abs = false;
bool satisfied_rel = false;
- if (_abstol) {
+ if (_abstol)
+ {
satisfied_abs = norm_res < *_abstol;
}
- if (_reltol && !_is_first_iteration) {
+ if (_reltol && !_is_first_iteration)
+ {
satisfied_rel =
checkRelativeTolerance(*_reltol, norm_res, _residual_norm_0);
}
diff --git a/NumLib/ODESolver/ConvergenceCriterionResidual.h b/NumLib/ODESolver/ConvergenceCriterionResidual.h
index b598cfbc83a3920dcbcdf4df61ba4473b5d57b41..b63af27e9b9af5010a7e7f297341b567b4db28d0 100644
--- a/NumLib/ODESolver/ConvergenceCriterionResidual.h
+++ b/NumLib/ODESolver/ConvergenceCriterionResidual.h
@@ -27,7 +27,7 @@ public:
ConvergenceCriterionResidual(
boost::optional<double>&& absolute_tolerance,
boost::optional<double>&& relative_tolerance,
- MathLib::VecNormType norm_type);
+ const MathLib::VecNormType norm_type);
bool hasDeltaXCheck() const override { return true; }
bool hasResidualCheck() const override { return true; }
@@ -41,7 +41,6 @@ public:
private:
const boost::optional<double> _abstol;
const boost::optional<double> _reltol;
- const MathLib::VecNormType _norm_type;
double _residual_norm_0;
};
diff --git a/NumLib/ODESolver/TimeDiscretization.cpp b/NumLib/ODESolver/TimeDiscretization.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e7cf70c71188bf60c952ccd3af9c1d07cf47ae3b
--- /dev/null
+++ b/NumLib/ODESolver/TimeDiscretization.cpp
@@ -0,0 +1,138 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file TimeDiscretization.cpp
+ * Created on March 31, 2017, 10:30 AM
+ */
+
+#include "TimeDiscretization.h"
+
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
+
+namespace NumLib
+{
+double TimeDiscretization::computeRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x,
+ GlobalVector const& x_old,
+ MathLib::VecNormType norm_type)
+{
+ if (norm_type == MathLib::VecNormType::INVALID)
+ {
+ OGS_FATAL("An invalid norm type has been passed");
+ }
+
+ if (!_dx)
+ _dx = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
+
+ auto& dx = *_dx;
+ MathLib::LinAlg::copy(x, dx); // copy x to dx.
+
+ // dx = x - x_old --> x - dx --> dx
+ MathLib::LinAlg::axpy(dx, -1.0, x_old);
+ const double norm_dx = MathLib::LinAlg::norm(dx, norm_type);
+
+ const double norm_x = MathLib::LinAlg::norm(x, norm_type);
+ if (norm_x > std::numeric_limits<double>::epsilon())
+ return norm_dx / norm_x;
+
+ // Both of norm_x and norm_dx are close to zero
+ if (norm_dx < std::numeric_limits<double>::epsilon())
+ return 1.0;
+
+ // Only norm_x is close to zero
+ return norm_dx / std::numeric_limits<double>::epsilon();
+}
+
+double BackwardEuler::getRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x, MathLib::VecNormType norm_type)
+{
+ return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
+}
+
+double ForwardEuler::getRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x, MathLib::VecNormType norm_type)
+{
+ return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
+}
+
+double CrankNicolson::getRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x, MathLib::VecNormType norm_type)
+{
+ return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
+}
+
+double BackwardDifferentiationFormula::getRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x, MathLib::VecNormType norm_type)
+{
+ return computeRelativeChangeFromPreviousTimestep(
+ x, *_xs_old[_offset], norm_type);
+}
+
+void BackwardDifferentiationFormula::pushState(const double,
+ GlobalVector const& x,
+ InternalMatrixStorage const&)
+{
+ namespace LinAlg = MathLib::LinAlg;
+ // TODO use boost circular buffer?
+
+ // until _xs_old is filled, lower-order BDF formulas are used.
+ if (_xs_old.size() < _num_steps)
+ {
+ _xs_old.push_back(&NumLib::GlobalVectorProvider::provider.getVector(x));
+ }
+ else
+ {
+ LinAlg::copy(x, *_xs_old[_offset]);
+ _offset =
+ (_offset + 1) % _num_steps; // treat _xs_old as a circular buffer
+ }
+}
+
+namespace detail
+{
+//! Coefficients used in the backward differentiation formulas.
+const double BDF_Coeffs[6][7] = {
+ // leftmost column: weight of the solution at the new timestep
+ // signs of columns > 1 are flipped compared to standard BDF tableaus
+ {1.0, 1.0},
+ {1.5, 2.0, -0.5},
+ {11.0 / 6.0, 3.0, -1.5, 1.0 / 3.0},
+ {25.0 / 12.0, 4.0, -3.0, 4.0 / 3.0, -0.25},
+ {137.0 / 60.0, 5.0, -5.0, 10.0 / 3.0, -1.25, 1.0 / 5.0},
+ {147.0 / 60.0, 6.0, -7.5, 20.0 / 3.0, -3.75, 6.0 / 5.0, -1.0 / 6.0}
+ // coefficient of (for BDF(6), the oldest state, x_n, is always rightmost)
+ // x_+6, x_+5, x_+4, x_+3, x_+2, x_+1, x_n
+};
+}
+
+double BackwardDifferentiationFormula::getNewXWeight() const
+{
+ auto const k = eff_num_steps();
+ return detail::BDF_Coeffs[k - 1][0] / _delta_t;
+}
+
+void BackwardDifferentiationFormula::getWeightedOldX(GlobalVector& y) const
+{
+ namespace LinAlg = MathLib::LinAlg;
+
+ auto const k = eff_num_steps();
+ auto const* const BDFk = detail::BDF_Coeffs[k - 1];
+
+ // compute linear combination \sum_{i=0}^{k-1} BDFk_{k-i} \cdot x_{n+i}
+ LinAlg::copy(*_xs_old[_offset], y); // _xs_old[offset] = x_n
+ LinAlg::scale(y, BDFk[k]);
+
+ for (unsigned i = 1; i < k; ++i)
+ {
+ auto const off = (_offset + i) % k;
+ LinAlg::axpy(y, BDFk[k - i], *_xs_old[off]);
+ }
+
+ LinAlg::scale(y, 1.0 / _delta_t);
+}
+
+} // end of namespace NumLib
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index 0db170d920e23a839d14025ff862284e4f07e7d6..ede512492ec2320005e834958b5e2d148584305a 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -114,9 +114,20 @@ public:
class TimeDiscretization
{
public:
+ TimeDiscretization() = default;
+
//! Sets the initial condition.
virtual void setInitialState(const double t0, GlobalVector const& x0) = 0;
+ /*! Get the relative change of solutions between two successive time steps
+ * by \f$ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| \f$.
+ *
+ * \param x The solution at the current timestep.
+ * \param norm_type The type of global vector norm.
+ */
+ virtual double getRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x, MathLib::VecNormType norm_type) = 0;
+
/*! Indicate that the current timestep is done and that you will proceed to
* the next one.
*
@@ -200,6 +211,25 @@ public:
*/
virtual bool needsPreload() const { return false; }
//! @}
+
+protected:
+ std::unique_ptr<GlobalVector> _dx; ///< Used to store \f$ u_{n+1}-u_{n}\f$.
+
+ /**
+ * Compute and return the relative change of solutions between two successive
+ * time steps by \f$ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| \f$.
+ *
+ * @param x The current solution
+ * @param x_old The previous solution
+ * @param norm_type The norm type of global vector
+ * @return \f$ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| \f$.
+ *
+ * @warning the value of x_old is changed to x - x_old after this computation.
+ */
+ double computeRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x,
+ GlobalVector const& x_old,
+ MathLib::VecNormType norm_type);
};
//! Backward Euler scheme.
@@ -222,6 +252,9 @@ public:
MathLib::LinAlg::copy(x0, _x_old);
}
+ double getRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x, MathLib::VecNormType norm_type) override;
+
void pushState(const double /*t*/, GlobalVector const& x,
InternalMatrixStorage const&) override
{
@@ -256,7 +289,7 @@ class ForwardEuler final : public TimeDiscretization
{
public:
ForwardEuler()
- : _x_old(NumLib::GlobalVectorProvider::provider.getVector())
+ : _x_old(NumLib::GlobalVectorProvider::provider.getVector())
{
}
@@ -272,6 +305,9 @@ public:
MathLib::LinAlg::copy(x0, _x_old);
}
+ double getRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x, MathLib::VecNormType norm_type) override;
+
void pushState(const double /*t*/, GlobalVector const& x,
InternalMatrixStorage const&) override
{
@@ -346,6 +382,9 @@ public:
MathLib::LinAlg::copy(x0, _x_old);
}
+ double getRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x, MathLib::VecNormType norm_type) override;
+
void pushState(const double, GlobalVector const& x,
InternalMatrixStorage const& strg) override
{
@@ -382,23 +421,6 @@ private:
GlobalVector& _x_old; //!< the solution from the preceding timestep
};
-namespace detail
-{
-//! Coefficients used in the backward differentiation formulas.
-const double BDF_Coeffs[6][7] = {
- // leftmost column: weight of the solution at the new timestep
- // signs of columns > 1 are flipped compared to standard BDF tableaus
- {1.0, 1.0},
- {1.5, 2.0, -0.5},
- {11.0 / 6.0, 3.0, -1.5, 1.0 / 3.0},
- {25.0 / 12.0, 4.0, -3.0, 4.0 / 3.0, -0.25},
- {137.0 / 60.0, 5.0, -5.0, 10.0 / 3.0, -1.25, 1.0 / 5.0},
- {147.0 / 60.0, 6.0, -7.5, 20.0 / 3.0, -3.75, 6.0 / 5.0, -1.0 / 6.0}
- // coefficient of (for BDF(6), the oldest state, x_n, is always rightmost)
- // x_+6, x_+5, x_+4, x_+3, x_+2, x_+1, x_n
-};
-}
-
//! Backward differentiation formula.
class BackwardDifferentiationFormula final : public TimeDiscretization
{
@@ -416,7 +438,7 @@ public:
* the first timesteps.
*/
explicit BackwardDifferentiationFormula(const unsigned num_steps)
- : _num_steps(num_steps)
+ : _num_steps(num_steps)
{
assert(1 <= num_steps && num_steps <= 6);
_xs_old.reserve(num_steps);
@@ -435,25 +457,11 @@ public:
&NumLib::GlobalVectorProvider::provider.getVector(x0));
}
+ double getRelativeChangeFromPreviousTimestep(
+ GlobalVector const& x, MathLib::VecNormType norm_type) override;
+
void pushState(const double, GlobalVector const& x,
- InternalMatrixStorage const&) override
- {
- namespace LinAlg = MathLib::LinAlg;
- // TODO use boost cirular buffer?
-
- // until _xs_old is filled, lower-order BDF formulas are used.
- if (_xs_old.size() < _num_steps)
- {
- _xs_old.push_back(
- &NumLib::GlobalVectorProvider::provider.getVector(x));
- }
- else
- {
- LinAlg::copy(x, *_xs_old[_offset]);
- _offset = (_offset + 1) %
- _num_steps; // treat _xs_old as a circular buffer
- }
- }
+ InternalMatrixStorage const&) override;
void nextTimestep(const double t, const double delta_t) override
{
@@ -462,31 +470,10 @@ public:
}
double getCurrentTime() const override { return _t; }
- double getNewXWeight() const override
- {
- auto const k = eff_num_steps();
- return detail::BDF_Coeffs[k - 1][0] / _delta_t;
- }
-
- void getWeightedOldX(GlobalVector& y) const override
- {
- namespace LinAlg = MathLib::LinAlg;
-
- auto const k = eff_num_steps();
- auto const* const BDFk = detail::BDF_Coeffs[k - 1];
- // compute linear combination \sum_{i=0}^{k-1} BDFk_{k-i} \cdot x_{n+i}
- LinAlg::copy(*_xs_old[_offset], y); // _xs_old[offset] = x_n
- LinAlg::scale(y, BDFk[k]);
+ double getNewXWeight() const override;
- for (unsigned i = 1; i < k; ++i)
- {
- auto const off = (_offset + i) % k;
- LinAlg::axpy(y, BDFk[k - i], *_xs_old[off]);
- }
-
- LinAlg::scale(y, 1.0 / _delta_t);
- }
+ void getWeightedOldX(GlobalVector& y) const override;
private:
std::size_t eff_num_steps() const { return _xs_old.size(); }
diff --git a/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp b/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..84b63208c16ec97563017fcaeec4d48e414b0abc
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp
@@ -0,0 +1,66 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateEvolutionaryPIDcontroller.cpp
+ * Created on June 26, 2017, 4:43 PM
+ */
+
+#include "CreateEvolutionaryPIDcontroller.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/makeVectorUnique.h"
+
+#include "EvolutionaryPIDcontroller.h"
+#include "TimeStepAlgorithm.h"
+
+namespace NumLib
+{
+class TimeStepAlgorithm;
+std::unique_ptr<TimeStepAlgorithm> createEvolutionaryPIDcontroller(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{prj__time_loop__time_stepping__type}
+ config.checkConfigParameter("type", "EvolutionaryPIDcontroller");
+
+ //! \ogs_file_param{prj__time_loop__time_stepping__EvolutionaryPIDcontroller__t_initial}
+ auto const t0 = config.getConfigParameter<double>("t_initial");
+ //! \ogs_file_param{prj__time_loop__time_stepping__EvolutionaryPIDcontroller__t_end}
+ auto const t_end = config.getConfigParameter<double>("t_end");
+ //! \ogs_file_param{prj__time_loop__time_stepping__EvolutionaryPIDcontroller__dt_guess}
+ auto const h0 = config.getConfigParameter<double>("dt_guess");
+
+ //! \ogs_file_param{prj__time_loop__time_stepping__EvolutionaryPIDcontroller__dt_min}
+ auto const h_min = config.getConfigParameter<double>("dt_min");
+ //! \ogs_file_param{prj__time_loop__time_stepping__EvolutionaryPIDcontroller__dt_max}
+ auto const h_max = config.getConfigParameter<double>("dt_max");
+ //! \ogs_file_param{prj__time_loop__time_stepping__EvolutionaryPIDcontroller__rel_dt_min}
+ auto const rel_h_min = config.getConfigParameter<double>("rel_dt_min");
+ //! \ogs_file_param{prj__time_loop__time_stepping__EvolutionaryPIDcontroller__rel_dt_max}
+ auto const rel_h_max = config.getConfigParameter<double>("rel_dt_max");
+
+ auto specific_times =
+ //! \ogs_file_param{prj__time_loop__time_stepping__EvolutionaryPIDcontroller__specific_times}
+ config.getConfigParameter<std::vector<double>>(
+ "specific_times", std::vector<double>{});
+ if (!specific_times.empty())
+ {
+ // Sort in descending order.
+ std::sort(specific_times.begin(),
+ specific_times.end(),
+ std::greater<double>());
+ // Remove possible duplicated elements.
+ BaseLib::makeVectorUnique(specific_times);
+ }
+
+ //! \ogs_file_param{prj__time_loop__time_stepping__EvolutionaryPIDcontroller__tol}
+ auto const tol = config.getConfigParameter<double>("tol");
+
+ return std::make_unique<EvolutionaryPIDcontroller>(
+ t0, t_end, h0, h_min, h_max, rel_h_min, rel_h_max,
+ std::move(specific_times), tol);
+}
+} // end of namespace NumLib
diff --git a/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.h b/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.h
new file mode 100644
index 0000000000000000000000000000000000000000..84c35364026d0440df1925d00daaac17012c5ca2
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.h
@@ -0,0 +1,29 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateEvolutionaryPIDcontroller.h
+ * Created on June 26, 2017, 4:43 PM
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace NumLib
+{
+class TimeStepAlgorithm;
+
+/// Create an EvolutionaryPIDcontroller time stepper from the given
+/// configuration
+std::unique_ptr<TimeStepAlgorithm> createEvolutionaryPIDcontroller(
+ BaseLib::ConfigTree const& config);
+} // end of namespace NumLib
diff --git a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..edc2085f479ec4d408b3cce8f888cfa0d6dcd58e
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp
@@ -0,0 +1,82 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateFixedTimeStepping.cpp
+ * Created on June 26, 2017, 5:03 PM
+ */
+
+#include "CreateFixedTimeStepping.h"
+#include <string>
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/Error.h"
+
+#include "FixedTimeStepping.h"
+#include "TimeStepAlgorithm.h"
+
+namespace NumLib
+{
+class TimeStepAlgorithm;
+std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{prj__time_loop__time_stepping__type}
+ config.checkConfigParameter("type", "FixedTimeStepping");
+
+ //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__t_initial}
+ auto const t_initial = config.getConfigParameter<double>("t_initial");
+ //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__t_end}
+ auto const t_end = config.getConfigParameter<double>("t_end");
+ //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__timesteps}
+ auto const delta_ts = config.getConfigSubtree("timesteps");
+
+ std::vector<double> timesteps;
+ double t_curr = t_initial;
+ double delta_t = 0.0;
+
+ // TODO: consider adding call "listNonEmpty" to config tree
+ //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__timesteps__pair}
+ auto const range = delta_ts.getConfigSubtreeList("pair");
+ if (range.begin() == range.end())
+ {
+ OGS_FATAL("no timesteps have been given");
+ }
+ for (auto const pair : range)
+ {
+ //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__timesteps__pair__repeat}
+ auto const repeat = pair.getConfigParameter<std::size_t>("repeat");
+ //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__timesteps__pair__delta_t}
+ delta_t = pair.getConfigParameter<double>("delta_t");
+
+ if (repeat == 0)
+ {
+ OGS_FATAL("<repeat> is zero.");
+ }
+ if (delta_t <= 0.0)
+ {
+ OGS_FATAL("timestep <delta_t> is <= 0.0.");
+ }
+
+ if (t_curr <= t_end)
+ {
+ timesteps.resize(timesteps.size() + repeat, delta_t);
+
+ t_curr += repeat * delta_t;
+ }
+ }
+
+ // append last delta_t until t_end is reached
+ if (t_curr <= t_end)
+ {
+ auto const repeat =
+ static_cast<std::size_t>(std::ceil((t_end - t_curr) / delta_t));
+ timesteps.resize(timesteps.size() + repeat, delta_t);
+ }
+
+ return std::make_unique<FixedTimeStepping>(t_initial, t_end, timesteps);
+}
+} // end of namespace NumLib
diff --git a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h
new file mode 100644
index 0000000000000000000000000000000000000000..ff7e066a04390579652187a4f41ab99db3dd9955
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h
@@ -0,0 +1,29 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateFixedTimeStepping.h
+ * Created on June 26, 2017, 5:02 PM
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace NumLib
+{
+class TimeStepAlgorithm;
+
+/// Create a FixedTimeStepping time stepper from the given
+/// configuration
+std::unique_ptr<TimeStepAlgorithm> createFixedTimeStepping(
+ BaseLib::ConfigTree const& config);
+} // end of namespace NumLib
diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f3317e6e1e11a564b3904e3312621f64a18f67f8
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
@@ -0,0 +1,112 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file EvolutionaryPIDcontroller.cpp
+ * Created on March 31, 2017, 4:13 PM
+ */
+
+#include <algorithm>
+#include <functional>
+#include <limits>
+#include <vector>
+
+#include "EvolutionaryPIDcontroller.h"
+
+namespace NumLib
+{
+bool EvolutionaryPIDcontroller::next(const double solution_error)
+{
+ const double e_n = solution_error;
+ const double zero_threshlod = std::numeric_limits<double>::epsilon();
+ // step rejected.
+ if (e_n > _tol) // e_n < TOL
+ {
+ _is_accepted = false;
+
+ 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 = checkSpecificTimeReached(h_new);
+
+ _ts_current = _ts_prev;
+ _ts_current += h_new;
+ return false;
+ }
+
+ // step accepted.
+ _is_accepted = true;
+
+ if (_ts_current.steps() == 0)
+ {
+ _ts_prev = _ts_current;
+ _ts_current += _h0;
+ _e_n_minus1 = e_n;
+
+ _dt_vector.push_back(_h0);
+ }
+ else
+ {
+ const double h_n = _ts_current.dt();
+ double h_new = h_n;
+
+ if (e_n > zero_threshlod)
+ {
+ if (_e_n_minus1 > zero_threshlod)
+ {
+ if (_e_n_minus2 > zero_threshlod)
+ {
+ 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, h_n);
+ h_new = checkSpecificTimeReached(h_new);
+ _dt_vector.push_back(h_new);
+
+ _ts_prev = _ts_current;
+ _ts_current += h_new;
+
+ _e_n_minus2 = _e_n_minus1;
+ _e_n_minus1 = e_n;
+ }
+
+ 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_threshold = std::numeric_limits<double>::epsilon();
+ if ((specific_time > _ts_current.current()) &&
+ (_ts_current.current() + h_new - specific_time > zero_threshold))
+ {
+ _specific_times.pop_back();
+ return specific_time - _ts_current.current();
+ }
+
+ return h_new;
+}
+} // end of namespace NumLib
diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h
new file mode 100644
index 0000000000000000000000000000000000000000..5e5226667f64bc5624bbbd2968be39d5247cd14a
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h
@@ -0,0 +1,127 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file EvolutionaryPIDcontroller.h
+ * Created on March 31, 2017, 4:13 PM
+ */
+
+#pragma once
+
+#include <vector>
+#include <memory>
+
+#include "TimeStepAlgorithm.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace NumLib
+{
+/**
+ * This class gives an adaptive algorithm whose time step control is
+ * evolutionary PID controller. With an definition of relative solution change
+ * \f$e_n=\frac{\|u^{n+1}-u^n\|}{\|u^{n+1}\|}\f$, the algorithm gives a time
+ * step size estimation as
+ * \f[
+ * 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}
+ * \f]
+ * where \f$k_P=0.075\f$, \f$k_I=0.175\f$, \f$k_D=0.01\f$ are empirical PID
+ * parameters.
+ *
+ * In the computation, \f$ e_n\f$ is calculated firstly. If \f$e_n>TOL\f$, the
+ * current time step is rejected and repeated with a new time step size of
+ * \f$h=\frac{TOL}{e_n} h_n\f$.
+ *
+ * Limits of the time step size are given as
+ * \f[
+ * h_{\mbox{min}} \leq h_{n+1} \leq h_{\mbox{max}},
+ * l \leq \frac{h_{n+1}}{h_n} \leq L
+ * \f]
+ *
+ * Similar algorithm can be found in \cite ahmed2015adaptive .
+ */
+class EvolutionaryPIDcontroller final : public TimeStepAlgorithm
+{
+public:
+ 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 std::vector<double>&& specific_times,
+ const double tol)
+ : 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),
+ _specific_times(std::move(specific_times)),
+ _tol(tol),
+ _e_n_minus1(0.),
+ _e_n_minus2(0.),
+ _is_accepted(true)
+ {
+ }
+
+ /**
+ * move to the next time step
+ * @param solution_error \f$e_n\f$, solution error between two successive
+ time steps.
+ * @return true if the next step exists
+ */
+ bool next(const double solution_error) override;
+
+ /// return if current time step is accepted
+ bool accepted() const override { return _is_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
+ const double _kD = 0.01; ///< Parameter. \see EvolutionaryPIDcontroller
+
+ const double _h0; ///< initial time step size.
+ const double _h_min; ///< minimum step size.
+ const double _h_max; ///< maximum step size.
+
+ /// \f$l\f$ in \f$ h_{\mbox{min}} \leq h_{n+1} \leq h_{\mbox{max}},\f$
+ const double _rel_h_min;
+ /// \f$L\f$ in \f$ h_{\mbox{min}} \leq h_{n+1} \leq h_{\mbox{max}},\f$
+ const double _rel_h_max;
+
+ // Given times that steps have to reach.
+ std::vector<double> _specific_times;
+
+ const double _tol;
+
+ double _e_n_minus1; ///< \f$e_{n-1}\f$.
+ double _e_n_minus2; ///< \f$e_{n-2}\f$.
+
+ 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 checkSpecificTimeReached(const double h_new);
+};
+
+/// Create an EvolutionaryPIDcontroller time stepper from the given
+/// configuration
+std::unique_ptr<TimeStepAlgorithm> createEvolutionaryPIDcontroller(
+ BaseLib::ConfigTree const& config);
+
+} // end of namespace NumLib
diff --git a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
index 7ee39bc34fcf0d0ec5c275984514edf94f58a404..63fb6da19973272e151e3460eb632286906ee2fc 100644
--- a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
@@ -12,92 +12,30 @@
#include "FixedTimeStepping.h"
#include <algorithm>
-#include <cmath>
#include <numeric>
#include <limits>
#include <cassert>
-#include <logog/include/logog.hpp>
-
-#include "BaseLib/ConfigTree.h"
-#include "BaseLib/Error.h"
-
namespace NumLib
{
-
-FixedTimeStepping::FixedTimeStepping(double t0, double tn, const std::vector<double> &vec_all_dt)
-: _t_initial(t0), _t_end(computeEnd(t0, tn, vec_all_dt)), _dt_vector(vec_all_dt), _ts_prev(t0), _ts_current(t0)
-{}
-
-FixedTimeStepping::FixedTimeStepping(double t0, double tn, double dt)
-: _t_initial(t0), _t_end(tn), _dt_vector(static_cast<std::size_t>(std::ceil((tn-t0)/dt)), dt), _ts_prev(t0), _ts_current(t0)
-{}
-
-std::unique_ptr<ITimeStepAlgorithm>
-FixedTimeStepping::newInstance(BaseLib::ConfigTree const& config)
+FixedTimeStepping::FixedTimeStepping(double t0,
+ double tn,
+ const std::vector<double>& vec_all_dt)
+ : TimeStepAlgorithm(t0, computeEnd(t0, tn, vec_all_dt), vec_all_dt)
{
- //! \ogs_file_param{prj__time_loop__time_stepping__type}
- config.checkConfigParameter("type", "FixedTimeStepping");
-
- //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__t_initial}
- auto const t_initial = config.getConfigParameter<double>("t_initial");
- //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__t_end}
- auto const t_end = config.getConfigParameter<double>("t_end");
- //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__timesteps}
- auto const delta_ts = config.getConfigSubtree("timesteps");
-
- std::vector<double> timesteps;
- double t_curr = t_initial;
- double delta_t = 0.0;
-
- // TODO: consider adding call "listNonEmpty" to config tree
- //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__timesteps__pair}
- auto const range = delta_ts.getConfigSubtreeList("pair");
- if (range.begin() == range.end()) {
- OGS_FATAL("no timesteps have been given");
- }
- for (auto const pair : range)
- {
- //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__timesteps__pair__repeat}
- auto const repeat = pair.getConfigParameter<std::size_t>("repeat");
- //! \ogs_file_param{prj__time_loop__time_stepping__FixedTimeStepping__timesteps__pair__delta_t}
- delta_t = pair.getConfigParameter<double>("delta_t");
-
- if (repeat == 0) {
- OGS_FATAL("<repeat> is zero.");
- }
- if (delta_t <= 0.0) {
- OGS_FATAL("timestep <delta_t> is <= 0.0.");
- }
-
- if (t_curr <= t_end) {
- timesteps.resize(timesteps.size() + repeat, delta_t);
-
- t_curr += repeat * delta_t;
- }
- }
-
- // append last delta_t until t_end is reached
- if (t_curr <= t_end)
- {
- auto const repeat =
- static_cast<std::size_t>(std::ceil((t_end - t_curr) / delta_t));
- timesteps.resize(timesteps.size() + repeat, delta_t);
- }
-
- return std::make_unique<FixedTimeStepping>(t_initial, t_end, timesteps);
}
-const TimeStep FixedTimeStepping::getTimeStep() const
+FixedTimeStepping::FixedTimeStepping(double t0, double tn, double dt)
+ : TimeStepAlgorithm(t0, tn, dt)
{
- return _ts_current;
}
-bool FixedTimeStepping::next()
+bool FixedTimeStepping::next(const double /*solution_error*/)
{
// check if last time step
- if (_ts_current.steps() == _dt_vector.size()
- || std::abs(_ts_current.current()-_t_end) < std::numeric_limits<double>::epsilon())
+ if (_ts_current.steps() == _dt_vector.size() ||
+ std::abs(_ts_current.current() - _t_end) <
+ std::numeric_limits<double>::epsilon())
return false;
// confirm current time and move to the next if accepted
@@ -107,17 +45,20 @@ bool FixedTimeStepping::next()
// prepare the next time step info
_ts_current = _ts_prev;
double dt = _dt_vector[_ts_prev.steps()];
- if (_ts_prev.current() + dt > _t_end) // upper bound by t_end
+ if (_ts_prev.current() + dt > _t_end) // upper bound by t_end
dt = _t_end - _ts_prev.current();
_ts_current += dt;
return true;
}
-double FixedTimeStepping::computeEnd(double t_initial, double t_end, const std::vector<double> &dt_vector)
+double FixedTimeStepping::computeEnd(double t_initial,
+ double t_end,
+ const std::vector<double>& dt_vector)
{
- double t_sum = t_initial + std::accumulate(dt_vector.begin(), dt_vector.end(), 0.);
+ double t_sum =
+ t_initial + std::accumulate(dt_vector.begin(), dt_vector.end(), 0.);
return std::min(t_end, t_sum);
}
-} //NumLib
+} // NumLib
diff --git a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h
index af57f4f9ddb65137060d2ece618019513101ea81..8cdffd4250737599ca62bd4635841099594bf8fb 100644
--- a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h
+++ b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.h
@@ -14,20 +14,16 @@
#include <memory>
#include <vector>
-#include "ITimeStepAlgorithm.h"
-
-namespace BaseLib { class ConfigTree; }
+#include "TimeStepAlgorithm.h"
namespace NumLib
{
-
/**
* \brief Fixed time stepping algorithm
*
* This algorithm returns time step size defined by a user priori.
*/
-class FixedTimeStepping final
- : public ITimeStepAlgorithm
+class FixedTimeStepping final : public TimeStepAlgorithm
{
public:
/**
@@ -58,43 +54,18 @@ public:
* @param t_end finish time
* @param vec_all_dt a vector of all time steps
*/
- FixedTimeStepping(double t_initial, double t_end, const std::vector<double> &vec_all_dt);
-
- /// Create timestepper from the given configuration
- static std::unique_ptr<ITimeStepAlgorithm> newInstance(BaseLib::ConfigTree const& config);
-
- /// return the beginning of time steps
- double begin() const override { return _t_initial; }
-
- /// return the end of time steps
- double end() const override { return _t_end; }
-
- /// return current time step
- const TimeStep getTimeStep() const override;
+ FixedTimeStepping(double t_initial, double t_end,
+ const std::vector<double>& vec_all_dt);
/// move to the next time step
- bool next() override;
+ bool next(const double solution_error) override;
/// return if current time step is accepted
bool accepted() const override { return true; }
-
- /// return a history of time step sizes
- const std::vector<double>& getTimeStepSizeHistory() const override { return _dt_vector; }
-
private:
/// determine true end time
- static double computeEnd(double t_initial, double t_end, const std::vector<double> &dt_vector);
-
- /// initial time
- const double _t_initial;
- /// end time
- const double _t_end;
- /// a vector of time step sizes
- const std::vector<double> _dt_vector;
- /// previous time step information
- TimeStep _ts_prev;
- /// current time step information
- TimeStep _ts_current;
+ static double computeEnd(double t_initial, double t_end,
+ const std::vector<double>& dt_vector);
};
-} //NumLib
+} // NumLib
diff --git a/NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h b/NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h
deleted file mode 100644
index 2dda10b913d87b0b37b5d4df2923b3128405b9af..0000000000000000000000000000000000000000
--- a/NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/**
- * \author Norihiro Watanabe
- * \date 2012-08-03
- *
- * \copyright
- * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- */
-
-#pragma once
-
-#include <vector>
-
-#include "NumLib/TimeStepping/TimeStep.h"
-
-namespace NumLib
-{
-
-/**
- * \brief Interface of time stepping algorithms
- *
- */
-class ITimeStepAlgorithm
-{
-public:
- /// return the beginning of time steps
- virtual double begin() const = 0;
-
- /// return the end of time steps
- virtual double end() const = 0;
-
- /// move to the next time step
- /// \return true if the next step exists
- virtual bool next() = 0;
-
- /// return current time step
- virtual const TimeStep getTimeStep() const = 0;
-
- /// return if current time step is accepted or not
- virtual bool accepted() const = 0;
-
- /// return a history of time step sizes
- virtual const std::vector<double>& getTimeStepSizeHistory() const = 0;
-
- virtual ~ITimeStepAlgorithm() = default;
-};
-
-} //NumLib
diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp
index d6e27d960e357f54a1d690ce9fe0ea669dcfb62b..db913264a371c83e2bf344ffe0b4045b7303cd26 100644
--- a/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.cpp
@@ -22,57 +22,54 @@ namespace NumLib
IterationNumberBasedAdaptiveTimeStepping::
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)
- : _t_initial(t0),
- _t_end(tn),
- _iter_times_vector(std::move(iter_times_vector)),
- _multiplier_vector(std::move(multiplier_vector)),
+ std::vector<std::size_t>& iter_times_vector,
+ std::vector<double>& multiplier_vector)
+ : TimeStepAlgorithm(t0, tn),
+ _iter_times_vector(iter_times_vector),
+ _multiplier_vector(multiplier_vector),
_min_ts(min_ts),
_max_ts(max_ts),
_initial_ts(initial_ts),
_max_iter(_iter_times_vector.empty() ? 0 : _iter_times_vector.back()),
_iter_times(0),
- _ts_pre(t0),
- _ts_current(t0),
_n_rejected_steps(0)
{
assert(iter_times_vector.size() == multiplier_vector.size());
}
-bool IterationNumberBasedAdaptiveTimeStepping::next()
+bool IterationNumberBasedAdaptiveTimeStepping::next(
+ const double /*solution_error*/)
{
// check current time step
- if (std::abs(_ts_current.current()-_t_end) < std::numeric_limits<double>::epsilon())
+ if (std::abs(_ts_current.current() - _t_end) <
+ std::numeric_limits<double>::epsilon())
return false;
// confirm current time and move to the next if accepted
- if (accepted()) {
- _ts_pre = _ts_current;
+ if (accepted())
+ {
+ _ts_prev = _ts_current;
_dt_vector.push_back(_ts_current.dt());
- } else {
+ }
+ else
+ {
++_n_rejected_steps;
}
// prepare the next time step info
- _ts_current = _ts_pre;
+ _ts_current = _ts_prev;
_ts_current += getNextTimeStepSize();
return true;
}
-const TimeStep IterationNumberBasedAdaptiveTimeStepping::getTimeStep() const
-{
- return _ts_current;
-}
-
double IterationNumberBasedAdaptiveTimeStepping::getNextTimeStepSize() const
{
double dt = 0.0;
// if this is the first time step
// then we use initial guess provided by a user
- if ( _ts_pre.steps() == 0 )
+ if (_ts_prev.steps() == 0)
{
dt = _initial_ts;
}
@@ -83,29 +80,29 @@ double IterationNumberBasedAdaptiveTimeStepping::getNextTimeStepSize() const
if (!_multiplier_vector.empty())
tmp_multiplier = _multiplier_vector[0];
// finding the right multiplier
- for (std::size_t i=0; i<_iter_times_vector.size(); i++ )
- if ( this->_iter_times > _iter_times_vector[i] )
+ for (std::size_t i = 0; i < _iter_times_vector.size(); i++)
+ if (this->_iter_times > _iter_times_vector[i])
tmp_multiplier = _multiplier_vector[i];
// multiply the the multiplier
- dt = _ts_pre.dt() * tmp_multiplier;
+ dt = _ts_prev.dt() * tmp_multiplier;
}
// check whether out of the boundary
- if ( dt < _min_ts )
+ if (dt < _min_ts)
dt = _min_ts;
- else if ( dt > _max_ts )
+ else if (dt > _max_ts)
dt = _max_ts;
- double t_next = dt + _ts_pre.current();
+ double t_next = dt + _ts_prev.current();
if (t_next > end())
- dt = end() - _ts_pre.current();
+ dt = end() - _ts_prev.current();
return dt;
}
bool IterationNumberBasedAdaptiveTimeStepping::accepted() const
{
- return ( this->_iter_times <= this->_max_iter );
+ return (_iter_times <= _max_iter);
}
-} // NumLib
+} // NumLib
diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.h b/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.h
index 38058e29d314b10b0cacf9fa6a4ccabddbc7b936..0cdf734d5caf27de8dd378dd8ac8c5f2bcefade6 100644
--- a/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.h
+++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedAdaptiveTimeStepping.h
@@ -13,24 +13,26 @@
#include <vector>
-#include "ITimeStepAlgorithm.h"
+#include "TimeStepAlgorithm.h"
namespace NumLib
{
-
/**
* \brief Iteration number based adaptive time stepping
*
* This algorithm estimates a time step size depending on the number of
* iterations (e.g. of iterative linear solvers, nonlinear methods, partitioned
- * coupling) needed in the last time step (see Hoffmann (2010) for Newton-Raphson case).
+ * coupling) needed in the last time step (see Hoffmann (2010) for
+ * Newton-Raphson case).
* The new time step \f$\Delta t_{n+1}\f$ size is calculated as
* \f[
* \Delta t_{n+1} = \alpha \Delta t_n
* \f]
- * with the previous time step size \f$\Delta t_{n}\f$ and a multiplier coefficient
+ * with the previous time step size \f$\Delta t_{n}\f$ and a multiplier
+ * coefficient
* \f$\alpha\f$ depending on the iteration number.
- * Note that a time step size is always bounded by the minimum and maximum allowed
+ * Note that a time step size is always bounded by the minimum and maximum
+ * allowed
* value.
* \f[
* \Delta t_{\min} \le \Delta t \le \Delta t_{\max}
@@ -39,24 +41,30 @@ namespace NumLib
* For example, users can setup the following time stepping strategy based on
* the iteration number of the Newton-Raphson method in the last time step.
* <table border="1">
- * <tr><th>Num. of Newton steps</th><th>0-2</th><th>3-6</th><th>7-8</th><th>9<</th></tr>
- * <tr><th>Time step size multiplier</th><th>1.6</th><th>1.</th><th>0.5</th><th>0.25 (repeat time step)</th></tr>
- * <tr><th>Upper and lower bound</th><th colspan="4"> \f$ 1. \le \Delta t \le 10.\f$ </th></tr>
+ * <tr><th>Num. of Newton
+ * steps</th><th>0-2</th><th>3-6</th><th>7-8</th><th>9<</th></tr>
+ * <tr><th>Time step size
+ * multiplier</th><th>1.6</th><th>1.</th><th>0.5</th><th>0.25 (repeat time
+ * step)</th></tr>
+ * <tr><th>Upper and lower bound</th><th colspan="4"> \f$ 1. \le \Delta t \le
+ * 10.\f$ </th></tr>
* </table>
- * A time step size is increased for the small iteration number, and decreased for the
- * large iteration number. If the iteration number exceeds a user-defined threshold (e.g. 9),
+ * A time step size is increased for the small iteration number, and decreased
+ * for the
+ * large iteration number. If the iteration number exceeds a user-defined
+ * threshold (e.g. 9),
* a time step is repeated with a smaller time step size.
*
* Reference
* - Hoffmann J (2010) Reactive Transport and Mineral Dissolution/Precipitation
* in Porous Media:Efficient Solution Algorithms, Benchmark Computations and
- * Existence of Global Solutions. PhD thesis. pp82. Friedrich-Alexander-Universität Erlangen-Nürnberg.
+ * Existence of Global Solutions. PhD thesis. pp82.
+ * Friedrich-Alexander-Universität Erlangen-Nürnberg.
*
*/
-class IterationNumberBasedAdaptiveTimeStepping : public ITimeStepAlgorithm
+class IterationNumberBasedAdaptiveTimeStepping final : public TimeStepAlgorithm
{
public:
-
/**
* Constructor
*
@@ -71,55 +79,41 @@ public:
* If an iteration number is larger than \f$i_n\f$, current time step is
* repeated with the new time step size.
* @param multiplier_vector a vector of multiplier coefficients
- * (\f$a_1\f$, \f$a_2\f$, ..., \f$a_n\f$) corresponding to the intervals given by iter_times_vector.
+ * (\f$a_1\f$, \f$a_2\f$, ..., \f$a_n\f$) corresponding to the intervals
+ * given by iter_times_vector.
* A time step size is calculated by \f$\Delta t_{n+1} = a * \Delta t_{n}\f$
*/
- IterationNumberBasedAdaptiveTimeStepping(double t_initial,
- double t_end,
- double min_ts,
- double max_ts,
- double initial_ts,
- std::vector<std::size_t>
- iter_times_vector,
- std::vector<double>
- multiplier_vector);
+ IterationNumberBasedAdaptiveTimeStepping(
+ double t_initial,
+ double t_end,
+ double min_ts,
+ double max_ts,
+ double initial_ts,
+ std::vector<std::size_t>& iter_times_vector,
+ std::vector<double>& multiplier_vector);
~IterationNumberBasedAdaptiveTimeStepping() override = default;
- /// return the beginning of time steps
- double begin() const override { return _t_initial; }
- /// return the end of time steps
- double end() const override { return _t_end; }
- /// return current time step
- const TimeStep getTimeStep() const override;
-
/// move to the next time step
- bool next() override;
+ bool next(const double solution_error) override;
/// return if the current step is accepted
bool accepted() const override;
- /// return a history of time step sizes
- const std::vector<double>& getTimeStepSizeHistory() const override
- {
- return this->_dt_vector;
- }
-
/// set the number of iterations
- void setNIterations(std::size_t n_itr) {this->_iter_times = n_itr;}
-
+ void setIterationNumber(std::size_t n_itr) { this->_iter_times = n_itr; }
/// return the number of repeated steps
- std::size_t getNumberOfRepeatedSteps() const {return this->_n_rejected_steps;}
+ std::size_t getNumberOfRepeatedSteps() const
+ {
+ return this->_n_rejected_steps;
+ }
private:
/// calculate the next time step size
double getNextTimeStepSize() const;
- /// initial time
- const double _t_initial;
- /// end time
- const double _t_end;
- /// this vector stores the number of iterations to which the respective multiplier coefficient will be applied
+ /// this vector stores the number of iterations to which the respective
+ /// multiplier coefficient will be applied
const std::vector<std::size_t> _iter_times_vector;
/// this vector stores the multiplier coefficients
const std::vector<double> _multiplier_vector;
@@ -133,14 +127,8 @@ private:
const std::size_t _max_iter;
/// the number of nonlinear iterations
std::size_t _iter_times;
- /// previous time step
- TimeStep _ts_pre;
- /// current time step
- TimeStep _ts_current;
- /// history of time step sizes
- std::vector<double> _dt_vector;
/// the number of rejected steps
std::size_t _n_rejected_steps;
};
-} // NumLib
+} // NumLib
diff --git a/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h b/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h
new file mode 100644
index 0000000000000000000000000000000000000000..1c459e5b7972963f64f992988d2e963cd53a32c1
--- /dev/null
+++ b/NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h
@@ -0,0 +1,99 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-08-03
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+
+#include "NumLib/TimeStepping/TimeStep.h"
+
+namespace NumLib
+{
+/**
+ * \brief Interface of time stepping algorithms
+ *
+ */
+class TimeStepAlgorithm
+{
+public:
+ TimeStepAlgorithm(const double t0, const double t_end)
+ : _t_initial(t0), _t_end(t_end), _ts_prev(t0), _ts_current(t0)
+ {
+ }
+
+ TimeStepAlgorithm(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)
+ {
+ }
+
+ TimeStepAlgorithm(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)
+ {
+ }
+
+ virtual ~TimeStepAlgorithm() = default;
+
+ /// return the beginning of time steps
+ double begin() const { return _t_initial; }
+ /// return the end of time steps
+ double end() const { return _t_end; }
+ /// return current time step
+ const TimeStep getTimeStep() const { return _ts_current; }
+ /// reset the current step size from the previous time
+ void resetCurrentTimeStep(const double dt)
+ {
+ _ts_current = _ts_prev;
+ _ts_current += dt;
+ }
+
+ /// move to the next time step
+ /// \param solution_error Solution error between two successive time steps.
+ /// \return true if the next step exists
+ virtual bool next(const double solution_error) = 0;
+
+ /// return if current time step is accepted or not
+ virtual bool accepted() const = 0;
+
+ /// return a history of time step sizes
+ const std::vector<double>& getTimeStepSizeHistory() const
+ {
+ return _dt_vector;
+ }
+
+ /// 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;
+ /// end time
+ const double _t_end;
+
+ /// previous time step information
+ TimeStep _ts_prev;
+ /// current time step information
+ TimeStep _ts_current;
+
+ /// a vector of time step sizes
+ std::vector<double> _dt_vector;
+};
+
+} // NumLib
diff --git a/NumLib/TimeStepping/CreateTimeStepper.cpp b/NumLib/TimeStepping/CreateTimeStepper.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3326b397756a75fb8836eb8cf5136e1875562d02
--- /dev/null
+++ b/NumLib/TimeStepping/CreateTimeStepper.cpp
@@ -0,0 +1,61 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateTimeStepper.cpp
+ * Created on May 2, 2017, 12:18 PM
+ */
+
+#include "CreateTimeStepper.h"
+
+#include <memory>
+#include <string>
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/Error.h"
+
+#include "NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h"
+#include "NumLib/TimeStepping/Algorithms/FixedTimeStepping.h"
+#include "NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.h"
+
+namespace NumLib
+{
+std::unique_ptr<TimeStepAlgorithm> createTimeStepper(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{prj__time_loop__time_stepping__type}
+ auto const type = config.peekConfigParameter<std::string>("type");
+
+ std::unique_ptr<NumLib::TimeStepAlgorithm> timestepper;
+
+ if (type == "SingleStep")
+ {
+ //! \ogs_file_param_special{prj__time_loop__time_stepping__SingleStep}
+ config.ignoreConfigParameter("type");
+ timestepper =
+ std::make_unique<NumLib::FixedTimeStepping>(0.0, 1.0, 1.0);
+ }
+ else if (type == "FixedTimeStepping")
+ {
+ timestepper = NumLib::createFixedTimeStepping(config);
+ }
+ else if (type == "EvolutionaryPIDcontroller")
+ {
+ timestepper = NumLib::createEvolutionaryPIDcontroller(config);
+ }
+ else
+ {
+ OGS_FATAL(
+ "Unknown time stepping type: `%s'. "
+ "The available types are \n\tSingleStep, \n\tFixedTimeStepping"
+ "\n\tEvolutionaryPIDcontroller\n",
+ type.data());
+ }
+
+ return timestepper;
+}
+
+} // end of namespace NumLib
\ No newline at end of file
diff --git a/NumLib/TimeStepping/CreateTimeStepper.h b/NumLib/TimeStepping/CreateTimeStepper.h
new file mode 100644
index 0000000000000000000000000000000000000000..25434a686b83cac992c5e122e2ec23e2922050bc
--- /dev/null
+++ b/NumLib/TimeStepping/CreateTimeStepper.h
@@ -0,0 +1,26 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateTimeStepper.h
+ * Created on May 2, 2017, 12:18 PM
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace NumLib
+{
+class TimeStepAlgorithm;
+std::unique_ptr<TimeStepAlgorithm> createTimeStepper(
+ BaseLib::ConfigTree const& config);
+};
diff --git a/NumLib/TimeStepping/TimeStep.h b/NumLib/TimeStepping/TimeStep.h
index 1a593a7035b26c1923639d1552c76132ab4e9e1e..a85ebb4827fecf13d5e4c2feed6b5c0d171ad5b9 100644
--- a/NumLib/TimeStepping/TimeStep.h
+++ b/NumLib/TimeStepping/TimeStep.h
@@ -15,12 +15,13 @@
namespace NumLib
{
-
/**
* \brief Time step object
*
- * TimeStep class contains previous time(\f$t_{n}\f$), current time (\f$t_{n+1}\f$),
- * time step length (\f$\Delta t_{n+1}\f$), and the number of time steps (\f$n+1\f$).
+ * TimeStep class contains previous time(\f$t_{n}\f$), current time
+ * (\f$t_{n+1}\f$),
+ * time step length (\f$\Delta t_{n+1}\f$), and the number of time steps
+ * (\f$n+1\f$).
*/
class TimeStep
{
@@ -30,7 +31,12 @@ public:
* @param current_time current time
*/
explicit TimeStep(double current_time)
- : _previous(current_time), _current(current_time), _dt(_current-_previous), _steps(0) {}
+ : _previous(current_time),
+ _current(current_time),
+ _dt(_current - _previous),
+ _steps(0)
+ {
+ }
/**
* Initialize a time step
@@ -39,11 +45,21 @@ public:
* @param n the number of time steps
*/
TimeStep(double previous_time, double current_time, std::size_t n)
- : _previous(previous_time), _current(current_time), _dt(_current-_previous), _steps(n) {}
+ : _previous(previous_time),
+ _current(current_time),
+ _dt(_current - _previous),
+ _steps(n)
+ {
+ }
/// copy a time step
- TimeStep(const TimeStep &src)
- : _previous(src._previous), _current(src._current), _dt(_current-_previous), _steps(src._steps) {}
+ TimeStep(const TimeStep& src)
+ : _previous(src._previous),
+ _current(src._current),
+ _dt(_current - _previous),
+ _steps(src._steps)
+ {
+ }
/// copy a time step
TimeStep& operator=(const TimeStep& src) = default;
@@ -57,7 +73,7 @@ public:
}
/// increment time step
- TimeStep &operator+=(const double dt)
+ TimeStep& operator+=(const double dt)
{
_previous = _current;
_current += dt;
@@ -67,35 +83,33 @@ public:
}
/// compare current time
- bool operator<(const TimeStep &t) const {return (_current < t._current);}
-
+ bool operator<(const TimeStep& t) const { return (_current < t._current); }
/// compare current time
- bool operator<(const double &t) const {return (_current < t);}
-
+ bool operator<(const double& t) const { return (_current < t); }
/// compare current time
- bool operator<=(const TimeStep &t) const {return (_current <= t._current);}
+ bool operator<=(const TimeStep& t) const
+ {
+ return (_current <= t._current);
+ }
/// compare current time
- bool operator<=(const double &t) const {return (_current <= t);}
-
+ bool operator<=(const double& t) const { return (_current <= t); }
/// compare current time
- bool operator==(const TimeStep &t) const {return (_current == t._current);}
+ bool operator==(const TimeStep& t) const
+ {
+ return (_current == t._current);
+ }
/// compare current time
- bool operator==(const double &t) const {return (_current == t);}
-
+ bool operator==(const double& t) const { return (_current == t); }
/// return previous time step
- double previous() const {return _previous;}
-
+ double previous() const { return _previous; }
/// return current time step
- double current() const {return _current;}
-
+ double current() const { return _current; }
/// time step size from _previous
- double dt() const {return _dt;}
-
+ double dt() const { return _dt; }
/// the number of time _steps
- std::size_t steps() const {return _steps;}
-
+ std::size_t steps() const { return _steps; }
private:
/// previous time step
double _previous;
@@ -107,4 +121,4 @@ private:
std::size_t _steps;
};
-} //NumLib
+} // NumLib
diff --git a/ProcessLib/AbstractJacobianAssembler.h b/ProcessLib/AbstractJacobianAssembler.h
index 4505fc8c6f95e6ab94892d853cbf14c35e0c262b..d081da64d4f45b66c986fd9feb01dfc19a2fcebe 100644
--- a/ProcessLib/AbstractJacobianAssembler.h
+++ b/ProcessLib/AbstractJacobianAssembler.h
@@ -32,7 +32,7 @@ public:
//! Assembles the Jacobian, the matrices \f$M\f$ and \f$K\f$, and the vector
//! \f$b\f$ with coupling.
- virtual void assembleWithJacobianAndCouping(
+ virtual void assembleWithJacobianAndCoupling(
LocalAssemblerInterface& /*local_assembler*/, double const /*t*/,
std::vector<double> const& /*local_x*/,
std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
diff --git a/ProcessLib/LiquidFlow/Tests.cmake b/ProcessLib/LiquidFlow/Tests.cmake
index cb3d9dfc9df5fb0070f63b2cddb2070723144eda..d2933143d52615305e4acdd59e0ac982addc226f 100644
--- a/ProcessLib/LiquidFlow/Tests.cmake
+++ b/ProcessLib/LiquidFlow/Tests.cmake
@@ -86,12 +86,29 @@ AddTest(
WRAPPER time
TESTER vtkdiff
REQUIREMENTS NOT OGS_USE_MPI
- ABSTOL 1e-16 RELTOL 1e-10
+ ABSTOL 1.1 RELTOL 1e-7
DIFF_DATA
- mesh2D.vtu gravity_driven_pcs_1_ts_5_t_300.000000.vtu pressure_ref pressure
- mesh2D.vtu gravity_driven_pcs_1_ts_5_t_300.000000.vtu temperature_ref temperature
- mesh2D.vtu gravity_driven_pcs_1_ts_5_t_300.000000.vtu v_x_ref v_x
- mesh2D.vtu gravity_driven_pcs_1_ts_5_t_300.000000.vtu v_y_ref v_y
+ mesh2D.vtu gravity_driven_pcs_1_ts_10_t_300.000000.vtu OGS5_PRESSURE1 pressure
+ mesh2D.vtu gravity_driven_pcs_1_ts_10_t_300.000000.vtu OGS5_TEMPERATURE1 temperature
+ mesh2D.vtu gravity_driven_pcs_1_ts_10_t_300.000000.vtu v_x_ref v_x
+ mesh2D.vtu gravity_driven_pcs_1_ts_10_t_300.000000.vtu v_y_ref v_y
+)
+
+# Coupling
+AddTest(
+ NAME Adaptive_dt_StaggeredTH_ThermalDensityDrivenFlow2D
+ PATH StaggeredCoupledProcesses/TH/ThermalDensityDrivenFlow2D
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS thermal_gravity_driven_flow2d_adaptive_dt.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ ABSTOL 1.1 RELTOL 1e-6
+ DIFF_DATA
+ mesh2D.vtu gravity_driven_adaptive_dt_pcs_1_ts_11_t_300.000000.vtu OGS5_PRESSURE1 pressure
+ mesh2D.vtu gravity_driven_adaptive_dt_pcs_1_ts_11_t_300.000000.vtu OGS5_TEMPERATURE1 temperature
+ mesh2D.vtu gravity_driven_adaptive_dt_pcs_1_ts_11_t_300.000000.vtu v_x_ref v_x
+ mesh2D.vtu gravity_driven_adaptive_dt_pcs_1_ts_11_t_300.000000.vtu v_y_ref v_y
)
# PETSc/MPI
@@ -181,10 +198,28 @@ AddTest(
WRAPPER_ARGS -np 1
TESTER vtkdiff
REQUIREMENTS OGS_USE_MPI
- ABSTOL 1e-10 RELTOL 1e-10
+ ABSTOL 1.1 RELTOL 1e-7
+ DIFF_DATA
+ mesh2D.vtu gravity_driven_pcs_1_ts_10_t_300_000000_0.vtu OGS5_PRESSURE1 pressure
+ mesh2D.vtu gravity_driven_pcs_1_ts_10_t_300_000000_0.vtu OGS5_TEMPERATURE1 temperature
+# To be activated when the output of velocity is correct under PETSc version.
+# mesh2D.vtu gravity_driven_pcs_1_ts_10_t_300_000000_0.vtu v_x_ref v_x
+# mesh2D.vtu gravity_driven_pcs_1_ts_10_t_300_000000_0.vtu v_y_ref v_y
+
+)
+
+AddTest(
+ NAME Adaptive_dt_StaggeredTH_ThermalDensityDrivenFlow2D
+ PATH StaggeredCoupledProcesses/TH/ThermalDensityDrivenFlow2D
+ EXECUTABLE_ARGS thermal_gravity_driven_flow2d_adaptive_dt.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ REQUIREMENTS OGS_USE_MPI
+ ABSTOL 1.1 RELTOL 1e-7
DIFF_DATA
- mesh2D.vtu gravity_driven_pcs_1_ts_5_t_300_000000_0.vtu pressure_ref pressure
- mesh2D.vtu gravity_driven_pcs_1_ts_5_t_300_000000_0.vtu temperature_ref temperature
+ mesh2D.vtu gravity_driven_adaptive_dt_pcs_1_ts_11_t_300_000000_0.vtu OGS5_PRESSURE1 pressure
+ mesh2D.vtu gravity_driven_adaptive_dt_pcs_1_ts_11_t_300_000000_0.vtu OGS5_TEMPERATURE1 temperature
# To be activated when the output of velocity is correct under PETSc version.
# mesh2D.vtu gravity_driven_pcs_1_ts_5_t_300_000000_0.vtu v_x_ref v_x
# mesh2D.vtu gravity_driven_pcs_1_ts_5_t_300_000000_0.vtu v_y_ref v_y
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index c626e7cb71cd2b915e6f5fad759817bc2318929a..0732fd810248f3a57ea15a6d24cd4e8c71352049 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -13,7 +13,6 @@
namespace ProcessLib
{
-
void LocalAssemblerInterface::assembleWithCoupledTerm(
double const /*t*/, std::vector<double> const& /*local_x*/,
std::vector<double>& /*local_M_data*/,
@@ -22,8 +21,8 @@ void LocalAssemblerInterface::assembleWithCoupledTerm(
LocalCouplingTerm const& /*coupling_term*/)
{
OGS_FATAL(
- "The assembleWithCoupledTerm() function is not implemented in the local "
- "assembler.");
+ "The assembleWithCoupledTerm() function is not implemented in the "
+ "local assembler.");
}
void LocalAssemblerInterface::assembleWithJacobian(
@@ -39,7 +38,7 @@ void LocalAssemblerInterface::assembleWithJacobian(
"assembler.");
}
-void LocalAssemblerInterface::assembleWithJacobianAndCouping(
+void LocalAssemblerInterface::assembleWithJacobianAndCoupling(
double const /*t*/, std::vector<double> const& /*local_x*/,
std::vector<double> const& /*local_xdot*/, const double /*dxdot_dx*/,
const double /*dx_dx*/, std::vector<double>& /*local_M_data*/,
@@ -49,15 +48,14 @@ void LocalAssemblerInterface::assembleWithJacobianAndCouping(
LocalCouplingTerm const& /*coupling_term*/)
{
OGS_FATAL(
- "The assembleWithJacobianAndCouping() function is not implemented in"
+ "The assembleWithJacobianAndCoupling() function is not implemented in"
" the local assembler.");
}
void LocalAssemblerInterface::computeSecondaryVariable(
- std::size_t const mesh_item_id,
- NumLib::LocalToGlobalIndexMap const& dof_table,
- double const t, GlobalVector const& x,
- StaggeredCouplingTerm const& coupled_term)
+ std::size_t const mesh_item_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table, double const t,
+ GlobalVector const& x, StaggeredCouplingTerm const& coupled_term)
{
auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
auto const local_x = x.get(indices);
@@ -68,11 +66,14 @@ void LocalAssemblerInterface::computeSecondaryVariable(
}
else
{
- auto const local_coupled_xs
- = getCurrentLocalSolutionsOfCoupledProcesses(
- coupled_term.coupled_xs, indices);
- computeSecondaryVariableWithCoupledProcessConcrete(t, local_x,
- local_coupled_xs);
+ auto const local_coupled_xs =
+ getCurrentLocalSolutionsOfCoupledProcesses(coupled_term.coupled_xs,
+ indices);
+ if (!local_coupled_xs.empty())
+ computeSecondaryVariableWithCoupledProcessConcrete(
+ t, local_x, local_coupled_xs);
+ else
+ computeSecondaryVariableConcrete(t, local_x);
}
}
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 040e6123f39f15b892bde0074900094143fdd00f..c9ad2950358f366fa3cf994720d5bd4c1b7c2d49 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -56,7 +56,7 @@ public:
std::vector<double>& local_b_data,
std::vector<double>& local_Jac_data);
- virtual void assembleWithJacobianAndCouping(double const t,
+ virtual void assembleWithJacobianAndCoupling(double const t,
std::vector<double> const& local_x,
std::vector<double> const& local_xdot,
const double dxdot_dx, const double dx_dx,
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index f31970e6894ce354fe4bdacf8ee519e8e72746aa..e8b122ca6a79706d77a373a1f1d0657cff094ff9 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -8,42 +8,16 @@
*/
#include "UncoupledProcessesTimeLoop.h"
+
#include "BaseLib/uniqueInsert.h"
#include "BaseLib/RunTime.h"
#include "NumLib/ODESolver/TimeDiscretizationBuilder.h"
#include "NumLib/ODESolver/TimeDiscretizedODESystem.h"
#include "NumLib/ODESolver/ConvergenceCriterionPerComponent.h"
-#include "NumLib/TimeStepping/Algorithms/FixedTimeStepping.h"
+#include "NumLib/TimeStepping/CreateTimeStepper.h"
#include "MathLib/LinAlg/LinAlg.h"
-std::unique_ptr<NumLib::ITimeStepAlgorithm> createTimeStepper(
- BaseLib::ConfigTree const& config)
-{
- //! \ogs_file_param{prj__time_loop__time_stepping__type}
- auto const type = config.peekConfigParameter<std::string>("type");
-
- std::unique_ptr<NumLib::ITimeStepAlgorithm> timestepper;
-
- if (type == "SingleStep")
- {
- //! \ogs_file_param_special{prj__time_loop__time_stepping__SingleStep}
- config.ignoreConfigParameter("type");
- timestepper =
- std::make_unique<NumLib::FixedTimeStepping>(0.0, 1.0, 1.0);
- }
- else if (type == "FixedTimeStepping")
- {
- timestepper = NumLib::FixedTimeStepping::newInstance(config);
- }
- else
- {
- OGS_FATAL("Unknown timestepper type: `%s'.", type.c_str());
- }
-
- return timestepper;
-}
-
std::unique_ptr<ProcessLib::Output> createOutput(
BaseLib::ConfigTree const& config, std::string const& output_directory)
{
@@ -129,6 +103,7 @@ struct SingleProcessData
{
template <NumLib::NonlinearSolverTag NLTag>
SingleProcessData(
+ std::unique_ptr<NumLib::TimeStepAlgorithm>&& timestepper_,
NumLib::NonlinearSolver<NLTag>& nonlinear_solver,
std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit_,
std::unique_ptr<NumLib::TimeDiscretization>&& time_disc_,
@@ -139,6 +114,13 @@ struct SingleProcessData
SingleProcessData(SingleProcessData&& spd);
+ std::unique_ptr<NumLib::TimeStepAlgorithm> timestepper;
+
+ //! Flag of skiping time stepping. It is used in the modelling of
+ //! coupled processes. If the stepping of any process reaches a steady state
+ //! or the ending time, the flag is set to true.
+ bool skip_time_stepping = false;
+
//! Tag containing the missing type information necessary to cast the
//! other members of this struct to their concrety types.
NumLib::NonlinearSolverTag const nonlinear_solver_tag;
@@ -159,13 +141,15 @@ struct SingleProcessData
template <NumLib::NonlinearSolverTag NLTag>
SingleProcessData::SingleProcessData(
+ std::unique_ptr<NumLib::TimeStepAlgorithm>&& timestepper_,
NumLib::NonlinearSolver<NLTag>& nonlinear_solver,
std::unique_ptr<NumLib::ConvergenceCriterion>&& conv_crit_,
std::unique_ptr<NumLib::TimeDiscretization>&& time_disc_,
Process& process_,
std::unordered_map<std::type_index, Process const&>&& coupled_processes_,
ProcessOutput&& process_output_)
- : nonlinear_solver_tag(NLTag),
+ : timestepper(std::move(timestepper_)),
+ nonlinear_solver_tag(NLTag),
nonlinear_solver(nonlinear_solver),
conv_crit(std::move(conv_crit_)),
time_disc(std::move(time_disc_)),
@@ -176,7 +160,8 @@ SingleProcessData::SingleProcessData(
}
SingleProcessData::SingleProcessData(SingleProcessData&& spd)
- : nonlinear_solver_tag(spd.nonlinear_solver_tag),
+ : timestepper(std::move(spd.timestepper)),
+ nonlinear_solver_tag(spd.nonlinear_solver_tag),
nonlinear_solver(spd.nonlinear_solver),
conv_crit(std::move(spd.conv_crit)),
time_disc(std::move(spd.time_disc)),
@@ -243,6 +228,7 @@ void setTimeDiscretizedODESystem(SingleProcessData& spd)
}
std::unique_ptr<SingleProcessData> makeSingleProcessData(
+ std::unique_ptr<NumLib::TimeStepAlgorithm>&& timestepper,
NumLib::NonlinearSolverBase& nonlinear_solver,
Process& process,
std::unique_ptr<NumLib::TimeDiscretization>&& time_disc,
@@ -257,18 +243,18 @@ std::unique_ptr<SingleProcessData> makeSingleProcessData(
&nonlinear_solver))
{
return std::make_unique<SingleProcessData>(
- *nonlinear_solver_picard, std::move(conv_crit),
- std::move(time_disc), process, std::move(coupled_processes),
- std::move(process_output));
+ std::move(timestepper), *nonlinear_solver_picard,
+ std::move(conv_crit), std::move(time_disc), process,
+ std::move(coupled_processes), std::move(process_output));
}
if (auto* nonlinear_solver_newton =
dynamic_cast<NumLib::NonlinearSolver<Tag::Newton>*>(
&nonlinear_solver))
{
return std::make_unique<SingleProcessData>(
- *nonlinear_solver_newton, std::move(conv_crit),
- std::move(time_disc), process, std::move(coupled_processes),
- std::move(process_output));
+ std::move(timestepper), *nonlinear_solver_newton,
+ std::move(conv_crit), std::move(time_disc), process,
+ std::move(coupled_processes), std::move(process_output));
}
OGS_FATAL("Encountered unknown nonlinear solver type. Aborting");
@@ -302,6 +288,10 @@ std::vector<std::unique_ptr<SingleProcessData>> createPerProcessData(
//! \ogs_file_param{prj__time_loop__processes__process__time_discretization}
pcs_config.getConfigSubtree("time_discretization"));
+ auto timestepper = NumLib::createTimeStepper(
+ //! \ogs_file_param{prj__time_loop__processes__process__time_stepping}
+ pcs_config.getConfigSubtree("time_stepping"));
+
auto conv_crit = NumLib::createConvergenceCriterion(
//! \ogs_file_param{prj__time_loop__processes__process__convergence_criterion}
pcs_config.getConfigSubtree("convergence_criterion"));
@@ -336,8 +326,9 @@ std::vector<std::unique_ptr<SingleProcessData>> createPerProcessData(
ProcessOutput process_output{pcs_config.getConfigSubtree("output")};
per_process_data.emplace_back(makeSingleProcessData(
- nl_slv, pcs, std::move(time_disc), std::move(conv_crit),
- std::move(coupled_processes), std::move(process_output)));
+ std::move(timestepper), nl_slv, pcs, std::move(time_disc),
+ std::move(conv_crit), std::move(coupled_processes),
+ std::move(process_output)));
}
if (per_process_data.size() != processes.size())
@@ -370,10 +361,6 @@ std::unique_ptr<UncoupledProcessesTimeLoop> createUncoupledProcessesTimeLoop(
coupling_config->getConfigSubtree("convergence_criterion"));
}
- auto timestepper =
- //! \ogs_file_param{prj__time_loop__time_stepping}
- createTimeStepper(config.getConfigSubtree("time_stepping"));
-
auto output =
//! \ogs_file_param{prj__time_loop__output}
createOutput(config.getConfigSubtree("output"), output_directory);
@@ -382,14 +369,27 @@ std::unique_ptr<UncoupledProcessesTimeLoop> createUncoupledProcessesTimeLoop(
//! \ogs_file_param{prj__time_loop__processes}
config.getConfigSubtree("processes"), processes, nonlinear_solvers);
+ const auto minmax_iter = std::minmax_element(
+ per_process_data.begin(),
+ per_process_data.end(),
+ [](std::unique_ptr<SingleProcessData> const& a,
+ std::unique_ptr<SingleProcessData> const& b) {
+ return (a->timestepper->end() < b->timestepper->end());
+ });
+ const double start_time =
+ per_process_data[minmax_iter.first - per_process_data.begin()]
+ ->timestepper->begin();
+ const double end_time =
+ per_process_data[minmax_iter.second - per_process_data.begin()]
+ ->timestepper->end();
+
return std::make_unique<UncoupledProcessesTimeLoop>(
- std::move(timestepper), std::move(output), std::move(per_process_data),
- max_coupling_iterations, std::move(coupling_conv_crit));
+ std::move(output), std::move(per_process_data), max_coupling_iterations,
+ std::move(coupling_conv_crit), start_time, end_time);
}
std::vector<GlobalVector*> setInitialConditions(
double const t0,
- double const delta_t0,
std::vector<std::unique_ptr<SingleProcessData>> const& per_process_data)
{
std::vector<GlobalVector*> process_solutions;
@@ -419,10 +419,8 @@ std::vector<GlobalVector*> setInitialConditions(
auto& nonlinear_solver = spd->nonlinear_solver;
auto& mat_strg = *spd->mat_strg;
auto& conv_crit = *spd->conv_crit;
- time_disc.nextTimestep(t0, delta_t0);
setEquationSystem(nonlinear_solver, ode_sys, conv_crit, nl_tag);
-
nonlinear_solver.assemble(
x0, ProcessLib::createVoidStaggeredCouplingTerm());
time_disc.pushState(
@@ -468,21 +466,19 @@ bool solveOneTimeStepOneProcess(GlobalVector& x, std::size_t const timestep,
bool nonlinear_solver_succeeded =
nonlinear_solver.solve(x, coupling_term, post_iteration_callback);
- auto& mat_strg = *process_data.mat_strg;
- time_disc.pushState(t, x, mat_strg);
-
return nonlinear_solver_succeeded;
}
UncoupledProcessesTimeLoop::UncoupledProcessesTimeLoop(
- std::unique_ptr<NumLib::ITimeStepAlgorithm>&& timestepper,
std::unique_ptr<Output>&& output,
std::vector<std::unique_ptr<SingleProcessData>>&& per_process_data,
const unsigned global_coupling_max_iterations,
- std::unique_ptr<NumLib::ConvergenceCriterion>&& global_coupling_conv_crit)
- : _timestepper{std::move(timestepper)},
- _output(std::move(output)),
+ std::unique_ptr<NumLib::ConvergenceCriterion>&& global_coupling_conv_crit,
+ const double start_time, const double end_time)
+ : _output(std::move(output)),
_per_process_data(std::move(per_process_data)),
+ _start_time(start_time),
+ _end_time(end_time),
_global_coupling_max_iterations(global_coupling_max_iterations),
_global_coupling_conv_crit(std::move(global_coupling_conv_crit))
{
@@ -491,7 +487,7 @@ UncoupledProcessesTimeLoop::UncoupledProcessesTimeLoop(
bool UncoupledProcessesTimeLoop::setCoupledSolutions()
{
// Do nothing if process are not coupled
- if (!_global_coupling_conv_crit && _global_coupling_max_iterations == 1)
+ if ((!_global_coupling_conv_crit) || _global_coupling_max_iterations == 1)
return false;
unsigned pcs_idx = 0;
@@ -542,6 +538,129 @@ bool UncoupledProcessesTimeLoop::setCoupledSolutions()
return true;
}
+double UncoupledProcessesTimeLoop::computeTimeStepping(
+ const double prev_dt, double& t, std::size_t& accepted_steps,
+ std::size_t& rejected_steps)
+{
+ bool all_process_steps_accepted = true;
+ // Get minimum time step size among step sizes of all processes.
+ double dt = std::numeric_limits<double>::max();
+ for (std::size_t i = 0; i < _per_process_data.size(); i++)
+ {
+ auto& ppd = *_per_process_data[i];
+ const auto& timestepper = ppd.timestepper;
+ if (t > timestepper->end())
+ {
+ // skip the process that already reaches the ending time.
+ ppd.skip_time_stepping = true;
+ continue;
+ }
+
+ auto& time_disc = ppd.time_disc;
+ auto const& x = *_process_solutions[i];
+
+ auto const& conv_crit = ppd.conv_crit;
+ const MathLib::VecNormType norm_type =
+ (conv_crit) ? conv_crit->getVectorNormType()
+ : MathLib::VecNormType::NORM2;
+
+ const double solution_error =
+ (timestepper->isSolutionErrorComputationNeeded())
+ ? ((t == timestepper->begin())
+ ? 0. // Always accepts the zeroth step
+ : time_disc->getRelativeChangeFromPreviousTimestep(
+ x, norm_type))
+ : 0.;
+ if (!timestepper->next(solution_error) &&
+ // In case of FixedTimeStepping, which makes timestepper->next(...)
+ // return false when the ending time is reached.
+ t + std::numeric_limits<double>::epsilon() < timestepper->end())
+ {
+ // Not all processes have accepted steps.
+ all_process_steps_accepted = false;
+ }
+
+ if (timestepper->getTimeStep().dt() >
+ std::numeric_limits<double>::min() ||
+ std::abs(t - timestepper->end()) <
+ std::numeric_limits<double>::epsilon())
+ {
+ if (timestepper->getTimeStep().dt() < dt)
+ {
+ dt = timestepper->getTimeStep().dt();
+ }
+ }
+ else
+ {
+ // dt being close to 0 only happens when
+ // t_n + dt > t_s, and dt is forced to be zero. Where t_n the time
+ // of previous time step, and t_s is the specified time taken from
+ // input or the end time. Under this condition, the time stepping
+ // is skipped.
+ ppd.skip_time_stepping = true;
+ }
+ }
+
+ bool is_initial_step = false;
+ // Reset the time step with the minimum step size, dt
+ // Update the solution of the previous time step in time_disc.
+ for (std::size_t i = 0; i < _per_process_data.size(); i++)
+ {
+ const auto& ppd = *_per_process_data[i];
+ auto& timestepper = ppd.timestepper;
+ timestepper->resetCurrentTimeStep(dt);
+
+ if (ppd.skip_time_stepping)
+ continue;
+
+ if (t == timestepper->begin())
+ {
+ is_initial_step = true;
+ continue;
+ }
+
+ if (all_process_steps_accepted)
+ {
+ auto& time_disc = ppd.time_disc;
+ auto& mat_strg = *ppd.mat_strg;
+ auto const& x = *_process_solutions[i];
+ time_disc->pushState(t, x, mat_strg);
+ }
+ else
+ {
+ if (t < _end_time)
+ {
+ WARN(
+ "Time step %d is rejected. "
+ "The computation is back to the previous time.",
+ accepted_steps + rejected_steps);
+ }
+ }
+ }
+
+ if (!is_initial_step)
+ {
+ if (all_process_steps_accepted)
+ accepted_steps++;
+ else
+ {
+ if (t < _end_time)
+ {
+ t -= prev_dt;
+ rejected_steps++;
+ }
+ }
+ }
+
+ // Adjust step size if t < _end_time, while t+dt exceeds the end time
+ if (t < _end_time && t + dt > _end_time)
+ {
+ dt = _end_time - t;
+ }
+
+ return dt;
+}
+
/*
* TODO:
* Now we have a structure inside the time loop which is very similar to the
@@ -574,12 +693,8 @@ bool UncoupledProcessesTimeLoop::loop()
}
}
- auto const t0 = _timestepper->getTimeStep().current(); // time of the IC
- auto const delta_t0 =
- _timestepper->getTimeStep().dt(); // initial time increment
-
// init solution storage
- _process_solutions = setInitialConditions(t0, delta_t0, _per_process_data);
+ _process_solutions = setInitialConditions(_start_time, _per_process_data);
// output initial conditions
{
@@ -589,46 +704,80 @@ bool UncoupledProcessesTimeLoop::loop()
auto& pcs = spd->process;
auto const& x0 = *_process_solutions[pcs_idx];
- pcs.preTimestep(x0, t0, _timestepper->getTimeStep().dt());
- _output->doOutput(pcs, spd->process_output, 0, t0, x0);
+ pcs.preTimestep(x0, _start_time,
+ spd->timestepper->getTimeStep().dt());
+ _output->doOutput(pcs, spd->process_output, 0, _start_time, x0);
++pcs_idx;
}
}
const bool is_staggered_coupling = setCoupledSolutions();
- double t = t0;
- std::size_t timestep = 1; // the first timestep really is number one
+ double t = _start_time;
+ std::size_t accepted_steps = 0;
+ std::size_t rejected_steps = 0;
bool nonlinear_solver_succeeded = true;
- while (_timestepper->next())
+ double dt = computeTimeStepping(0.0, t, accepted_steps, rejected_steps);
+
+ while (t < _end_time)
{
BaseLib::RunTime time_timestep;
time_timestep.start();
- auto const ts = _timestepper->getTimeStep();
- auto const delta_t = ts.dt();
- t = ts.current();
- timestep = ts.steps();
+ t += dt;
+ const double prev_dt = dt;
- INFO("=== timestep #%u (t=%gs, dt=%gs) ==============================",
- timestep, t, delta_t);
+ const std::size_t timesteps = accepted_steps + rejected_steps + 1;
+ // TODO, input option for time unit.
+ INFO("=== Time stepping at step #%u and time %g with step size %g",
+ timesteps, t, dt);
if (is_staggered_coupling)
nonlinear_solver_succeeded =
- solveCoupledEquationSystemsByStaggeredScheme(t, delta_t,
- timestep);
+ solveCoupledEquationSystemsByStaggeredScheme(t, dt, timesteps);
else
nonlinear_solver_succeeded =
- solveUncoupledEquationSystems(t, delta_t, timestep);
+ solveUncoupledEquationSystems(t, dt, timesteps);
- INFO("[time] Time step #%u took %g s.", timestep,
+ INFO("[time] Time step #%u took %g s.", timesteps,
time_timestep.elapsed());
if (!nonlinear_solver_succeeded)
+ {
+ WARN(
+ "Time step %d is rejected due to "
+ "the divergence of the non-linear solver.\n"
+ "\tThe time stepping steps back to the previous time\n"
+ "\tand starts again with the half of the current step size.",
+ timesteps);
+ t -= prev_dt;
+ dt *= 0.5;
+ rejected_steps++;
+ continue;
+ }
+
+ dt = computeTimeStepping(prev_dt, t, accepted_steps, rejected_steps);
+
+ if (t + dt > _end_time ||
+ t + std::numeric_limits<double>::epsilon() > _end_time)
+ break;
+
+ if (dt < std::numeric_limits<double>::epsilon())
+ {
+ WARN(
+ "Time step size of %g is too small.\n"
+ "Time stepping stops at step %u and at time of %g.",
+ dt, timesteps, t);
break;
+ }
}
+ INFO(
+ "The whole computation of the time stepping took %u steps, in which\n"
+ "\t the accepted steps are %u, and the rejected steps are %u.\n",
+ accepted_steps + rejected_steps, accepted_steps, rejected_steps)
+
// output last time step
if (nonlinear_solver_succeeded)
{
@@ -637,7 +786,8 @@ bool UncoupledProcessesTimeLoop::loop()
{
auto& pcs = spd->process;
auto const& x = *_process_solutions[pcs_idx];
- _output->doOutputLastTimestep(pcs, spd->process_output, timestep, t,
+ _output->doOutputLastTimestep(pcs, spd->process_output,
+ accepted_steps + rejected_steps, t,
x);
++pcs_idx;
@@ -654,11 +804,18 @@ bool UncoupledProcessesTimeLoop::solveUncoupledEquationSystems(
unsigned pcs_idx = 0;
for (auto& spd : _per_process_data)
{
- auto& pcs = spd->process;
+ if (spd->skip_time_stepping)
+ {
+ INFO("Process %u is skipped in the time stepping.", pcs_idx);
+ ++pcs_idx;
+ continue;
+ }
+
BaseLib::RunTime time_timestep_process;
time_timestep_process.start();
auto& x = *_process_solutions[pcs_idx];
+ auto& pcs = spd->process;
pcs.preTimestep(x, t, dt);
const auto void_staggered_coupling_term =
@@ -708,7 +865,13 @@ bool UncoupledProcessesTimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
unsigned pcs_idx = 0;
for (auto& spd : _per_process_data)
{
- auto& pcs = spd->process;
+ if (spd->skip_time_stepping)
+ {
+ INFO("Process %u is skipped in the time stepping.", pcs_idx);
+ ++pcs_idx;
+ continue;
+ }
+
BaseLib::RunTime time_timestep_process;
time_timestep_process.start();
@@ -719,7 +882,7 @@ bool UncoupledProcessesTimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
// belongs to process. For some problems, both of the current
// solution and the solution of the previous time step are
// required for the coupling computation.
- pcs.preTimestep(x, t, dt);
+ spd->process.preTimestep(x, t, dt);
// Set the flag of the first iteration be true.
_global_coupling_conv_crit->preFirstIteration();
@@ -750,8 +913,8 @@ bool UncoupledProcessesTimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
timestep_id, t, pcs_idx);
// save unsuccessful solution
- _output->doOutputAlways(pcs, spd->process_output, timestep_id,
- t, x);
+ _output->doOutputAlways(spd->process, spd->process_output,
+ timestep_id, t, x);
break;
}
@@ -791,6 +954,12 @@ bool UncoupledProcessesTimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
unsigned pcs_idx = 0;
for (auto& spd : _per_process_data)
{
+ if (spd->skip_time_stepping)
+ {
+ ++pcs_idx;
+ continue;
+ }
+
auto& pcs = spd->process;
auto& x = *_process_solutions[pcs_idx];
pcs.postTimestep(x);
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.h b/ProcessLib/UncoupledProcessesTimeLoop.h
index 497e1276a9238786284f560d9a33c42aab076800..df0c9e058b3fc1ba66587f458fbe733ee990053c 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.h
+++ b/ProcessLib/UncoupledProcessesTimeLoop.h
@@ -17,7 +17,7 @@
#include <logog/include/logog.hpp>
#include "NumLib/ODESolver/NonlinearSolver.h"
-#include "NumLib/TimeStepping/Algorithms/ITimeStepAlgorithm.h"
+#include "NumLib/TimeStepping/Algorithms/TimeStepAlgorithm.h"
#include "Output.h"
#include "Process.h"
@@ -38,12 +38,12 @@ class UncoupledProcessesTimeLoop
{
public:
explicit UncoupledProcessesTimeLoop(
- std::unique_ptr<NumLib::ITimeStepAlgorithm>&& timestepper,
std::unique_ptr<Output>&& output,
std::vector<std::unique_ptr<SingleProcessData>>&& per_process_data,
const unsigned global_coupling_max_iterations,
std::unique_ptr<NumLib::ConvergenceCriterion>&&
- global_coupling_conv_crit);
+ global_coupling_conv_crit,
+ const double start_time, const double end_time);
bool loop();
@@ -63,10 +63,12 @@ public:
private:
std::vector<GlobalVector*> _process_solutions;
- std::unique_ptr<NumLib::ITimeStepAlgorithm> _timestepper;
std::unique_ptr<Output> _output;
std::vector<std::unique_ptr<SingleProcessData>> _per_process_data;
+ const double _start_time;
+ const double _end_time;
+
/// Maximum iterations of the global coupling.
const unsigned _global_coupling_max_iterations;
/// Convergence criteria of the global coupling iterations.
@@ -110,6 +112,21 @@ private:
*/
bool solveCoupledEquationSystemsByStaggeredScheme(
const double t, const double dt, const std::size_t timestep_id);
+
+ /**
+ * Find the minimum time step size among the predicted step sizes of
+ * processes and step it as common time step size.
+ *
+ * @param prev_dt Previous time step size.
+ * @param t Current time.
+ * @param accepted_steps Accepted time steps that are counted in this
+ * function.
+ * @param rejected_steps Rejected time steps that are counted in this
+ * function.
+ */
+ double computeTimeStepping(const double prev_dt, double& t,
+ std::size_t& accepted_steps,
+ std::size_t& rejected_steps);
};
//! Builds an UncoupledProcessesTimeLoop from the given configuration.
diff --git a/ProcessLib/VectorMatrixAssembler.cpp b/ProcessLib/VectorMatrixAssembler.cpp
index 0420af82b719d319bcb3ab2e163b44d9836d59df..04b9e6ba71bccdd05b26d795828d6e436e9905ff 100644
--- a/ProcessLib/VectorMatrixAssembler.cpp
+++ b/ProcessLib/VectorMatrixAssembler.cpp
@@ -79,13 +79,22 @@ void VectorMatrixAssembler::assemble(
getPreviousLocalSolutionsOfCoupledProcesses(coupling_term, indices);
auto local_coupled_xs = getCurrentLocalSolutionsOfCoupledProcesses(
coupling_term.coupled_xs, indices);
- ProcessLib::LocalCouplingTerm local_coupling_term(
- coupling_term.dt, coupling_term.coupled_processes,
- std::move(local_coupled_xs0), std::move(local_coupled_xs));
- local_assembler.assembleWithCoupledTerm(t, local_x, _local_M_data,
- _local_K_data, _local_b_data,
- local_coupling_term);
+ if (local_coupled_xs0.empty() || local_coupled_xs.empty())
+ {
+ local_assembler.assemble(t, local_x, _local_M_data, _local_K_data,
+ _local_b_data);
+ }
+ else
+ {
+ ProcessLib::LocalCouplingTerm local_coupling_term(
+ coupling_term.dt, coupling_term.coupled_processes,
+ std::move(local_coupled_xs0), std::move(local_coupled_xs));
+
+ local_assembler.assembleWithCoupledTerm(
+ t, local_x, _local_M_data, _local_K_data, _local_b_data,
+ local_coupling_term);
+ }
}
auto const num_r_c = indices.size();
@@ -137,14 +146,23 @@ void VectorMatrixAssembler::assembleWithJacobian(
getPreviousLocalSolutionsOfCoupledProcesses(coupling_term, indices);
auto local_coupled_xs = getCurrentLocalSolutionsOfCoupledProcesses(
coupling_term.coupled_xs, indices);
- ProcessLib::LocalCouplingTerm local_coupling_term(
- coupling_term.dt, coupling_term.coupled_processes,
- std::move(local_coupled_xs0), std::move(local_coupled_xs));
-
- _jacobian_assembler->assembleWithJacobianAndCouping(
- local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
- _local_M_data, _local_K_data, _local_b_data, _local_Jac_data,
- local_coupling_term);
+ if (local_coupled_xs0.empty() || local_coupled_xs.empty())
+ {
+ _jacobian_assembler->assembleWithJacobian(
+ local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
+ _local_M_data, _local_K_data, _local_b_data, _local_Jac_data);
+ }
+ else
+ {
+ ProcessLib::LocalCouplingTerm local_coupling_term(
+ coupling_term.dt, coupling_term.coupled_processes,
+ std::move(local_coupled_xs0), std::move(local_coupled_xs));
+
+ _jacobian_assembler->assembleWithJacobianAndCoupling(
+ local_assembler, t, local_x, local_xdot, dxdot_dx, dx_dx,
+ _local_M_data, _local_K_data, _local_b_data, _local_Jac_data,
+ local_coupling_term);
+ }
}
auto const num_r_c = indices.size();
diff --git a/Tests/Data b/Tests/Data
index 277b6617f11ab54c72f3e6c27560674722869b5e..ca1db9c602b485f66cead0e82924181c26b86539 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 277b6617f11ab54c72f3e6c27560674722869b5e
+Subproject commit ca1db9c602b485f66cead0e82924181c26b86539
diff --git a/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp b/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1e509148fcf6b0f8545a63f5141ba51edc79e4e8
--- /dev/null
+++ b/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp
@@ -0,0 +1,110 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file TestTimeSteppingEvolutionaryPIDcontroller.cpp
+ * Created on April 5, 2017, 12:09 PM
+ */
+
+#include <gtest/gtest.h>
+
+#include <vector>
+#include <memory>
+
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/ConfigTree.h"
+
+#include "NumLib/TimeStepping/TimeStep.h"
+#include "NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.h"
+
+#include "Tests/TestTools.h"
+
+std::unique_ptr<NumLib::TimeStepAlgorithm> createTestTimeStepper(
+ const char xml[])
+{
+ auto const ptree = readXml(xml);
+ BaseLib::ConfigTree conf(ptree, "", BaseLib::ConfigTree::onerror,
+ BaseLib::ConfigTree::onwarning);
+ auto const& sub_config = conf.getConfigSubtree("time_stepping");
+ return NumLib::createEvolutionaryPIDcontroller(sub_config);
+}
+
+TEST(NumLibTimeStepping, testEvolutionaryPIDcontroller)
+{
+ const char xml[] =
+ "<time_stepping>"
+ " <type>EvolutionaryPIDcontroller</type>"
+ " <t_initial> 0.0 </t_initial>"
+ " <t_end> 10 </t_end>"
+ " <dt_guess> 0.01 </dt_guess>"
+ " <dt_min> 0.001 </dt_min>"
+ " <dt_max> 1 </dt_max>"
+ " <rel_dt_min> 0.01 </rel_dt_min>"
+ " <rel_dt_max> 5 </rel_dt_max>"
+ " <tol> 1.e-3 </tol>"
+ "</time_stepping>";
+ auto const PIDStepper = createTestTimeStepper(xml);
+
+ double solution_error = 0.;
+ // 1st step
+ ASSERT_TRUE(PIDStepper->next(solution_error));
+ NumLib::TimeStep ts = PIDStepper->getTimeStep();
+ double h_new = 0.01;
+ double t_previous = 0.;
+ ASSERT_EQ(1u, ts.steps());
+ ASSERT_EQ(t_previous, ts.previous());
+ ASSERT_EQ(t_previous + h_new, ts.current());
+ ASSERT_EQ(h_new, ts.dt());
+ ASSERT_TRUE(PIDStepper->accepted());
+ t_previous += h_new;
+
+ // e_n_minus1 is filled.
+ solution_error = 1.0e-4;
+ PIDStepper->next(solution_error);
+ ts = PIDStepper->getTimeStep();
+ h_new = ts.dt();
+ ASSERT_EQ(2u, ts.steps());
+ const double tol = 1.e-16;
+ ASSERT_NEAR(t_previous, ts.previous(), tol);
+ ASSERT_NEAR(t_previous + h_new, ts.current(), tol);
+ ASSERT_TRUE(PIDStepper->accepted());
+ t_previous += h_new;
+
+ // e_n_minus2 is filled.
+ solution_error = 0.5e-3;
+ PIDStepper->next(solution_error);
+ ts = PIDStepper->getTimeStep();
+ h_new = ts.dt();
+ ASSERT_EQ(3u, ts.steps());
+ ASSERT_NEAR(t_previous, ts.previous(), tol);
+ ASSERT_NEAR(t_previous + h_new, ts.current(), tol);
+ ASSERT_TRUE(PIDStepper->accepted());
+
+ // error > TOL=1.3-3, step rejected and new step size estimated.
+ solution_error = 0.01;
+ PIDStepper->next(solution_error);
+ ts = PIDStepper->getTimeStep();
+ h_new = ts.dt();
+ // No change in ts.steps
+ ASSERT_EQ(3u, ts.steps());
+ // No change in ts.previous(), which is the same as that of the previous
+ // step.
+ ASSERT_NEAR(t_previous, ts.previous(), tol);
+ ASSERT_NEAR(t_previous + h_new, ts.current(), tol);
+ ASSERT_FALSE(PIDStepper->accepted());
+ t_previous += h_new;
+
+ // With e_n, e_n_minus1, e_n_minus2
+ solution_error = 0.4e-3;
+ PIDStepper->next(solution_error);
+ ts = PIDStepper->getTimeStep();
+ h_new = ts.dt();
+ ASSERT_EQ(4u, ts.steps());
+ ASSERT_NEAR(t_previous, ts.previous(), tol);
+ ASSERT_NEAR(t_previous + h_new, ts.current(), tol);
+ ASSERT_TRUE(PIDStepper->accepted());
+}
diff --git a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
index baf12ca923216bb2022bfc3c8720c6cf22bf1363..00c82b9813715a989d08b59050b1e14c26e07c8f 100644
--- a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
+++ b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
@@ -29,7 +29,9 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
std::vector<double> multiplier_vector = {2.0, 1.0, 0.5, 0.25};
NumLib::IterationNumberBasedAdaptiveTimeStepping alg(1, 31, 1, 10, 1, iter_times_vector, multiplier_vector);
- ASSERT_TRUE(alg.next()); // t=2, dt=1
+ const double solution_error = 0.;
+
+ ASSERT_TRUE(alg.next(solution_error)); // t=2, dt=1
NumLib::TimeStep ts = alg.getTimeStep();
ASSERT_EQ(1u, ts.steps());
ASSERT_EQ(1., ts.previous());
@@ -37,11 +39,11 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
ASSERT_EQ(1., ts.dt());
ASSERT_TRUE(alg.accepted());
- ASSERT_TRUE(alg.next()); // t=4, dt=2
+ ASSERT_TRUE(alg.next(solution_error)); // t=4, dt=2
// dt*=2
- alg.setNIterations(3);
- ASSERT_TRUE(alg.next()); // t=8, dt=4
+ alg.setIterationNumber(3);
+ ASSERT_TRUE(alg.next(solution_error)); // t=8, dt=4
ts = alg.getTimeStep();
ASSERT_EQ(3u, ts.steps());
ASSERT_EQ(4., ts.previous());
@@ -50,8 +52,8 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
ASSERT_TRUE(alg.accepted());
// dt*=1
- alg.setNIterations(5);
- ASSERT_TRUE(alg.next()); // t=12, dt=4
+ alg.setIterationNumber(5);
+ ASSERT_TRUE(alg.next(solution_error)); // t=12, dt=4
ts = alg.getTimeStep();
ASSERT_EQ(4u, ts.steps());
ASSERT_EQ(8., ts.previous());
@@ -60,8 +62,8 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
ASSERT_TRUE(alg.accepted());
// dt*=0.5
- alg.setNIterations(7);
- ASSERT_TRUE(alg.next()); // t=14, dt=2
+ alg.setIterationNumber(7);
+ ASSERT_TRUE(alg.next(solution_error)); // t=14, dt=2
ts = alg.getTimeStep();
ASSERT_EQ(5u, ts.steps());
ASSERT_EQ(12., ts.previous());
@@ -70,8 +72,8 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
ASSERT_TRUE(alg.accepted());
// dt*=0.25 but dt_min = 1
- alg.setNIterations(8); // exceed max
- ASSERT_TRUE(alg.next()); // t=13, dt=1
+ alg.setIterationNumber(8); // exceed max
+ ASSERT_TRUE(alg.next(solution_error)); // t=13, dt=1
ts = alg.getTimeStep();
ASSERT_EQ(5u, ts.steps());
ASSERT_EQ(12., ts.previous());
@@ -80,8 +82,8 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
ASSERT_FALSE(alg.accepted());
// restart, dt*=1
- alg.setNIterations(4);
- ASSERT_TRUE(alg.next()); // t=14, dt=1
+ alg.setIterationNumber(4);
+ ASSERT_TRUE(alg.next(solution_error)); // t=14, dt=1
ts = alg.getTimeStep();
ASSERT_EQ(6u, ts.steps());
ASSERT_EQ(13., ts.previous());
@@ -114,7 +116,7 @@ TEST(NumLib, TimeSteppingIterationNumberBased2)
{
std::size_t n = (i<_nr_iterations.size()) ? _nr_iterations[i++] : 0;
//INFO("-> NR-iterations=%d", n);
- obj.setNIterations(n);
+ obj.setIterationNumber(n);
}
};
diff --git a/Tests/NumLib/TimeSteppingTestingTools.h b/Tests/NumLib/TimeSteppingTestingTools.h
index 0922d246883df3aa5556ba540dd5a3fa77b5a535..eaeca4df6100f4c60c50f0b6f338e2ae19872537 100644
--- a/Tests/NumLib/TimeSteppingTestingTools.h
+++ b/Tests/NumLib/TimeSteppingTestingTools.h
@@ -31,7 +31,7 @@ std::vector<double> timeStepping(T_TIME_STEPPING &algorithm, T* obj=nullptr)
std::vector<double> vec_t;
vec_t.push_back(algorithm.begin());
- while (algorithm.next()) {
+ while (algorithm.next(0.0)) {
NumLib::TimeStep t = algorithm.getTimeStep();
//INFO("t: n=%d,t=%g,dt=%g", t.steps(), t.current(), t.dt());
if (obj)