Merge branch 'master' of ssh://git.edgarpierre.fr:39529/m2cce/internship

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

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],
+    x0[:-1] * 1e-2,
-    z0[:-1],
+    z0[:-1] * 1e-2,
-    np.diff(x0)[:],
+    np.diff(x0)[:] * 1e-2,
-    np.diff(z0)[:],
+    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")
 

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

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")
@@ -47,7 +48,37 @@ 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))
 
-fig, ax = plt.subplots()
+fig, ax = plt.subplots(figsize=(8 / 2.54, 2 / 3 * 10 / 2.54), constrained_layout=True)
-ax.plot(t, raw_ts["z"])
+tp = np.datetime64("2017-02-28T17:00:00") + t.astype(np.timedelta64)[-(t.size // 3) :]
-ax.set(xlabel="t (s)", ylabel="z (cm)")
+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)

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);
+
+// ************************************************************************* //

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;
+}
+
+// ************************************************************************* //

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
+(
+);
+
+// ************************************************************************* //

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;
+            }
+        );
+    }
+}
+*/
+// ************************************************************************* //

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"
+
+// ************************************************************************* //

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"
+
+// ************************************************************************* //

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,9 +68,17 @@ 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()
+fig = plt.figure(
-gs = GridSpec(3, 1, figure=fig, height_ratios=[1, 0.05, 0.05])
+    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])
+if not args.initial:
     cax1 = fig.add_subplot(gs[1])
     cax2 = fig.add_subplot(gs[2])
 aw_m = ax.imshow(
@@ -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")
+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)

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

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)),
+    figsize=figsize,
-    dpi=100,
+    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],
+                    _data[::-1],
-                    vmin=0,
+                    cmap="PiYG",
-                    vmax=20,
+                    alpha=np.clip(_model.post_fields[args.func]["alpha.water"], 0, 1).T[::-1],
-                    cmap="inferno_r",
                     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.set(xlabel="t (s)", ylabel="z (m)")
-        _ax.grid(color="k", alpha=0.2)
+        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",
                 )

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

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

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"))
+    x_ = np.load(inp_comp.joinpath("x.npy"))
-    t_ = np.load(inp_comp.joinpath("tsec.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_eta_ = sgl.welch(eta_, f_, nperseg=nperseg)
-    phi_u_ = sgl.welch(u_, f, nperseg=nperseg)
+    phi_u_ = sgl.welch(u_, f_, nperseg=nperseg)
-    phi_eta_u_ = sgl.csd(eta_, 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

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

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

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)
+fig2, ax2 = plt.subplots(
-ax2.plot(t[cr0[i0 - 5] : cr0[i0 + 7]] * 1e-3, w0[cr0[i0 - 5] : cr0[i0 + 7]], color="k", label="Case 1")
+    figsize=(10 / 2.54, 2 / 3 * 10 / 2.54), constrained_layout=True
-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")
+)
+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():%}"
+)