Repository: SophT-Team/SophT

Issue:
h5 field dumping functionality for eulerian grid fields
For more info:
<img width="726" alt="Screen Shot 2022-11-04 at 12 11 41 AM" src="https://user-images.githubusercontent.com/33580851/199894267-dc64d41f-ca6d-479f-862b-ac2f2bbc7319.png">

@armantekinalp do you want to have a shot at this?

Discussion:

Patch:
diff --git a/examples/2d_examples/FlowPastRodCase/flow_past_rod.py b/examples/2d_examples/FlowPastRodCase/flow_past_rod.py
index c97b976..12b5085 100644
--- a/examples/2d_examples/FlowPastRodCase/flow_past_rod.py
+++ b/examples/2d_examples/FlowPastRodCase/flow_past_rod.py
@@ -131,20 +131,13 @@ def flow_past_rod_case(
     )
     # ==================FLOW-ROD COMMUNICATOR SETUP END======
     if save_flow_data:
-        # setup IO
-        # TODO internalise this in flow simulator as dump_fields
-        io_origin = np.array(
-            [
-                flow_sim.position_field[y_axis_idx].min(),
-                flow_sim.position_field[x_axis_idx].min(),
-            ]
-        )
-        io_dx = flow_sim.dx * np.ones(grid_dim)
-        io_grid_size = np.array(grid_size)
-        io = spu.IO(dim=grid_dim, real_dtype=real_t)
-        io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
-        io.add_as_eulerian_fields_for_io(
-            vorticity=flow_sim.vorticity_field, velocity=flow_sim.velocity_field
+        # setup flow IO
+        io = spu.EulerianFieldIO(
+            position_field=flow_sim.position_field,
+            eulerian_fields_dict={
+                "vorticity": flow_sim.vorticity_field,
+                "velocity": flow_sim.velocity_field,
+            },
         )
         # Initialize rod IO
         rod_io = spu.CosseratRodIO(
diff --git a/examples/2d_examples/VortexInducedFlagFlappingCase/turek_and_hron_benchmark.py b/examples/2d_examples/VortexInducedFlagFlappingCase/turek_and_hron_benchmark.py
index df0467c..aeb1e50 100644
--- a/examples/2d_examples/VortexInducedFlagFlappingCase/turek_and_hron_benchmark.py
+++ b/examples/2d_examples/VortexInducedFlagFlappingCase/turek_and_hron_benchmark.py
@@ -238,20 +238,13 @@ def flow_past_rod_case(
     # ==================FLOW-ROD COMMUNICATOR SETUP END======
     # =================FLOW H5 DUMPING====================
     if save_flow_data:
-        # setup IO
-        # TODO internalise this in flow simulator as dump_fields
-        io_origin = np.array(
-            [
-                flow_sim.position_field[y_axis_idx].min(),
-                flow_sim.position_field[x_axis_idx].min(),
-            ]
-        )
-        io_dx = flow_sim.dx * np.ones(grid_dim)
-        io_grid_size = np.array(grid_size)
-        io = spu.IO(dim=grid_dim, real_dtype=real_t)
-        io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
-        io.add_as_eulerian_fields_for_io(
-            vorticity=flow_sim.vorticity_field, velocity=flow_sim.velocity_field
+        # setup flow IO
+        io = spu.EulerianFieldIO(
+            position_field=flow_sim.position_field,
+            eulerian_fields_dict={
+                "vorticity": flow_sim.vorticity_field,
+                "velocity": flow_sim.velocity_field,
+            },
         )
         # Initialize rod IO
         rod_io = spu.CosseratRodIO(
diff --git a/examples/3d_examples/ElasticNetCase/immersed_elastic_net.py b/examples/3d_examples/ElasticNetCase/immersed_elastic_net.py
index e697ec4..56e1714 100644
--- a/examples/3d_examples/ElasticNetCase/immersed_elastic_net.py
+++ b/examples/3d_examples/ElasticNetCase/immersed_elastic_net.py
@@ -21,7 +21,6 @@ def immersed_elastic_net_case(
     grid_dim = 3
     real_t = spu.get_real_t(precision)
     x_axis_idx = spu.VectorField.x_axis_idx()
-    y_axis_idx = spu.VectorField.y_axis_idx()
     z_axis_idx = spu.VectorField.z_axis_idx()
     z_range, y_range, x_range = domain_range
     grid_size_y = round(y_range / x_range * grid_size_x)
@@ -72,21 +71,13 @@ def immersed_elastic_net_case(
         )
     # ==================FLOW-ROD COMMUNICATOR SETUP END======
     if save_flow_data:
-        # setup IO
-        # TODO internalise this in flow simulator as dump_fields
-        io_origin = np.array(
-            [
-                flow_sim.position_field[z_axis_idx].min(),
-                flow_sim.position_field[y_axis_idx].min(),
-                flow_sim.position_field[x_axis_idx].min(),
-            ]
-        )
-        io_dx = flow_sim.dx * np.ones(grid_dim)
-        io_grid_size = np.array(grid_size)
-        io = spu.IO(dim=grid_dim, real_dtype=real_t)
-        io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
-        io.add_as_eulerian_fields_for_io(
-            vorticity=flow_sim.vorticity_field, velocity=flow_sim.velocity_field
+        # setup flow IO
+        io = spu.EulerianFieldIO(
+            position_field=flow_sim.position_field,
+            eulerian_fields_dict={
+                "vorticity": flow_sim.vorticity_field,
+                "velocity": flow_sim.velocity_field,
+            },
         )
         # Initialize carpet IO
         elastic_net_io = spu.IO(dim=grid_dim, real_dtype=real_t)
diff --git a/examples/3d_examples/FlowPastRodCase/flow_past_rod_case.py b/examples/3d_examples/FlowPastRodCase/flow_past_rod_case.py
index 8335628..b3ffff8 100644
--- a/examples/3d_examples/FlowPastRodCase/flow_past_rod_case.py
+++ b/examples/3d_examples/FlowPastRodCase/flow_past_rod_case.py
@@ -28,7 +28,6 @@ def flow_past_rod_case(
     grid_size_z, grid_size_y, grid_size_x = grid_size
     real_t = spu.get_real_t(precision)
     x_axis_idx = spu.VectorField.x_axis_idx()
-    y_axis_idx = spu.VectorField.y_axis_idx()
     z_axis_idx = spu.VectorField.z_axis_idx()
     rho_f = 1.0
     U_free_stream = 1.0
@@ -139,21 +138,13 @@ def flow_past_rod_case(
     )
 
     if save_data:
-        # setup IO
-        # TODO internalise this in flow simulator as dump_fields
-        io_origin = np.array(
-            [
-                flow_sim.position_field[z_axis_idx].min(),
-                flow_sim.position_field[y_axis_idx].min(),
-                flow_sim.position_field[x_axis_idx].min(),
-            ]
-        )
-        io_dx = flow_sim.dx * np.ones(grid_dim)
-        io_grid_size = np.array(grid_size)
-        io = spu.IO(dim=grid_dim, real_dtype=real_t)
-        io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
-        io.add_as_eulerian_fields_for_io(
-            vorticity=flow_sim.vorticity_field, velocity=flow_sim.velocity_field
+        # setup flow IO
+        io = spu.EulerianFieldIO(
+            position_field=flow_sim.position_field,
+            eulerian_fields_dict={
+                "vorticity": flow_sim.vorticity_field,
+                "velocity": flow_sim.velocity_field,
+            },
         )
         # Initialize sphere IO
         rod_io = spu.CosseratRodIO(
diff --git a/examples/3d_examples/FlowPastSphereCase/flow_past_sphere_case.py b/examples/3d_examples/FlowPastSphereCase/flow_past_sphere_case.py
index 3ad57cf..540408b 100644
--- a/examples/3d_examples/FlowPastSphereCase/flow_past_sphere_case.py
+++ b/examples/3d_examples/FlowPastSphereCase/flow_past_sphere_case.py
@@ -22,7 +22,6 @@ def flow_past_sphere_case(
     grid_size_z, grid_size_y, grid_size_x = grid_size
     real_t = spu.get_real_t(precision)
     x_axis_idx: int = spu.VectorField.x_axis_idx()
-    y_axis_idx = spu.VectorField.y_axis_idx()
     z_axis_idx = spu.VectorField.z_axis_idx()
     x_range = 1.0
     far_field_velocity = 1.0
@@ -79,21 +78,13 @@ def flow_past_sphere_case(
     # ==================FLOW-BODY COMMUNICATOR SETUP END======
 
     if save_flow_data:
-        # setup IO
-        # TODO internalise this in flow simulator as dump_fields
-        io_origin = np.array(
-            [
-                flow_sim.position_field[z_axis_idx].min(),
-                flow_sim.position_field[y_axis_idx].min(),
-                flow_sim.position_field[x_axis_idx].min(),
-            ]
-        )
-        io_dx = flow_sim.dx * np.ones(grid_dim)
-        io_grid_size = np.array(grid_size)
-        io = spu.IO(dim=grid_dim, real_dtype=real_t)
-        io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
-        io.add_as_eulerian_fields_for_io(
-            vorticity=flow_sim.vorticity_field, velocity=flow_sim.velocity_field
+        # setup flow IO
+        io = spu.EulerianFieldIO(
+            position_field=flow_sim.position_field,
+            eulerian_fields_dict={
+                "vorticity": flow_sim.vorticity_field,
+                "velocity": flow_sim.velocity_field,
+            },
         )
         # Initialize sphere IO
         sphere_io = spu.IO(dim=grid_dim, real_dtype=real_t)
diff --git a/examples/3d_examples/FlowThroughPipeCase/flow_through_circular_pipe_case.py b/examples/3d_examples/FlowThroughPipeCase/flow_through_circular_pipe_case.py
index 343fdde..430a4bf 100644
--- a/examples/3d_examples/FlowThroughPipeCase/flow_through_circular_pipe_case.py
+++ b/examples/3d_examples/FlowThroughPipeCase/flow_through_circular_pipe_case.py
@@ -26,7 +26,6 @@ def flow_through_circular_pipe_case(
     real_t = spu.get_real_t(precision)
     x_axis_idx = spu.VectorField.x_axis_idx()
     y_axis_idx = spu.VectorField.y_axis_idx()
-    z_axis_idx = spu.VectorField.z_axis_idx()
     x_range = 1.0
     nu = 1e-2
     flow_sim = sps.UnboundedNavierStokesFlowSimulator3D(
@@ -72,21 +71,13 @@ def flow_through_circular_pipe_case(
     # ==================FLOW-BODY COMMUNICATOR SETUP END======
 
     if save_data:
-        # setup IO
-        # TODO internalise this in flow simulator as dump_fields
-        io_origin = np.array(
-            [
-                flow_sim.position_field[z_axis_idx].min(),
-                flow_sim.position_field[y_axis_idx].min(),
-                flow_sim.position_field[x_axis_idx].min(),
-            ]
-        )
-        io_dx = flow_sim.dx * np.ones(grid_dim)
-        io_grid_size = np.array(grid_size)
-        io = spu.IO(dim=grid_dim, real_dtype=real_t)
-        io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
-        io.add_as_eulerian_fields_for_io(
-            vorticity=flow_sim.vorticity_field, velocity=flow_sim.velocity_field
+        # setup flow IO
+        io = spu.EulerianFieldIO(
+            position_field=flow_sim.position_field,
+            eulerian_fields_dict={
+                "vorticity": flow_sim.vorticity_field,
+                "velocity": flow_sim.velocity_field,
+            },
         )
 
     t_end = 1.0
diff --git a/examples/3d_examples/HillSphericalVortexCase/hill_sphere_vortex_case.py b/examples/3d_examples/HillSphericalVortexCase/hill_sphere_vortex_case.py
index 812940e..cb8c287 100644
--- a/examples/3d_examples/HillSphericalVortexCase/hill_sphere_vortex_case.py
+++ b/examples/3d_examples/HillSphericalVortexCase/hill_sphere_vortex_case.py
@@ -15,7 +15,6 @@ def hill_sphere_vortex_case(
     the velocity recovery and vortex stretching steps, for the solver,
     by comparing against analytical expressions.
     """
-    grid_dim = 3
     grid_size_z, grid_size_y, grid_size_x = grid_size
     real_t = spu.get_real_t(precision)
     x_axis_idx = spu.VectorField.x_axis_idx()
@@ -116,23 +115,14 @@ def hill_sphere_vortex_case(
     fig2.savefig("midplane_radial_velocity.png")
 
     if save_data:
-        # setup IO
-        # TODO internalise this in flow simulator as dump_fields
-        io_origin = np.array(
-            [
-                flow_sim.position_field[z_axis_idx].min(),
-                flow_sim.position_field[y_axis_idx].min(),
-                flow_sim.position_field[x_axis_idx].min(),
-            ]
-        )
-        io_dx = flow_sim.dx * np.ones(grid_dim)
-        io_grid_size = np.array(grid_size)
-        io = spu.IO(dim=grid_dim, real_dtype=real_t)
-        io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
-        io.add_as_eulerian_fields_for_io(
-            vorticity=flow_sim.vorticity_field,
-            stream_func=flow_sim.stream_func_field,
-            velocity=flow_sim.velocity_field,
+        # setup flow IO
+        io = spu.EulerianFieldIO(
+            position_field=flow_sim.position_field,
+            eulerian_fields_dict={
+                "vorticity": flow_sim.vorticity_field,
+                "stream_func": flow_sim.stream_func_field,
+                "velocity": flow_sim.velocity_field,
+            },
         )
         io.save(h5_file_name="sopht.h5")
 
diff --git a/examples/3d_examples/MagneticCiliaCarpetCase/immersed_magnetic_cilia_carpet.py b/examples/3d_examples/MagneticCiliaCarpetCase/immersed_magnetic_cilia_carpet.py
index 7ff549a..ecac8b4 100644
--- a/examples/3d_examples/MagneticCiliaCarpetCase/immersed_magnetic_cilia_carpet.py
+++ b/examples/3d_examples/MagneticCiliaCarpetCase/immersed_magnetic_cilia_carpet.py
@@ -19,7 +19,6 @@ def immersed_magnetic_cilia_carpet_case(
     grid_dim = 3
     real_t = spu.get_real_t(precision)
     x_axis_idx = spu.VectorField.x_axis_idx()
-    y_axis_idx = spu.VectorField.y_axis_idx()
     z_axis_idx = spu.VectorField.z_axis_idx()
     z_range, y_range, x_range = domain_range
     grid_size_y = round(y_range / x_range * grid_size_x)
@@ -65,21 +64,13 @@ def immersed_magnetic_cilia_carpet_case(
         )
     # ==================FLOW-ROD COMMUNICATOR SETUP END======
     if save_data:
-        # setup IO
-        # TODO internalise this in flow simulator as dump_fields
-        io_origin = np.array(
-            [
-                flow_sim.position_field[z_axis_idx].min(),
-                flow_sim.position_field[y_axis_idx].min(),
-                flow_sim.position_field[x_axis_idx].min(),
-            ]
-        )
-        io_dx = flow_sim.dx * np.ones(grid_dim)
-        io_grid_size = np.array(grid_size)
-        io = spu.IO(dim=grid_dim, real_dtype=real_t)
-        io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
-        io.add_as_eulerian_fields_for_io(
-            vorticity=flow_sim.vorticity_field, velocity=flow_sim.velocity_field
+        # setup flow IO
+        io = spu.EulerianFieldIO(
+            position_field=flow_sim.position_field,
+            eulerian_fields_dict={
+                "vorticity": flow_sim.vorticity_field,
+                "velocity": flow_sim.velocity_field,
+            },
         )
         # Initialize carpet IO
         carpet_io = spu.IO(dim=grid_dim, real_dtype=real_t)
diff --git a/examples/3d_examples/OctopusArmCase3D/run_tapered_arm_and_sphere_with_flow.py b/examples/3d_examples/OctopusArmCase3D/run_tapered_arm_and_sphere_with_flow.py
index d44ee98..64307c4 100644
--- a/examples/3d_examples/OctopusArmCase3D/run_tapered_arm_and_sphere_with_flow.py
+++ b/examples/3d_examples/OctopusArmCase3D/run_tapered_arm_and_sphere_with_flow.py
@@ -183,22 +183,13 @@ def tapered_arm_and_cylinder_flow_coupling(
     # ==================FLOW-ROD COMMUNICATOR SETUP END======
 
     if save_data:
-        # setup IO
-        # TODO internalise this in flow simulator as dump_fields
-        io_origin = np.array(
-            [
-                flow_sim.position_field[z_axis_idx].min(),
-                flow_sim.position_field[y_axis_idx].min(),
-                flow_sim.position_field[x_axis_idx].min(),
-            ]
-        )
-        io_dx = flow_sim.dx * np.ones(grid_dim)
-        io_grid_size = np.array(grid_size)
-        # Initialize flow eulerian grid IO
-        io = spu.IO(dim=grid_dim, real_dtype=real_t)
-        io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
-        io.add_as_eulerian_fields_for_io(
-            vorticity=flow_sim.vorticity_field, velocity=flow_sim.velocity_field
+        # setup flow IO
+        io = spu.EulerianFieldIO(
+            position_field=flow_sim.position_field,
+            eulerian_fields_dict={
+                "vorticity": flow_sim.vorticity_field,
+                "velocity": flow_sim.velocity_field,
+            },
         )
         # Initialize rod io
         rod_io = spu.CosseratRodIO(
diff --git a/examples/3d_examples/PointSourceAdvectAndDiffuseCase/point_source_advection_diffusion.py b/examples/3d_examples/PointSourceAdvectAndDiffuseCase/point_source_advection_diffusion.py
index f97b37b..ee400ef 100644
--- a/examples/3d_examples/PointSourceAdvectAndDiffuseCase/point_source_advection_diffusion.py
+++ b/examples/3d_examples/PointSourceAdvectAndDiffuseCase/point_source_advection_diffusion.py
@@ -15,7 +15,6 @@ def point_source_advection_diffusion_case(
     constant velocity in 3D, while it diffuses in time.
     """
     grid_dim = 3
-    grid_size_z, grid_size_y, grid_size_x = grid_size
     real_t = spu.get_real_t(precision)
     x_axis_idx = spu.VectorField.x_axis_idx()
     y_axis_idx = spu.VectorField.y_axis_idx()
@@ -62,20 +61,13 @@ def point_source_advection_diffusion_case(
     flow_sim.velocity_field[...] = velocity_free_stream
 
     if save_data:
-        # setup IO
-        # TODO internalise this in flow simulator as dump_fields
-        io_origin = np.array(
-            [
-                flow_sim.position_field[z_axis_idx].min(),
-                flow_sim.position_field[y_axis_idx].min(),
-                flow_sim.position_field[x_axis_idx].min(),
-            ]
+        # setup flow IO
+        io = spu.EulerianFieldIO(
+            position_field=flow_sim.position_field,
+            eulerian_fields_dict={
+                "vorticity": vorticity_field,
+            },
         )
-        io_dx = flow_sim.dx * np.ones(grid_dim)
-        io_grid_size = np.array([grid_size_z, grid_size_y, grid_size_x])
-        io = spu.IO(dim=grid_dim, real_dtype=real_t)
-        io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
-        io.add_as_eulerian_fields_for_io(vorticity=vorticity_field)
 
     foto_timer = 0.0
     foto_timer_limit = (t_end - t_start) / 20
diff --git a/sopht/utils/__init__.py b/sopht/utils/__init__.py
index 8c506ef..8340ac4 100644
--- a/sopht/utils/__init__.py
+++ b/sopht/utils/__init__.py
@@ -2,7 +2,7 @@ from .lab_cmap import get_lab_cmap
 from .plot_field import create_figure_and_axes, save_and_clear_fig
 from .post_process import make_video_from_image_series, make_dir_and_transfer_h5_data
 from .field import VectorField
-from .io import IO, CosseratRodIO
+from .io import IO, CosseratRodIO, EulerianFieldIO
 from .rod_viz import plot_video_of_rod_surface
 from .precision import get_real_t, get_test_tol
 from .pyst_kernel_config import get_pyst_dtype, get_pyst_kernel_config
diff --git a/sopht/utils/io.py b/sopht/utils/io.py
index 6de0adb..d757e3a 100644
--- a/sopht/utils/io.py
+++ b/sopht/utils/io.py
@@ -2,6 +2,7 @@
 import h5py
 import numpy as np
 from elastica.rod.cosserat_rod import CosseratRod
+import sopht.utils as spu
 
 
 class IO:
@@ -301,7 +302,7 @@ class IO:
         if self.lagrangian_fields:
             self.generate_xdmf_lagrangian(h5_file_name=h5_file_name, time=time)
 
-    def load(self, h5_file_name: str) -> None:  # noqa: C901
+    def load(self, h5_file_name: str) -> float:  # noqa: C901
         """Load fields from hdf5 file.
 
         Field arrays need to be allocated and added to `eulerian_fields` and/or
@@ -647,3 +648,52 @@ class CosseratRodIO(IO):
             self.cosserat_rod.position_collection[: self.dim, 1:]
             + self.cosserat_rod.position_collection[: self.dim, :-1]
         )
+
+
+class EulerianFieldIO(IO):
+    """
+    Derived IO class for saving Eulerian grid fields.
+    """
+
+    def __init__(
+        self, position_field: np.ndarray, eulerian_fields_dict: dict[str, np.ndarray]
+    ) -> None:
+        """Class initializer
+
+        :param position_field: Array with position coordinates, of shape (grid_dim, grid_size)
+        :param eulerian_fields_dict: Dictionary with keys as names of Eulerian fields
+        to be saved, with corresponding values representing the field array names
+
+        """
+        grid_dim = position_field.shape[0]
+        grid_size = position_field.shape[1:]
+        real_dtype = position_field.dtype
+        super().__init__(dim=grid_dim, real_dtype=real_dtype)
+
+        match grid_dim:
+            case 2:
+                io_origin = np.array(
+                    [
+                        position_field[spu.VectorField.y_axis_idx()].min(),
+                        position_field[spu.VectorField.x_axis_idx()].min(),
+                    ]
+                )
+            case 3:
+                io_origin = np.array(
+                    [
+                        position_field[spu.VectorField.z_axis_idx()].min(),
+                        position_field[spu.VectorField.y_axis_idx()].min(),
+                        position_field[spu.VectorField.x_axis_idx()].min(),
+                    ]
+                )
+            case _:
+                raise ValueError("Position field of invalid shape.")
+        # get grid spacing from X coodinate field
+        x_grid_flattened_field = position_field[spu.VectorField.x_axis_idx()].reshape(
+            -1,
+        )
+        dx = x_grid_flattened_field[1] - x_grid_flattened_field[0]
+        io_dx = dx * np.ones(grid_dim)
+        io_grid_size = np.array(grid_size)
+        self.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
+        self.add_as_eulerian_fields_for_io(**eulerian_fields_dict)


Test patch:
diff --git a/tests/test_utils/test_io.py b/tests/test_utils/test_io.py
index 3c1e752..20fcae2 100644
--- a/tests/test_utils/test_io.py
+++ b/tests/test_utils/test_io.py
@@ -441,3 +441,72 @@ def test_cosserat_rod_io(dim, precision):
     )
     np.testing.assert_allclose(time, time_loaded, atol=testing_atol)
     os.system("rm -f *h5 *xmf")
+
+
+@pytest.mark.parametrize("precision", ["single", "double"])
+@pytest.mark.parametrize("grid_dim", [2, 3])
+@pytest.mark.parametrize("grid_size_x", [8, 16])
+def test_eulerian_field_io(grid_dim, precision, grid_size_x):
+    real_t = spu.get_real_t(precision)
+    grid_size = (grid_size_x,) * grid_dim
+    scalar_field = np.random.rand(*grid_size).astype(real_t)
+    vector_field = np.random.rand(grid_dim, *grid_size).astype(real_t)
+
+    x_range = 2.0
+    dx = real_t(x_range / grid_size_x)
+    eul_grid_shift = dx / 2.0
+    x = np.linspace(eul_grid_shift, x_range - eul_grid_shift, grid_size_x).astype(
+        real_t
+    )
+    position_field = np.flipud(np.array(np.meshgrid(*((x,) * grid_dim), indexing="ij")))
+    h5_file_name = "eulerian_field.h5"
+    time = 2.0
+    eulerian_field_dict = {"scalar_field": scalar_field, "vector_field": vector_field}
+    test_io = spu.EulerianFieldIO(
+        position_field=position_field, eulerian_fields_dict=eulerian_field_dict
+    )
+    test_io.save(h5_file_name=h5_file_name, time=time)
+    del test_io
+
+    # Load saved HDF5 file for checking
+    scalar_field_loaded = np.zeros_like(scalar_field)
+    vector_field_loaded = np.zeros_like(vector_field)
+    io = spu.IO(dim=grid_dim, real_dtype=real_t)
+    match grid_dim:
+        case 2:
+            io_origin = np.array(
+                [
+                    position_field[spu.VectorField.y_axis_idx()].min(),
+                    position_field[spu.VectorField.x_axis_idx()].min(),
+                ]
+            )
+        case 3:
+            io_origin = np.array(
+                [
+                    position_field[spu.VectorField.z_axis_idx()].min(),
+                    position_field[spu.VectorField.y_axis_idx()].min(),
+                    position_field[spu.VectorField.x_axis_idx()].min(),
+                ]
+            )
+        case _:
+            raise ValueError("Position field of invalid shape.")
+    io_dx = dx * np.ones(grid_dim)
+    io_grid_size = np.array(grid_size)
+    io.define_eulerian_grid(origin=io_origin, dx=io_dx, grid_size=io_grid_size)
+    io.add_as_eulerian_fields_for_io(scalar_field=scalar_field_loaded)
+    io.add_as_eulerian_fields_for_io(vector_field=vector_field_loaded)
+    time_loaded = io.load(h5_file_name=h5_file_name)
+
+    # Check values
+    np.testing.assert_allclose(
+        scalar_field,
+        scalar_field_loaded,
+        atol=spu.get_test_tol(precision),
+    )
+    np.testing.assert_allclose(
+        vector_field,
+        vector_field_loaded,
+        atol=spu.get_test_tol(precision),
+    )
+    np.testing.assert_allclose(time, time_loaded, atol=spu.get_test_tol(precision))
+    os.system("rm -f *h5 *xmf")