From 3e2de12a9292b153481d52b4f14ddfbdbf9fa1a7 Mon Sep 17 00:00:00 2001 From: "Edgar P. Burkhart" Date: Mon, 7 Mar 2022 10:32:36 +0100 Subject: [PATCH 1/7] Swash buoy spectrum input --- swash/sws/SPEC_buoy.sws | 58 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 swash/sws/SPEC_buoy.sws diff --git a/swash/sws/SPEC_buoy.sws b/swash/sws/SPEC_buoy.sws new file mode 100644 index 0000000..96ae124 --- /dev/null +++ b/swash/sws/SPEC_buoy.sws @@ -0,0 +1,58 @@ +$************************* GENERAL *************************************** +PROJ 'GW' 'T1' +SET NAUT +SET LEVEL 0.94 +SET MAXERR 1 +SET DEPMIN 0.001 +MODE DYN ONED +$************************ GRIDS *************************************** +CGRID REG -1251 0 0 1251 0 1250 0 +INPGRID BOT REG -1251 0 0 1250 0 1 0 $x0 y0 theta nx-1 ny-1 dx dy +VERT 10 $nb couches +READ BOTTOM -1 'bathy.dat' 3 0 FREE + +INPGRID PORO REG -1251 0 0 1250 0 1 0 +INPGRID PSIZ REG -1251 0 0 1250 0 1 0 +INPGRID HSTRUC REG -1251 0 0 1250 0 1 0 + +READINP PORO 1 'poro.dat' 3 0 FREE +READINP PSIZ 1 'psize.dat' 3 0 FREE +READINP HSTRUC 1 'hstru.dat' 3 0 FREE + + +$*********************** BOUNDARIES **************************************** +INIT ZERO +BOUN SHAP JON SIG PEAK DSPR DEGR +BOUN SIDE W BTYPE WEAK SMOOT 10 SEC ADDBOUNDWAVE CON SPECT 7.31 16.9 0 0 CYCLE 20 MIN +BOUN SIDE E BTYPE SOMMERFELD + + +$*********************** PHYSICS ******************************************* +BRE 0.6 0.2 +FRICTION MANN 0.08 +POROsity 4 200 1.1 13 +VISC V KEPS + +$*********************** NUMERICS ****************************************** +NONHYD BOX 1.0 PREC ILU +DISCRET UPW MOM +TIMEI 0.2 0.6 + +$*********************** OUTPUTS ******************************************* +SET OUTLEV 1 +$BLOCK 'COMPGRID' HEAD 'test.txt' LAY 4 TSEC XP DEP BOTL WATL PRESS DISCH USTAR VEL VZ VELK ZK BRKP OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'tsec.dat' TSEC OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'xp.dat' XP OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'dep.dat' DEP OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'botl.dat' BOTL OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'watl.dat' WATL OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'press.dat' PRESS OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'disch.dat' DISCH OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'ustar.dat' USTAR OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'vel.dat' VEL OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'vz.dat' VZ OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'velk.dat' VELK OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'zk.dat' ZK OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'brkp.dat' BRKP OUTPUT 000000.00 0.25 SEC +COMPUTE 000000.000 0.015 SEC 003000.000 $timeini dtini unit timefin +STOP From d94577045342ce192f74ab3b562a4ae9ef6ed25b Mon Sep 17 00:00:00 2001 From: "Edgar P. Burkhart" Date: Mon, 7 Mar 2022 10:52:22 +0100 Subject: [PATCH 2/7] Added spectrum to reflection plot --- swash/processing/post.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/swash/processing/post.py b/swash/processing/post.py index 4ad36b1..ec4977c 100644 --- a/swash/processing/post.py +++ b/swash/processing/post.py @@ -7,6 +7,7 @@ import matplotlib.pyplot as plt import numpy as np import pandas as pd import scipy.signal as sgl +import scipy.fft as fft from .read_swash import * @@ -68,9 +69,18 @@ ax_vel.axvline(t0, c="k", alpha=0.2) fig.tight_layout() fig_r, ax_r = plt.subplots() +ax_fft = ax_r.twinx() -ax_r.plot(1 / phi_eta[0, 1:], R[1:], marker="+", lw=1) -ax_r.set(xlim=(1, 30), ylim=(0, 1), xlabel="t (s)", ylabel="R") +ax_fft.plot( + fft.rfftfreq(eta.size, dt), + np.abs(fft.rfft(eta)), + lw=1, + c="k", + alpha=0.2, +) +ax_r.plot(phi_eta[0], R, marker="+") +ax_r.set(xlim=(0, 0.3), ylim=(0, 1), xlabel="f (Hz)", ylabel="R") +ax_fft.set(ylim=0, ylabel="FFT") ax_r.grid() fig_x, ax_x = plt.subplots() @@ -88,8 +98,8 @@ out = pathlib.Path(config.get("post", "out")).joinpath(f"t{t0}x{x0}") log.info(f"Saving plots in '{out}'") out.mkdir(parents=True, exist_ok=True) -fig.savefig(out.joinpath(f"t.png")) -fig_r.savefig(out.joinpath(f"R.png")) -fig_x.savefig(out.joinpath(f"x.png")) +fig.savefig(out.joinpath("t.png")) +fig_r.savefig(out.joinpath("R.png")) +fig_x.savefig(out.joinpath("x.png")) log.info("Finished post-processing") From fa2e30deec3c5c095dbdc885cc64ea69d9defd37 Mon Sep 17 00:00:00 2001 From: "Edgar P. Burkhart" Date: Mon, 7 Mar 2022 12:25:38 +0100 Subject: [PATCH 3/7] Switched fft to psd --- swash/config.ini | 2 ++ swash/processing/post.py | 20 ++++++++++++-------- 2 files changed, 14 insertions(+), 8 deletions(-) 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") From 371acdf9d13ebe8957dab7e4966482a2750f796b Mon Sep 17 00:00:00 2001 From: "Edgar P. Burkhart" Date: Mon, 7 Mar 2022 12:30:50 +0100 Subject: [PATCH 4/7] Changed styling --- swash/processing/post.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/swash/processing/post.py b/swash/processing/post.py index 8a1ea30..fbe32ff 100644 --- a/swash/processing/post.py +++ b/swash/processing/post.py @@ -55,13 +55,13 @@ R = np.sqrt( log.info("Plotting results") fig, (ax_watl, ax_vel) = plt.subplots(2) -ax_watl.plot(t, data["watl"][:, arg_x0], label="watl") +ax_watl.plot(t, data["watl"][:, arg_x0], lw=1, label="watl") ax_watl.set(xlabel="t (s)", ylabel="z (m)") ax_watl.autoscale(axis="x", tight=True) ax_watl.grid() ax_watl.axvline(t0, c="k", alpha=0.2) -ax_vel.plot(t, data["vel"][:, 0, arg_x0], label="vel") +ax_vel.plot(t, data["vel"][:, 0, arg_x0], lw=1, label="vel") ax_vel.set(xlabel="t (s)", ylabel="U (m/s)") ax_vel.autoscale(axis="x", tight=True) ax_vel.grid() @@ -77,7 +77,7 @@ ax_fft.plot( lw=1, c="k", alpha=0.2, - label="PSD", + label="PSD ($\\eta$)", ) 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") @@ -87,7 +87,7 @@ ax_r.legend(loc="upper left") ax_fft.legend(loc="upper right") fig_r.tight_layout() -fig_x, ax_x = plt.subplots() +fig_x, ax_x = plt.subplots(figsize=(10, 1)) ax_x.plot(data["x"], -data["botl"], color="k") ax_x.fill_between( data["x"], @@ -97,6 +97,8 @@ ax_x.fill_between( ax_x.axvline(x0, c="k", alpha=0.2) ax_x.set(xlabel="x (m)", ylabel="z (m)") ax_x.autoscale(axis="x", tight=True) +ax_x.set(aspect="equal") +fig_x.tight_layout() out = pathlib.Path(config.get("post", "out")).joinpath(f"t{t0}x{x0}") log.info(f"Saving plots in '{out}'") From ff7d5380a8957f3001892a3a18ad69b66071b189 Mon Sep 17 00:00:00 2001 From: "Edgar P. Burkhart" Date: Mon, 7 Mar 2022 12:39:00 +0100 Subject: [PATCH 5/7] Add bathy without breakwater --- swash/config.ini | 1 + swash/processing/bathy_nb.py | 94 ++++++++++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+) create mode 100644 swash/processing/bathy_nb.py diff --git a/swash/config.ini b/swash/config.ini index f815c62..560cfa6 100644 --- a/swash/config.ini +++ b/swash/config.ini @@ -9,6 +9,7 @@ hstru=Hstru.dat poro=Poro.dat psize=Psize.dat out=out_data +out_nb=out_data_nb [swash] input=sws/INPUT.sws diff --git a/swash/processing/bathy_nb.py b/swash/processing/bathy_nb.py new file mode 100644 index 0000000..dace50a --- /dev/null +++ b/swash/processing/bathy_nb.py @@ -0,0 +1,94 @@ +import argparse +import configparser +import logging +import pathlib +import sys + +import numpy as np +import pandas as pd + +try: + import matplotlib.pyplot as plt +except ImportError: + plt = None + + +parser = argparse.ArgumentParser(description="Pre-process bathymetry") +parser.add_argument("-v", "--verbose", action="count", default=0) +args = parser.parse_args() + +logging.basicConfig(level=max((10, 20 - 10 * args.verbose))) +log = logging.getLogger("bathy") + +log.info("Starting bathymetry pre-processing") +config = configparser.ConfigParser() +config.read("config.ini") + +root = pathlib.Path(config.get("data", "root")) + +log.info(f"Reading input data from '{root}'") +bathy_hires = np.loadtxt(root.joinpath(config.get("data", "hires"))) +bathy_lores = np.loadtxt(root.joinpath(config.get("data", "bathy"))) +hstru = np.loadtxt(root.joinpath(config.get("data", "hstru"))) +poro = np.loadtxt(root.joinpath(config.get("data", "poro"))) +psize = np.loadtxt(root.joinpath(config.get("data", "psize"))) + +log.info("Generating grid") +x_hires = -np.arange(0, 0.5 * bathy_hires.size, 0.5)[::-1] +x_lores = -np.arange(0, 1 * bathy_lores.size, 1)[::-1] +x_hstru = -np.arange(0, 0.5 * hstru.size, 0.5)[::-1] + +log.info("Generating output data") +bathy_hires_pd = pd.Series(bathy_hires.copy(), index=x_hires) +bathy_lores_pd = pd.Series(bathy_lores.copy(), index=x_lores) + +bathy = pd.DataFrame( + index=bathy_lores_pd.index.union(bathy_hires_pd.index), + columns=("z", "hstru", "poro", "psize"), + data=0, +) +bathy.z[bathy_lores_pd.index] = bathy_lores_pd +bathy.z[bathy_hires_pd.index] = bathy_hires_pd +bathy.z = np.minimum(bathy.z, -15) + +#bathy.loc[x_hstru, ("hstru", "poro", "psize")] = np.array( +# (hstru, poro, psize) +#).T + +bathy = bathy.reindex(bathy_lores_pd.index) +log.debug(f"Bathymetry:\n{bathy}") +log.info( + f"xmin: {bathy.index.min()}, " + f"xmax: {bathy.index.max()}, " + f"n: {bathy.index.size}" +) + +if config.has_option("data", "out_nb"): + out = pathlib.Path(config.get("data", "out_nb")) + log.info(f"Writing output data to '{out}'") + out.mkdir(exist_ok=True) + np.savetxt(out.joinpath("bathy.dat"), bathy.z, newline=" ") + np.savetxt(out.joinpath("hstru.dat"), bathy.hstru, newline=" ") + np.savetxt(out.joinpath("poro.dat"), bathy.poro, newline=" ") + np.savetxt(out.joinpath("psize.dat"), bathy.psize, newline=" ") + + bathy.to_hdf(out.joinpath("bathy.h5"), "bathy", mode="w") + +if config.getboolean("proc", "plot", fallback=False): + if plt is None: + log.error("Could not import PyPlot") + sys.exit(1) + + log.info("Plotting data") + + fig, ax = plt.subplots() + ax.plot(x_hires, bathy_hires, label="High-res") + ax.plot(x_lores, bathy_lores, label="Low-res") + ax.plot(bathy.index, bathy.z, ls="-.", c="k", label="Combined") + ax.plot(bathy.index, bathy.z + bathy.hstru, label="Hstru") + + ax.grid() + ax.legend() + plt.show(block=True) + +log.info("Processing finished") From c4c751721b5ca2dba614baeb80b45eed16bb48f6 Mon Sep 17 00:00:00 2001 From: "Edgar P. Burkhart" Date: Mon, 7 Mar 2022 12:42:11 +0100 Subject: [PATCH 6/7] Add swash model nobreakwater --- swash/config.ini | 1 + swash/processing/swash.py | 6 +++++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/swash/config.ini b/swash/config.ini index 560cfa6..657a41b 100644 --- a/swash/config.ini +++ b/swash/config.ini @@ -12,6 +12,7 @@ out=out_data out_nb=out_data_nb [swash] +nb=True input=sws/INPUT.sws path=/opt/swash out=out_sws diff --git a/swash/processing/swash.py b/swash/processing/swash.py index d73cce8..c86a769 100644 --- a/swash/processing/swash.py +++ b/swash/processing/swash.py @@ -30,8 +30,12 @@ with tempfile.TemporaryDirectory(prefix="swash_", dir=".") as tmp_raw: log.info(f"Copying files to '{tmpdir}'") shutil.copy2(inp, tmpdir) + if config.getboolean("swash", "nb", fallback=False): + path = pathlib.Path(config.get("data", "out_nb")) + else: + path = pathlib.Path(config.get("data", "out")) shutil.copytree( - pathlib.Path(config.get("data", "out")), tmpdir, dirs_exist_ok=True + path, tmpdir, dirs_exist_ok=True ) if config.has_option("swash", "mpi"): From 6fd766f51593bb35b970414ee23876f7e0d6fd97 Mon Sep 17 00:00:00 2001 From: "Edgar P. Burkhart" Date: Mon, 7 Mar 2022 13:24:24 +0100 Subject: [PATCH 7/7] Add swash input with sponge layer --- swash/sws/SPEC_buoy_nb.sws | 59 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100644 swash/sws/SPEC_buoy_nb.sws diff --git a/swash/sws/SPEC_buoy_nb.sws b/swash/sws/SPEC_buoy_nb.sws new file mode 100644 index 0000000..b077878 --- /dev/null +++ b/swash/sws/SPEC_buoy_nb.sws @@ -0,0 +1,59 @@ +$************************* GENERAL *************************************** +PROJ 'GW' 'T1' +SET NAUT +SET LEVEL 0.94 +SET MAXERR 1 +SET DEPMIN 0.001 +MODE DYN ONED +$************************ GRIDS *************************************** +CGRID REG -1251 0 0 1251 0 1250 0 +INPGRID BOT REG -1251 0 0 1250 0 1 0 $x0 y0 theta nx-1 ny-1 dx dy +VERT 10 $nb couches +READ BOTTOM -1 'bathy.dat' 3 0 FREE + +INPGRID PORO REG -1251 0 0 1250 0 1 0 +INPGRID PSIZ REG -1251 0 0 1250 0 1 0 +INPGRID HSTRUC REG -1251 0 0 1250 0 1 0 + +READINP PORO 1 'poro.dat' 3 0 FREE +READINP PSIZ 1 'psize.dat' 3 0 FREE +READINP HSTRUC 1 'hstru.dat' 3 0 FREE + + +$*********************** BOUNDARIES **************************************** +INIT ZERO +BOUN SHAP JON SIG PEAK DSPR DEGR +BOUN SIDE W BTYPE WEAK SMOOT 10 SEC ADDBOUNDWAVE CON SPECT 7.31 16.9 0 0 CYCLE 20 MIN +$BOUN SIDE E BTYPE SOMMERFELD +SPON E 50 + + +$*********************** PHYSICS ******************************************* +BRE 0.6 0.2 +FRICTION MANN 0.08 +POROsity 4 200 1.1 13 +VISC V KEPS + +$*********************** NUMERICS ****************************************** +NONHYD BOX 1.0 PREC ILU +DISCRET UPW MOM +TIMEI 0.2 0.6 + +$*********************** OUTPUTS ******************************************* +SET OUTLEV 1 +$BLOCK 'COMPGRID' HEAD 'test.txt' LAY 4 TSEC XP DEP BOTL WATL PRESS DISCH USTAR VEL VZ VELK ZK BRKP OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'tsec.dat' TSEC OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'xp.dat' XP OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'dep.dat' DEP OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'botl.dat' BOTL OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'watl.dat' WATL OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'press.dat' PRESS OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'disch.dat' DISCH OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'ustar.dat' USTAR OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'vel.dat' VEL OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'vz.dat' VZ OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'velk.dat' VELK OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'zk.dat' ZK OUTPUT 000000.00 0.25 SEC +BLOCK 'COMPGRID' NOHEAD 'brkp.dat' BRKP OUTPUT 000000.00 0.25 SEC +COMPUTE 000000.000 0.015 SEC 003000.000 $timeini dtini unit timefin +STOP