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():%}" )