Dedalus setting to handle machine precision when dumping files

[tlunet]

Hi dear Dedalus dev team,

there is currently a minor issue when writing files for analysis, for instance using the Spherical Shallow Water test problem given in examples. Because of some rounding error when adding dt to solver.sim_time, the fields are not written in file every 1 hour, as planned, but some time after a given number of hour + dt.

On my computer, I got the following simulation time for the first written files :

[  0.        ,   1.16666667,   2.        ,   3.16666667, 4.16666667,   5.        ,   6.        ,   7.        , ...]

This is due to the particular value for the time-step and the computation of sim_div in evaluate_scheduled. For instance, here we consider dt=5/30 which corresponds to 10 minutes. The first field should be written after 6 times step, be since dt+dt+dt+dt+dt+dt=0.9999999999999999, then the integer division in evaluate_scheduled does not call the handler, even if 0.9999 is way closer to 1 than 1.166666.

This is actually problematic for my team, since we are doing some time-convergence analysis with different time-steps, and even if we choose a time-step value that is a divider of 1 hour, machine precision produces some offsets in the compared data. Since this could be problematic for other users, here is a proposal for an easy fix :

• new dedalus parameter : FILEHANDLER_SIMTIME_EPSILON (with default value = 0), that is added to sim_time when computing sim_div in evaluate_scheduled
• if needed, one can simply set a different value for this epsilon in a personal dedalus.cfg file. For the SWE test case, it appears that FILEHANDLER_SIMTIME_EPSILON=1e-10 was adapted.

So here is a first idea, don't know if that is of interest for you (eventually improved ...).

Cheers !

[kburns]

Hi Thibaut, yes this is a problem we've encountered before, and it's just an unfortunate reality of floating point. I think adding a tolerance just sort of kicks the can down the road -- eventually the deviation will build up past any epsilon.

But there's a simple solution -- you can set the outputs to occur on a cadence set by the iteration count, rather than the simulation time, and since that's an integer everything is fine. When you create the file handler, you can just pass iter=6 for instance, instead of sim_dt. If you want the same outputs as you change the timestep (and they evenly divide), you still have to be a bit careful, but can do iter=int(np.round(1/dt)), for instance, to output at integer times.

[tlunet]

Thanks for the tip. The iter argument seems indeed a better option for the convergence analysis than the sim_dt for the evaluator. So should I interpret from your answer that you are not interested in introducing epsilon when comparing simulation time (at least, to handle situations when the floating point error is not catastrophically high ...) ?

[kburns]

I think it would be good to have this functionality in some way, but we're generally trying to keep the config options as minimal as possible in favor of other interfaces. This approach does seem pretty robust but there are some edge cases I can think of, including:

• the floating point error could eventually accumulate past epsilon
• if the minimum timestep in a simulation with e.g. adaptive timestepping drops below epsilon, you might end up taking an output too early, i.e. not on the timestep that's actually closest to the cadence.

But these both seem relatively minor. Maybe a good interface would be adding this via sim_dt_rtol and sim_dt_atol, for relative and/or absolute tolerances, as keyword arguments in the FileHandler? Any other thoughts @jsoishi or @lecoanet?

[kburns]

Actually why not just always use a relative tolerance of 0.5, or set that as the default? The outputs would then always be within half a timestep (from below) of the target time, and it seems like this would only be a suboptimal choice when the timestep is decreasing exactly 1 iteration before output, which seems fairly innocuous . Anyone see a problem with that approach? I feel like we talked about this years ago...

[bpbrown]

That seems like a good approach to me. Clean and simple.

On Fri, Nov 18, 2022 at 2:21 PM Keaton J. Burns @.***> wrote:

Actually why not just always use a relative tolerance of 0.5, or set that as the default? The outputs would then always be within half a timestep (from below) of the target time, and it seems like this would only be a suboptimal cohice when the timestep is decreasing exactly 1 iteration before output, which seems fairly innocuous . Anyone see a problem with that approach? I feel like we talked about this years ago...

— Reply to this email directly, view it on GitHub https://github.com/DedalusProject/dedalus/pull/232#issuecomment-1320044883, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABN7BYMINZLDBTYXSRD2SVDWI6FBZANCNFSM6AAAAAASDOY2WU . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[tlunet]

Actually why not just always use a relative tolerance of 0.5, or set that as the default? The outputs would then always be within half a timestep (from below) of the target time, and it seems like this would only be a suboptimal choice when the timestep is decreasing exactly 1 iteration before output, which seems fairly innocuous . Anyone see a problem with that approach? I feel like we talked about this years ago...

Do you mean, computing sim_div like this ?

sim_div  = (sim_time + 0.5*sim_time) // handler.sim_dt

[kburns]

Thinking about this a bit more, as the IVP progresses we do output at the beginning of each step, and so we always know the current time as well as the next iteration time, since we have the current timestep on hand. This means we could setup logic so that we do the output if either the current or next iteration meet the new divisor criteria, AND the current timestep is closer to the output time than the next one will be. I think this requires no free parameters and will always output at the nearest timestep to the output cadence.

[kburns]

We've implemented some logic similar to this in commit 0598602123aab89592df3c204da24dc3192e63db which will be merged with #238. If you're doing constant timestepping, it should always give the right output times, and it seems to do ok in other cases too, including when the timestep is longer than the output cadence. We've tested with your specific cadence and timestep (thanks for including those), and now it outputs on the near-integer steps as expected. I think this is a good way to go since it doesn't require any extra parameters. How does that sound?

[tlunet]

Sounds great, thanks ! Guess this PR can be closed now :smiley:

[kburns]

Thanks!