m2cce/internship
m2cce/internship
1
Fork 0
Code Issues Activity

Server side scripts & all

Browse source
This commit is contained in:
Edgar P. Burkhart 2022-06-24 16:50:38 +02:00
parent a000c67e93
commit b92e52ecbb
Signed by: edpibu
GPG key ID: 9833D3C5A25BD227
20 changed files with 629 additions and 58 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

29
data/config2.ini Normal file
View file

@ -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

13
data/processing/orbitals.py
View file

@ -4,6 +4,7 @@ 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")
@ -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) 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], 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", 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 (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.pdf")
fig2dv.savefig("out_orbitals.jpg") fig2dv.savefig("out_orbitals.jpg")

8
data/processing/projection.py
View file

@ -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("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() fig, ax = plt.subplots(figsize=(16 / 2.54, 2 / 3 * 10 / 2.54), constrained_layout=True)
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"))

43
data/processing/ts.py
View file

@ -4,6 +4,7 @@ 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")
@ -25,12 +26,14 @@ out_ts = out_root.joinpath("ts.dat")
raw_ts = [] raw_ts = []
for tsi in config.get("inp", "raw_ts").split(","): for tsi in config.get("inp", "raw_ts").split(","):
raw_ts.append(np.loadtxt( raw_ts.append(
np.loadtxt(
inp_root.joinpath(tsi), inp_root.joinpath(tsi),
dtype=[("state", int), ("z", float), ("y", float), ("x", float)], dtype=[("state", int), ("z", float), ("y", float), ("x", float)],
delimiter=",", delimiter=",",
max_rows=2304, max_rows=2304,
)) )
)
n = len(raw_ts) n = len(raw_ts)
raw_ts = np.concatenate(raw_ts) raw_ts = np.concatenate(raw_ts)
log.debug(f"{raw_ts=}") log.debug(f"{raw_ts=}")
@ -45,7 +48,37 @@ 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() 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.pdf"))
fig.savefig(out_root.joinpath("ts.jpg"), dpi=200)

29
olaflow/of_ts_fine_1/constant/porosityDict Normal file
View file

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

29
olaflow/of_ts_fine_1/constant/turbulenceProperties Normal file
View file

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

81
olaflow/of_ts_fine_1/system/blockMeshDict Normal file
View file

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

138
olaflow/of_ts_fine_1/system/controlDict Normal file
View file

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

25
olaflow/of_ts_fine_1/system/graphUniform Normal file
View file

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

25
olaflow/of_ts_fine_1/system/graphUniform2 Normal file
View file

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

33
olaflow/processing/animate.py
View file

@ -8,6 +8,7 @@ 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
@ -28,6 +29,12 @@ 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)))
@ -61,9 +68,17 @@ 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(
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]) 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(
@ -86,10 +101,15 @@ 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)
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") 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()
@ -137,6 +157,13 @@ 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)

2
olaflow/processing/animateU.py
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, np.nanmax(UU)), clim=(0, 20),
) )
# (wat_p,) = axU.plot(x0, watl[0]) # (wat_p,) = axU.plot(x0, watl[0])

62
olaflow/processing/diff.py
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
import matplotlib.gridspec as gridspec from matplotlib.ticker import MultipleLocator
import numpy as np import numpy as np
from scipy import interpolate from scipy import interpolate
@ -67,11 +67,12 @@ 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=(6, 1.5 * len(models)), figsize=figsize,
dpi=100, dpi=200,
constrained_layout=True, constrained_layout=True,
squeeze=False, squeeze=False,
) )
@ -90,32 +91,71 @@ 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(
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=( 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)", 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( 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":
@ -123,7 +163,9 @@ else:
_ax.tricontour( _ax.tricontour(
_model.x, _model.x,
_model.z, _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,), levels=(0.5,),
colors="k", colors="k",
) )

2
swash/config-ts.ini
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=-1250 x0=-50
t0=180 t0=180

2
swash/config.ini
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

22
swash/processing/post.py
View file

@ -7,8 +7,6 @@ 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")
@ -30,7 +28,8 @@ 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")
@ -65,22 +64,25 @@ 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("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")) 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]))
@ -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)
/ (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

7
swash/processing/transmission.py
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, 5) t0 = np.linspace(23 * 60 + 8, 23 * 60 + 8 + 100, 6)
# 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(
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])) 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") 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,5 +67,6 @@ 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")

87
swash/processing/wavelet.py Normal file
View file

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

30
swash/processing/zero_cross.py
View file

@ -38,6 +38,9 @@ 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(
@ -72,9 +75,22 @@ 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(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.set(xlabel="t (s)", ylabel="z (m)")
ax2.autoscale(True, "x", True) ax2.autoscale(True, "x", True)
ax2.grid() ax2.grid()
@ -82,7 +98,11 @@ 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(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"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(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():%}"
)