Olaflow processing for article

data/processing/projection.py

@@ -60,6 +60,12 @@ x = np.arange(
 )
 theta = np.angle(D.dot((1, 1j)))
 
+log.info(f"N points: {bathy.size:e}")
+S = bathy[:, 0].ptp() * bathy[:, 1].ptp()
+log.info(f"Surface: {S*1e-6:.2f}km^2")
+res = np.sqrt(S / bathy.size)
+log.info(f"Resolution: {res:.2f}m")
+
 coords = artha + (x * np.stack((np.cos(theta), np.sin(theta)))).T
 
 log.info("Interpolating bathymetry in 1D")

data/processing/zero_cross.py

@@ -78,6 +78,7 @@ dj = 0.5
 J = 1 / dj * np.log2(N * dt / s0)
 j = np.arange(0, J)
 sj = s0 * 2 ** (j * dj)
+Tj = 2 * sj * np.pi / 5
 # sj = s0 * np.arange(1, 2 ** (J * dj))
 Mw = sj / dt
 Mlims = sj[[0, -1]]
@@ -99,7 +100,7 @@ ax.set(zorder=1, frame_on=False)
 ax.semilogy()
 a = [t0[0], t0[-1], *Mlims]
 # c = ax.imshow(M, extent=a, aspect="auto", cmap="plasma", vmin=0)
-c = ax.contourf(t, sj, M, cmap="Greys", vmin=0, vmax=v)
+c = ax.contourf(t, Tj, M, cmap="Greys", vmin=0, vmax=v)
 fig.colorbar(c)
 
 H13 = np.quantile(wave, 2 / 3)
@@ -131,7 +132,7 @@ for w, ax2, ax in zip(bigw, ax_[::2].flatten(), ax_[1::2].flatten()):
     log.info(f"Wave [{w}] size: {ws:.2f}m")
     c = ax2.contourf(
         t[i0:i1],
-        sj,
+        Tj,
         M[:, i0:i1],
         cmap="Greys",
         vmin=0,

olaflow/processing/flow_velocity.py

@@ -0,0 +1,137 @@
+import argparse
+import gzip
+from itertools import starmap
+import logging
+from multiprocessing import pool
+import pathlib
+import pickle
+import sys
+
+from cycler import cycler
+import matplotlib.pyplot as plt
+from matplotlib.ticker import MultipleLocator
+import numpy as np
+from scipy import interpolate
+
+from .olaflow import OFModel
+
+parser = argparse.ArgumentParser(description="Post-process olaflow results")
+parser.add_argument("-v", "--verbose", action="count", default=0)
+parser.add_argument(
+    "-o",
+    "--output",
+    type=pathlib.Path,
+    help="Output directory for pickled data",
+    required=True,
+)
+parser.add_argument(
+    "-f",
+    "--func",
+    type=str,
+    help="Post-process function to compare",
+    default="graphUniform",
+    choices=("graphUniform", "graphUniform2"),
+)
+
+args = parser.parse_args()
+
+logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
+log = logging.getLogger("ola_post")
+
+log.info("Plotting comparison of model output")
+
+
+def get_pickle(out):
+    with (
+        path.open("rb")
+        if (path := out.joinpath("pickle")).exists()
+        else gzip.open(path.with_suffix(".gz"), "rb")
+    ) as f:
+        return pickle.load(f)
+
+
+model = get_pickle(args.output)
+
+figsize = 15 / 2.54, 6 / 2.54
+
+fig, ax_ = plt.subplots(
+    2,
+    figsize=figsize,
+    dpi=200,
+    constrained_layout=True,
+)
+
+_ax = ax_[0]
+v = np.nanmax(np.abs(np.where(
+    model.post_fields[args.func]["alpha.water"] > 0.5,
+    #np.linalg.norm(model.post_fields[args.func]["U"], axis=2),
+    model.post_fields[args.func]["U"][..., 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]["U"], axis=2),
+    model.post_fields[args.func]["U"][..., 0],
+    np.nan,
+)))
+_data = model.post_fields[args.func]["U"][..., 0].T
+#_c = _ax.contourf(
+#    model.t,
+#    model.post_fields[args.func]["x_U"],
+#    _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(
+    _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]["x_U"].min(),
+        model.post_fields[args.func]["x_U"].max(),
+    ),
+    vmin=-v150,
+    vmax=v150,
+    aspect="auto",
+)
+_ax.set(xlim=(0, 400))
+_ax.set(facecolor="k")
+_ax.xaxis.set_minor_locator(MultipleLocator(10))
+_ax.yaxis.set_minor_locator(MultipleLocator(1))
+_ax.set(ylabel="z (m)")
+_ax.axes.set_xticklabels([])
+fig.colorbar(_c, label=f"U (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:%}")
+
+
+x = model.post_fields[args.func]["x_U"]
+i0 = np.argmin(np.abs(x - 5))
+_data = model.post_fields[args.func]["U"][..., i0, 0]
+_alpha = model.post_fields[args.func]["alpha.water"][..., i0]
+
+ax = ax_[1]
+ax.fill_between(model.t, np.where(_alpha > 0.5, _data, 0), lw=1, color="#898989", edgecolor="k")
+#ax.autoscale(True, "x", True)
+ax.set(xlim=(0, 400))
+ax.set(xlabel="t (s)", ylabel="U (m/s)")
+ax.grid(c="k", alpha=.2)
+ax.xaxis.set_minor_locator(MultipleLocator(10))
+ax.yaxis.set_minor_locator(MultipleLocator(2))
+
+fig.savefig(
+    args.output.joinpath(
+        f"U_{args.func}.pdf"
+    )
+)
+fig.savefig(
+    args.output.joinpath(
+        f"U_{args.func}.jpg"
+    )
+)

olaflow/processing/snap.py

@@ -0,0 +1,79 @@
+import argparse
+import gzip
+import logging
+import multiprocessing as mp
+import pathlib
+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
+
+from .olaflow import OFModel
+
+parser = argparse.ArgumentParser(description="Post-process olaflow results")
+parser.add_argument("-v", "--verbose", action="count", default=0)
+parser.add_argument(
+    "-o",
+    "--output",
+    type=pathlib.Path,
+    help="Output directory for pickled data",
+    required=True,
+)
+parser.add_argument(
+    "-m",
+    "--max",
+    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)))
+log = logging.getLogger("ola_post")
+
+log.info("Animating olaFlow output")
+out = args.output
+out.mkdir(parents=True, exist_ok=True)
+
+with (
+    path.open("rb")
+    if (path := out.joinpath("pickle")).exists()
+    else gzip.open(path.with_suffix(".gz"), "rb")
+) as f:
+    model = pickle.load(f)
+
+x0, idx0 = np.unique(model.x.astype(np.half), return_inverse=True)
+z0, idz0 = np.unique(model.z.astype(np.half), return_inverse=True)
+
+ix0 = np.argsort(x0)
+iz0 = np.argsort(z0)[::-1]
+
+X, Z = np.meshgrid(x0, z0)
+
+P = np.full((model.t.size, *X.shape), np.nan)
+P[:, iz0[idz0], ix0[idx0]] = model.fields["porosity"]
+
+AW = np.full((model.t.size, *X.shape), np.nan)
+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)
+
+i0 = np.argmin(np.abs(model.t[:, None] - np.asarray((102, 118, 144.5, 176.5))[None, :]), axis=0) 
+
+fig, ax_ = plt.subplots(
+    2, 2, figsize=(15 / 2.54, 4 / 2.54), dpi=200, constrained_layout=True
+)
+for ax, i in zip(ax_.flatten(), i0):
+    ax.imshow(AW[i], cmap="Blues", vmin=0, vmax=1)
+
+fig.savefig(out.joinpath("snap.pdf"))

swash/processing/wavelet.py

@@ -4,6 +4,7 @@ import logging
 import pathlib
 
 import matplotlib.pyplot as plt
+from matplotlib.ticker import MultipleLocator, LogLocator, NullFormatter
 import numpy as np
 import scipy.signal as sgl
 
@@ -49,12 +50,13 @@ 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])
+sig = np.std(watl[:, i0])
 M = np.stack([(np.abs(sgl.cwt(watl[:, i], sgl.morlet2, Mw))/sig)**2 for i in i0])
 v = np.max(M)
+T = 2 * sj * np.pi / 5
 
 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")
+    c = ax_x.contourf(t, T, 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(
@@ -71,6 +73,10 @@ for ax_x, M_, x_ in zip(ax.reshape(-1), M, x[i0]):
     ax_x.set_rasterization_zorder(1.5)
     ax_x.set(ylabel="T (s)", ylim=(sj[0], sj[-1]))
 
+    ax_x.xaxis.set_minor_locator(MultipleLocator(100))
+    ax_x.yaxis.set_major_locator(LogLocator(10, numticks=2**10))
+    ax_x.yaxis.set_minor_locator(LogLocator(10, subs=np.arange(10), numticks=2**10))
+    ax_x.yaxis.set_minor_formatter(NullFormatter())
     if ax_x != ax.reshape(-1)[-1]:
         ax_x.axes.set_xticklabels([])
     else:

swash/processing/zero_cross.py

@@ -76,7 +76,7 @@ ax.grid()
 fig.savefig(out.joinpath("wsize.pdf"))
 
 fig2, ax2 = plt.subplots(
-    figsize=(10 / 2.54, 2 / 3 * 10 / 2.54), constrained_layout=True
+    figsize=(10 / 2.54, 4 / 2.54), constrained_layout=True
 )
 ax2.plot(
     t[cr0[i0 - 5] : cr0[i0 + 7]] * 1e-3,