131 lines
4.6 KiB
Python
131 lines
4.6 KiB
Python
"""CPU-only cycle extraction. All criteria frozen before any Re100 GPU result."""
|
|
|
|
import numpy as np
|
|
|
|
GATES = dict(
|
|
min_tstar=250.0,
|
|
cycles=5,
|
|
min_period=2.0,
|
|
max_period=5.0,
|
|
min_lift_amplitude=0.2,
|
|
period_relative_range=0.005,
|
|
extrema_absolute_range=0.005,
|
|
mass_relative_range=1e-4,
|
|
cycle_mean_flux_relative=0.005,
|
|
)
|
|
REFERENCE_GATES = dict(
|
|
cd_relative=0.02, cl_absolute=0.03, dp_relative=0.02, st_relative=0.02
|
|
)
|
|
|
|
|
|
def minima(t, y):
|
|
out = []
|
|
for i in range(1, len(y) - 1):
|
|
if y[i] < y[i - 1] and y[i] <= y[i + 1]:
|
|
# Quadratic sub-sample location; only separated physical cycles accepted.
|
|
coef = np.polyfit(t[i - 1 : i + 2] - t[i], y[i - 1 : i + 2], 2)
|
|
v = t[i] - coef[1] / (2 * coef[0]) if coef[0] > 0 else t[i]
|
|
if t[i - 1] <= v <= t[i + 1] and (not out or v - out[-1] > 1.0):
|
|
out.append(float(v))
|
|
return out
|
|
|
|
|
|
def extract(rows, min_time=250.0):
|
|
if len(rows) < 3:
|
|
return []
|
|
t = np.array([r["tstar_force"] for r in rows])
|
|
cl = np.array([r["cl"] for r in rows])
|
|
ends = minima(t, cl)
|
|
cycles = []
|
|
for a, b in zip(ends, ends[1:]):
|
|
if a < min_time or not GATES["min_period"] < b - a < GATES["max_period"]:
|
|
continue
|
|
grid = np.linspace(a, b, 257)
|
|
c = dict(start=a, end=b, period=b - a, st=1 / (b - a))
|
|
for key in ("cd", "cl", "dp", "mass", "flux_in", "flux_out"):
|
|
times = np.array(
|
|
[r["tstar" if key not in ("cd", "cl") else "tstar_force"] for r in rows]
|
|
)
|
|
# No endpoint extrapolation: pressure/mass must cover the whole cycle.
|
|
if a < times[0] or b > times[-1]:
|
|
break
|
|
v = np.interp(grid, times, [r[key] for r in rows])
|
|
c[key + "_min"] = float(v.min())
|
|
c[key + "_max"] = float(v.max())
|
|
c[key + "_midrange"] = float((v.min() + v.max()) / 2)
|
|
c[key + "_amplitude"] = float(v.max() - v.min())
|
|
c[key + "_time_mean"] = float(
|
|
np.sum((v[1:] + v[:-1]) * np.diff(grid) / 2) / (b - a)
|
|
)
|
|
if key in ("cd", "cl", "dp"):
|
|
c[key + "_waveform"] = v.tolist()
|
|
else:
|
|
c["flux_relative"] = abs(
|
|
c["flux_in_time_mean"] - c["flux_out_time_mean"]
|
|
) / max(abs(c["flux_in_time_mean"]), 1e-30)
|
|
cycles.append(c)
|
|
return cycles
|
|
|
|
|
|
def assess(rows):
|
|
cycles = extract(rows, min_time=GATES["min_tstar"])
|
|
tail = cycles[-GATES["cycles"] :]
|
|
checks = {}
|
|
checks["enough_cycles"] = len(tail) == GATES["cycles"]
|
|
if checks["enough_cycles"]:
|
|
checks["consecutive"] = all(
|
|
abs(a["end"] - b["start"]) < 1e-8 for a, b in zip(tail, tail[1:])
|
|
)
|
|
checks["recent"] = (
|
|
0 <= rows[-1]["tstar_force"] - tail[-1]["end"] <= GATES["max_period"]
|
|
)
|
|
checks["shedding"] = (
|
|
min(c["cl_amplitude"] for c in tail) >= GATES["min_lift_amplitude"]
|
|
)
|
|
periods = [c["period"] for c in tail]
|
|
checks["period"] = (
|
|
np.ptp(periods) / np.mean(periods) <= GATES["period_relative_range"]
|
|
)
|
|
checks["extrema"] = all(
|
|
np.ptp([c[k + "_" + e] for c in tail]) <= GATES["extrema_absolute_range"]
|
|
for k in ("cd", "cl", "dp")
|
|
for e in ("min", "max")
|
|
)
|
|
mass = [c["mass_time_mean"] for c in tail]
|
|
checks["mass"] = (
|
|
np.ptp(mass) / abs(np.mean(mass)) <= GATES["mass_relative_range"]
|
|
)
|
|
checks["flux"] = (
|
|
max(c["flux_relative"] for c in tail) <= GATES["cycle_mean_flux_relative"]
|
|
)
|
|
return dict(
|
|
periodic=bool(checks["enough_cycles"] and all(checks.values())),
|
|
checks={k: bool(v) for k, v in checks.items()},
|
|
cycles=cycles,
|
|
)
|
|
|
|
|
|
def compare(cycle, ref):
|
|
errors = {}
|
|
passed = True
|
|
for key in ("cd", "cl", "dp"):
|
|
for ext in ("min", "max"):
|
|
k = key + "_" + ext
|
|
e = abs(cycle[k] - ref[k])
|
|
e = e if key == "cl" else e / abs(ref[k])
|
|
errors[k] = e
|
|
passed &= (
|
|
e
|
|
<= REFERENCE_GATES[key + ("_absolute" if key == "cl" else "_relative")]
|
|
)
|
|
errors[key + "_waveform_rms"] = float(
|
|
np.sqrt(
|
|
np.mean(
|
|
(np.array(cycle[key + "_waveform"]) - ref[key + "_waveform"]) ** 2
|
|
)
|
|
)
|
|
)
|
|
errors["st_relative"] = abs(cycle["st"] - ref["st"]) / ref["st"]
|
|
passed &= errors["st_relative"] <= REFERENCE_GATES["st_relative"]
|
|
return dict(reference_agreement=bool(passed), errors=errors)
|