Merge branch 'swash'

swash/config.ini

@@ -9,8 +9,10 @@ hstru=Hstru.dat
 poro=Poro.dat
 psize=Psize.dat
 out=out_data
+out_nb=out_data_nb
 
 [swash]
+nb=True
 input=sws/INPUT.sws
 path=/opt/swash
 out=out_sws
@@ -18,7 +20,9 @@ mpi=4
 
 [post]
 inp=inp_post
+case=sws_spec_buoy.npz
 out=out_post
+#nperseg=1024
 dt=0.25
 x0=-1000
 t0=180

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")

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 *
 
@@ -25,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
@@ -35,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_eta = np.abs(sgl.welch(eta, f, nperseg=nperseg))
-phi_u = np.abs(sgl.csd(u, u, f))
+phi_u = np.abs(sgl.welch(u, f, nperseg=nperseg))
-phi_eta_u = np.abs(sgl.csd(eta, u, f))
+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])
@@ -53,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()
@@ -68,12 +70,24 @@ 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_fft.plot(
-ax_r.set(xlim=(1, 30), ylim=(0, 1), xlabel="t (s)", ylabel="R")
+    *sgl.welch(eta, 1/dt, nperseg=nperseg),
+    lw=1,
+    c="k",
+    alpha=0.2,
+    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")
+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()
+fig_x, ax_x = plt.subplots(figsize=(10, 1))
 ax_x.plot(data["x"], -data["botl"], color="k")
 ax_x.fill_between(
     data["x"],
@@ -83,13 +97,15 @@ 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}'")
 out.mkdir(parents=True, exist_ok=True)
 
-fig.savefig(out.joinpath(f"t.png"))
+fig.savefig(out.joinpath("t.png"))
-fig_r.savefig(out.joinpath(f"R.png"))
+fig_r.savefig(out.joinpath("R.png"))
-fig_x.savefig(out.joinpath(f"x.png"))
+fig_x.savefig(out.joinpath("x.png"))
 
 log.info("Finished post-processing")

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"):

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

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