Create DynamicsData from Solution (#61)

* Create DynamicsData from Solution

* Update test_solution.py

tdgl/device/device.py

@@ -71,7 +71,7 @@ def __init__(
         film: Polygon,
         holes: Union[List[Polygon], None] = None,
         terminals: Union[List[Polygon], None] = None,
-        probe_points: Sequence[float] = None,
+        probe_points: Optional[Sequence[Tuple[float, float]]] = None,
         length_units: str = "um",
     ):
         self.name = name

tdgl/solution/data.py

@@ -442,6 +442,57 @@ def to_hdf5(self, h5group: h5py.Group) -> None:
         if self.screening_iterations is not None:
             h5group["screening_iterations"] = self.screening_iterations
 
+    @staticmethod
+    def from_solution(
+        solution_path: str,
+        probe_points: Optional[Sequence[Tuple[float, float]]] = None,
+    ) -> "DynamicsData":
+        """Load :class:`DynamicsData` from the saved time steps of a :class:`tdgl.Solution`.
+
+        Args:
+            solution_path: Path to the :class:`tdgl.Solution`
+            probe_points: The probe coordinates for which to extract dynamics.
+                If ``None``, defaults to ``solution.device.probe_points``.
+
+        Returns:
+            A new :class:`DynamicsData` instance
+        """
+        from .solution import Solution
+
+        solution = Solution.from_hdf5(solution_path)
+        device = solution.device
+        mesh = device.mesh
+        if probe_points is None:
+            probe_points = device.probe_points
+        if probe_points is None:
+            raise ValueError("No probe points were provided.")
+        probe_points = np.asarray(probe_points).squeeze()
+        if probe_points.ndim != 2 or probe_points.shape[1] != 2:
+            raise ValueError(
+                f"Probe points must have shape (n, 2), got {probe_points.shape}."
+            )
+        if not device.contains_points(probe_points).all():
+            raise ValueError("All probe points must lie within the film.")
+
+        xi = device.coherence_length.magnitude
+        probe_point_indices = [mesh.closest_site(xy) for xy in probe_points / xi]
+        step_min, step_max = solution.data_range
+
+        num_probes = len(probe_points)
+        num_steps = step_max - step_min + 1
+        times = np.zeros(num_steps)
+        mus = np.zeros((num_probes, num_steps))
+        thetas = np.zeros((num_probes, num_steps))
+
+        with h5py.File(solution_path, "r") as h5file:
+            for i in range(step_min, step_max + 1):
+                grp = h5file[f"data/{i}"]
+                times[i] = float(grp.attrs["time"])
+                mus[:, i] = np.array(grp["mu"])[probe_point_indices]
+                thetas[:, i] = np.angle(np.array(grp["psi"]))[probe_point_indices]
+
+        return DynamicsData(dt=np.diff(times), mu=mus, theta=thetas)
+
     def __eq__(self, other: Any) -> bool:
         return dataclass_equals(self, other)
 

tdgl/test/test_solution.py

@@ -121,6 +121,11 @@ def test_dynamics(solution: tdgl.Solution):
     assert len(time) == (solution.data_range[1] + 1)
     assert solution.closest_solve_step(0) == 0
 
+    _ = DynamicsData.from_solution(solution.path, probe_points=None)
+    _ = DynamicsData.from_solution(
+        solution.path, probe_points=solution.device.probe_points
+    )
+
 
 @pytest.mark.parametrize("dataset", [None, "supercurrent", "normal_current", "invalid"])
 @pytest.mark.parametrize("interp_method", ["linear", "cubic", "invalid"])