V1

docs/source/api_reference/datatypes.rst

@@ -7,4 +7,7 @@ lizzy.datatypes
 .. autoclass:: lizzy.datatypes.SimulationParameters
 
 
+.. autoclass:: lizzy.datatypes.Solution
+
+
 

examples/scripts/benchmark_solver_comparison.py

@@ -24,8 +24,7 @@ def test_solver_config(solver_type, use_masked, name, **kwargs):
         solution = model.solve()
         elapsed = time.time() - start
         print(f"  ✓ Success: {elapsed:.2f}s")
-        print(f"    - Fill time: {solution['time'][-1]:.5f}s")
-        print(f"    - Time steps: {solution['time_steps']}")
+        print(f"    - Fill time: {solution.time[-1]:.5f}s")
         return True, elapsed, solution
     except Exception as e:
         print(f"  ✗ Failed: {e}")

src/lizzy/_core/datatypes/__init__.py

@@ -1,3 +1,2 @@
 from .solution import Solution
-from .timestep import TimeStep
 from .simparams import SimulationParameters

src/lizzy/_core/datatypes/solution.py

@@ -1,9 +1,33 @@
 from dataclasses import dataclass
 import numpy as np
 
-@dataclass(slots=True)
+@dataclass(slots=True, frozen=True)
 class Solution:
-    time_steps : int
+    """A data class that stores the solution of a simulation.
+    It stores a number of time steps (the ones that were flagged for write-out), up to the instant of its creation.
+
+    Attributes
+    ----------
+
+    time_steps_in_solution : int
+        The number of time steps stored in the solution.
+    time_step_idx : ndarray of int, shape (time_steps_in_solution,)
+        The indices of the time steps stored in the solution. The last index corresponds to the time step number at which this solution was saved.
+    p : np.ndarray of float, shape (time_steps_in_solution, N_nodes)
+        The pressure values at each step.
+    v : np.ndarray of float, shape (time_steps_in_solution, N_elements, 3)
+        The velocity values at each step.
+    v_nodal : np.ndarray of float, shape (time_steps_in_solution, N_nodes, 3)
+        The nodal velocity values at each step.
+    time : np.ndarray of float, shape (time_steps_in_solution,)
+        The simulation time values at each step.
+    fill_factor : np.ndarray of float, shape (time_steps_in_solution, N_nodes)
+        The fill factor values at each step.
+    free_surface : np.ndarray of int, shape (time_steps_in_solution, N_nodes)
+        The free surface values at each step.
+    """
+    time_steps_in_solution : int
+    time_step_idx : np.ndarray
     p : np.ndarray
     v : np.ndarray
     v_nodal : np.ndarray

src/lizzy/_core/io/writer.py

@@ -11,6 +11,7 @@
 import numpy as np
 import meshio
 import textwrap
+from lizzy._core.datatypes import Solution
 
 
 class Writer:
@@ -32,7 +33,7 @@ def __init__(self, mesh):
         """
         self.mesh = mesh
 
-    def save_results(self, solution:dict, result_name:str, **kwargs):
+    def save_results(self, solution:Solution, result_name:str, **kwargs):
         """Save the results contained in the solution dictionary into an XDMF file.
 
         Parameters
@@ -53,28 +54,6 @@ def save_results(self, solution:dict, result_name:str, **kwargs):
         cells_list = []
         for i in range(len(cells)) :
             cells_list.append(cells[i])
-        # Iterate over time steps
-        for i in range(solution["time_steps"]):
-            # Create point data and cell data for the current time step
-            point_data = {
-                "FillFactor": solution["fill_factor"][i],  # Point data
-                "Pressure": solution["p"][i],  # Point data
-                "Time" : solution["time"],
-                "FreeSurface" : solution["free_surface"][i],
-                "Velocity" : solution["v_nodal"][i],
-            }
-            cell_data = {
-                "Velocity": [solution["v"][i]],  # Cell data for velocity
-            }
-            if _format == "vtk":# Create the meshio object with correct point and cell data
-                mesh_res = meshio.Mesh(
-                    points=points,
-                    cells=[("triangle", cells_list)],  # Triangle connectivity
-                    point_data=point_data,
-                    cell_data=cell_data,
-                )
-                # write the time step
-                mesh_res.write(destination_path / f"{result_name}_RES_{i}.vtk")
 
         if save_cv_mesh:
             mesh_cv = meshio.Mesh(
@@ -84,19 +63,19 @@ def save_results(self, solution:dict, result_name:str, **kwargs):
             mesh_cv.write(destination_path / f"{result_name}_CV.vtk")
 
         if _format == "xdmf":
-            filename = f"{result_name}_RES.xdmf"
+            filename = f"{result_name}.xdmf"
             with meshio.xdmf.TimeSeriesWriter(filename) as writer:
                 writer.write_points_cells(points, [("triangle", cells_list)])
-                for j in range(solution["time_steps"]):
-                    time = solution["time"][j]
-                    point_data = {  "Pressure" : np.array(solution["p"][j]),
-                                    "FillFactor" : np.array(solution["fill_factor"][j]),
-                                    "FreeSurface" : np.array(solution["free_surface"][j]),
-                                    "Velocity" : np.array(solution["v_nodal"][j])
+                for i in range(solution.time_steps_in_solution):
+                    time = solution.time[i]
+                    point_data = {  "Pressure" : solution.p[i],
+                                    "FillFactor" : solution.fill_factor[i],
+                                    "FreeSurface" : solution.free_surface[i],
+                                    "Velocity" : solution.v_nodal[i]
                                  }
-                    cell_data = { "Velocity" : np.array(solution["v"][j]) }
+                    cell_data = { "Velocity" : solution.v[i] }
                     writer.write_data(time, point_data=point_data, cell_data=cell_data)
             shutil.move(filename, destination_path / filename)
-            shutil.move(f"{result_name}_RES.h5", destination_path / f"{result_name}_RES.h5")
+            shutil.move(f"{result_name}.h5", destination_path / f"{result_name}.h5")
 
         print(f"Results saved in {destination_path}")

src/lizzy/_core/solver/solver.py

@@ -6,6 +6,8 @@
 
 from __future__ import annotations
 from typing import TYPE_CHECKING
+
+from lizzy._core.cvmesh import mesh
 if TYPE_CHECKING:
     from lizzy._core.sensors import SensorManager
     from lizzy._core.bcond import BCManager
@@ -33,7 +35,7 @@ def __init__(self, mesh, bc_manager, simulation_parameters, material_manager, se
         self.bc_manager : BCManager = bc_manager 
         self.simulation_parameters = simulation_parameters
         self.material_manager = material_manager
-        self.time_step_manager = TimeStepManager()
+        self.time_step_manager = TimeStepManager(mesh.nodes.N, mesh.triangles.N)
         self._sensor_manager = sensor_manager
         self.bcs = SolverBCs()
         self.vsolver = None
@@ -139,6 +141,23 @@ def update_n_empty_cvs(self):
         Must be called AFTER calling "update_empty_nodes_idx()"
         """
         self.n_empty_cvs = len(self.bcs.p0_idx)
+
+    def generate_initial_time_step(self):
+        time_0 = 0
+        dt_0 = 0
+        p_0 = np.zeros(self.mesh.nodes.N)
+        fill_factor_0 = np.zeros(self.mesh.nodes.N)
+        flow_front_0 = np.zeros(self.mesh.nodes.N)
+        for idx, val in zip(self.bcs.dirichlet_idx, self.bcs.dirichlet_vals):
+            p_0[idx] = val
+            fill_factor_0[idx] = 1
+            flow_front_0[idx] = 1
+        v_0 = np.zeros((self.mesh.triangles.N, 3))
+        v_nodal_0 = np.zeros((self.mesh.nodes.N, 3))
+        write_out_0 = True
+        initial_time_step = (time_0, dt_0, p_0, v_0, v_nodal_0, fill_factor_0, flow_front_0, write_out_0)
+        return initial_time_step
+
 
     def initialise_new_solution(self):
         """
@@ -160,10 +179,11 @@ def initialise_new_solution(self):
         active_cvs_ids, self.solver_vars["free_surface_array"] = self.fill_solver.find_free_surface_cvs(
             self.solver_vars["fill_factor_array"], self.cv_support_cvs_array)
         self.time_step_manager.reset()
-        self.time_step_manager.save_initial_timestep(self.mesh, self.bcs)
+        initial_time_step = self.generate_initial_time_step()
+        self.time_step_manager.save_timestep(*initial_time_step)
         self._sensor_manager.reset_sensors()
         # TODO: this first probe is temporary and should be cleaner
-        self._sensor_manager.probe_current_solution(self.time_step_manager.time_steps[0].P, self.time_step_manager.time_steps[0].V_nodal, self.time_step_manager.time_steps[0].fill_factor, 0.0)
+        self._sensor_manager.probe_current_solution(self.time_step_manager.p_buffer[0], self.time_step_manager.v_nodal_buffer[0], self.time_step_manager.fill_factor_buffer[0], 0.0)
 
     def handle_wo_criterion(self, dt):
         write_out = False

src/lizzy/_core/solver/timestep_manager.py

@@ -6,66 +6,89 @@
 
 
 import numpy as np
-from lizzy._core.datatypes import TimeStep, Solution
+from lizzy._core.datatypes import Solution
 
 
 class TimeStepManager:
-    def __init__(self):
-        self.time_steps = []
-        self.time_step_count = 0
+    def __init__(self, n_nodes, n_elements):
+        self.n_nodes : int = n_nodes
+        self.n_elements : int = n_elements
+        self.time_step_buffer_size : int = None
+        self.time_step_count : int = None
+        self.time_buffer : np.ndarray = None
+        self.dt_buffer : np.ndarray = None
+        self.p_buffer : np.ndarray = None
+        self.v_buffer : np.ndarray = None
+        self.v_nodal_buffer : np.ndarray = None
+        self.fill_factor_buffer : np.ndarray = None
+        self.flow_front_buffer : np.ndarray = None
+        self.write_out_buffer : np.ndarray = None
+        self.reset()
 
     def save_timestep(self, time, dt, P, v_array, v_nodal_array, fill_factor, flow_front, write_out):
-        if(v_array.shape[1]==3):
-            v_full = v_array
-            v_nodal_full = v_nodal_array
-        else:
-            v3_nul = np.zeros((np.size(v_array,0), 1))
-            v3_nodal_nul = np.zeros((np.size(v_nodal_array,0), 1))
-            v_full = np.hstack((v_array, v3_nul))
-            v_nodal_full = np.hstack((v_nodal_array, v3_nodal_nul))
-        timestep = TimeStep(self.time_step_count, time, dt, P, v_full, v_nodal_full, np.clip(fill_factor, 0, 1), flow_front, write_out)
-        self.time_steps.append(timestep)
+        if self.time_step_count >= self.time_step_buffer_size:
+            self.grow_buffers()
+        self.time_buffer[self.time_step_count] = time
+        self.dt_buffer[self.time_step_count] = dt
+        self.p_buffer[self.time_step_count, :] = P
+        self.v_buffer[self.time_step_count, :, :] = v_array
+        self.v_nodal_buffer[self.time_step_count, :, :] = v_nodal_array
+        self.fill_factor_buffer[self.time_step_count, :] = fill_factor
+        self.flow_front_buffer[self.time_step_count, :] = flow_front
+        self.write_out_buffer[self.time_step_count] = write_out
         self.time_step_count += 1
 
-    def get_write_out_steps(self):
-        return [step for step in self.time_steps if step.write_out == True]
+
+    def grow_buffers(self):
+            new_size = self.time_step_buffer_size * 2
+            self.time_buffer = np.resize(self.time_buffer, new_size)
+            self.dt_buffer = np.resize(self.dt_buffer, new_size)
+            self.p_buffer = np.resize(self.p_buffer, (new_size, self.p_buffer.shape[1]))
+            self.v_buffer = np.resize(self.v_buffer, (new_size, self.v_buffer.shape[1], self.v_buffer.shape[2]))
+            self.v_nodal_buffer = np.resize(self.v_nodal_buffer, (new_size, self.v_nodal_buffer.shape[1], self.v_nodal_buffer.shape[2]))
+            self.fill_factor_buffer = np.resize(self.fill_factor_buffer, (new_size, self.fill_factor_buffer.shape[1]))
+            self.flow_front_buffer = np.resize(self.flow_front_buffer, (new_size, self.flow_front_buffer.shape[1]))
+            self.write_out_buffer = np.resize(self.write_out_buffer, new_size)
+            self.time_step_buffer_size = new_size
+
 
-    def save_initial_timestep(self, mesh, bcs):
-        time_0 = 0
-        p_0 = [0] * mesh.nodes.N
-        fill_factor_0 = [0] * mesh.nodes.N
-        flow_front_0 = [0] * mesh.nodes.N
-        for i, val in enumerate(bcs.dirichlet_idx):
-            p_0[val] = bcs.dirichlet_vals[i]
-            fill_factor_0[val] = 1
-            flow_front_0[val] = 1
-        v_0 = np.zeros((mesh.triangles.N, 2))
-        v_nodal_0 = np.zeros((mesh.nodes.N, 2))
-        self.save_timestep(time_0, 0, p_0, v_0, v_nodal_0, fill_factor_0, flow_front_0, True)
+    def get_write_out_indices(self):
+        write_out_array = self.write_out_buffer[:self.time_step_count]
+        return np.nonzero(write_out_array)[0]
 
     def pack_solution(self):
         # flag the last time step as write-out regardless of its setting:
-        self.time_steps[-1].write_out = True
+        if self.time_step_count > 0:
+            self.write_out_buffer[self.time_step_count - 1] = True
         # populate solution with write-out time steps:
-        wo_time_steps = self.get_write_out_steps()
-        solution_obj = Solution(len(wo_time_steps),
-                                    np.array([step.P for step in wo_time_steps]),
-                                    np.array([step.V.tolist() for step in wo_time_steps]),
-                                    np.array([step.V_nodal for step in wo_time_steps]),
-                                    np.array([step.time for step in wo_time_steps]),
-                                    np.array([step.fill_factor for step in wo_time_steps]),
-                                    np.array([step.flow_front for step in wo_time_steps]),
-                                    )
-        solution = {"time_steps" : len(wo_time_steps),
-                    "p" : [step.P for step in wo_time_steps],
-                    "v" : [step.V.tolist() for step in wo_time_steps],
-                    "v_nodal" : [step.V_nodal for step in wo_time_steps],
-                    "time" : [step.time for step in wo_time_steps],
-                    "fill_factor" : [step.fill_factor for step in wo_time_steps],
-                    "free_surface" : [step.flow_front for step in wo_time_steps],
-                    }
-        return solution
+        wo_idx = self.get_write_out_indices()
+        solution_obj = Solution(len(wo_idx),
+                                wo_idx,
+                                self.p_buffer[wo_idx, :],
+                                self.v_buffer[wo_idx, :, :],
+                                self.v_nodal_buffer[wo_idx, :, :],
+                                self.time_buffer[wo_idx],
+                                self.fill_factor_buffer[wo_idx, :],
+                                self.flow_front_buffer[wo_idx, :],
+                                )
+        # solution = {"time_steps" : len(wo_idx),
+        #             "p" : [step.P for step in wo_time_steps],
+        #             "v" : [step.V.tolist() for step in wo_time_steps],
+        #             "v_nodal" : [step.V_nodal for step in wo_time_steps],
+        #             "time" : [step.time for step in wo_time_steps],
+        #             "fill_factor" : [step.fill_factor for step in wo_time_steps],
+        #             "free_surface" : [step.flow_front for step in wo_time_steps],
+        #             }
+        return solution_obj
 
     def reset(self):
-        self.time_steps = []
-        self.time_step_count = 0
+        self.time_step_buffer_size = 1000
+        self.time_step_count = 0
+        self.time_buffer = np.empty(self.time_step_buffer_size, dtype=float)
+        self.dt_buffer = np.empty(self.time_step_buffer_size, dtype=float)
+        self.p_buffer = np.empty((self.time_step_buffer_size, self.n_nodes), dtype=float)
+        self.v_buffer = np.empty((self.time_step_buffer_size, self.n_elements, 3), dtype=float)
+        self.v_nodal_buffer = np.empty((self.time_step_buffer_size, self.n_nodes, 3), dtype=float)
+        self.fill_factor_buffer = np.empty((self.time_step_buffer_size, self.n_nodes), dtype=float)
+        self.flow_front_buffer = np.empty((self.time_step_buffer_size, self.n_nodes), dtype=int)
+        self.write_out_buffer = np.empty(self.time_step_buffer_size, dtype=bool)

src/lizzy/datatypes/__init__.py

@@ -1 +1,7 @@
-from lizzy._core.datatypes import SimulationParameters
+#  Copyright 2025-2026 Simone Bancora, Paris Mulye
+#
+#  This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
+#  This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+#  You should have received a copy of the GNU General Public License along with this program. If not, see <https://www.gnu.org/licenses/>.
+
+from lizzy._core.datatypes import SimulationParameters, Solution