Use consistent time point tolerance - #133
Open
jpeloquin wants to merge 1 commit into
Open
Conversation
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.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Without this patch,
FETimeStepController::CheckMustPointsuses a relative tolerance = 1e-12 for time, butFEAnalysis::Solveuses 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::Solvemay 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 byFEBioModel::WritePlotifpstep->m_timeController->m_nmust >= 0, andFETimeStepController::CheckMustPointsleavesm_nmust = -1unless it needs to shift the time step chosen byFETimeStepController::AutoTimeStepearlier in time. IfSolveends an analysis step early,CheckMustPointsdoes not adjust the time step, because according to its 1e-12 relative tolerance the end of the analysis step has not yet been reached. InWritePlotm_nmusttherefore 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::Solveends a time step "early" (according to the 1e12 relative tolerance used byCheckMustPoints), it callsfem.SetStartTime(fem.GetCurrentTime()), which setsm_ftime0to the early termination time. When the next analysis step starts,FEAnalysis::Activatesetsm_tstart = fem.GetStartTime(), andfem.GetStartTimereturnsm_ftime0, which is now early ( vs.CheckMustPoints) by a relative error up to 1e-7.FEAnalysis::Activatefurthermore setsm_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,
CheckMustPointswill adjust the time step to pass through the step's final must point. However, due to the error inm_tend, the time point it is attempting to pass through is now later thanm_tend. This triggers the "make sure we are not exceeding the final time" condition inFETimeStepController::AutoTimeStep, becausetold + dtn > m_step->m_tend.AutoTimeStepdoes not adjustm_nmust, so it is still has the value ≠ -1 that was set byCheckMustPoints. Therefore, the time step is written to the .xplt byFEBioModel::WritePloteven if it does not equal the next must point time according to the 1e-12 tolerance inCheckMustPoints. If there is a subsequent analysis step,CheckMustPointswill 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::Solvenow uses a relative tolerance = 1e-12 for time, consistent withFETimeStepController::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 theCheckMustPoints, which attempts to choose a time stepdtthat 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, andCheckMustPointscan only returndt = tmust - t, not the desiredt.t + (tmust - t) == t_mustis not always true.A complete fix would require that
FETimeStepControllerprovide both dt and t_next (so they are never recalculated and roundoff error can be controlled in one place) or forFETimeStepControllerto keep track of drift and periodically adjust dt to compensate for it. Callingfem.SetStartTime(endtime)at the end ofFEAnalysis::Solvewas 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.