Expose tolerance for milestone checks

[ee-nn]

Description

Proposed fix for #933 (and my first attempt at a FESTIM codebase contribution!)

Summary

Adds a milestone_tolerance argument that allows the user to manually tune the relative tolerance (rtol) in the np.isclose check of whether the time is near a milestone. If none is specified, it passes the default of 1e-5.

Motivation and Context

For large times, the default relative tolerance can miss specified milestones. Exposing milestone_tolerance to the user provides a mitigation strategy.

Type of Change

• 🐛 Bug fix (non-breaking change which fixes an issue)
• ✨ New feature (non-breaking change which adds functionality)
• 💥 Breaking change (fix or feature that would cause existing functionality to not work as expected)
• 🔨 Code refactoring (no functional changes, no API changes)
• 📝 Documentation update
• ✅ Test update (adding missing tests or correcting existing tests)
• 🔧 Build/CI configuration change

Testing

• All existing tests pass locally (pytest)
• I have added new tests that prove my fix is effective or that my feature works

Code Quality Checklist

• My code follows the code style of this project (Ruff formatted: ruff format .)
• My code passes linting checks (ruff check .)
• I have performed a self-review of my own code
• I have commented my code, particularly in hard-to-understand areas

I kept comments to a bare minimum to keep with the rest of the code style, I can add more explanation if needed

Documentation

• I have updated the documentation accordingly (if applicable)
• I have added docstrings to new functions/classes following the project conventions

Screenshots/Examples

Running the "Simple Simulation" example (slightly modified, see below) with initial_value=1e-3 and milestone_tolerance=1e-2 reproduces the bug described in #933 ; the t=1 s milestone is missed:

Omitting milestone_tolerance causes fallback to the default 1e-5, successfully reproducing the graph on the simple simulation page.

import matplotlib.pyplot as plt
import numpy as np

import festim as F

model = F.HydrogenTransportProblem()
model.mesh = F.Mesh1D(vertices=np.linspace(0, 1, num=1001))

mat = F.Material(D_0=1, E_D=0.0)

volume_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, 1], material=mat)
boundary_left = F.SurfaceSubdomain1D(id=1, x=0)
boundary_right = F.SurfaceSubdomain1D(id=2, x=1)
model.subdomains = [volume_subdomain, boundary_left, boundary_right]

H = F.Species("H")
model.species = [H]

model.temperature = 300

model.boundary_conditions = [
    F.FixedConcentrationBC(subdomain=boundary_left, value=1, species=H),
    F.FixedConcentrationBC(subdomain=boundary_right, value=0, species=H),
]

model.settings = F.Settings(atol=1e-10, rtol=1e-10, final_time=2)

model.settings.stepsize = F.Stepsize(
    initial_value=1.0e-3,
    growth_factor=1.1,
    cutback_factor=0.9,
    target_nb_iterations=4,
    milestones=[0.05, 0.1, 0.2, 0.5, 1],
)


profile = F.Profile1DExport(
    field=H, subdomain=volume_subdomain, times=model.settings.stepsize.milestones
)

model.exports = [
    F.VTXSpeciesExport(field=H, filename="h_concentration.bp"),
    profile,
]

model.initialise()
model.run()
x = model.mesh.mesh.geometry.x[:, 0]
for time, data in zip(profile.t, profile.data):
    plt.plot(x, data, label=f"{time} s")

plt.xlabel("x (m)")
plt.ylabel("Mobile concentration (H/m3)")
plt.legend()
plt.show()

Additional Notes

[codecov]

Codecov Report

❌ Patch coverage is 88.88889% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 94.76%. Comparing base (ba4f051) to head (81c62ac).

Files with missing lines	Patch %	Lines
src/festim/stepsize.py	88.88%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1184      +/-   ##
==========================================
- Coverage   94.77%   94.76%   -0.02%     
==========================================
  Files          46       46              
  Lines        3428     3437       +9     
==========================================
+ Hits         3249     3257       +8     
- Misses        179      180       +1     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[RemDelaporteMathurin]

Thanks @ee-nn and congratulations on your first PR!

A few comments on my end.

• it wuold be nice to add a test in the test suite. You proved your fix is effective in the PR description but ideally we would catch this bug programmatically in the test suite. Let us know if you don't know how to write tests

[RemDelaporteMathurin]

Here the actual argument milestone_tolerance defaults to None. If this is the expected behaviour then we should change line 49. And perhaps removing the setter.

I think we should be more accurate in the description:

Suggested change

	milestone_tolerance (float, optional): relative tolerance for how closely
	a time must align with a milestone to be triggered. Defaults to 1e-5.
	milestone_tolerance (float, optional): relative tolerance passed to numpy.isclose (rtol) Defaults to 1e-5.

[RemDelaporteMathurin]

@ee-nn keeping the setter discussion in a thread for simplicity:

When you said "perhaps removing the setter" did you envision moving the error check up to the variable assignment?

I initially thought removing the setter/getter methods completely but I see you have a check for rtol < 0 that I think is good to keep.

[ee-nn]

Sorry, I added this comment in the above suggestion chain that I resolved, so re-adding it here...

Hi @RemDelaporteMathurin , great, thank you for the quick review! Setting the default directly in the init makes sense to me, I'll go ahead and apply your suggested changes.

I haven't written test suites so if you have any basic pointers I'd love to hear them.

[ee-nn]

When you said "perhaps removing the setter" did you envision moving the error check up to the variable assignment?

 def __init__(
    self,
    initial_value,
    growth_factor=None,
    cutback_factor=None,
    target_nb_iterations=None,
    max_stepsize=None,
    milestones=None,
    milestone_tolerance=1e-5,
 ) -> None:
    self.initial_value = initial_value
    self.growth_factor = growth_factor
    self.cutback_factor = cutback_factor
    self.target_nb_iterations = target_nb_iterations
    self.max_stepsize = max_stepsize
    self.milestones = milestones or []

    if milestone_tolerance <=0:  # <--- Put the check here 
        raise ValueError("milestone tolerance should be greater than zero")

    self.milestone_tolerance = milestone_tolerance

    if milestones and (growth_factor is None or cutback_factor is None):
        raise ValueError(
            "Milestones are only relevant if the stepsize is adaptive. "
            "Please provide growth and cutback factors."
        )

[RemDelaporteMathurin]

Sorry, I added this comment in the above suggestion chain that I resolved, so re-adding it here...

You can always "unresolve" conversations!

[RemDelaporteMathurin]

I haven't written test suites so if you have any basic pointers I'd love to hear them.

@ee-nn of course. I would suggest adding a test function in this file https://github.com/festim-dev/FESTIM/blob/main/test/test_stepsize.py

def test_milestone_tolerance_miss_tolerance():
    """
    This tests confirms that with a high rtol the milestone is missed
    See issue #number of issue
    """
   
   ....
   
   
   my_model.run()

   # check
   assert timestep_is_not_in_export

   ...
   # lower tolerance and rerun

   # check
   assert timestep_is_in_export