diff --git a/data/config2.ini b/data/config2.ini new file mode 100644 index 0000000..ebda150 --- /dev/null +++ b/data/config2.ini @@ -0,0 +1,29 @@ +[inp] +root=data +base=Database_20220224.xyz +hires=bathyhires.dat +hstru=Hstru.dat +poro=Poro.dat +psize=Psize.dat +raw_ts=cerema/raw/201702281700.raw +raw_spec=cerema/spt/201702281715.spt +hires_step=0.5 +cycle=14400 + +[out] +margin=0.005 +#no_breakwater=True +root=out2 +sub=bathy_sub.npy +out=bathy.npy +step=1 +left=-300 +right=150 + +[artha] +lat=43.398450 +lon=-1.673097 + +[buoy] +lat=43.408333 +lon=-1.681667 diff --git a/data/processing/orbitals.py b/data/processing/orbitals.py index 116231f..e800c7c 100644 --- a/data/processing/orbitals.py +++ b/data/processing/orbitals.py @@ -4,6 +4,7 @@ import logging import pathlib import matplotlib.pyplot as plt +from matplotlib.ticker import MultipleLocator import numpy as np parser = argparse.ArgumentParser(description="Plot orbitals") @@ -85,16 +86,18 @@ fig2dv, ax2dv = plt.subplots(figsize=(5/2.54, 2/3*10/2.54), dpi=200, constrained x0 = ts_flt["x"] * np.cos(theta) + ts_flt["y"] * np.sin(theta) #ax2dv.plot(x0, z0, c="#0066ff", lw=1) ax2dv.quiver( - x0[:-1], - z0[:-1], - np.diff(x0)[:], - np.diff(z0)[:], + x0[:-1] * 1e-2, + z0[:-1] * 1e-2, + np.diff(x0)[:] * 1e-2, + np.diff(z0)[:] * 1e-2, color="k", scale_units="xy", scale=1, ) ax2dv.grid(c="k", alpha=.2) -ax2dv.set(aspect="equal", xlabel="x (cm)", ylabel="z (cm)") +ax2dv.set(aspect="equal", xlabel="x (m)", ylabel="z (m)") +ax2dv.set(ylim=(-10, 10)) +ax2dv.yaxis.set_minor_locator(MultipleLocator(1)) fig2dv.savefig("out_orbitals.pdf") fig2dv.savefig("out_orbitals.jpg") diff --git a/data/processing/projection.py b/data/processing/projection.py index 0aa1e16..2c22998 100644 --- a/data/processing/projection.py +++ b/data/processing/projection.py @@ -131,8 +131,12 @@ np.savetxt(out_root.joinpath("hstru.dat"), hstru[::-1], newline=" ") np.savetxt(out_root.joinpath("poro.dat"), poro[::-1], newline=" ") np.savetxt(out_root.joinpath("psize.dat"), psize[::-1], newline=" ") -fig, ax = plt.subplots() +fig, ax = plt.subplots(figsize=(16 / 2.54, 2 / 3 * 10 / 2.54), constrained_layout=True) ax.plot(-x, z, color="k") ax.fill_between(-x, z + hstru, z, color="k", alpha=0.2) -ax.set_title(f"N={z.size-1}, x=[{-x.max()};{-x.min()}]") +#ax.set_title(f"N={z.size-1}, x=[{-x.max()};{-x.min()}]") +ax.set(ylim=(-30, 15)) +ax.set(xlabel="x (m)", ylabel="z (m)") +ax.autoscale(True, "x", True) +ax.grid() fig.savefig(out_root.joinpath("bathy.pdf")) diff --git a/data/processing/spec.py b/data/processing/spec.py index b8a4b03..a6b08bd 100644 --- a/data/processing/spec.py +++ b/data/processing/spec.py @@ -45,7 +45,7 @@ if cycle is None: f = inp["f"] S = inp["S"] * Sm else: - f = np.arange(inp["f"].min(), inp["f"].max() + 1/cycle, 1/cycle) + f = np.arange(inp["f"].min(), inp["f"].max() + 1 / cycle, 1 / cycle) S = griddata(inp["f"], inp["S"] * Sm, f) with out_spec.open("w") as out: diff --git a/data/processing/ts.py b/data/processing/ts.py index 1a9af04..aacbb81 100644 --- a/data/processing/ts.py +++ b/data/processing/ts.py @@ -4,6 +4,7 @@ import logging import pathlib import matplotlib.pyplot as plt +from matplotlib.ticker import MultipleLocator import numpy as np parser = argparse.ArgumentParser(description="Pre-process time-series") @@ -25,12 +26,14 @@ out_ts = out_root.joinpath("ts.dat") raw_ts = [] for tsi in config.get("inp", "raw_ts").split(","): - raw_ts.append(np.loadtxt( - inp_root.joinpath(tsi), - dtype=[("state", int), ("z", float), ("y", float), ("x", float)], - delimiter=",", - max_rows=2304, - )) + raw_ts.append( + np.loadtxt( + inp_root.joinpath(tsi), + dtype=[("state", int), ("z", float), ("y", float), ("x", float)], + delimiter=",", + max_rows=2304, + ) + ) n = len(raw_ts) raw_ts = np.concatenate(raw_ts) log.debug(f"{raw_ts=}") @@ -39,13 +42,43 @@ if (errs := np.count_nonzero(raw_ts["state"])) != 0: log.warning(f"{errs} transmission errors!") log.debug(f"{dict(zip(*np.unique(raw_ts['state'], return_counts=True)))}") -t = np.linspace(0, 30 * 60 * n, 2304*n+1)[:-1] +t = np.linspace(0, 30 * 60 * n, 2304 * n + 1)[:-1] log.debug(f"{t=}") log.info(f"Saving timeseries to '{out_ts}'") -np.savetxt(out_ts, np.stack((t, raw_ts["z"]/100), axis=1)) +np.savetxt(out_ts, np.stack((t, raw_ts["z"] / 100), axis=1)) -fig, ax = plt.subplots() -ax.plot(t, raw_ts["z"]) -ax.set(xlabel="t (s)", ylabel="z (cm)") +fig, ax = plt.subplots(figsize=(8 / 2.54, 2 / 3 * 10 / 2.54), constrained_layout=True) +tp = np.datetime64("2017-02-28T17:00:00") + t.astype(np.timedelta64)[-(t.size // 3) :] +ax.plot( + tp, + raw_ts["z"][-(t.size // 3) :] * 1e-2, + color="k", + lw=1, +) +ax.axvline( + np.datetime64("2017-02-28T17:00:00") + np.timedelta64(23 * 60 + 8), + color="k", + alpha=0.1, + lw=20, +) +ax.autoscale(True, "x", True) +ax.set(xlabel="t (s)", ylabel="z (m)") +yabs_max = abs(max(ax.get_ylim(), key=abs)) +ax.set(ylim=(-10, 10)) +ax.set( + xticks=( + np.datetime64("2017-02-28T17:20:00"), + np.datetime64("2017-02-28T17:25:00"), + np.datetime64("2017-02-28T17:30:00"), + ), + xticklabels=( + "17:20", + "17:25", + "17:30", + ), +) +ax.yaxis.set_minor_locator(MultipleLocator(1)) +ax.grid(color="k", alpha=0.2) fig.savefig(out_root.joinpath("ts.pdf")) +fig.savefig(out_root.joinpath("ts.jpg"), dpi=200) diff --git a/olaflow/of_ts_fine_1/constant/porosityDict b/olaflow/of_ts_fine_1/constant/porosityDict new file mode 100644 index 0000000..1d10f00 --- /dev/null +++ b/olaflow/of_ts_fine_1/constant/porosityDict @@ -0,0 +1,29 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: 2.1.0 | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object porosityDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Materials: clear region, core, secondary armour layer, primary armour layer +// a,b,c: tuning parameters +a 2(0 50); +b 2(0 1.2); +c 2(0 0.34); + +// D50: mean nominal diameter +D50 2(1 4); +// porosity (phi) +porosity 2(1 0.4); + +// ************************************************************************* // diff --git a/olaflow/of_ts_fine_1/constant/turbulenceProperties b/olaflow/of_ts_fine_1/constant/turbulenceProperties new file mode 100644 index 0000000..187f845 --- /dev/null +++ b/olaflow/of_ts_fine_1/constant/turbulenceProperties @@ -0,0 +1,29 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: 2.1.0 | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object turbulenceProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType RAS; + +RAS +{ + RASModel kOmegaSST; + + turbulence on; + + printCoeffs on; +} + +// ************************************************************************* // diff --git a/olaflow/of_ts_fine_1/system/blockMeshDict b/olaflow/of_ts_fine_1/system/blockMeshDict new file mode 100644 index 0000000..c365fbc --- /dev/null +++ b/olaflow/of_ts_fine_1/system/blockMeshDict @@ -0,0 +1,81 @@ +/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: 1.7.1 | +| \\ / A nd | Web: www.OpenFOAM.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +scale 1; + +vertices +( + (-150 0 -30) + (0 0 -30) + (0 0 30) + (-150 0 30) + (-150 1 -30) + (0 1 -30) + (0 1 30) + (-150 1 30) +); + +blocks +( + hex (0 1 5 4 3 2 6 7) (750 1 300) simpleGrading (1 1 1) +); + +edges +( +); + +boundary +( + inlet + { + type patch; + faces + ( + (0 4 7 3) + ); + } + /*outlet + { + type patch; + faces + ( + (1 5 6 2) + ); + }*/ + wall1 + { + type wall; + faces + ( + (0 1 5 4) + ); + } + atmosphere + { + type patch; + faces + ( + (3 2 6 7) + (1 5 6 2) + ); + } +); + +mergePatchPairs +( +); + +// ************************************************************************* // diff --git a/olaflow/of_ts_fine_1/system/controlDict b/olaflow/of_ts_fine_1/system/controlDict new file mode 100644 index 0000000..ca393bb --- /dev/null +++ b/olaflow/of_ts_fine_1/system/controlDict @@ -0,0 +1,138 @@ +/*---------------------------------------------------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: 1.3 | +| \\ / A nd | Web: http://www.openfoam.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + location "system"; + class dictionary; + object controlDict; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application olaFlow; + +startFrom latestTime; + +startTime 0; + +stopAt endTime; + +endTime 400; + +deltaT 0.1; + +writeControl adjustableRunTime; + +writeInterval 0.5; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +compression on; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable yes; + +adjustTimeStep yes; + +maxCo 0.45; +maxAlphaCo 0.45; + +maxDeltaT 0.5; +/* +functions +{ + gaugesVOF + { + type sets; + libs ("libsampling.so"); + writeControl outputTime; + writeInterval 1; + setFormat raw; + surfaceFormat raw; + interpolationScheme cell; + fields ( alpha.water ); + sets + ( + GaugeVOF01 + { + type lineCellFace; + axis xyz; + start ( 0.5 0.001 0 ); + end ( 0.5 0.001 1.2 ); + } + GaugeVOF02 + { + type lineCellFace; + axis xyz; + start ( 9.25 0.001 0 ); + end ( 9.25 0.001 1.2 ); + } + GaugeVOF03 + { + type lineCellFace; + axis xyz; + start ( 15.75 0.001 0 ); + end ( 15.75 0.001 1.2 ); + } + GaugeVOF04 + { + type lineCellFace; + axis xyz; + start ( 17.75 0.001 0 ); + end ( 17.75 0.001 1.2 ); + } + GaugeVOF05 + { + type lineCellFace; + axis xyz; + start ( 21.1 0.001 0 ); + end ( 21.1 0.001 1.2 ); + } + ); + } + gaugesP + { + type sets; + libs ("libsampling.so"); + writeControl outputTime; + writeInterval 1; + setFormat raw; + surfaceFormat raw; + interpolationScheme cellPointFace; + fields ( p ); + sets + ( + GaugesP + { + type boundaryPoints; + axis xyz; + patches 1(caisson); + points ((18.0 0.01 0.75) + (18.00 0.01 0.80) + (18.00 0.01 0.85) + (18.00 0.01 0.95) + (18.01 0.01 0.70) + (18.25 0.01 0.70) + (18.50 0.01 0.70) + (18.75 0.01 0.70)); + maxDistance 0.01; + } + ); + } +} +*/ +// ************************************************************************* // diff --git a/olaflow/of_ts_fine_1/system/graphUniform b/olaflow/of_ts_fine_1/system/graphUniform new file mode 100644 index 0000000..cb575db --- /dev/null +++ b/olaflow/of_ts_fine_1/system/graphUniform @@ -0,0 +1,25 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 9 + \\/ M anipulation | +------------------------------------------------------------------------------- +Description + Writes graph data for specified fields along a line, specified by start and + end points. A specified number of graph points are used, distributed + uniformly along the line. + +\*---------------------------------------------------------------------------*/ + +start (-50 0.5 -15); +end (-50 0.5 15); +nPoints 100; + +fields (alpha.water U); + +axis z; + +#includeEtc "caseDicts/postProcessing/graphs/graphUniform.cfg" + +// ************************************************************************* // diff --git a/olaflow/of_ts_fine_1/system/graphUniform2 b/olaflow/of_ts_fine_1/system/graphUniform2 new file mode 100644 index 0000000..185937c --- /dev/null +++ b/olaflow/of_ts_fine_1/system/graphUniform2 @@ -0,0 +1,25 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 9 + \\/ M anipulation | +------------------------------------------------------------------------------- +Description + Writes graph data for specified fields along a line, specified by start and + end points. A specified number of graph points are used, distributed + uniformly along the line. + +\*---------------------------------------------------------------------------*/ + +start (-20 0.5 -15); +end (-20 0.5 15); +nPoints 100; + +fields (alpha.water U); + +axis z; + +#includeEtc "caseDicts/postProcessing/graphs/graphUniform.cfg" + +// ************************************************************************* // diff --git a/olaflow/processing/animate.py b/olaflow/processing/animate.py index bf89e4d..165a8d0 100644 --- a/olaflow/processing/animate.py +++ b/olaflow/processing/animate.py @@ -8,6 +8,7 @@ import pickle import matplotlib.pyplot as plt import matplotlib.animation as animation from matplotlib.gridspec import GridSpec +from matplotlib.ticker import MultipleLocator import numpy as np from scipy import interpolate @@ -28,6 +29,12 @@ parser.add_argument( help="Only compute maximum rather than animation", action="store_true", ) +parser.add_argument( + "-i", + "--initial", + help="Only compute initial domain", + action="store_true", +) args = parser.parse_args() logging.basicConfig(level=max((10, 20 - 10 * args.verbose))) @@ -61,11 +68,19 @@ AW[:, iz0[idz0], ix0[idx0]] = model.fields["alpha.water"] U = np.full((model.t.size, *X.shape), np.nan) U[:, iz0[idz0], ix0[idx0]] = np.linalg.norm(model.fields["U"], axis=1) -fig = plt.figure() -gs = GridSpec(3, 1, figure=fig, height_ratios=[1, 0.05, 0.05]) +fig = plt.figure( + figsize=(15 / 2.54, 4 / 2.54), dpi=200, constrained_layout=True +) +gs = GridSpec( + 3 if not args.initial else 1, + 1, + figure=fig, + height_ratios=[1, 0.1, 0.1] if not args.initial else [1], +) ax = fig.add_subplot(gs[0]) -cax1 = fig.add_subplot(gs[1]) -cax2 = fig.add_subplot(gs[2]) +if not args.initial: + cax1 = fig.add_subplot(gs[1]) + cax2 = fig.add_subplot(gs[2]) aw_m = ax.imshow( AW[0], vmin=0, @@ -86,10 +101,15 @@ p_m = ax.imshow( ax.axhline(4.5, ls="-.", lw=1, c="k", alpha=0.2, zorder=1.2) -fig.colorbar(aw_m, label=r"$\alpha_w$", cax=cax1, shrink=0.6, orientation="horizontal") -fig.colorbar(p_m, label=r"Porosity", cax=cax2, shrink=0.6, orientation="horizontal") +if not args.initial: + fig.colorbar( + aw_m, label=r"$\alpha_w$", cax=cax1, shrink=0.6, orientation="horizontal" + ) + fig.colorbar(p_m, label=r"Porosity", cax=cax2, shrink=0.6, orientation="horizontal") ax.set(xlabel="x (m)", ylabel="z (m)", aspect="equal", facecolor="#000000") ax.grid(c="k", alpha=0.2) +ax.xaxis.set_minor_locator(MultipleLocator(5)) +ax.yaxis.set_minor_locator(MultipleLocator(5)) figU = plt.figure() @@ -137,6 +157,13 @@ if args.max: fig.savefig(out.joinpath("max_aw.pdf")) figU.savefig(out.joinpath("max_U.pdf")) +elif args.initial: + aw_m.set_array(AW[0]) + ax.vlines(-20, -15, 15, color="k", lw=1, ls="--", label="Measurements") + ax.text(-20, 15, "Measurements", ha="right", va="bottom") + + fig.savefig(out.joinpath("aw_t0.pdf")) + fig.savefig(out.joinpath("aw_t0.jpg"), dpi=200) else: fig.set(figwidth=19.2, figheight=10.8, dpi=100) figU.set(figwidth=19.2, figheight=10.8, dpi=100) diff --git a/olaflow/processing/animateU.py b/olaflow/processing/animateU.py index 35a9eb5..915cc2e 100644 --- a/olaflow/processing/animateU.py +++ b/olaflow/processing/animateU.py @@ -74,7 +74,7 @@ u_m = axU.quiver( UU[0], alpha=alp[0], cmap="inferno_r", - clim=(0, np.nanmax(UU)), + clim=(0, 20), ) # (wat_p,) = axU.plot(x0, watl[0]) diff --git a/olaflow/processing/diff.py b/olaflow/processing/diff.py index 6263143..321f264 100644 --- a/olaflow/processing/diff.py +++ b/olaflow/processing/diff.py @@ -9,7 +9,7 @@ import sys from cycler import cycler import matplotlib.pyplot as plt -import matplotlib.gridspec as gridspec +from matplotlib.ticker import MultipleLocator import numpy as np from scipy import interpolate @@ -67,11 +67,12 @@ def get_pickle(out): models = list(map(get_pickle, args.output)) +figsize = 15 / 2.54, 4 / 2.54 * len(models) fig, ax_ = plt.subplots( len(models), - figsize=(6, 1.5 * len(models)), - dpi=100, + figsize=figsize, + dpi=200, constrained_layout=True, squeeze=False, ) @@ -90,32 +91,71 @@ if args.timestep is None: ) case "U": for i, (_ax, _model) in enumerate(zip(ax, models)): + v = np.nanmax(np.abs(np.where( + _model.post_fields[args.func]["alpha.water"] > 0.5, + #np.linalg.norm(_model.post_fields[args.func][args.field], axis=2), + _model.post_fields[args.func][args.field][..., 0], + np.nan, + ))) + v150 = np.nanmax(np.abs(np.where( + (_model.post_fields[args.func]["alpha.water"] > 0.5) & (_model.t[:, None] > 170) & (_model.t[:, None] < 200), + #np.linalg.norm(_model.post_fields[args.func][args.field], axis=2), + _model.post_fields[args.func][args.field][..., 0], + np.nan, + ))) + _data = _model.post_fields[args.func][args.field][..., 0].T + #_c = _ax.contourf( + # _model.t, + # _model.post_fields[args.func][f"x_{args.field}"], + # _data, + # cmap="PiYG", + # #levels=[-15, -10, -5, -2, -1, 0, 1, 2, 5, 10, 15], + # vmin=-np.nanmax(np.abs(_data)), + # vmax=np.nanmax(np.abs(_data)), + # extend="both", + #) _c = _ax.imshow( - np.where(_model.post_fields[args.func]["alpha.water"] > 0.5, np.linalg.norm(_model.post_fields[args.func][args.field], axis=2), np.nan).T[::-1], - vmin=0, - vmax=20, - cmap="inferno_r", + _data[::-1], + cmap="PiYG", + alpha=np.clip(_model.post_fields[args.func]["alpha.water"], 0, 1).T[::-1], extent=( _model.t.min(), _model.t.max(), _model.post_fields[args.func][f"x_{args.field}"].min(), _model.post_fields[args.func][f"x_{args.field}"].max(), ), + vmin=-v150, + vmax=v150, + aspect="auto", ) + _ax.set(xlim=(100, 300)) + _ax.set(facecolor="k") + _ax.xaxis.set_minor_locator(MultipleLocator(5)) + _ax.yaxis.set_minor_locator(MultipleLocator(1)) fig.colorbar(_c, label=f"{args.field} (m/s)", ax=_ax) + log.info(f"Vitesse max: {v}m/s") + log.info(f"Vitesse max [170,200]: {v150}m/s") + log.info(f"Écart: {abs(np.nanmax(_data)-17.7)/17.7:%}") case _: log.error(f"Cannot plot field {args.field} from {args.func}") sys.exit(1) for i, (_ax, _model) in enumerate(zip(ax, models)): - _ax.set(xlabel="t (s)", ylabel="z (m)", title=f"Case {i}") - _ax.grid(color="k", alpha=0.2) + _ax.set(xlabel="t (s)", ylabel="z (m)") + if len(models) > 1: + _ax.set(title=f"Case {i}") + #_ax.grid(color="#797979", alpha=0.5) fig.savefig( args.output[0].joinpath( f"diff_{args.func}_{args.field}_{'_'.join([o.name for o in args.output])}.pdf" ) ) + fig.savefig( + args.output[0].joinpath( + f"diff_{args.func}_{args.field}_{'_'.join([o.name for o in args.output])}.jpg" + ) + ) else: match args.field: case "alpha.water": @@ -123,7 +163,9 @@ else: _ax.tricontour( _model.x, _model.z, - _model.fields[args.field][np.where(_model.t == args.timestep)[0]][0], + _model.fields[args.field][np.where(_model.t == args.timestep)[0]][ + 0 + ], levels=(0.5,), colors="k", ) diff --git a/swash/config-ts.ini b/swash/config-ts.ini index f0a3752..3371d49 100644 --- a/swash/config-ts.ini +++ b/swash/config-ts.ini @@ -10,5 +10,5 @@ mpi=8 [post] inp=inp_post/ts_4lay_1h out=out_post/ts_4lay_1h -x0=-1250 +x0=-50 t0=180 diff --git a/swash/config.ini b/swash/config.ini index 647c623..3e21397 100644 --- a/swash/config.ini +++ b/swash/config.ini @@ -9,7 +9,7 @@ mpi=8 [post] inp=inp_post/real_spec_interp -#compare=inp_post/real_spec_interp_nb +compare=inp_post/real_spec_interp_nb out=out_post/real_spec_interp x0=-1250 t0=180 diff --git a/swash/processing/post.py b/swash/processing/post.py index 92f31eb..cc24b73 100644 --- a/swash/processing/post.py +++ b/swash/processing/post.py @@ -7,8 +7,6 @@ import matplotlib.pyplot as plt import numpy as np import scipy.signal as sgl -from .read_swash import * - 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") @@ -30,7 +28,8 @@ t = np.load(inp.joinpath("t.npy")) botl = np.load(inp.joinpath("botl.npy")) watl = np.load(inp.joinpath("watl.npy")) -vel = np.load(inp.joinpath("vel_x.npy")) +#vel = np.load(inp.joinpath("vel_x.npy")) +vel = np.load(inp.joinpath("vel.npy"))[0] # Cospectral calculations x0 = config.getint("post", "x0") @@ -65,22 +64,25 @@ R = np.sqrt( if config.has_option("post", "compare"): inp_comp = pathlib.Path(config.get("post", "compare")) - x_ = np.load(inp_comp.joinpath("xp.npy")) - t_ = np.load(inp_comp.joinpath("tsec.npy")) + x_ = np.load(inp_comp.joinpath("x.npy")) + t_ = np.load(inp_comp.joinpath("t.npy")) botl_ = np.load(inp_comp.joinpath("botl.npy")) watl_ = np.load(inp_comp.joinpath("watl.npy")) vel_ = np.load(inp_comp.joinpath("vel_x.npy")) + # Cospectral calculations arg_x0_ = np.abs(x_ - x0).argmin() arg_t0_ = np.abs(t_ - t0).argmin() + dt_ = np.diff(t_).mean() * 1e-3 + f_ = 1 / dt_ eta_ = watl_[t_ > t0, arg_x0_] u_ = vel_[t_ > t0, arg_x0_] - phi_eta_ = sgl.welch(eta_, f, nperseg=nperseg) - phi_u_ = sgl.welch(u_, f, nperseg=nperseg) - phi_eta_u_ = sgl.csd(eta_, u_, f, nperseg=nperseg) + phi_eta_ = sgl.welch(eta_, f_, nperseg=nperseg) + phi_u_ = sgl.welch(u_, f_, nperseg=nperseg) + phi_eta_u_ = sgl.csd(eta_, u_, f_, nperseg=nperseg) H_ = np.sqrt(np.abs(phi_eta_[1])) U_ = np.sqrt(np.abs(phi_u_[1])) @@ -91,10 +93,6 @@ if config.has_option("post", "compare"): (np.abs(phi_eta_[1]) + np.abs(phi_u_[1]) - 2 * phi_eta_u_[1].real) / (np.abs(phi_eta_[1]) + np.abs(phi_u_[1]) + 2 * phi_eta_u_[1].real) ) - # R_ = np.sqrt( - # (1 + G_**2 - 2 * G_ * np.cos(th_eta_u_)) - # / (1 + G_**2 + 2 * G_ * np.cos(th_eta_u_)) - # ) # Plotting diff --git a/swash/processing/transmission.py b/swash/processing/transmission.py index b272a2d..bbd47ba 100644 --- a/swash/processing/transmission.py +++ b/swash/processing/transmission.py @@ -30,7 +30,7 @@ botl = np.load(inp.joinpath("botl.npy")) watl = np.load(inp.joinpath("watl.npy")) vel = np.load(inp.joinpath("vel.npy"))[0] -t0 = np.linspace(23 * 60 + 8, 23 * 60 + 8 + 100, 5) +t0 = np.linspace(23 * 60 + 8, 23 * 60 + 8 + 100, 6) # Plotting log.info("Plotting results") @@ -38,10 +38,10 @@ log.info("Plotting results") vlim = np.nanmin(np.maximum(watl, -botl)), np.nanmax(np.maximum(watl, -botl)) fig_x, ax = plt.subplots( - len(t0), figsize=(10 / 2.54, 4 / 3 * 10 / 2.54), constrained_layout=True + 3, 2, figsize=(15 / 2.54, 2 / 3 * 10 / 2.54), constrained_layout=True ) i0 = np.argmin(np.abs(t * 1e-3 - t0[0])) -for ax_x, t0_x in zip(ax, t0): +for ax_x, t0_x in zip(ax.reshape(-1), t0): ax_x.plot(x, -botl, color="k") i = np.argmin(np.abs(t * 1e-3 - t0_x)) ax_x.plot( @@ -67,5 +67,6 @@ log.info(f"Saving plots in '{out}'") out.mkdir(parents=True, exist_ok=True) fig_x.savefig(out.joinpath("x.pdf")) +fig_x.savefig(out.joinpath("x.jpg"), dpi=200) log.info("Finished post-processing") diff --git a/swash/processing/wavelet.py b/swash/processing/wavelet.py new file mode 100644 index 0000000..d4fe66b --- /dev/null +++ b/swash/processing/wavelet.py @@ -0,0 +1,87 @@ +import argparse +import configparser +import logging +import pathlib + +import matplotlib.pyplot as plt +import numpy as np +import scipy.signal as sgl + +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")) +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")) * 1e-3 + +botl = np.load(inp.joinpath("botl.npy")) +watl = np.load(inp.joinpath("watl.npy")) +vel = np.load(inp.joinpath("vel.npy"))[0] + +# Plotting +log.info("Plotting results") + +vlim = np.nanmin(np.maximum(watl, -botl)), np.nanmax(np.maximum(watl, -botl)) + +x0 = np.linspace(-600, -200, 5) +i0 = np.argmin(np.abs(x[:, None] - x0), axis=0) + +fig_x, ax = plt.subplots( + 5, 1, figsize=(15 / 2.54, 15/ 2.54), constrained_layout=True +) +dt = np.mean(np.diff(t)) +N = t.size +s0 = 2 * dt +dj = 0.5 +J = 1 / dj * np.log2(N * dt / s0) +j = np.arange(0, J) +sj = s0 * 2 ** (j * dj) +Mw = sj / dt +sig = np.var(watl[:, i0]) +M = np.stack([(np.abs(sgl.cwt(watl[:, i], sgl.morlet2, Mw))/sig)**2 for i in i0]) +v = np.max(M) + +for ax_x, M_, x_ in zip(ax.reshape(-1), M, x[i0]): + c = ax_x.contourf(t, sj, M_, cmap="Greys", vmin=0, levels=[1, 2.5, 5, 10, 20, 40], extend="both") + fig_x.colorbar(c, ax=ax_x, label="NWPS") + ax_x.grid(color="k", alpha=0.2) + ax_x.text( + 0.95, + 0.95, + f"x={x_:.0f}m", + horizontalalignment="right", + verticalalignment="top", + transform=ax_x.transAxes, + #c="w", + ) + ax_x.semilogy() + ax_x.autoscale(True, "both", True) + ax_x.set_rasterization_zorder(1.5) + ax_x.set(ylabel="T (s)", ylim=(sj[0], sj[-1])) + + if ax_x != ax.reshape(-1)[-1]: + ax_x.axes.set_xticklabels([]) + else: + ax_x.set(xlabel="t (s)") + +out = pathlib.Path(config.get("post", "out")).joinpath(f"trans") +log.info(f"Saving plots in '{out}'") +out.mkdir(parents=True, exist_ok=True) + +fig_x.savefig(out.joinpath("wavelet.pdf"), dpi=300) +fig_x.savefig(out.joinpath("wavelet.jpg"), dpi=200) +fig_x.savefig(out.joinpath("wavelet.png"), dpi=200) + +log.info("Finished post-processing") diff --git a/swash/processing/zero_cross.py b/swash/processing/zero_cross.py index 5d5c0e5..2a4e45c 100644 --- a/swash/processing/zero_cross.py +++ b/swash/processing/zero_cross.py @@ -38,6 +38,9 @@ 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( @@ -72,9 +75,22 @@ ax.autoscale(True, "x", True) ax.grid() fig.savefig(out.joinpath("wsize.pdf")) -fig2, ax2 = plt.subplots(figsize=(10/2.54, 2/3*10/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="Case 1") -ax2.plot(t[cr0[i0 - 5] : cr0[i0 + 7]] * 1e-3, w0_comp[cr0[i0 - 5] : cr0[i0 + 7]], color="k", ls="-.", label="Case 2") +fig2, ax2 = plt.subplots( + figsize=(10 / 2.54, 2 / 3 * 10 / 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() @@ -82,7 +98,11 @@ 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"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():%}") +log.info( + f"Maximum wave size difference: {abs(wave_comp.max()-wave.max())/wave.max():%}" +)