1
Fork 0
internship/swash/processing/zero_cross.py

109 lines
2.8 KiB
Python

import argparse
import configparser
import logging
import pathlib
import matplotlib.pyplot as plt
import numpy as np
parser = argparse.ArgumentParser(description="Post-process swash output")
parser.add_argument("-v", "--verbose", action="count", default=0)
parser.add_argument("-c", "--config", default="config.ini")
args = parser.parse_args()
logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
log = logging.getLogger("post")
log.info("Starting post-processing")
config = configparser.ConfigParser()
config.read(args.config)
inp = pathlib.Path(config.get("post", "inp"))
inp_comp = pathlib.Path(config.get("post", "compare"))
root = pathlib.Path(config.get("swash", "out"))
log.info(f"Reading data from '{inp}'")
x = np.load(inp.joinpath("x.npy"))
t = np.load(inp.joinpath("t.npy"))
watl = np.load(inp.joinpath("watl.npy"))
watl_comp = np.load(inp_comp.joinpath("watl.npy"))
# Cospectral calculations
x0 = config.getint("post", "x0")
arg_x0 = np.abs(x - x0).argmin()
w0 = watl[:, arg_x0]
w0_comp = watl_comp[:, arg_x0]
cr0 = np.where(np.diff(np.sign(w0)))[0]
cr0_comp = np.where(np.diff(np.sign(w0_comp)))[0]
log.info(f"1: {cr0.size} waves")
log.info(f"2: {cr0_comp.size} waves")
wave = np.fromiter(
(
np.abs(
np.max(np.abs(w0[cr0[i - 1] : cr0[i]]))
+ np.max(np.abs(w0[cr0[i] : cr0[i + 1]]))
)
for i in range(1, len(cr0) - 1)
),
dtype=np.single,
)
wave_comp = np.fromiter(
(
np.abs(
np.max(np.abs(w0_comp[cr0_comp[i - 1] : cr0_comp[i]]))
+ np.max(np.abs(w0_comp[cr0_comp[i] : cr0_comp[i + 1]]))
)
for i in range(1, len(cr0) - 1)
),
dtype=np.single,
)
i0 = np.argmax(wave)
out = pathlib.Path(config.get("post", "out")).joinpath(f"x{x0}")
log.info(f"Saving plots in '{out}'")
out.mkdir(parents=True, exist_ok=True)
fig, ax = plt.subplots()
ax.plot(t[cr0[1:-1]] * 1e-3, wave)
ax.set(xlabel="t (s)", ylabel="z (m)")
ax.autoscale(True, "x", True)
ax.grid()
fig.savefig(out.joinpath("wsize.pdf"))
fig2, ax2 = plt.subplots(
figsize=(10 / 2.54, 4 / 2.54), constrained_layout=True
)
ax2.plot(
t[cr0[i0 - 5] : cr0[i0 + 7]] * 1e-3,
w0[cr0[i0 - 5] : cr0[i0 + 7]],
color="k",
label="Cas 1",
)
ax2.plot(
t[cr0[i0 - 5] : cr0[i0 + 7]] * 1e-3,
w0_comp[cr0[i0 - 5] : cr0[i0 + 7]],
color="k",
ls="-.",
label="Cas 2",
)
ax2.set(xlabel="t (s)", ylabel="z (m)")
ax2.autoscale(True, "x", True)
ax2.grid()
ax2.legend()
fig2.savefig(out.joinpath("maxw.pdf"))
fig2.savefig(out.joinpath("maxw.jpg"), dpi=200)
log.info(
f"RMS difference: {np.sqrt(np.mean((w0_comp-w0)**2))}m ; {np.sqrt(np.mean((w0_comp-w0)**2))/(w0.max()-w0.min()):%}"
)
log.info(f"Bias: {np.mean(w0_comp-w0)}m")
log.info(f"Maximum wave size: {wave.max()}m ; {wave_comp.max()}m")
log.info(
f"Maximum wave size difference: {abs(wave_comp.max()-wave.max())/wave.max():%}"
)