From 33ae49bc2bee39ba0d247bc670d289f1efcb759e Mon Sep 17 00:00:00 2001 From: Reza Rastak Date: Wed, 25 Mar 2020 22:52:43 -0700 Subject: [PATCH] New const member functions added to DiscreteTime --- include/deal.II/base/discrete_time.h | 266 ++++++++++++++++++++++++++- source/base/discrete_time.cc | 22 ++- tests/base/discrete_time_1.cc | 11 +- tests/base/discrete_time_1.output | 47 ++--- 4 files changed, 307 insertions(+), 39 deletions(-) diff --git a/include/deal.II/base/discrete_time.h b/include/deal.II/base/discrete_time.h index 3427ed7f03..de7fa92bc2 100644 --- a/include/deal.II/base/discrete_time.h +++ b/include/deal.II/base/discrete_time.h @@ -24,9 +24,16 @@ DEAL_II_NAMESPACE_OPEN * This class provides a means to keep track of the simulation time in a * time-dependent simulation. It manages stepping forward from a start time * $T_{\text{start}}$ to an end time $T_{\text{end}}$. It also allows adjusting - * the time step size during the simulation. It is guaranteed that at all times - * the current simulation time is in the closed interval between the start time - * and the end time. + * the time step size during the simulation. + * + * This class provides a number of invariants that are guaranteed to be + * true at all times. + * + * * The current simulation time is within the closed interval between the + * start time and the end time ($T_{\text{start}} \le t \le T_{\text{end}}$). + * * Whenever time is incremented, the step size is positive ($dt > 0$). + * In other words, time advances in strictly ascending order + * ($m < n \Leftrightarrow t_m < t_n$). * * The model this class follows is that one sets a *desired* time step length * either through the constructor or using set_desired_next_step_size() @@ -38,7 +45,7 @@ DEAL_II_NAMESPACE_OPEN * Let's say that you loop over all of the time steps by using a for loop * @code * for (DiscreteTime time(0., 1., 0.3); - * time.get_current_time() != time.get_end_time(); + * time.is_at_end() == false; * time.advance_time()) * { * // Insert simulation code here @@ -63,7 +70,7 @@ DEAL_II_NAMESPACE_OPEN * end time: * @code * for (DiscreteTime time(0., 1.21, 0.3); - * time.get_current_time() != time.get_end_time(); + * time.is_at_end() == false; * time.advance_time()) * { * // Insert simulation code here @@ -78,13 +85,143 @@ DEAL_II_NAMESPACE_OPEN * only a *desired* step size. You can query the actual time step size using the * get_next_step_size() function. * + * + * ### Details of time-stepping + * + * Since time is marched forward in a discrete manner in our simulations, we + * need to discuss how we increment time. During time stepping we enter two + * separate alternating regimes in every step. + * + * * The **snapshot** stage (the **current** stage, the **consistent** + * stage): In this part of the algorithm, we are at $t = t_n$ and all + * quantities of the simulation (displacements, strains, temperatures, etc.) + * are up-to-date for $t = t_n$. In this stage, *current time* refers to + * $t_n$, *next time* refers to $t_{n+1}$, *previous time* refers to + * $t_{n-1}$. The other useful notation quantities are the *next* time step + * size $t_{n+1} - t_n$ and *previous* time step size $t_n - t_{n-1}$. In + * this stage, it is a perfect occasion to generate text output using print + * commands within the user's code. Additionally, post-processed outputs can + * be prepared here which can be viewed later by visualization programs such + * as `Tecplot`, `Paraview`, and `VisIt`. Additionally, during the snapshot + * stage, the code can assess the quality of the previous step and decide + * whether it wants to increase or decrease the time step size. The step + * size for the next time step can be modified here. + * * The **update** stage (the **transition** stage, the **inconsistent** + * stage): In this section of the program, the internal state of the + * simulation is getting updated from $t_n$ to $t_{n+1}$. All of the + * variables need to be updated one by one, the step number is incremented, + * the time is incremented by $dt = t_{n+1} - t_n$, and time-integration + * algorithms are used to update the other simulation quantities. In the + * middle of this stage, some variables have been updated to $t_{n+1}$ but + * other variables still represent their value at $t_n$. Thus, we call this + * the inconsistent stage, requiring that no post-processing output related + * to the state variables take place within it. The state variables, namely + * those related to time, the solution field and any internal variables, are + * not synchronized and then get updated one by one. In general, the order of + * updating variables is arbitrary, but some care should be taken if there + * are interdependencies between them. For example, if some variable such as + * $x$ depends on the calculation of another variable such as $y$, then $y$ + * must be updated before $x$ can be updated. + * + * The question arises whether time should be incremented before updating + * state quantities. Multiple possibilities exist, depending on program and + * formulation requirements, and possibly the programmer's preferences: + * * Time is incremented before the rest of the updates. In this case, even + * though time is incremented to $t_{n+1}$, not all variables are updated + * yet. During this update phase, $dt$ equals the *previous* time step + * size. *Previous* means that it is referring to the $dt$ of the + * `advance_time()` command that was performed previously. In the + * following example code, we are assuming that `a` and `b` are two state + * variables that need to be updated in this time step. + * @code + * time.advance_time(); + * new_a = update_a(a, b, time.get_previous_step_size()); + * b = update_b(a, b, time.get_previous_step_size()); + * a = new_a; + * @endcode + * * Time is incremented from $t_n$ to $t_{n+1}$ after all variables have + * already been updated for $t_{n+1}$. During the update stage, $dt$ is + * denoted as the *next* time step size. *Next* means that $dt$ of the + * step corresponds to the `advance_time()` command that will happen + * subsequently. + * @code + * new_a = update_a(a, b, time.get_next_step_size()); + * b = update_b(a, b, time.get_next_step_size()); + * a = new_a; + * time.advance_time(); + * @endcode + * * Time is incremented in the middle of the other updates: In this case + * $dt$ would correspond to *next* or *previous* depending of whether it + * is used before or after the call to `advance_time()`. + * @code + * new_a = update_a(a, b, time.get_next_step_size()); + * time.advance_time(); + * b = update_b(a, b, time.get_previous_step_size()); + * a = new_a; + * @endcode + * + * One thing to note is that, during the update phase, $dt$ is referred to + * either **next** or **previous** time step size, depending on whether the + * command `advance_time()` has been called yet. The notion of *current* time + * step size is ill-defined. In fact, in the update stage the definition of + * every variable depends on whether it has been updated yet or not, hence the + * name **the inconsistent stage**. + * + * The following code snippet shows the code sections for the snapshot stage + * and the update stage in the context of a complete time-dependent + * simulation. This code follows the coding conventions incorporated in the + * tutorial examples. Note that even though this example is written in the + * format of a `for` loop, it can equivalently be written as a `while` or + * `do while` loop (as shown in step-21). + * @code + * // pre-processing/setup stage { + * make_grid(); + * setup_system(); + * for (DiscreteTime time(0., 1., 0.1); // } end pre-processing/setup stage + * time.is_at_end() == false; + * time.advance_time()) // part of the update stage, runs at + * // the end of loop body + * { + * // snapshot stage { + * const double time_of_simulation = time.get_next_time(); + * const double timestep_size = time.get_next_step_size(); + * + * std::cout + * << "Timestep: " << time.get_step_number() << " -- " + * << "Solving for the solution at " + * << "t = " << time_of_simulation << " with " + * << "dt = " << timestep_size << "." << std::endl; + * // } end snapshot stage + * + * // update stage { + * assemble_system(time_of_simulation, timestep_size); + * solve(); + * update_solutions(); + * // } end update stage + * + * // snapshot stage { + * output_results(time_of_solution); + * + * // propose a new timestep size if need be + * // time.set_desired_next_step_size(...); + * // } end snapshot stage + * } + * @endcode + * * @author Reza Rastak, 2019 */ class DiscreteTime { public: /** - * Constructor + * Constructor. + * + * @pre @p start_step_size must be non-negative. + * + * @note If @p start_step_size is specified as zero, it indicates that the + * desired size for the time step will be calculated at a different location + * in the code. In this case, the created object cannot increment time until + * the step size is changed by calling set_desired_next_step_size(). */ DiscreteTime(const double start_time, const double end_time, @@ -96,6 +233,25 @@ public: double get_current_time() const; + /** + * Return the next time that we would reach if we were to advance the time + * by one step. + * + * @note If the simulation is at the end time, this method returns the + * end time. + */ + double + get_next_time() const; + + /** + * Return the time we were at before `advance_time()` was called last time. + * + * @note If the simulation is at the start time, this method returns the + * start time. + */ + double + get_previous_time() const; + /** * Return the start time. */ @@ -105,23 +261,52 @@ public: /** * Return the end of the time interval. * The final time step ends exactly at this point. This exact floating-point - * equality is very important because it allows us to use the expression - * time.get_current_time() != time.get_end_time() as the - * conditional statement in a for loop to check if the end time is reached. + * equality is very important because it allows us to equality-compare + * current time with end time and decide whether we have reached the end of + * the simulation. */ double get_end_time() const; + /** + * Return whether no step has taken place yet. + */ + bool + is_at_start() const; + + /** + * Return whether time has reached the end time. + */ + bool + is_at_end() const; + /** * Return the size of the step from current time step to the * next. As discussed in the introduction to the class, this is the * *actual* time step, and may differ from the *desired* time step * set in the constructor or through the * set_desired_next_step_size() function. + * + * @note If the simulation is at the end time, this method returns zero. */ double get_next_step_size() const; + /** + * Return the step size of the previous step. + * + * @note If the simulation is at the start time, this method returns zero. + */ + double + get_previous_step_size() const; + + /** + * Return the number of times the simulation time has been incremented. + * Return zero when the simulation is at the start time. + */ + unsigned int + get_step_number() const; + /** * Set the value of the next time step size. The next time advance_time() * is called, the newly set @p time_step_size will be used to advance @@ -148,6 +333,10 @@ public: * advance time if it is already at the end time. This rule is created to * avoid the creation of an infinite loop when advance_time() is called * inside a loop. + * + * @pre The time step size must be nonzero. If the step size is currently + * zero, change it by calling set_desired_next_step_size() before calling + * advance_time(). */ void advance_time(); @@ -192,6 +381,17 @@ private: * floating-point value of the time exactly matches the end time. */ double next_time; + + /** + * The previous time. + */ + double previous_time; + + /** + * The step number i.e. the number of times the simulation time ha been + * incremented. + */ + unsigned int step_number; }; @@ -214,6 +414,22 @@ DiscreteTime::get_end_time() const +inline bool +DiscreteTime::is_at_start() const +{ + return step_number == 0; +} + + + +inline bool +DiscreteTime::is_at_end() const +{ + return current_time == end_time; +} + + + inline double DiscreteTime::get_next_step_size() const { @@ -222,6 +438,14 @@ DiscreteTime::get_next_step_size() const +inline double +DiscreteTime::get_previous_step_size() const +{ + return current_time - previous_time; +} + + + inline double DiscreteTime::get_current_time() const { @@ -229,6 +453,30 @@ DiscreteTime::get_current_time() const } + +inline double +DiscreteTime::get_next_time() const +{ + return next_time; +} + + + +inline double +DiscreteTime::get_previous_time() const +{ + return previous_time; +} + + + +inline unsigned int +DiscreteTime::get_step_number() const +{ + return step_number; +} + + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/source/base/discrete_time.cc b/source/base/discrete_time.cc index d45141f69b..f32d6eee4e 100644 --- a/source/base/discrete_time.cc +++ b/source/base/discrete_time.cc @@ -25,9 +25,9 @@ namespace // - The next time exceeds the end time. // - The next time is smaller but very close to the end time. double - get_next_time(const double current_time, - const double step_size, - const double end_time) + calculate_next_time(const double current_time, + const double step_size, + const double end_time) { Assert(step_size >= 0., ExcMessage("Time step size must be non-negative")); Assert(end_time >= current_time, ExcInternalError()); @@ -49,7 +49,9 @@ DiscreteTime::DiscreteTime(const double start_time, , end_time{end_time} , start_step_size{start_step_size} , current_time{start_time} - , next_time{get_next_time(start_time, start_step_size, end_time)} + , next_time{calculate_next_time(start_time, start_step_size, end_time)} + , previous_time{start_time} + , step_number{0} {} @@ -57,7 +59,7 @@ DiscreteTime::DiscreteTime(const double start_time, void DiscreteTime::set_next_step_size(const double next_step_size) { - next_time = get_next_time(current_time, next_step_size, end_time); + next_time = calculate_next_time(current_time, next_step_size, end_time); } @@ -70,8 +72,10 @@ DiscreteTime::advance_time() "You can't advance time further." "Either dt == 0 or you are at the end of the simulation time.")); const double step_size = get_next_step_size(); + previous_time = current_time; current_time = next_time; - next_time = get_next_time(current_time, step_size, end_time); + ++step_number; + next_time = calculate_next_time(current_time, step_size, end_time); } @@ -79,8 +83,10 @@ DiscreteTime::advance_time() void DiscreteTime::restart() { - current_time = start_time; - next_time = get_next_time(current_time, start_step_size, end_time); + previous_time = start_time; + current_time = start_time; + next_time = calculate_next_time(current_time, start_step_size, end_time); + step_number = 0; } DEAL_II_NAMESPACE_CLOSE diff --git a/tests/base/discrete_time_1.cc b/tests/base/discrete_time_1.cc index 363fbd4a06..39d3e13f1f 100644 --- a/tests/base/discrete_time_1.cc +++ b/tests/base/discrete_time_1.cc @@ -22,8 +22,17 @@ void print_time(const DiscreteTime &time) { + if (time.is_at_start()) + deallog << "Simulation started" << std::endl; deallog << "Current time = " << time.get_current_time() - << ", next step size = " << time.get_next_step_size() << std::endl; + << ", next = " << time.get_next_time() + << ", previous = " << time.get_previous_time() + << ", step number = " << time.get_step_number() + << ", next step size = " << time.get_next_step_size() + << ", previous step size = " << time.get_previous_step_size() + << std::endl; + if (time.is_at_end()) + deallog << "Simulation ended" << std::endl; } void diff --git a/tests/base/discrete_time_1.output b/tests/base/discrete_time_1.output index 24806bc48e..6110561f23 100644 --- a/tests/base/discrete_time_1.output +++ b/tests/base/discrete_time_1.output @@ -1,27 +1,32 @@ DEAL:Start to end::Start time = 0.00000 DEAL:Start to end::End time = 1.50000 -DEAL:Start to end::Current time = 0.00000, next step size = 0.123000 -DEAL:Start to end::Current time = 0.123000, next step size = 0.123000 -DEAL:Start to end::Current time = 0.246000, next step size = 0.123000 -DEAL:Start to end::Current time = 0.369000, next step size = 0.123000 -DEAL:Start to end::Current time = 0.492000, next step size = 0.123000 -DEAL:Start to end::Current time = 0.615000, next step size = 0.123000 -DEAL:Start to end::Current time = 0.738000, next step size = 0.123000 -DEAL:Start to end::Current time = 0.861000, next step size = 0.123000 -DEAL:Start to end::Current time = 0.984000, next step size = 0.123000 -DEAL:Start to end::Current time = 1.10700, next step size = 0.123000 -DEAL:Start to end::Current time = 1.23000, next step size = 0.123000 -DEAL:Start to end::Current time = 1.35300, next step size = 0.123000 -DEAL:Start to end::Current time = 1.47600, next step size = 0.0240000 -DEAL:Start to end::Current time = 1.50000, next step size = 0.00000 +DEAL:Start to end::Simulation started +DEAL:Start to end::Current time = 0.00000, next = 0.123000, previous = 0.00000, step number = 0, next step size = 0.123000, previous step size = 0.00000 +DEAL:Start to end::Current time = 0.123000, next = 0.246000, previous = 0.00000, step number = 1, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 0.246000, next = 0.369000, previous = 0.123000, step number = 2, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 0.369000, next = 0.492000, previous = 0.246000, step number = 3, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 0.492000, next = 0.615000, previous = 0.369000, step number = 4, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 0.615000, next = 0.738000, previous = 0.492000, step number = 5, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 0.738000, next = 0.861000, previous = 0.615000, step number = 6, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 0.861000, next = 0.984000, previous = 0.738000, step number = 7, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 0.984000, next = 1.10700, previous = 0.861000, step number = 8, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 1.10700, next = 1.23000, previous = 0.984000, step number = 9, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 1.23000, next = 1.35300, previous = 1.10700, step number = 10, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 1.35300, next = 1.47600, previous = 1.23000, step number = 11, next step size = 0.123000, previous step size = 0.123000 +DEAL:Start to end::Current time = 1.47600, next = 1.50000, previous = 1.35300, step number = 12, next step size = 0.0240000, previous step size = 0.123000 +DEAL:Start to end::Current time = 1.50000, next = 1.50000, previous = 1.47600, step number = 13, next step size = 0.00000, previous step size = 0.0240000 +DEAL:Start to end::Simulation ended DEAL:Start to end::Restarted -DEAL:Start to end::Current time = 0.00000, next step size = 0.123000 +DEAL:Start to end::Simulation started +DEAL:Start to end::Current time = 0.00000, next = 0.123000, previous = 0.00000, step number = 0, next step size = 0.123000, previous step size = 0.00000 DEAL:Start to end::OK -DEAL:Adjust time step size::Current time = 0.400000, next step size = 0.150000 -DEAL:Adjust time step size::Current time = 0.550000, next step size = 0.150000 -DEAL:Adjust time step size::Current time = 0.910000, next step size = 0.360000 -DEAL:Adjust time step size::Current time = 1.27000, next step size = 0.360000 -DEAL:Adjust time step size::Current time = 1.88000, next step size = 0.220000 -DEAL:Adjust time step size::Current time = 2.10000, next step size = 0.00000 +DEAL:Adjust time step size::Simulation started +DEAL:Adjust time step size::Current time = 0.400000, next = 0.550000, previous = 0.400000, step number = 0, next step size = 0.150000, previous step size = 0.00000 +DEAL:Adjust time step size::Current time = 0.550000, next = 0.700000, previous = 0.400000, step number = 1, next step size = 0.150000, previous step size = 0.150000 +DEAL:Adjust time step size::Current time = 0.910000, next = 1.27000, previous = 0.550000, step number = 2, next step size = 0.360000, previous step size = 0.360000 +DEAL:Adjust time step size::Current time = 1.27000, next = 1.63000, previous = 0.910000, step number = 3, next step size = 0.360000, previous step size = 0.360000 +DEAL:Adjust time step size::Current time = 1.88000, next = 2.10000, previous = 1.27000, step number = 4, next step size = 0.220000, previous step size = 0.610000 +DEAL:Adjust time step size::Current time = 2.10000, next = 2.10000, previous = 1.88000, step number = 5, next step size = 0.00000, previous step size = 0.220000 +DEAL:Adjust time step size::Simulation ended DEAL:Adjust time step size::OK -- 2.39.5