jorgensd / jorgensd/adios4dolfinx

Float precision results in invalid time stamps

Open
#186 0 comments 0 reactions 0 assignees View on GitHub
Dominant language
Python
Stars
33
Forks
8
Avg merge
22d 12h
Merged PRs (30d)
1

Description

Consider the following example where we want to save some function after a time loop and then later we would like to read in the same function (this could also be in a different postprocessing script)
```python
from pathlib import Path
from mpi4py import MPI
import adios4dolfinx
import dolfinx

N = 1
comm = MPI.COMM_WORLD
domain = dolfinx.mesh.create_unit_square(comm, N, N)
V = dolfinx.fem.functionspace(domain, ("CG", 1))
u = dolfinx.fem.Function(V, name="u")
t = 0.0
dt = 0.01 # Time step size
T_final = 1.0 # Final time

checkpoint_fname = Path("checkpoint.bp")
adios4dolfinx.write_mesh(checkpoint_fname, domain)
while t < T_final:
t += dt
print("Save at time:", t)
print(abs(t - T_final))

adios4dolfinx.write_function(checkpoint_fname, u, time=t, name="u")

new_mesh = adios4dolfinx.read_mesh(filename=checkpoint_fname, comm=MPI.COMM_WORLD)
new_V = dolfinx.fem.functionspace(new_mesh, ("CG", 1))
new_u = dolfinx.fem.Function(new_V, name="u")
adios4dolfinx.read_function(checkpoint_fname, new_u, name="u", time=1.0)
print("u at time 1.0:", new_u.x.array)
```
The output of the program is as follows
```bash
$ python mwe.py
Save at time: 1.0000000000000007
6.661338147750939e-16
Traceback (most recent call last):
File "/home/shared/examples-fenicsx/mms/mwe.py", line 28, in
adios4dolfinx.read_function(checkpoint_fname, new_u, name="u", time=1.0)
File "/dolfinx-env/lib/python3.12/site-packages/adios4dolfinx/checkpointing.py", line 446, in read_function
input_array, starting_pos = read_array(
^^^^^^^^^^^
File "/dolfinx-env/lib/python3.12/site-packages/adios4dolfinx/adios2_helpers.py", line 308, in read_array
raise KeyError(
KeyError: 'No data associated with u_time=1.0 found in checkpoint.bp'
```
Now, one way to mitigate this is to use time_stamps stored in the file i.e
```python
import numpy as np

time_stamps = adios4dolfinx.read_timestamps(
checkpoint_fname, comm=comm, function_name="u"
)
time = 1.0
index = np.argmin(np.abs(time_stamps - time))
assert np.isclose(time_stamps[index], time)
adios4dolfinx.read_function(checkpoint_fname, new_u, name="u", time=time_stamps[index])
print("u at time 1.0:", new_u.x.array)
```
which works. However, I wonder whether we should do this within `adios4dolfinx` instead, since this would be a very common usecase which might trip off many users.

There are two other options to solve this which I think are also reasonable, so one of the following fixes should be implemented I think (ordered by preference)
1. Handle float precision within `adios4dolfinx` and only raise KeyError if a time stamp is not within a given tolerance of the available time stamps (this tolerance should also be possible to specify). In this case we need to update the API for `read_function` to also take an optional tolerance which could have a default value close to the float precision.
2. Print a more informative error message saying that the KeyError might be due to a float precision error and link to documentation of how to solve this (i.e the fix above). Then this needs to also be documented
3. Only allow integer time stamps to avoid floating number precision errors.

Contributor guide

Open the contributing guide

Research direction

Start with adios4dolfinx/checkpointing.py around read_function and adios4dolfinx/adios2_helpers.py around read_array, then reproduce the provided MWE. Compare the requested time with read_timestamps output and determine the intended tolerance or error behavior; done means the chosen approach is implemented, tested, and documented for floating-point timestamps.

Written by the indexing model from the issue text.

Assessment

Tech stack
python
Domain
data
Issue type
Bug
Difficulty
5/5
Estimated time
Over a week
Activity status
Stale
Clarity
Mostly clear
Newbie friendliness
35/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.