Skip to content

Use consistent time point tolerance - #133

Open
jpeloquin wants to merge 1 commit into
febiosoftware:developfrom
jpeloquin:develop
Open

Use consistent time point tolerance#133
jpeloquin wants to merge 1 commit into
febiosoftware:developfrom
jpeloquin:develop

Conversation

@jpeloquin

Copy link
Copy Markdown

Without this patch, FETimeStepController::CheckMustPoints uses a relative tolerance = 1e-12 for time, but FEAnalysis::Solve uses a relative tolerance = 1e-7. This can occasionally cause a must point to be omitted or an extra must point to be inserted.

Consider a file with with <plot_level>PLOT_MUST_POINTS</plot_level> and a must point at the begging & end of each analysis step.

FEAnalysis::Solve may decide an analysis step is done "early" (from the perspective of CheckMustPoints) while it is on a time point chosen by the auto-stepper. If this is the last analysis step, the last must point is consequently not written to the .xplt. A must point is only written by FEBioModel::WritePlot if pstep->m_timeController->m_nmust >= 0, and FETimeStepController::CheckMustPoints leaves m_nmust = -1 unless it needs to shift the time step chosen by FETimeStepController::AutoTimeStep earlier in time. If Solve ends an analysis step early, CheckMustPoints does not adjust the time step, because according to its 1e-12 relative tolerance the end of the analysis step has not yet been reached. In WritePlot m_nmust therefore is equal to -1 and the time step is not written.

This discrepancy can also cause an extra must point to be inserted. Consider the same file, but with multiple analysis steps. When FEAnalysis::Solve ends a time step "early" (according to the 1e12 relative tolerance used by CheckMustPoints), it calls fem.SetStartTime(fem.GetCurrentTime()), which sets m_ftime0 to the early termination time. When the next analysis step starts, FEAnalysis::Activate sets m_tstart = fem.GetStartTime(), and fem.GetStartTime returns m_ftime0, which is now early ( vs. CheckMustPoints) by a relative error up to 1e-7. FEAnalysis::Activate furthermore sets m_tend = m_tstart + Dt, so the analysis step's termination time is now early by up to 1e-7 relative error.

As the time steps approach the end of the analysis step, CheckMustPoints will adjust the time step to pass through the step's final must point. However, due to the error in m_tend, the time point it is attempting to pass through is now later than m_tend. This triggers the "make sure we are not exceeding the final time" condition in FETimeStepController::AutoTimeStep, because told + dtn > m_step->m_tend. AutoTimeStep does not adjust m_nmust, so it is still has the value ≠ -1 that was set by CheckMustPoints. Therefore, the time step is written to the .xplt by FEBioModel::WritePlot even if it does not equal the next must point time according to the 1e-12 tolerance in CheckMustPoints. If there is a subsequent analysis step, CheckMustPoints will once again attempt to pass through the same must point, this time successfully, and it will be written to the .xplt, resulting in one more time point than is specified by the file.

With this patch, FEAnalysis::Solve now uses a relative tolerance = 1e-12 for time, consistent with FETimeStepController::CheckMustPoints. This should almost entirely prevent both issues provided the user sets their must points exactly equal (equal 64-bit float representations) to the end times FEBio will calculate from dt * nsteps. Then each analysis step will end with a time step chosen by the CheckMustPoints, which attempts to choose a time step dt that will result in a final time equal to the analysis step's final must point. It is not a perfect fix because in IEEE-754 floating point subtraction and addition are rounded separately, and CheckMustPoints can only return dt = tmust - t, not the desired t. t + (tmust - t) == t_must is not always true.

A complete fix would require that FETimeStepController provide both dt and t_next (so they are never recalculated and roundoff error can be controlled in one place) or for FETimeStepController to keep track of drift and periodically adjust dt to compensate for it. Calling fem.SetStartTime(endtime) at the end of FEAnalysis::Solve was tempting but desynchronizing reported time from the simulated physical duration, however slightly, seems risky. Nevertheless, making the tolerances equal will greatly reduce the likelihood of must point errors.

Without this patch, `FETimeStepController::CheckMustPoints` uses a relative tolerance = 1e-12 for time, but `FEAnalysis::Solve` uses a relative tolerance = 1e-7.  This can occasionally cause a must point to be omitted or an extra must point to be inserted.

Consider a file with with `<plot_level>PLOT_MUST_POINTS</plot_level>` and a must point at the begging & end of each analysis step.

`FEAnalysis::Solve` may decide an analysis step is done "early" (from the perspective of CheckMustPoints) while it is on a time point chosen by the auto-stepper.  If this is the last analysis step, the last must point is consequently not written to the .xplt.  A must point is only written by `FEBioModel::WritePlot` if `pstep->m_timeController->m_nmust >= 0`, and `FETimeStepController::CheckMustPoints` leaves `m_nmust = -1` unless it needs to shift the time step chosen by `FETimeStepController::AutoTimeStep` earlier in time.  If `Solve` ends an analysis step early, `CheckMustPoints` does not adjust the time step, because according to its 1e-12 relative tolerance the end of the analysis step has not yet been reached. In `WritePlot` `m_nmust` therefore is equal to -1 and the time step is not written.

This discrepancy can also cause an _extra_ must point to be inserted.  Consider the same file, but with multiple analysis steps.  When `FEAnalysis::Solve` ends a time step "early" (according to the 1e12 relative tolerance used by `CheckMustPoints`), it calls `fem.SetStartTime(fem.GetCurrentTime())`, which sets `m_ftime0` to _the early termination time_.  When the next analysis step starts, `FEAnalysis::Activate` sets `m_tstart = fem.GetStartTime()`, and `fem.GetStartTime` returns `m_ftime0`, which is now early ( vs. `CheckMustPoints`) by a relative error up to 1e-7. `FEAnalysis::Activate` furthermore sets `m_tend = m_tstart + Dt`, so the analysis step's termination time is now early by up to 1e-7 relative error.

As the time steps approach the end of the analysis step, `CheckMustPoints` will adjust the time step to pass through the step's final must point.  However, due to the error in `m_tend`, the time point it is attempting to pass through is now _later than_ `m_tend`.  This triggers the "make sure we are not exceeding the final time" condition in `FETimeStepController::AutoTimeStep`, because `told + dtn > m_step->m_tend`.  `AutoTimeStep` does _not_ adjust `m_nmust`, so it is still has the value ≠ -1 that was set by `CheckMustPoints`.  Therefore, the time step is written to the .xplt by `FEBioModel::WritePlot` even if it _does not_ equal the next must point time according to the 1e-12 tolerance in `CheckMustPoints`.  If there is a subsequent analysis step, `CheckMustPoints` will once again attempt to pass through the same must point, this time successfully, and it will be written to the .xplt, resulting in one more time point than is specified by the file.

With this patch, `FEAnalysis::Solve` now uses a relative tolerance = 1e-12 for time, consistent with `FETimeStepController::CheckMustPoints`.  This should almost entirely prevent both issues provided the user sets their must points _exactly_ equal (equal 64-bit float representations) to the end times FEBio will calculate from dt * nsteps. Then each analysis step will end with a time step chosen by the `CheckMustPoints`, which attempts to choose a time step `dt` that will result in a final time equal to the analysis step's final must point.  It is not a perfect fix because in IEEE-754 floating point subtraction and addition are rounded separately, and `CheckMustPoints` can only return `dt = tmust - t`, not the desired `t`. `t + (tmust - t) == t_must` is not always true.

A complete fix would require that `FETimeStepController` provide both dt and t_next (so they are never recalculated and roundoff error can be controlled in one place) or for `FETimeStepController` to keep track of drift and periodically adjust dt to compensate for it.  Calling `fem.SetStartTime(endtime)` at the end of `FEAnalysis::Solve` was tempting but desynchronizing reported time from the simulated physical duration, however slightly, seems risky. Nevertheless, making the tolerances equal will greatly reduce the likelihood of must point errors.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant