Skills › Scientific Agent Skills › openpiv › analyze.py
analyze.py
153 lines · skills/openpiv/scripts/analyze.py at 5f2bc03
1"""Post-processing helpers for a params.npz written by runner.py.
3Verified against openpiv 0.26.1. Pure numpy -- no OpenPIV import needed here.
4"""
6from pathlib import Path
7from typing import Dict, Optional, Tuple
9import numpy as np
12class PIVAnalyzer:
13 """Derived quantities from a saved PIV velocity field.
15 The fields in params.npz are already scaled to physical units by runner.py, so
16 x and y are in the same unit as the scaling factor and u and v are that unit per
17 second. Gradient methods infer spacing and orientation from saved coordinates.
18 Derivatives use interpolated fluid vectors; permanently masked regions stay NaN.
19 """
21 def __init__(self, params_file: str):
22 self.params_file = Path(params_file)
24 with np.load(self.params_file, allow_pickle=False) as data:
25 self.x = data["x"]
26 self.y = data["y"]
27 self.u = data["u"]
28 self.v = data["v"]
29 self.flags = data["flags"].astype(bool)
30 self.mask = data["mask"].astype(bool) if "mask" in data else np.zeros_like(self.flags)
31 self.u = np.where(self.mask, np.nan, self.u)
32 self.v = np.where(self.mask, np.nan, self.v)
34 @property
35 def grid_spacing(self) -> Tuple[float, float]:
36 """(dx, dy) inferred from the coordinate arrays, in physical units."""
37 dx = float(np.abs(np.diff(self.x, axis=1)).mean()) if self.x.shape[1] > 1 else 1.0
38 dy = float(np.abs(np.diff(self.y, axis=0)).mean()) if self.y.shape[0] > 1 else 1.0
39 return dx, dy
41 @property
42 def axis_signs(self) -> Tuple[float, float]:
43 """(+-1, +-1): whether x and y increase or decrease along the array axes.
45 runner.py converts to a physical right-handed (y-up) frame while leaving rows in
46 image order -- so physical y *decreases* as the row index grows. np.gradient
47 only sees the array, so differentiating with a positive spacing would return
48 -du/dy there and change the vorticity and shear-strain contributions.
49 These signs put the derivatives back on the physical axes.
50 """
51 x_sign = -1.0 if self.x.shape[1] > 1 and self.x[0, 1] < self.x[0, 0] else 1.0
52 y_sign = -1.0 if self.y.shape[0] > 1 and self.y[1, 0] < self.y[0, 0] else 1.0
53 return x_sign, y_sign
55 def _steps(
56 self, dx: Optional[float], dy: Optional[float]
57 ) -> Tuple[float, float, float, float]:
58 """(dx, dy, x_sign, y_sign) for a gradient call, defaulting to the grid."""
59 gx, gy = self.grid_spacing
60 x_sign, y_sign = self.axis_signs
61 return (
62 gx if dx is None else abs(float(dx)),
63 gy if dy is None else abs(float(dy)),
64 x_sign,
65 y_sign,
66 )
68 def plot_vector_field(
69 self,
70 scale: int = 50,
71 width: float = 0.0035,
72 save_path: Optional[str] = None,
73 ):
74 """Quiver plot of the valid vectors. Saves to save_path, or shows interactively."""
75 import matplotlib.pyplot as plt
77 valid = ~(self.flags | self.mask) & np.isfinite(self.u) & np.isfinite(self.v)
78 fig, ax = plt.subplots(figsize=(10, 8))
79 ax.quiver(
80 self.x[valid],
81 self.y[valid],
82 self.u[valid],
83 self.v[valid],
84 scale=scale,
85 width=width,
86 )
87 ax.set_xlabel("x")
88 ax.set_ylabel("y")
89 ax.set_title("Velocity Field")
90 ax.set_aspect("equal")
91 if save_path:
92 fig.savefig(save_path, dpi=150, bbox_inches="tight")
93 plt.close(fig)
94 else:
95 plt.show()
96 return fig
98 def get_velocity_magnitude(self) -> np.ndarray:
99 return np.sqrt(self.u**2 + self.v**2)
101 def compute_vorticity(
102 self, dx: Optional[float] = None, dy: Optional[float] = None
103 ) -> np.ndarray:
104 """Out-of-plane vorticity, dv/dx - du/dy. Defaults to the inferred grid spacing.
106 dx and dy are spacing magnitudes; the axis orientation comes from the saved
107 coordinates (see axis_signs), so a counter-clockwise flow gives positive
108 vorticity whichever way the rows run.
109 """
110 dx, dy, x_sign, y_sign = self._steps(dx, dy)
111 dv_dx = x_sign * np.gradient(self.v, dx, axis=1)
112 du_dy = y_sign * np.gradient(self.u, dy, axis=0)
113 return dv_dx - du_dy
115 def compute_strain(
116 self, dx: Optional[float] = None, dy: Optional[float] = None
117 ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
118 """Return (exx, eyy, exy) of the 2D strain-rate tensor."""
119 dx, dy, x_sign, y_sign = self._steps(dx, dy)
120 du_dx = x_sign * np.gradient(self.u, dx, axis=1)
121 du_dy = y_sign * np.gradient(self.u, dy, axis=0)
122 dv_dx = x_sign * np.gradient(self.v, dx, axis=1)
123 dv_dy = y_sign * np.gradient(self.v, dy, axis=0)
124 return du_dx, dv_dy, 0.5 * (du_dy + dv_dx)
126 def compute_statistics(self, include_interpolated: bool = False) -> Dict[str, float]:
127 """Spatial mean and RMS over this single frame.
129 This is NOT Reynolds decomposition: subtracting one frame's spatial mean
130 measures spatial variance, which equals turbulent intensity only for a
131 homogeneous field. True turbulence statistics need an ensemble of pairs --
132 average over the time axis, then subtract that mean field from each frame.
133 By default exclude flagged/interpolated vectors and masked obstacles.
134 The legacy key 'tke' is only two-component spatial variance energy, not
135 three-component turbulent kinetic energy.
136 """
137 valid = ~self.mask & np.isfinite(self.u) & np.isfinite(self.v)
138 if not include_interpolated:
139 valid &= ~self.flags
140 if not valid.any():
141 raise ValueError("No valid fluid vectors remain for statistics")
142 u, v = self.u[valid], self.v[valid]
143 u_prime = u - np.mean(u)
144 v_prime = v - np.mean(v)
145 rms_u = float(np.nanstd(u_prime))
146 rms_v = float(np.nanstd(v_prime))
147 return {
148 "u_mean": float(np.mean(u)),
149 "v_mean": float(np.mean(v)),
150 "rms_u": rms_u,
151 "rms_v": rms_v,
152 "tke": 0.5 * (rms_u**2 + rms_v**2),
153 }