1
Fork 0

Compare commits

..

No commits in common. "d4e9c3fd7840b9b627fbf9427d37a00145af8995" and "25e0f91bf0acce0ed5d16359bf50d6c864fc139a" have entirely different histories.

19 changed files with 49 additions and 618 deletions

View file

@ -1,29 +0,0 @@
[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

View file

@ -4,7 +4,6 @@ import logging
import pathlib import pathlib
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
import numpy as np import numpy as np
parser = argparse.ArgumentParser(description="Plot orbitals") parser = argparse.ArgumentParser(description="Plot orbitals")
@ -86,18 +85,16 @@ 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) x0 = ts_flt["x"] * np.cos(theta) + ts_flt["y"] * np.sin(theta)
#ax2dv.plot(x0, z0, c="#0066ff", lw=1) #ax2dv.plot(x0, z0, c="#0066ff", lw=1)
ax2dv.quiver( ax2dv.quiver(
x0[:-1] * 1e-2, x0[:-1],
z0[:-1] * 1e-2, z0[:-1],
np.diff(x0)[:] * 1e-2, np.diff(x0)[:],
np.diff(z0)[:] * 1e-2, np.diff(z0)[:],
color="k", color="k",
scale_units="xy", scale_units="xy",
scale=1, scale=1,
) )
ax2dv.grid(c="k", alpha=.2) ax2dv.grid(c="k", alpha=.2)
ax2dv.set(aspect="equal", xlabel="x (m)", ylabel="z (m)") ax2dv.set(aspect="equal", xlabel="x (cm)", ylabel="z (cm)")
ax2dv.set(ylim=(-10, 10))
ax2dv.yaxis.set_minor_locator(MultipleLocator(1))
fig2dv.savefig("out_orbitals.pdf") fig2dv.savefig("out_orbitals.pdf")
fig2dv.savefig("out_orbitals.jpg") fig2dv.savefig("out_orbitals.jpg")

View file

@ -131,12 +131,8 @@ 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("poro.dat"), poro[::-1], newline=" ")
np.savetxt(out_root.joinpath("psize.dat"), psize[::-1], newline=" ") np.savetxt(out_root.joinpath("psize.dat"), psize[::-1], newline=" ")
fig, ax = plt.subplots(figsize=(16 / 2.54, 2 / 3 * 10 / 2.54), constrained_layout=True) fig, ax = plt.subplots()
ax.plot(-x, z, color="k") ax.plot(-x, z, color="k")
ax.fill_between(-x, z + hstru, z, color="k", alpha=0.2) 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")) fig.savefig(out_root.joinpath("bathy.pdf"))

View file

@ -4,7 +4,6 @@ import logging
import pathlib import pathlib
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
import numpy as np import numpy as np
parser = argparse.ArgumentParser(description="Pre-process time-series") parser = argparse.ArgumentParser(description="Pre-process time-series")
@ -48,37 +47,7 @@ log.debug(f"{t=}")
log.info(f"Saving timeseries to '{out_ts}'") 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(figsize=(8 / 2.54, 2 / 3 * 10 / 2.54), constrained_layout=True) fig, ax = plt.subplots()
tp = np.datetime64("2017-02-28T17:00:00") + t.astype(np.timedelta64)[-(t.size // 3) :] ax.plot(t, raw_ts["z"])
ax.plot( ax.set(xlabel="t (s)", ylabel="z (cm)")
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.pdf"))
fig.savefig(out_root.joinpath("ts.jpg"), dpi=200)

View file

@ -1,29 +0,0 @@
/*--------------------------------*- 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);
// ************************************************************************* //

View file

@ -1,29 +0,0 @@
/*--------------------------------*- 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;
}
// ************************************************************************* //

View file

@ -1,81 +0,0 @@
/*--------------------------------*- 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
(
);
// ************************************************************************* //

View file

@ -1,138 +0,0 @@
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / 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;
}
);
}
}
*/
// ************************************************************************* //

View file

@ -1,25 +0,0 @@
/*--------------------------------*- 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"
// ************************************************************************* //

View file

@ -1,25 +0,0 @@
/*--------------------------------*- 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"
// ************************************************************************* //

View file

@ -8,7 +8,6 @@ import pickle
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.animation as animation import matplotlib.animation as animation
from matplotlib.gridspec import GridSpec from matplotlib.gridspec import GridSpec
from matplotlib.ticker import MultipleLocator
import numpy as np import numpy as np
from scipy import interpolate from scipy import interpolate
@ -29,12 +28,6 @@ parser.add_argument(
help="Only compute maximum rather than animation", help="Only compute maximum rather than animation",
action="store_true", action="store_true",
) )
parser.add_argument(
"-i",
"--initial",
help="Only compute initial domain",
action="store_true",
)
args = parser.parse_args() args = parser.parse_args()
logging.basicConfig(level=max((10, 20 - 10 * args.verbose))) logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
@ -68,19 +61,11 @@ AW[:, iz0[idz0], ix0[idx0]] = model.fields["alpha.water"]
U = np.full((model.t.size, *X.shape), np.nan) U = np.full((model.t.size, *X.shape), np.nan)
U[:, iz0[idz0], ix0[idx0]] = np.linalg.norm(model.fields["U"], axis=1) U[:, iz0[idz0], ix0[idx0]] = np.linalg.norm(model.fields["U"], axis=1)
fig = plt.figure( fig = plt.figure()
figsize=(15 / 2.54, 4 / 2.54), dpi=200, constrained_layout=True gs = GridSpec(3, 1, figure=fig, height_ratios=[1, 0.05, 0.05])
)
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]) ax = fig.add_subplot(gs[0])
if not args.initial: cax1 = fig.add_subplot(gs[1])
cax1 = fig.add_subplot(gs[1]) cax2 = fig.add_subplot(gs[2])
cax2 = fig.add_subplot(gs[2])
aw_m = ax.imshow( aw_m = ax.imshow(
AW[0], AW[0],
vmin=0, vmin=0,
@ -101,15 +86,10 @@ p_m = ax.imshow(
ax.axhline(4.5, ls="-.", lw=1, c="k", alpha=0.2, zorder=1.2) ax.axhline(4.5, ls="-.", lw=1, c="k", alpha=0.2, zorder=1.2)
if not args.initial: fig.colorbar(aw_m, label=r"$\alpha_w$", cax=cax1, shrink=0.6, orientation="horizontal")
fig.colorbar( fig.colorbar(p_m, label=r"Porosity", cax=cax2, shrink=0.6, orientation="horizontal")
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.set(xlabel="x (m)", ylabel="z (m)", aspect="equal", facecolor="#000000")
ax.grid(c="k", alpha=0.2) ax.grid(c="k", alpha=0.2)
ax.xaxis.set_minor_locator(MultipleLocator(5))
ax.yaxis.set_minor_locator(MultipleLocator(5))
figU = plt.figure() figU = plt.figure()
@ -157,13 +137,6 @@ if args.max:
fig.savefig(out.joinpath("max_aw.pdf")) fig.savefig(out.joinpath("max_aw.pdf"))
figU.savefig(out.joinpath("max_U.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: else:
fig.set(figwidth=19.2, figheight=10.8, dpi=100) fig.set(figwidth=19.2, figheight=10.8, dpi=100)
figU.set(figwidth=19.2, figheight=10.8, dpi=100) figU.set(figwidth=19.2, figheight=10.8, dpi=100)

View file

@ -74,7 +74,7 @@ u_m = axU.quiver(
UU[0], UU[0],
alpha=alp[0], alpha=alp[0],
cmap="inferno_r", cmap="inferno_r",
clim=(0, 20), clim=(0, np.nanmax(UU)),
) )
# (wat_p,) = axU.plot(x0, watl[0]) # (wat_p,) = axU.plot(x0, watl[0])

View file

@ -9,7 +9,7 @@ import sys
from cycler import cycler from cycler import cycler
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator import matplotlib.gridspec as gridspec
import numpy as np import numpy as np
from scipy import interpolate from scipy import interpolate
@ -67,12 +67,11 @@ def get_pickle(out):
models = list(map(get_pickle, args.output)) models = list(map(get_pickle, args.output))
figsize = 15 / 2.54, 4 / 2.54 * len(models)
fig, ax_ = plt.subplots( fig, ax_ = plt.subplots(
len(models), len(models),
figsize=figsize, figsize=(6, 1.5 * len(models)),
dpi=200, dpi=100,
constrained_layout=True, constrained_layout=True,
squeeze=False, squeeze=False,
) )
@ -91,71 +90,32 @@ if args.timestep is None:
) )
case "U": case "U":
for i, (_ax, _model) in enumerate(zip(ax, models)): 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( _c = _ax.imshow(
_data[::-1], 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],
cmap="PiYG", vmin=0,
alpha=np.clip(_model.post_fields[args.func]["alpha.water"], 0, 1).T[::-1], vmax=20,
cmap="inferno_r",
extent=( extent=(
_model.t.min(), _model.t.min(),
_model.t.max(), _model.t.max(),
_model.post_fields[args.func][f"x_{args.field}"].min(), _model.post_fields[args.func][f"x_{args.field}"].min(),
_model.post_fields[args.func][f"x_{args.field}"].max(), _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) 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 _: case _:
log.error(f"Cannot plot field {args.field} from {args.func}") log.error(f"Cannot plot field {args.field} from {args.func}")
sys.exit(1) sys.exit(1)
for i, (_ax, _model) in enumerate(zip(ax, models)): for i, (_ax, _model) in enumerate(zip(ax, models)):
_ax.set(xlabel="t (s)", ylabel="z (m)") _ax.set(xlabel="t (s)", ylabel="z (m)", title=f"Case {i}")
if len(models) > 1: _ax.grid(color="k", alpha=0.2)
_ax.set(title=f"Case {i}")
#_ax.grid(color="#797979", alpha=0.5)
fig.savefig( fig.savefig(
args.output[0].joinpath( args.output[0].joinpath(
f"diff_{args.func}_{args.field}_{'_'.join([o.name for o in args.output])}.pdf" 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: else:
match args.field: match args.field:
case "alpha.water": case "alpha.water":
@ -163,9 +123,7 @@ else:
_ax.tricontour( _ax.tricontour(
_model.x, _model.x,
_model.z, _model.z,
_model.fields[args.field][np.where(_model.t == args.timestep)[0]][ _model.fields[args.field][np.where(_model.t == args.timestep)[0]][0],
0
],
levels=(0.5,), levels=(0.5,),
colors="k", colors="k",
) )

View file

@ -10,5 +10,5 @@ mpi=8
[post] [post]
inp=inp_post/ts_4lay_1h inp=inp_post/ts_4lay_1h
out=out_post/ts_4lay_1h out=out_post/ts_4lay_1h
x0=-50 x0=-1250
t0=180 t0=180

View file

@ -9,7 +9,7 @@ mpi=8
[post] [post]
inp=inp_post/real_spec_interp 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 out=out_post/real_spec_interp
x0=-1250 x0=-1250
t0=180 t0=180

View file

@ -7,6 +7,8 @@ import matplotlib.pyplot as plt
import numpy as np import numpy as np
import scipy.signal as sgl import scipy.signal as sgl
from .read_swash import *
parser = argparse.ArgumentParser(description="Post-process swash output") parser = argparse.ArgumentParser(description="Post-process swash output")
parser.add_argument("-v", "--verbose", action="count", default=0) parser.add_argument("-v", "--verbose", action="count", default=0)
parser.add_argument("-c", "--config", default="config.ini") parser.add_argument("-c", "--config", default="config.ini")
@ -28,8 +30,7 @@ t = np.load(inp.joinpath("t.npy"))
botl = np.load(inp.joinpath("botl.npy")) botl = np.load(inp.joinpath("botl.npy"))
watl = np.load(inp.joinpath("watl.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 # Cospectral calculations
x0 = config.getint("post", "x0") x0 = config.getint("post", "x0")
@ -64,25 +65,22 @@ R = np.sqrt(
if config.has_option("post", "compare"): if config.has_option("post", "compare"):
inp_comp = pathlib.Path(config.get("post", "compare")) inp_comp = pathlib.Path(config.get("post", "compare"))
x_ = np.load(inp_comp.joinpath("x.npy")) x_ = np.load(inp_comp.joinpath("xp.npy"))
t_ = np.load(inp_comp.joinpath("t.npy")) t_ = np.load(inp_comp.joinpath("tsec.npy"))
botl_ = np.load(inp_comp.joinpath("botl.npy")) botl_ = np.load(inp_comp.joinpath("botl.npy"))
watl_ = np.load(inp_comp.joinpath("watl.npy")) watl_ = np.load(inp_comp.joinpath("watl.npy"))
vel_ = np.load(inp_comp.joinpath("vel_x.npy")) vel_ = np.load(inp_comp.joinpath("vel_x.npy"))
# Cospectral calculations
arg_x0_ = np.abs(x_ - x0).argmin() arg_x0_ = np.abs(x_ - x0).argmin()
arg_t0_ = np.abs(t_ - t0).argmin() arg_t0_ = np.abs(t_ - t0).argmin()
dt_ = np.diff(t_).mean() * 1e-3
f_ = 1 / dt_
eta_ = watl_[t_ > t0, arg_x0_] eta_ = watl_[t_ > t0, arg_x0_]
u_ = vel_[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])) H_ = np.sqrt(np.abs(phi_eta_[1]))
U_ = np.sqrt(np.abs(phi_u_[1])) U_ = np.sqrt(np.abs(phi_u_[1]))
@ -93,6 +91,10 @@ 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)
/ (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 # Plotting

View file

@ -30,7 +30,7 @@ botl = np.load(inp.joinpath("botl.npy"))
watl = np.load(inp.joinpath("watl.npy")) watl = np.load(inp.joinpath("watl.npy"))
vel = np.load(inp.joinpath("vel.npy"))[0] vel = np.load(inp.joinpath("vel.npy"))[0]
t0 = np.linspace(23 * 60 + 8, 23 * 60 + 8 + 100, 6) t0 = np.linspace(23 * 60 + 8, 23 * 60 + 8 + 100, 5)
# Plotting # Plotting
log.info("Plotting results") 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)) vlim = np.nanmin(np.maximum(watl, -botl)), np.nanmax(np.maximum(watl, -botl))
fig_x, ax = plt.subplots( fig_x, ax = plt.subplots(
3, 2, figsize=(15 / 2.54, 2 / 3 * 10 / 2.54), constrained_layout=True len(t0), figsize=(10 / 2.54, 4 / 3 * 10 / 2.54), constrained_layout=True
) )
i0 = np.argmin(np.abs(t * 1e-3 - t0[0])) i0 = np.argmin(np.abs(t * 1e-3 - t0[0]))
for ax_x, t0_x in zip(ax.reshape(-1), t0): for ax_x, t0_x in zip(ax, t0):
ax_x.plot(x, -botl, color="k") ax_x.plot(x, -botl, color="k")
i = np.argmin(np.abs(t * 1e-3 - t0_x)) i = np.argmin(np.abs(t * 1e-3 - t0_x))
ax_x.plot( ax_x.plot(
@ -67,6 +67,5 @@ log.info(f"Saving plots in '{out}'")
out.mkdir(parents=True, exist_ok=True) out.mkdir(parents=True, exist_ok=True)
fig_x.savefig(out.joinpath("x.pdf")) fig_x.savefig(out.joinpath("x.pdf"))
fig_x.savefig(out.joinpath("x.jpg"), dpi=200)
log.info("Finished post-processing") log.info("Finished post-processing")

View file

@ -1,87 +0,0 @@
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")

View file

@ -38,9 +38,6 @@ w0_comp = watl_comp[:, arg_x0]
cr0 = np.where(np.diff(np.sign(w0)))[0] cr0 = np.where(np.diff(np.sign(w0)))[0]
cr0_comp = np.where(np.diff(np.sign(w0_comp)))[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( wave = np.fromiter(
( (
np.abs( np.abs(
@ -75,22 +72,9 @@ ax.autoscale(True, "x", True)
ax.grid() ax.grid()
fig.savefig(out.joinpath("wsize.pdf")) fig.savefig(out.joinpath("wsize.pdf"))
fig2, ax2 = plt.subplots( fig2, ax2 = plt.subplots(figsize=(10/2.54, 2/3*10/2.54), constrained_layout=True)
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")
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.set(xlabel="t (s)", ylabel="z (m)")
ax2.autoscale(True, "x", True) ax2.autoscale(True, "x", True)
ax2.grid() ax2.grid()
@ -98,11 +82,7 @@ ax2.legend()
fig2.savefig(out.joinpath("maxw.pdf")) fig2.savefig(out.joinpath("maxw.pdf"))
fig2.savefig(out.joinpath("maxw.jpg"), dpi=200) fig2.savefig(out.joinpath("maxw.jpg"), dpi=200)
log.info( 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()):%}")
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"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: {wave.max()}m ; {wave_comp.max()}m")
log.info( log.info(f"Maximum wave size difference: {abs(wave_comp.max()-wave.max())/wave.max():%}")
f"Maximum wave size difference: {abs(wave_comp.max()-wave.max())/wave.max():%}"
)