Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve DiscreteTime and use it in step 21 #9748

Merged
merged 4 commits into from Apr 23, 2020

Conversation

rezarastak
Copy link
Contributor

@rezarastak rezarastak commented Mar 27, 2020

Builds on #8824

This PR shows how the class DiscreteTime can be used to streamline tracking of time and time-stepping in our simulations. I have modified step-21 to use the basic features of this class. This PR has multiple new method definitions that were not implemented in #8824 .

The users of the library learn the concept of time-stepping first in step-18 using the traditional coding style for time incrementation.

  present_time += present_timestep;
  ++timestep_no;

Then in step-21, they learn how to use DiscreteTime to simplify these operations and to avoid hard-to-find mistakes.

time.advance_time();
std::cout << "   Now at t=" << time.get_current_time()
          << ", dt=" << time.get_previous_step_size() << '.'

Later, in some other example after step 21, I will show the advanced features of DiscreteTime, which I haven't implemented yet.

We can directly merge this PR and close the PR #8824. Alternatively, we can discuss #8824 further and merge it before incorporating this PR. Note that there is still some discussion in #8824 on the exact name for this class and the namespace that holds this class. I hope the content I present here helps us in deciding the naming aspects of this class and its methods.

@rezarastak

This comment has been minimized.

@rezarastak
Copy link
Contributor Author

Now that #8824 is merged, I updated this PR and rebased the commits so that it can be reviewed.

Copy link
Member

@bangerth bangerth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just minor stuff!

I have some changes I want to make to the documentation and suggest a different function name, but will wait until this has been merged.

@@ -0,0 +1,4 @@
Added: New methods added to the DiscreteTime class including get_step_number()
, get_previous_time(), is_at_start(), is_at_end(), get_previous_step_size().
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

put the comma at the end of the previous line

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed

@@ -250,6 +250,9 @@ cell term to get an equation as follows:
(S^n,\sigma)_\Omega +
\triangle t \sum_K \left(F(S^n) q^{n+1}, \sigma\right)_K.
@f}
In order to keep track of the current value of time and the time step in the code,
an object of type DiscreteTime is used. This class encapsulates many complexities
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe more as general advice for writing English language: Try to avoid using the passive voice ("is used"). There are just so many ways in which passive voice can be misunderstood. Just say "we use an object of type DiscreteTime"

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for the advice, I changed some of the sentences to use the active voice.

@@ -1,5 +1,11 @@
<h1>Results</h1>

The code as presented here does actually not compute the results
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
The code as presented here does actually not compute the results
The code as presented here does not actually compute the results

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed

// DiscreteTime::get_previous_step_size() instead of
// DiscreteTime::get_next_step_size(). After many steps, when the simulation
// reaches the end time, the last dt is chosen by the DiscreteTime class in
// way that the last step finishes exactly at the end time.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
// way that the last step finishes exactly at the end time.
// such a way that the last step finishes exactly at the end time.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed

@@ -62,6 +62,18 @@ class DiscreteTime
double
get_current_time() const;

/**
* Return the next time that we reach if we advance time by one step.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
* Return the next time that we reach if we advance time by one step.
* Return the next time that we would reach if we were to advance the time by one step.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you. I fixed it

get_next_time() const;

/**
* Return the previous time.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe say "Return the time we were at before advance_time() was called the last time" to make clear what "previous" corresponds to.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's a good way of saying it. I modified the doc.

@@ -267,9 +267,9 @@
* <td>step-21</td>
* <td> The time dependent two-phase flow in
* porous media. Extensions of mixed Laplace discretizations. More
* complicated block solvers. Simple time stepping.
* complicated block solvers. Simple time stepping using DiscreteTime.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that we have, for the most part, managed to keep the descriptions of the tutorials separated from the classes that have been used in them; we've kept the keywords section for that. I'd be keen to keep this line as it was.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, let's drop this here. But it's worth mentioning in the DiscreteTime class that step-21 uses it!

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I fixed it here and I added a mention to step-21 in discrete_time.h

@@ -0,0 +1,3 @@
Improved: Step-21 adapted the use of the new class DiscreteTime
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Improved: Step-21 adapted the use of the new class DiscreteTime
Improved: Step-21 has been adapted to use the new class DiscreteTime.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed

The code as presented here does actually not compute the results
found on the web page. The reason is, that even on a decent
computer it runs more than a day. If you want to reproduce these
results, modify the end time of the DiscreteTime obejct to 250 within the
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
results, modify the end time of the DiscreteTime obejct to 250 within the
results, modify the end time of the DiscreteTime object to `250` within the

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed

@@ -512,8 +517,7 @@ namespace Step21
1)
, dof_handler(triangulation)
, n_refinement_steps(5)
, time_step(0)
, timestep_number(1)
, time(/*start time*/ 0., /*end time*/ 1., /*time step*/ 0.)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't the start_step_size be non-zero?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How about using the constructor I'm proposing in #9913?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Setting the starting time step to zero is the idiomatic way (for this class) to say that we are going to provide the first time step size later in the code. I added more explanation and information about this in the documentation. I further modified the constructor such that the third parameter is optional and it is automatically set to zero if it is not provided by the user. Even though the class allows having dt ==0, it does not allow incrementing time when dt == 0, forcing the user to provide a positive (nonzero) step size before calling advance_time().

@jppelteret
Copy link
Member

I've checked out this branch at aeb7a6b221 and have written a small test case that tries to present the two time-stepping concepts used in step-21 and step-44 (as I understand them). I've also added an example that adopts the methodology that the documentation suggests (advancing time in the for loop). I've attached the test and output for inspection, and I'll extract a part of the output that highlights both the conceptual difficulties that I'm seeing.

==================================================
Solving with a step-21 type setup.
==================================================

Current time @0 with (forward) dt = 0.1
Solving for equilibrium solution at timestep 0... Done!
Current time @0.1 with (forward) dt = 0.1
Solving for equilibrium solution at timestep 1... Done!
...
Current time @0.9 with (forward) dt = 0.1
Solving for equilibrium solution at timestep 9... Done!

// Current time never hits the end time
// But (current time + forward dt) will hit the end time

==================================================
Solving with a step-44 type setup.
==================================================

Doing some post-processing at initial configuration... // i.e t=0

Current time @0.1 with (backward) dt = 0.1
Solving for equilibrium solution at timestep 1... Done!
Current time @0.2 with (backward) dt = 0.1
Solving for equilibrium solution at timestep 2... Done!
...
Current time @1 with (backward) dt = 0.1
Solving for equilibrium solution at timestep 10... Done!

// Current time hits end time

There also looks to be a discrepancy between what the current documentation explains should happen when the time advancement is done within the "loop-expression" of the for loop. If the time interval is not divided equally by the time increment, then you never actually increment to an end time (although the current time does remain within the designated interval). I think that this is because the is_at_end() check is too strict.

One possible solution would be to let the current time run right up to the end time, but set a flag that indicates that any further time increment would result in the current time exceeding the final time. Then the next time that is_at_end() is called, this flag would indicate that the time-stepping should not occur (i.e. you're done). I think, though, that this might invalidate the modifications proposed for step-21, as it would then want to take another time-step with get_next_time_step() == 0.

==========================================================================================
Investigating how time stepping works for different end points with a fixed step size.
==========================================================================================

--------------------------------------------------
Advance after solve (like step-21).
Settings:   t_begin = 0  dt = 5  t_end = 24
--------------------------------------------------

Current time @0 with (backward) dt = 0 and (forward) dt = 5
Current time @5 with (backward) dt = 5 and (forward) dt = 5
Current time @10 with (backward) dt = 5 and (forward) dt = 5
Current time @15 with (backward) dt = 5 and (forward) dt = 5
Current time @20 with (backward) dt = 5 and (forward) dt = 4

// Current time never hits end time

--------------------------------------------------
Advance before solve (like step-44).
Settings:   t_begin = 0  dt = 5  t_end = 24
--------------------------------------------------

Initial values:   Current time @0 with (backward) dt = 0 and (forward) dt = 5
Current time @5 with (backward) dt = 5 and (forward) dt = 5
Current time @10 with (backward) dt = 5 and (forward) dt = 5
Current time @15 with (backward) dt = 5 and (forward) dt = 5
Current time @20 with (backward) dt = 5 and (forward) dt = 4
Current time @24 with (backward) dt = 4 and (forward) dt = 0

// Current time always hits end time

--------------------------------------------------
Advance in for loop (like documentation).
Settings:   t_begin = 0  dt = 5  t_end = 24
--------------------------------------------------

Current time @0 with (backward) dt = 0 and (forward) dt = 5
Current time @5 with (backward) dt = 5 and (forward) dt = 5
Current time @10 with (backward) dt = 5 and (forward) dt = 5
Current time @15 with (backward) dt = 5 and (forward) dt = 5
Current time @20 with (backward) dt = 5 and (forward) dt = 4

// Current time never hits end time

discrete_time_02.cc.txt
discrete_time_02.output.txt
stdout.txt

@bangerth
Copy link
Member

bangerth commented Apr 9, 2020 via email

@rezarastak
Copy link
Contributor Author

Based on the questions that were asked here and in other places about the meaning of current, next, and previous in the context of incrementing time, I wrote an explanation in the documentation in a section called details of time-stepping and added it to this PR as a commit. Let me know if you agree with my analysis and please provide suggestions on how to improve the explanation.

@jppelteret
Copy link
Member

@rezarastak I went through the time-dependent tutorials to try to better understand why there's a mixture of schemes used to drive the incrementation of time. More often than not it boils down to whether or not the time-stepping strategy is explicit (e.g. a forward Euler method), implicit (e.g. a backward Euler method), or something more complicated. Here's a short summary of what I found:

  • The tutorials that seem share the same time incrementation strategy as step-21 (i.e. incrementing time after solving the timestep) are steps 31, 32, 33, 43 and 67. They are predominantly explicit methods, but steps 33 uses a fractional step approach.
  • The tutorials that seem share the same time incrementation strategy as step-44 (i.e. incrementing time before solving the timestep) are steps 18, 23, 24, 25, 26 and 48. These are all either fully implicit or fractional step methods, with the exception of step 48 which is an explicit method.
  • Steps 52 and 69 do something a little more complicated, as they evolve the timestep.

I've opened #9873 to better document the time-stepping strategies. Ultimately I think that what you've done here is correct, but best fits the "explicit" category of time integrators. My viewpoint has always been that that seems to better suit the second category of time integrators, namely the fully implicit or fractional step methods (although its clear that with some careful consideration one can theoretically setup your problem in either setting). So apologies for getting this wrong -- ideally I should have looked more carefully through the documentation to the problem in the first place.

So I guess it might be fair to say that we've both been approaching this with an incomplete view to how time might need to evolve for any time integration scheme in general. Currently, I see two plausible ways to accommodate both schemes:

  1. Introduce a flag to the constructor that specifies whether or not the last timestep should be one that gives current_time == end_time.
  2. Split this class up into a base class and then two other classes, one that's specialised for "forward" time integration (i.e. one for which you're trying to solve for the solution at current_time + get_next_time_step()), and one that is specialised for "backward" time integration (i.e. one for which you're trying to solve for the solution at current_time with the previous solution being at current_time - get_previous_time_step()).

What do you think?

@jppelteret
Copy link
Member

Let me know if you agree with my analysis and please provide suggestions on how to improve the explanation.

Thank you, I'll take a look at it as soon as I can :-)

Currently, I see two plausible ways to accommodate both schemes:

To qualify, I meant that we could do this in some future PR.

@bangerth
Copy link
Member

bangerth commented Apr 11, 2020 via email

@jppelteret
Copy link
Member

jppelteret commented Apr 11, 2020

I think that it's no more helpful to have a base and two derived classes.

I fully agree. I just wanted to supply what I felt were the plausible options to make this class work with the for-loop style of time-stepping.

What's your sense how hard it may be to unify the tutorial programs on one
style or another?

Just checking: Is this question posed to me?

@bangerth
Copy link
Member

bangerth commented Apr 11, 2020 via email

@jppelteret
Copy link
Member

jppelteret commented Apr 11, 2020

@bangerth If I were to put a number to it, I think about 80% of them seem to follow roughly the same recipe (of which there are variants for the explicit and implicit time stepping / solving strategies). But if we get a nice solution for both of these setups individually then I think that it'd be easy to apply to most of the tutorials. In some cases there's some complicated logic mixed into the time loops, so it might not be possible to do clean for(DiscreteTime time (...); !time.at_end(); time.advance()){...} without changing (at the very least) what's printed to screen. But the proposed modifications for step-21 pretty much apply directly to the other explicit cases, and I think that for(DiscreteTime time (...); !time.at_end(); ){... time.advance(); ...} can be used for the majority of the implicit and fractional step time discretisations. The Runga-Kutta tutorials are a bit more involved, so that might take some additional thought. But I didn't look into them too deeply.

Copy link
Member

@jppelteret jppelteret left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rezarastak I think that the "details of time-stepping" section that you've provided is very useful and insightful. Thanks very much for taking the time to think about it and make such a carefully structured write-up on the topic. I think that everyone will benefit from having it.

Apart from the final example, I'm quite happy with this. I've left some comments suggesting how the text might be improved, just adding detail where appropriate and restructuring a few sentences to prevent double word use or improve readability.

include/deal.II/base/discrete_time.h Outdated Show resolved Hide resolved
include/deal.II/base/discrete_time.h Outdated Show resolved Hide resolved
include/deal.II/base/discrete_time.h Outdated Show resolved Hide resolved
include/deal.II/base/discrete_time.h Outdated Show resolved Hide resolved
include/deal.II/base/discrete_time.h Outdated Show resolved Hide resolved
include/deal.II/base/discrete_time.h Outdated Show resolved Hide resolved
include/deal.II/base/discrete_time.h Outdated Show resolved Hide resolved
include/deal.II/base/discrete_time.h Outdated Show resolved Hide resolved
include/deal.II/base/discrete_time.h Outdated Show resolved Hide resolved
include/deal.II/base/discrete_time.h Show resolved Hide resolved
Copy link
Member

@bangerth bangerth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you address the comments/suggestions I and @jppelteret left, I'll be happy to take another look!

@@ -267,9 +267,9 @@
* <td>step-21</td>
* <td> The time dependent two-phase flow in
* porous media. Extensions of mixed Laplace discretizations. More
* complicated block solvers. Simple time stepping.
* complicated block solvers. Simple time stepping using DiscreteTime.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, let's drop this here. But it's worth mentioning in the DiscreteTime class that step-21 uses it!

@@ -512,8 +517,7 @@ namespace Step21
1)
, dof_handler(triangulation)
, n_refinement_steps(5)
, time_step(0)
, timestep_number(1)
, time(/*start time*/ 0., /*end time*/ 1., /*time step*/ 0.)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How about using the constructor I'm proposing in #9913?

@rezarastak
Copy link
Contributor Author

@jppelteret @bangerth I incorporated all of your comments. Thank you for the enormous amount of time you provided to review this code. Based on the comments, in addition to changes to the documentation, I made two important modifications.

  • As @bangerth previously suggested, I fully renamed set_next_step_size() to set_desired_next_step_size(). I thought about creating a separate PR for it but I decided to include it in this rather large PR. let me know if you prefer to separate it.
  • Modified the constructor to make its third parameter optional. This change effectively achieves the same purpose as Add a second constructor to DiscreteTime. #9913.

Sorry for having a large number of commits in this PR. I rebased and squashed them in a way that each commit is a distinct, self-contained change.

@rezarastak
Copy link
Contributor Author

I want to specifically thank @jppelteret for reading all related tutorial steps and providing very useful feedback with comprehensive example codes and analysis of the code results. I read all of your comments and it helped me write the section titled details of time-stepping in the documentation.

Here I want to talk about using loop structures (for loop, while loop) to control DiscreteTime in an elegant and correct way.

The simplest coding style is the for loop that is presented at the beginning of the documentation of DiscreteTime.

for (DiscreteTime time(0, 1, 0.1); time.is_at_end() == false; time.advance_time())
{
  do_all_stuff();
}

But as we all know, this simplified code does not work for many scenarios. In a real simulation, as I explained, each operation can be categorize to be part of the snapshot stage or the update stage. The problem is that the snapshot stage shouldn't be run for the same number of times as the update stage. The snapshot stage should be run one more time. For example, if our simulation goes from 0 to 1 with dt = 0.1, then we have 10 distinct update stages (0 to 0.1, 0.1 to 0.2, ..., 0.9 to 1). However, there are 11 snapshot stages (0, 0.1, 0.2, ..., 0.9, 1). In many of our tutorials, we ignore one of the snapshot stages, for example, we don't export vtu files for t = 0. In reality, I think we should always call all n+1 snapshot stages.

There are multiple ways to approach this using loops:
Code style 1: elegant but has unwanted code duplication

make_grid(); //setup
setup_system(); //setup
for (DiscreteTime time(0, 1, 0.1); //setup
     time.is_at_end() == false;
     time.advance_time()) // update
{
  print_results(); // snapshot
  output_results(); // snapshot
  assemble_system(); // update
  solve(); // update
  update_solutions(); // update
}
// Now we have reached end time, we have to run one more snapshot stage
print_results(); // snapshot
output_results(); // snapshot

Code style 2: Has to use code duplication or ignore the output at t = 0

make_grid(); //setup
setup_system(); //setup
DiscreteTime time(0, 1, 0.1); //setup
// Before entering the for loop call the output routines for t = 0
// These might be ignored by some programmers because they are not interested in t = 0 
print_results(); // snapshot
output_results(); // snapshot
for (; time.is_at_end() == false;)
{
  assemble_system(); // update
  solve(); // update
  update_solutions(); // update
  time.advance_time(); // update
  print_results(); // snapshot
  output_results(); // snapshot
}

As you saw, none of these styles are perfect. I think the best way is to use a for loop or while loop with a break statement in the middle. It does not have code duplication, it does not require ignoring the output at t = 0, but it is not elegant because it uses a break statement in the middle of the loop.

make_grid(); //setup
setup_system(); //setup
for (DiscreteTime time(0, 1, 0.1); //setup
    ;
    time.advance_time()) // update
{
  print_results(); // snapshot
  output_results(); // snapshot
  if (time.is_at_end() == false)
  {
    break;
  }
  assemble_system(); // update
  solve(); // update
  update_solutions(); // update
}

The last style is the best style in my opinion. But the tutorials don't follow it. Therefore I didn't write it in the documentation. I think for each tutorial we have to see how they can be used with DiscreteTime and choose an approach that fits that specific tutorial. Users can also freely choose their own style. My goal is to keep DiscretTime general enough in order that it can be used with any coding style.

@rezarastak
Copy link
Contributor Author

One separate note about the loop style: I believe it is important that the function calls for the snapshot stage does not interleave with the calls to update functions. The code snippet that I used in the documentation of DiscreteTime and is based on the use cases in the tutorial, unfortunately mixes these snapshot and update operations.

make_grid(); // setup
setup_system(); // setup
for (DiscreteTime time(0., 1., 0.1);  // setup
     time.is_at_end() == false;
     time.advance_time())             // update
{
  const double time_of_simulation = time.get_next_time();
  const double timestep_size      = time.get_next_step_size();

   std::cout // snapshot
     << "Timestep: " << time.get_step_number() << " -- "
     << "Solving for equilibrium solution at "
     << "t = " << time_of_solution << " with "
     << "dt = " << timestep_size << "." << std::endl;

  assemble_system(time_of_solution, timestep_size); // update
  solve(); // update
  update_solutions(); // update

   output_results(time_of_solution); // snapshot

   // propose a new timestep size if need be
   // time.set_desired_next_step_size(...); // snapshot
}

For example, output_results which is a snapshot operation is called after update_solutions i.e. an update operation. However, after output_results we go back and call time.advance_time() which is an update operation again and then we call std::cout << to do a snapshot operation. As you can see, it is interleaving the two types of operations. I put this code in the documentation and I think it is very useful for the users, but I wanted to comment on some of the shortcomings of this code snippet.

@bangerth
Copy link
Member

bangerth commented Apr 20, 2020 via email

Copy link
Member

@bangerth bangerth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like this a lot. I've only got a couple of comments and would have just merged but want to give @jppelteret time to also approve.

@@ -0,0 +1,4 @@
Modified: The method DiscreteTime::set_next_step_size() is renamed to
set_desired_step_size() to better reflect its purpose and behavior.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just drop this and the next changelog file. These changelog files are meant for users to peruse for each release, and so they only need to describe things that changed from one release to the next.

@@ -0,0 +1,3 @@
Improved: Step-21 has been adapted to use the new class DiscreteTime.
<br>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This one should stay, however.

// inside the loop. Since we are reporting the time and dt after we
// increment it, we have to call the method
// DiscreteTime::get_previous_step_size() instead of
// DiscreteTime::get_next_step_size(). After many steps, when the simulation
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's just not output the time step size at all -- then you don't need this explanation :-)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that removing the dt from text output helps with the explanation here. However, because this tutorial calculates a new value of dt every step and describes this adaptive time-stepping scheme, it seems very useful to keep printing dt values.

std::cout << " Now at t=" << time << ", dt=" << time_step << '.'
time.advance_time();
std::cout << " Now at t=" << time.get_current_time()
<< ", dt=" << time.get_previous_step_size() << '.'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

let's just drop this line

* 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 the get updated one by one. In general, the order of
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
* not synchronized and the get updated one by one. In general, the order of
* not synchronized and then get updated one by one. In general, the order of

* // pre-processing/setup stage {
* make_grid();
* setup_system();
* for (DiscreteTime time(0., 1., 0.1); // } end pre-processing/setup stage
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
* for (DiscreteTime time(0., 1., 0.1); // } end pre-processing/setup stage
* for (DiscreteTime time(0., 1., 0.1); // end pre-processing/setup stage

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems to be the end brace to pre-processing/setup stage {. So it's correct if all the bracing is kept.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since it is consistent with the way I denoted the start and end of each stage in // ... comments, I am going to keep it as it is now. I hope @bangerth is okay with it.

* 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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
* time.advance_time()) // part of the update stage
* time.advance_time()) // part of the update stage; happens at the end of the loop body's execution

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, I'm not too fond of writing the loop like this. I'd much rather the time.advance_time() call was explicitly put at the end of the loop body, and the initialization of the time variable moved before the for loop. The for loop would then simply become a while loop. I think that would be easier to read.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm... maybe you and @rezarastak need to negotiate on this. It seems to have been the purpose of the whole discussion, and is also aligned with the existing documentation in the class intro.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@bangerth What you've stated here is what I've always done timestepping in the past, but I'm trying to come around to @rezarastak's way of doing things. Clearly there's merit and a reason for using any of the permissible permutations of constructing a time-stepping loop. I like that the for loop is structured and self-contained (i.e. drives its own incrementation). In its simplest form it encapsulates, without the possibility for misinterpretation, the full set of operations required from the start of simulation to the end of simulation.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@bangerth @rezarastak Would you feel it to be of value if we write the equivalent loop using a while loop as an additional example? If so then I can add it as a follow up to this PR.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am happy with the current for loop style and I think users who read it can follow the order of operations happening within the loop. Using for loop is good for consistency with the beginning of the tutorial as @jppelteret said. Regardless, the users of this class can rewrite it in while loop style. The style preference is very subjective.

I can add to the documentation that this loop can also be written in while format and do { ... } while format as shown in step-21.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Leave it as you had it and I can follow up with a patch that shows an alternative. I should give you some artistic license :-)

*
* std::cout
* << "Timestep: " << time.get_step_number() << " -- "
* << "Solving for equilibrium solution at "
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
* << "Solving for equilibrium solution at "
* << "Solving for the solution at "

*
* @param[in] desired_start_step_size A desired step size for incrementing
* time for the first step. It is not guaranteed that this value will be
* actually used as the size of the first step.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
* actually used as the size of the first step.
* actually used as the size of the first step, as discussed in the introduction.

get_previous_step_size() const;

/**
* Return the number of times the simulation time is incremented. Return
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
* Return the number of times the simulation time is incremented. Return
* Return the number of times the simulation time has been incremented. Return

@jppelteret
Copy link
Member

jppelteret commented Apr 22, 2020

@rezarastak, thank you for giving my comments consideration. I'm struggling for time, and would only be able to give this another detailed look on Sunday. I don't really think that its fair for me to hold this up so I'm going to dismiss my review and let @bangerth merge it if he's happy with the changes. I trust his and your judgements on what's been done here. May I just ask that you look once more through the comments - from a cursory glance it looks to me like there are maybe one or two more that @bangerth made in the last couple of days that have not been addressed. (Maybe you were waiting for me. If so, then sorry!) If you can just indicate that you've seen all of them, even if you don't want to actually include / address the comment. Then we all know where we stand on all points by the time this is merged.

EDIT: I've staved off the need for sleep and have looked through it again in a very nonlinear manner. I've also marked all of my comments that have been attended to as resolved.

@jppelteret jppelteret dismissed their stale review April 22, 2020 20:38

In the interest of fairness, I should not block this PR any further.

* std::cout
* << "Timestep: " << time.get_step_number() << " -- "
* << "Solving for the solution at "
* << "t = " << time_of_solution << " with "
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think that time_of_solution is defined. I guess its time_of_simulation?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that is a typo. Thanks

@jppelteret
Copy link
Member

jppelteret commented Apr 22, 2020

My goal is to keep DiscretTime general enough in order that it can be used with any coding style.

I think that this is the best approach, and that we should not force our users to use the class in any particular way. You've certainly done a good job to make sure that this is possible. As I said before, I do though quite strongly believe that we should just be careful not to recommend any structure or scheme that's abnormal or confusing. That's why I made the recommendations that I did.

The code snippet that I used in the documentation of DiscreteTime and is based on the use cases in the tutorial, unfortunately mixes these snapshot and update operations... For example, output_results which is a snapshot operation is called after update_solutions i.e. an update operation. However, after output_results we go back and call time.advance_time() which is an update operation again and then we call std::cout << to do a snapshot operation. As you can see, it is interleaving the two types of operations.

The funny thing about programming is that it allows one to relax on correctness and inject a bit of pragmatism where required. One could write the loop like this

make_grid(); // setup
setup_system(); // setup
for (DiscreteTime time(0., 1., 0.1);  // setup
     time.is_at_end() == false;
     time.advance_time())             // update
{
  const double time_of_simulation = time.get_next_time();
  const double timestep_size      = time.get_next_step_size();

   std::cout // snapshot
     << "Timestep: " << time.get_step_number() << " -- "
     << "Solving for equilibrium solution at "
     << "t = " << time_of_solution << " with "
     << "dt = " << timestep_size << "." << std::endl;

  assemble_system(time_of_solution, timestep_size); // update
  solve(); // update

// ========================
   output_results(time_of_solution, solution_t1, delta_solution); // snapshot
   // propose a new timestep size if need be
   // time.set_desired_next_step_size(...); // snapshot
// ========================

   update_solutions(); // update
}

with

void output_results(time_of_solution, solution_t1, delta_solution)
{
  // Local to this function
  Vector solution = solution_t1 + delta_solution;
  ....
}

and

void update_solutions()
{
  // Class-owned vector
  solution = solution_t1 + delta_solution;
}

Then I believe there is no interleaving of the snapshot and update stages.

So you win somewhere and lose somewhere else. Here we've duplicated a vector accumulation, which may require MPI communication. In the other examples, one duplicates a single line of code to output at time zero / end time. One just needs to find the balance that suits oneself, irrespective of whether or not it fits into this delineation of updates, snapshots at the like. I'm certainly going to use concept that you've laid out as a benchmark for whether or not the loop structure is sensible, but I don't think that it is a hard and fast rule that has to be followed.

To put the postprocessing first seems like the wrong approach.

I fully agree. To me it feels like the order of operations are wrong, and then I'd end up spending a lot of time puzzling it over in my head to make sure that all the calls happen in the right sequence.

Copy link
Member

@jppelteret jppelteret left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that there are a couple of minor things outstanding, but it looks great otherwise. Thanks for all of the hard work @rezarastak!

@jppelteret
Copy link
Member

/rebuild

@rezarastak
Copy link
Contributor Author

May I just ask that you look once more through the comments - from a cursory glance it looks to me like there are maybe one or two more that @bangerth made in the last couple of days that have not been addressed. (Maybe you were waiting for me. If so, then sorry!) If you can just indicate that you've seen all of them, even if you don't want to actually include / address the comment. Then we all know where we stand on all points by the time this is merged.

Although I addressed some of @bangerth's comments already, I was waiting for your feedback before working on the rest of his comments. Thank you so much @jppelteret for the review. I will work on it and finalize this PR soon.

@bangerth
Copy link
Member

Let us know when you think you're ready and we can merge!

@rezarastak
Copy link
Contributor Author

Thank you. It is ready to me merged in my opinion. I think I have addressed all of the critical comments.

@bangerth
Copy link
Member

Let's go with it then! Thanks for all of your work on this, I think this is addressing a need that none of us has recognized in 20 years but that really helps with making things more systematic -- it's a collective failure of imagination you demonstrate with this patch, and I'm glad to be working in a project where we can have people come in with a different perspective and make things better! Thank you very much!

@bangerth bangerth merged commit dc8d05e into dealii:master Apr 23, 2020
@rezarastak
Copy link
Contributor Author

Thank you! I am glad this class is available in the next release of deal.II. Hopefully people can use it and provide feedback which will help us improve our approach to adaptive time-stepping.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants