diff --git a/swash/config.ini b/swash/config.ini index 08df32a..f815c62 100644 --- a/swash/config.ini +++ b/swash/config.ini @@ -18,7 +18,9 @@ mpi=4 [post] inp=inp_post +case=sws_spec_buoy.npz out=out_post +#nperseg=1024 dt=0.25 x0=-1000 t0=180 diff --git a/swash/processing/post.py b/swash/processing/post.py index ec4977c..8a1ea30 100644 --- a/swash/processing/post.py +++ b/swash/processing/post.py @@ -26,7 +26,7 @@ inp = pathlib.Path(config.get("post", "inp")) root = pathlib.Path(config.get("swash", "out")) log.info(f"Reading bathymetry from '{inp}'") -data = np.load(inp.joinpath("sws.npz")) +data = np.load(inp.joinpath(config.get("post", "case"))) x, t = data["x"], data["t"] # Cospectral calculations @@ -36,14 +36,15 @@ t0 = config.getfloat("post", "t0") arg_t0 = np.abs(t - t0).argmin() dt = config.getfloat("post", "dt") f = 1 / dt +nperseg = config.getint("post", "nperseg", fallback=None) log.info(f"Computing reflection coefficient at x={x0}") eta = data["watl"][t > t0, arg_x0] u = data["vel"][t > t0, 0, arg_x0] -phi_eta = np.abs(sgl.csd(eta, eta, f)) -phi_u = np.abs(sgl.csd(u, u, f)) -phi_eta_u = np.abs(sgl.csd(eta, u, f)) +phi_eta = np.abs(sgl.welch(eta, f, nperseg=nperseg)) +phi_u = np.abs(sgl.welch(u, f, nperseg=nperseg)) +phi_eta_u = np.abs(sgl.csd(eta, u, f, nperseg=nperseg)) R = np.sqrt( (phi_eta[1] + phi_u[1] - 2 * phi_eta_u[1]) @@ -72,16 +73,19 @@ fig_r, ax_r = plt.subplots() ax_fft = ax_r.twinx() ax_fft.plot( - fft.rfftfreq(eta.size, dt), - np.abs(fft.rfft(eta)), + *sgl.welch(eta, 1/dt, nperseg=nperseg), lw=1, c="k", alpha=0.2, + label="PSD", ) -ax_r.plot(phi_eta[0], R, marker="+") +ax_r.plot(phi_eta[0], R, marker="+", label="R") ax_r.set(xlim=(0, 0.3), ylim=(0, 1), xlabel="f (Hz)", ylabel="R") -ax_fft.set(ylim=0, ylabel="FFT") +ax_fft.set(ylim=0, ylabel="PSD (m²/Hz)") ax_r.grid() +ax_r.legend(loc="upper left") +ax_fft.legend(loc="upper right") +fig_r.tight_layout() fig_x, ax_x = plt.subplots() ax_x.plot(data["x"], -data["botl"], color="k")