From b44a17948569f0b28da6967098f40127e0b39e0c Mon Sep 17 00:00:00 2001 From: John Peloquin Date: Tue, 1 Sep 2026 17:40:01 -0400 Subject: [PATCH] Use consistent time point tolerance MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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_MUST_POINTS` 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. --- FECore/FEAnalysis.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/FECore/FEAnalysis.cpp b/FECore/FEAnalysis.cpp index 45a6e721e..ed35930ba 100644 --- a/FECore/FEAnalysis.cpp +++ b/FECore/FEAnalysis.cpp @@ -402,7 +402,11 @@ bool FEAnalysis::Solve() double starttime = fem.GetStartTime(); // double endtime = fem.m_ftime0 + m_ntime*m_dt0; double endtime = m_tend; - const double eps = endtime*1e-7; + // Use the same time tolerance as FETimeStepController::CheckMustPoints. Otherwise, + // we may terminate the analysis step before reaching the last must point, or get an + // extra must point at the transition between two analysis steps due to accumulated + // early-termination error in m_timeController->m_tstart and m_timeController->m_tend. + const double eps = endtime*1e-12; // if we restarted we need to update the timestep // before continuing