Switched to new storage convention for swash output in sws_ola

This commit is contained in:
Edgar P. Burkhart 2022-03-28 09:44:35 +02:00
parent b557712372
commit 71049d49ea
Signed by: edpibu
GPG key ID: 9833D3C5A25BD227
63 changed files with 174 additions and 5717 deletions

View file

@ -1,60 +1,3 @@
import argparse
import configparser
import logging
import pathlib
import matplotlib.animation as animation
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
parser = argparse.ArgumentParser(description="Animate swash output")
parser.add_argument("-v", "--verbose", action="count", default=0)
args = parser.parse_args()
logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
log = logging.getLogger("post")
log.info("Starting post-processing")
config = configparser.ConfigParser()
config.read("config.ini")
inp = pathlib.Path(config.get("post", "inp"))
root = pathlib.Path(config.get("swash", "out"))
bathy = pd.read_hdf(
pathlib.Path(config.get("data", "out")).joinpath("bathy.h5"), "bathy"
)
def data(var):
return np.load(inp.joinpath(f"{var}.npy"))
x = data("xp")
t = data("tsec")
watl = data("watl")
botl = data("botl")
wl = np.maximum(watl, -botl)
# print(x.size, -np.arange(0, 1 * bathy.hstru.size, 1)[::-1].size)
fig, ax = plt.subplots()
ax.plot(x, -botl, c="k")
# ax.fill_between(
# x, -botl, -data["botl"] + bathy.hstru, color="k", alpha=0.2
# )
(line,) = ax.plot(x, wl[0])
def animate(i):
line.set_ydata(wl[i])
return (line,)
ani = animation.FuncAnimation(
fig, animate, frames=wl[:, 0].size, interval=20, blit=True
)
plt.show(block=True)
version https://git-lfs.github.com/spec/v1
oid sha256:4ab57e0a813e6ea3086045717e251b2c44e00891a858d9f6365e3e604fefee03
size 1316

View file

@ -1,93 +1,3 @@
import argparse
import configparser
import logging
import pathlib
import sys
import numpy as np
import pandas as pd
try:
import matplotlib.pyplot as plt
except ImportError:
plt = None
parser = argparse.ArgumentParser(description="Pre-process bathymetry")
parser.add_argument("-v", "--verbose", action="count", default=0)
args = parser.parse_args()
logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
log = logging.getLogger("bathy")
log.info("Starting bathymetry pre-processing")
config = configparser.ConfigParser()
config.read("config.ini")
root = pathlib.Path(config.get("data", "root"))
log.info(f"Reading input data from '{root}'")
bathy_hires = np.loadtxt(root.joinpath(config.get("data", "hires")))
bathy_lores = np.loadtxt(root.joinpath(config.get("data", "bathy")))
hstru = np.loadtxt(root.joinpath(config.get("data", "hstru")))
poro = np.loadtxt(root.joinpath(config.get("data", "poro")))
psize = np.loadtxt(root.joinpath(config.get("data", "psize")))
log.info("Generating grid")
x_hires = -np.arange(0, 0.5 * bathy_hires.size, 0.5)[::-1]
x_lores = -np.arange(0, 1 * bathy_lores.size, 1)[::-1]
x_hstru = -np.arange(0, 0.5 * hstru.size, 0.5)[::-1]
log.info("Generating output data")
bathy_hires_pd = pd.Series(bathy_hires.copy(), index=x_hires)
bathy_lores_pd = pd.Series(bathy_lores.copy(), index=x_lores)
bathy = pd.DataFrame(
index=bathy_lores_pd.index.union(bathy_hires_pd.index),
columns=("z", "hstru", "poro", "psize"),
data=0,
)
bathy.z[bathy_lores_pd.index] = bathy_lores_pd
bathy.z[bathy_hires_pd.index] = bathy_hires_pd
bathy.loc[x_hstru, ("hstru", "poro", "psize")] = np.array(
(hstru, poro, psize)
).T
bathy = bathy.reindex(bathy_lores_pd.index)
log.debug(f"Bathymetry:\n{bathy}")
log.info(
f"xmin: {bathy.index.min()}, "
f"xmax: {bathy.index.max()}, "
f"n: {bathy.index.size}"
)
if config.has_option("data", "out"):
out = pathlib.Path(config.get("data", "out"))
log.info(f"Writing output data to '{out}'")
out.mkdir(exist_ok=True)
np.savetxt(out.joinpath("bathy.dat"), bathy.z, newline=" ")
np.savetxt(out.joinpath("hstru.dat"), bathy.hstru, newline=" ")
np.savetxt(out.joinpath("poro.dat"), bathy.poro, newline=" ")
np.savetxt(out.joinpath("psize.dat"), bathy.psize, newline=" ")
bathy.to_hdf(out.joinpath("bathy.h5"), "bathy", mode="w")
if config.getboolean("proc", "plot", fallback=False):
if plt is None:
log.error("Could not import PyPlot")
sys.exit(1)
log.info("Plotting data")
fig, ax = plt.subplots()
ax.plot(x_hires, bathy_hires, label="High-res")
ax.plot(x_lores, bathy_lores, label="Low-res")
ax.plot(bathy.index, bathy.z, ls="-.", c="k", label="Combined")
ax.plot(bathy.index, bathy.z + bathy.hstru, label="Hstru")
ax.grid()
ax.legend()
plt.show(block=True)
log.info("Processing finished")
version https://git-lfs.github.com/spec/v1
oid sha256:93e058c82e05d4d81cb1bb68459b978780f6c32fd0ba6c5d3662212542a0d118
size 2877

View file

@ -1,94 +1,3 @@
import argparse
import configparser
import logging
import pathlib
import sys
import numpy as np
import pandas as pd
try:
import matplotlib.pyplot as plt
except ImportError:
plt = None
parser = argparse.ArgumentParser(description="Pre-process bathymetry")
parser.add_argument("-v", "--verbose", action="count", default=0)
args = parser.parse_args()
logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
log = logging.getLogger("bathy")
log.info("Starting bathymetry pre-processing")
config = configparser.ConfigParser()
config.read("config.ini")
root = pathlib.Path(config.get("data", "root"))
log.info(f"Reading input data from '{root}'")
bathy_hires = np.loadtxt(root.joinpath(config.get("data", "hires")))
bathy_lores = np.loadtxt(root.joinpath(config.get("data", "bathy")))
hstru = np.loadtxt(root.joinpath(config.get("data", "hstru")))
poro = np.loadtxt(root.joinpath(config.get("data", "poro")))
psize = np.loadtxt(root.joinpath(config.get("data", "psize")))
log.info("Generating grid")
x_hires = -np.arange(0, 0.5 * bathy_hires.size, 0.5)[::-1]
x_lores = -np.arange(0, 1 * bathy_lores.size, 1)[::-1]
x_hstru = -np.arange(0, 0.5 * hstru.size, 0.5)[::-1]
log.info("Generating output data")
bathy_hires_pd = pd.Series(bathy_hires.copy(), index=x_hires)
bathy_lores_pd = pd.Series(bathy_lores.copy(), index=x_lores)
bathy = pd.DataFrame(
index=bathy_lores_pd.index.union(bathy_hires_pd.index),
columns=("z", "hstru", "poro", "psize"),
data=0,
)
bathy.z[bathy_lores_pd.index] = bathy_lores_pd
bathy.z[bathy_hires_pd.index] = bathy_hires_pd
bathy.z = np.minimum(bathy.z, -15)
# bathy.loc[x_hstru, ("hstru", "poro", "psize")] = np.array(
# (hstru, poro, psize)
# ).T
bathy = bathy.reindex(bathy_lores_pd.index)
log.debug(f"Bathymetry:\n{bathy}")
log.info(
f"xmin: {bathy.index.min()}, "
f"xmax: {bathy.index.max()}, "
f"n: {bathy.index.size}"
)
if config.has_option("data", "out_nb"):
out = pathlib.Path(config.get("data", "out_nb"))
log.info(f"Writing output data to '{out}'")
out.mkdir(exist_ok=True)
np.savetxt(out.joinpath("bathy.dat"), bathy.z, newline=" ")
np.savetxt(out.joinpath("hstru.dat"), bathy.hstru, newline=" ")
np.savetxt(out.joinpath("poro.dat"), bathy.poro, newline=" ")
np.savetxt(out.joinpath("psize.dat"), bathy.psize, newline=" ")
bathy.to_hdf(out.joinpath("bathy.h5"), "bathy", mode="w")
if config.getboolean("proc", "plot", fallback=False):
if plt is None:
log.error("Could not import PyPlot")
sys.exit(1)
log.info("Plotting data")
fig, ax = plt.subplots()
ax.plot(x_hires, bathy_hires, label="High-res")
ax.plot(x_lores, bathy_lores, label="Low-res")
ax.plot(bathy.index, bathy.z, ls="-.", c="k", label="Combined")
ax.plot(bathy.index, bathy.z + bathy.hstru, label="Hstru")
ax.grid()
ax.legend()
plt.show(block=True)
log.info("Processing finished")
version https://git-lfs.github.com/spec/v1
oid sha256:8e3a75a559d7bf52e3e1703c5e45c8e473cd813aa99a95169bd9a500285f5a2d
size 2923

View file

@ -1,92 +1,3 @@
import argparse
import configparser
import logging
import pathlib
import matplotlib.animation as animation
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
parser = argparse.ArgumentParser(description="Animate swash output")
parser.add_argument("-v", "--verbose", action="count", default=0)
args = parser.parse_args()
logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
log = logging.getLogger("post")
log.info("Starting post-processing")
config = configparser.ConfigParser()
config.read("config.ini")
inp = pathlib.Path(config.get("post", "inp"))
root = pathlib.Path(config.get("swash", "out"))
bathy = pd.read_hdf(
pathlib.Path(config.get("data", "out")).joinpath("bathy.h5"), "bathy"
)
def data(var):
return np.load(inp.joinpath(f"{var}.npy"))
x = data("xp")
t = data("tsec")
watl = data("watl")
botl = data("botl")
zk = data("zk")
velk = data("velk")
vz = data("vz")
wl = np.maximum(watl, -botl)
# print(x.size, -np.arange(0, 1 * bathy.hstru.size, 1)[::-1].size)
fig, ax = plt.subplots()
# ax.plot(x, -botl, c="k")
# ax.fill_between(
# x, -botl, -data["botl"] + bathy.hstru, color="k", alpha=0.2
# )
n = 0
vk = np.sqrt((velk[n] ** 2).sum(axis=1))
# print(vk.shape)
# plt.imshow(vk)
# plt.colorbar()
lines = ax.plot(x, zk[n].T, c="#0066cc")
quiv = []
for i in range(10):
quiv.append(
ax.quiver(
x[::50],
(zk[n, i, ::50] + zk[n, i + 1, ::50]) / 2,
velk[n, i, 0, ::50],
vz[n, i, ::50],
units="dots",
width=2,
scale=0.05,
)
)
ax.autoscale(True, "w", True)
ax.set_ylim(top=15)
def animate(k):
for i, q in enumerate(quiv):
q.set_UVC(
velk[k, i, 0, ::50],
vz[k, i, ::50],
)
for i, l in enumerate(lines):
l.set_ydata(zk[k, i])
return *quiv, *lines
ani = animation.FuncAnimation(
fig, animate, frames=wl[:, 0].size, interval=20, blit=True
)
plt.show(block=True)
version https://git-lfs.github.com/spec/v1
oid sha256:60aeefb03fe5f6ec669d64cd781e7254d7aeace3a2faab75098a4734343385a0
size 1992

View file

@ -1,169 +1,3 @@
import argparse
import configparser
import logging
import pathlib
import matplotlib.pyplot as plt
import numpy as np
import scipy.fft as fft
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)
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("config.ini")
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("xp.npy"))
t = np.load(inp.joinpath("tsec.npy"))
botl = np.load(inp.joinpath("botl.npy"))
watl = np.load(inp.joinpath("watl.npy"))
vel = np.load(inp.joinpath("vel.npy"))
# Cospectral calculations
x0 = config.getint("post", "x0")
arg_x0 = np.abs(x - x0).argmin()
t0 = config.getfloat("post", "t0")
arg_t0 = np.abs(t - t0).argmin()
dt = config.getfloat("post", "dt")
f = 1 / dt
nperseg = config.getint("post", "nperseg", fallback=None)
log.info(f"Computing reflection coefficient at x={x0}")
eta = watl[t > t0, arg_x0]
u = vel[t > t0, 0, arg_x0]
phi_eta = sgl.welch(eta, f, nperseg=nperseg)
phi_u = sgl.welch(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]))
G = H / U
th_eta_u = np.angle(phi_eta_u[1])
# R1 = np.sqrt(
# (np.abs(phi_eta[1]) + np.abs(phi_u[1]) - 2 * np.abs(phi_eta_u[1]))
# / (np.abs(phi_eta[1]) + np.abs(phi_u[1]) + 2 * np.abs(phi_eta_u[1]))
# )
R = np.sqrt(
(1 + G**2 - 2 * G * np.cos(th_eta_u))
/ (1 + G**2 + 2 * G * np.cos(th_eta_u))
)
if config.has_option("post", "compare"):
inp_comp = pathlib.Path(config.get("post", "compare"))
x_ = np.load(inp_comp.joinpath("xp.npy"))
t_ = np.load(inp_comp.joinpath("tsec.npy"))
botl_ = np.load(inp_comp.joinpath("botl.npy"))
watl_ = np.load(inp_comp.joinpath("watl.npy"))
vel_ = np.load(inp_comp.joinpath("vel.npy"))
arg_x0_ = np.abs(x_ - x0).argmin()
arg_t0_ = np.abs(t_ - t0).argmin()
eta_ = watl_[t_ > t0, arg_x0_]
u_ = vel_[t_ > t0, 0, arg_x0_]
phi_eta_ = sgl.welch(eta_, f, nperseg=nperseg)
phi_u_ = sgl.welch(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]))
G_ = H_ / U_
th_eta_u_ = np.angle(phi_eta_u_[1])
R_ = np.sqrt(
(1 + G_**2 - 2 * G_ * np.cos(th_eta_u_))
/ (1 + G_**2 + 2 * G_ * np.cos(th_eta_u_))
)
# Plotting
log.info("Plotting results")
fig, (ax_watl, ax_vel) = plt.subplots(2)
ax_watl.plot(t, watl[:, arg_x0], lw=1, label="watl")
ax_watl.set(xlabel="t (s)", ylabel="z (m)")
ax_watl.autoscale(axis="x", tight=True)
ax_watl.grid()
ax_watl.axvline(t0, c="k", alpha=0.2)
ax_vel.plot(t, vel[:, 0, arg_x0], lw=1, label="vel")
ax_vel.set(xlabel="t (s)", ylabel="U (m/s)")
ax_vel.autoscale(axis="x", tight=True)
ax_vel.grid()
ax_vel.axvline(t0, c="k", alpha=0.2)
fig.tight_layout()
fig_r, ax_r = plt.subplots()
ax_fft = ax_r.twinx()
ax_fft.plot(
*sgl.welch(eta, 1 / dt, nperseg=nperseg),
lw=1,
c="k",
alpha=0.2,
label="PSD ($\\eta$, cas 1)",
)
ax_r.plot(phi_eta[0], R, marker="+", label="R (cas 1)")
if config.has_option("post", "compare"):
ax_fft.plot(
*sgl.welch(eta_, 1 / dt, nperseg=nperseg),
lw=1,
c="k",
alpha=0.2,
label="PSD ($\\eta$, cas 2)",
)
ax_r.plot(phi_eta[0], R_, marker="+", label="R (cas 2)")
ax_r.set(xlim=(0, 0.3), ylim=(0, 1), xlabel="f (Hz)", ylabel="R")
ax_fft.set(ylim=0, ylabel="PSD (m²/Hz)")
ax_r.grid()
ax_r.legend(loc="upper left")
ax_fft.legend(loc="upper right")
fig_r.tight_layout()
fig_x, ax_x = plt.subplots(figsize=(10, 1))
ax_x.plot(x, -botl, color="k")
ax_x.plot(
x,
np.maximum(watl[arg_t0, :], -botl),
)
if config.has_option("post", "compare"):
ax_x.plot(x, -botl_, color="k", ls="-.")
ax_x.plot(
x,
np.maximum(watl_[arg_t0, :], -botl_),
ls="-.",
)
ax_x.axvline(x0, c="k", alpha=0.2)
ax_x.set(xlabel="x (m)", ylabel="z (m)")
ax_x.autoscale(axis="x", tight=True)
ax_x.set(aspect="equal")
fig_x.tight_layout()
out = pathlib.Path(config.get("post", "out")).joinpath(f"t{t0}x{x0}")
log.info(f"Saving plots in '{out}'")
out.mkdir(parents=True, exist_ok=True)
fig.savefig(out.joinpath("t.png"))
fig_r.savefig(out.joinpath("R.png"))
fig_x.savefig(out.joinpath("x.png"))
plt.show(block=True)
log.info("Finished post-processing")
version https://git-lfs.github.com/spec/v1
oid sha256:c8fefac1023d769a2a43c86f159cf33155f5462aed2dcd4ac968d782f79760d4
size 4665

View file

@ -1,52 +1,3 @@
import subprocess
import tempfile
import numpy as np
class ReadSwash:
def __init__(self):
self._n_x = None
self._n_t = None
self._t = None
self._x = None
@classmethod
def read_nohead(cls, path):
subprocess.run(("sed", "-i", r"s/ /\n/g", path))
return np.loadtxt(path)
def read_time(self, path):
self._t = np.unique(self.read_nohead(path))
self._n_t = self._t.size
return self.t
def read_x(self, path):
self._x = np.unique(self.read_nohead(path))
self._n_x = self._x.size
return self.x
def read_scalar(self, path, const=False):
if const:
return self.read_nohead(path).reshape((self._n_t, self._n_x))[0, :]
return self.read_nohead(path).reshape((self._n_t, self._n_x))
def read_const(self, path):
return self.read_scalar(path, const=True)
def read_vector(self, path):
return self.read_nohead(path).reshape((self._n_t, 2, self._n_x))
def read_scalar_lay(self, path):
return self.read_nohead(path).reshape((self._n_t, -1, self._n_x))
def read_vector_lay(self, path):
return self.read_nohead(path).reshape((self._n_t, -1, 2, self._n_x))
@property
def t(self):
return self._t
@property
def x(self):
return self._x
version https://git-lfs.github.com/spec/v1
oid sha256:81ddc1baff1a727e68dc2819fe3f39ca657c488ab8cbeb83c2d260ae7c9c13fd
size 1349

View file

@ -1,67 +1,3 @@
import argparse
import configparser
import logging
import os
import pathlib
import shutil
import subprocess
import sys
import tempfile
parser = argparse.ArgumentParser(description="Run swash model")
parser.add_argument("-v", "--verbose", action="count", default=0)
args = parser.parse_args()
logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
log = logging.getLogger("swash")
log.info("Starting swash model")
config = configparser.ConfigParser()
config.read("config.ini")
inp = pathlib.Path(config.get("swash", "input"))
out = pathlib.Path(config.get("swash", "out"))
if out.exists():
log.error(f"Swash output '{out}' already exists")
sys.exit(1)
with tempfile.TemporaryDirectory(prefix="swash_", dir=".") as tmp_raw:
tmpdir = pathlib.Path(tmp_raw)
log.info(f"Copying files to '{tmpdir}'")
shutil.copy2(inp, tmpdir)
if config.getboolean("swash", "nb", fallback=False):
path = pathlib.Path(config.get("data", "out_nb"))
else:
path = pathlib.Path(config.get("data", "out"))
shutil.copytree(path, tmpdir, dirs_exist_ok=True)
if config.has_option("swash", "mpi"):
mpi = ("-mpi", config.get("swash", "mpi"))
log.info(f"Using mpi with {mpi}")
else:
mpi = ()
with open(tmpdir.joinpath("sws.log"), "w") as logfile:
log.info(f"Runing swash in '{tmpdir}'")
path = pathlib.Path(config.get("swash", "path"))
cmd = (path.joinpath("swashrun"), "-input", inp.name, *mpi)
log.info(f"Running {cmd}")
swash_run = subprocess.Popen(
cmd,
cwd=tmpdir,
stdout=logfile,
stderr=logfile,
env={"PATH": f"{path}:{os.environ['PATH']}"},
)
code = swash_run.wait()
if code != 0:
log.error(f"Swash returned error code {code}")
log.info(f"Moving swash output to '{out}'")
shutil.move(tmpdir, out)
log.info(f"Swash model finished successfully")
version https://git-lfs.github.com/spec/v1
oid sha256:b17bd29f45896e07ef9e48ab27c905158de21bdedc7b9ecbf1a7943e4af7ae70
size 1962

View file

@ -1,51 +1,3 @@
import argparse
import configparser
import logging
import pathlib
from multiprocessing.pool import ThreadPool
import numpy as np
import pandas as pd
from .read_swash import ReadSwash
parser = argparse.ArgumentParser(description="Convert swash output to numpy")
parser.add_argument("-v", "--verbose", action="count", default=0)
args = parser.parse_args()
logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
log = logging.getLogger("sws_npz")
log.info("Starting sws -> npz converter")
config = configparser.ConfigParser()
config.read("config.ini")
sws_out = pathlib.Path(config.get("swash", "out"))
inp = pathlib.Path(config.get("post", "inp"))
log.info(f"Reading swash output from '{sws_out}'")
rsws = ReadSwash()
np.save(inp.joinpath("tsec"), rsws.read_time(sws_out.joinpath("tsec.dat")))
np.save(inp.joinpath("xp"), rsws.read_x(sws_out.joinpath("xp.dat")))
var = {
"dep": rsws.read_scalar,
"botl": rsws.read_const,
"watl": rsws.read_scalar,
"vel": rsws.read_vector,
"press": rsws.read_scalar,
"zk": rsws.read_scalar_lay,
"velk": rsws.read_vector_lay,
"vz": rsws.read_scalar_lay,
}
inp.mkdir(exist_ok=True)
with ThreadPool() as pool:
log.info("Converting all data")
pool.map(
lambda x: np.save(
inp.joinpath(x[0]),
x[1](sws_out.joinpath(x[0]).with_suffix(".dat")),
),
var.items(),
)
version https://git-lfs.github.com/spec/v1
oid sha256:f3b186f594089b65b3739d2f03855c691f45e0b2744e37196e07cdf60c896cfc
size 1393