8

I recently inherited a large body of legacy code that solves a very stiff, transient problem. I would like to demonstrate that the spatial and temporal step sizes are small enough that the qualitative nature of the computed solution would not change if they were decreased. In other words, I'd like to show that the solution is "converged" in a qualitative sense. Since I can explicitly set the spatial mesh size, that part is straightforward. However, as the code uses automatic time-step size control, I cannot set the time step size directly.

The algorithm changes the time step between two bounds based on the number of jacobian iterations needed to reach an error tolerance during the last $n$ time steps. The fact that it uses jacobian iteration makes me fairly certain that it is some sort of implicit scheme, but I cannot be absolutely certain. It does not account for the error it is seeing in the current time step, which leads to it running into the iteration limit on occasion (maybe a dozen times over the course of several thousand time steps, almost always during the most dynamic portions of the simulation). The current runs I am completing I am setting the time-step bounds two and a half orders of magnitude apart ($10^{-13}$ to $5 \cdot 10^{-11}$).

In the runs, I have control over the time-step bounds, the number of past time-steps it looks at to choose the current time step, the maximum change in the time step (ratios), the target number of jacobian iterations, the maximum number of iterations, and the error bound. I would like if someone could put me on the right path to analyzing the time-step independence, or at the very least figuring out what algorithm is used.

David Ketcheson
  • 16,522
  • 4
  • 54
  • 105
Godric Seer
  • 4,637
  • 2
  • 23
  • 38
  • Are you saying that you think it will be easier to reverse-engineer the time stepping algorithm than to just read the code? – David Ketcheson Sep 04 '12 at 01:41
  • the code is approximately 50k lines of fortran written over the past 20 years, so hunting down the details of the main loop is non-trivial. I believe that it is an implicit method, which is enough information for my purposes. I am more interested in what I need to change in separate runs to establish that my time-step is sufficiently small. – Godric Seer Sep 04 '12 at 02:26
  • I've tried to clarify what is being asked; please correct it if I've misinterpreted. Note that the solution cannot be literally "time-step independent", since the local errors will always depend on the time step. You can only hope that the errors are small enough for your purposes. – David Ketcheson Sep 04 '12 at 05:16

1 Answers1

5

The purpose of automatic error estimation and step size control is to free you from the problem of determining manually what a sufficiently small step size is. So your question is a bit like asking "somebody gave me this automatic transmission car; how can I tell what gear I'm in?" The point is that you shouldn't need to know. Of course, if the transmission is faulty, then you might need to take it apart and fix it, but that's a much bigger problem.

In your case, typically the right approach is to determine what kind of error is acceptable and impose that through the automatic step size control. It's imperfect because error control in this sense is usually only local error control, so you don't directly control the global error, which is what you probably care about.

One thing you could easily do if you're in doubt, is to run the simulation with a sequence of increasingly tight (i.e. small) error tolerances. Once the solution seems insensitive (in whatever your metric is) to decreasing the tolerance, you can stop.

Addendum: Regarding the issue of the maximum iteration limit being reached (which leads to a local error exceeding the specified tolerance), I suggest the following.

Apparently the code thinks that if it exceeds the maximum number of iterations, the right thing to do is to accept the step. I would say that is the wrong thing. A better approach is to reject the step and start that step again with a reduced step size. Of course, there is the danger of the step being reduced below the minimum step size. In that case, the right thing is to abort the simulation. But if you believe that a wrong solution is better than no solution, you could just accept the step and go on if both conditions are met: the minimum step size is reached and the maximum number of iterations is exceeded.

In a well-designed code, making these types of changes is trivial, but in an arbitrary code it may be arbitrarily difficult.

David Ketcheson
  • 16,522
  • 4
  • 54
  • 105
  • The last paragraph is the important one from a practical point of view. The asymptotic regime is reached once the plot of relevant solution characteristics against specified error tolerances looks regular enough. – Arnold Neumaier Sep 04 '12 at 07:17
  • That is good news, however I am still unsure how to handle those few occasions when it reaches the maximum iteration limit. I have seen it happily take the maximum time step, reaching $10^{-7}$ local error in half a dozen iterations, then the next time step take 100 iterations and only reach a local error of $10^2$. The time step after this it reduces the step size and reaches the error bound fine, but I still have that one large error propagating through the solution. – Godric Seer Sep 04 '12 at 11:13
  • I just noticed the addendum. Your description of what the code actually does agrees with what I see from logfiles. Hopefully the change you suggest is possible without too many headaches. – Godric Seer Sep 05 '12 at 16:17