1
Fork 0

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

This commit is contained in:
Edgar P. Burkhart 2022-03-30 14:18:00 +02:00
commit 18b6bea844
Signed by: edpibu
GPG key ID: 9833D3C5A25BD227
23 changed files with 111 additions and 311438 deletions

View file

@ -22,12 +22,61 @@ options:
## Bathymetry ## Bathymetry
Configuration:
```
[inp]
root: input directory
base: bathymetry database (.xyz) -- not provided in this repository
hires: 1d hires smoothed breakwater profile (.dat)
hstru: 1d hires smoothed breakwater porosity height (.dat)
poro: 1d breakwater porosity (.dat)
psize: 1d porosity size (.dat)
hires_step: mesh sizing of 1d profile
[out]
margin: margin around the buoy and breakwater for subdomain
no_breakwater: option to remove breakwater
root: output directory
sub: output file for subdomain
out: output file for 1d projection
step: mesh size
left: margin for 1d projection (breakwater side)
right: margin for 1d projection (buoy side)
[artha]
lat: latitude of the breakwater
lon: longitude of the breakwater
[buoy]
lat: latitude of the buoy
lon: longitude of the buoy
```
* Insert database in `data/data` * Insert database in `data/data`
* Run `processing.subdomain` to generate 2D data in a smaller domain * Run `processing.subdomain` to generate 2D data in a smaller domain
* Run `processing.projection` to generate a 1D bathymetry * Run `processing.projection` to generate a 1D bathymetry
## Swash ## Swash
Configuration:
```
[data]
out: input directory (output of bathymetry)
[swash]
input: swash input file (.sws)
path: path to swash directory (should contain swashrun and swash.exe)
out: output directory for swash (created at the end of the computation)
mpi: number of mpi threads (omit to use serial)
[post]
inp: input directory for post-processing (output of sws_npz)
compare: input directory for post-processing comparison (omit to display results from single case)
out: output directory for figures
x0: position of reflection calculations
t0: minimum time for reflection calculations (necessary because of boundary condition ramping)
```
* Run `processing.swash` to run the Swash model * Run `processing.swash` to run the Swash model
* Run `processing.sws_npz` to convert Swash output to Numpy files * Run `processing.sws_npz` to convert Swash output to Numpy files
* Run `processing.post` to plot wave height at a point, water level at an * Run `processing.post` to plot wave height at a point, water level at an
@ -37,6 +86,21 @@ options:
## OlaFlow ## OlaFlow
Configuration:
```
[swash]
np_out: input of olaflow from swash (output of sws_npz)
[bathy]
bathy: bathymetry to be used by olaflow
hstru: height of porous domain (see bathy)
scale: openscad scaling (mesh sizing in [x,y,z])
out: output folder for stl bathymetry
[olaflow]
root: olaflow computation directory
```
* Run `processing.bathy` to generate stl input files * Run `processing.bathy` to generate stl input files
* Run `blockMesh`, `snappyHexMesh -overwrite`, copy `0.org` to `0`, `setFields` * Run `blockMesh`, `snappyHexMesh -overwrite`, copy `0.org` to `0`, `setFields`
* Run `processing.sws_ola` to convert Swash output to OlaFlow input * Run `processing.sws_ola` to convert Swash output to OlaFlow input

View file

@ -8,12 +8,13 @@ psize=Psize.dat
hires_step=0.5 hires_step=0.5
[out] [out]
margin=0.005
no_breakwater=True no_breakwater=True
root=out_nb root=out_nb
sub=bathy_sub.npy sub=bathy_sub.npy
out=bathy.npy out=bathy.npy
step=1 step=1
left=0 left=-300
right=150 right=150
[artha] [artha]

View file

@ -14,7 +14,7 @@ root=out
sub=bathy_sub.npy sub=bathy_sub.npy
out=bathy.npy out=bathy.npy
step=1 step=1
left=0 left=-300
right=150 right=150
[artha] [artha]

View file

@ -91,13 +91,14 @@ x_min_hires = x[z_crossing] + (bathy_hires[:, 0].min() - bathy_hires[hires_cross
x_max_hires = x[z_crossing] + (bathy_hires[:, 0].max() - bathy_hires[hires_crossing, 0]) x_max_hires = x[z_crossing] + (bathy_hires[:, 0].max() - bathy_hires[hires_crossing, 0])
log.debug(f"Replacing range: [{x_min_hires},{x_max_hires}]") log.debug(f"Replacing range: [{x_min_hires},{x_max_hires}]")
flt_x = (x > x_min_hires) & (x < x_max_hires)
hstru = np.zeros(z.shape) hstru = np.zeros(z.shape)
poro = np.zeros(z.shape) poro = np.zeros(z.shape)
psize = np.zeros(z.shape) psize = np.zeros(z.shape)
if config.getboolean("out", "no_breakwater", fallback=False): if config.getboolean("out", "no_breakwater", fallback=False):
z[flt_x] = z[flt_x][-1] flt_x = np.abs(x) < 250
z[flt_x] = np.linspace(z[flt_x][0], z[flt_x][-1], flt_x.sum())
else: else:
flt_x = (x > x_min_hires) & (x < x_max_hires)
z[flt_x] = interpolate.griddata( z[flt_x] = interpolate.griddata(
(bathy_hires[:, 0],), (bathy_hires[:, 0],),
bathy_hires[:, 1], bathy_hires[:, 1],
@ -133,4 +134,5 @@ np.savetxt(out_root.joinpath("psize.dat"), psize[::-1], newline=" ")
fig, ax = plt.subplots() 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()}]")
fig.savefig(out_root.joinpath("bathy.pdf")) fig.savefig(out_root.joinpath("bathy.pdf"))

View file

@ -1,53 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class polyBoundaryMesh;
location "constant/polyMesh";
object boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
5
(
inlet
{
type patch;
nFaces 77;
startFace 38318;
}
outlet
{
type patch;
nFaces 34;
startFace 38395;
}
atmosphere
{
type patch;
nFaces 300;
startFace 38429;
}
defaultFaces
{
type empty;
inGroups List<word> 1(empty);
nFaces 38698;
startFace 38729;
}
bathy
{
type wall;
inGroups List<word> 1(wall);
nFaces 349;
startFace 77427;
}
)
// ************************************************************************* //

View file

@ -1,19 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class labelList;
location "constant/polyMesh";
object cellLevel;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
19349{0}
// ************************************************************************* //

View file

@ -1,20 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class regIOobject;
location "constant/polyMesh";
object cellZones;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
0
()
// ************************************************************************* //

View file

@ -1,20 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class regIOobject;
location "constant/polyMesh";
object faceZones;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
0
()
// ************************************************************************* //

File diff suppressed because it is too large Load diff

View file

@ -1,22 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class uniformDimensionedScalarField;
location "constant/polyMesh";
object level0Edge;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 0 0 0 0 0];
value 0.5;
// ************************************************************************* //

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

View file

@ -1,19 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class labelList;
location "constant/polyMesh";
object pointLevel;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39460{0}
// ************************************************************************* //

View file

@ -1,20 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class regIOobject;
location "constant/polyMesh";
object pointZones;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
0
()
// ************************************************************************* //

File diff suppressed because it is too large Load diff

File diff suppressed because it is too large Load diff

View file

@ -2,16 +2,14 @@
out=../data/out_nb out=../data/out_nb
[swash] [swash]
input=sws/SPEC_buoy_nb.sws input=sws/SPEC_buoy.sws
path=/data/code/swash path=/data/code/swash
out=out/spec_nb out=out/spec_nb_2lay
mpi=8 mpi=8
[post] [post]
inp=inp_post/spec_nb inp=inp_post/spec_nb_2lay
compare=inp_post_nb #compare=
out=out_post out=out_post/spec_nb_2lay
#nperseg=1024
dt=0.25
x0=-1250 x0=-1250
t0=180 t0=180

View file

@ -4,14 +4,12 @@ out=../data/out
[swash] [swash]
input=sws/SPEC_buoy.sws input=sws/SPEC_buoy.sws
path=/data/code/swash path=/data/code/swash
out=out/spec out=out/spec_2lay
mpi=8 mpi=8
[post] [post]
inp=inp_post/spec inp=inp_post/spec_2lay
compare=inp_post_nb compare=inp_post/spec_nb_2lay
out=out_post out=out_post/spec_2lay_bvsnb
#nperseg=1024
dt=0.25
x0=-1250 x0=-1250
t0=180 t0=180

View file

@ -55,7 +55,7 @@ vk = np.sqrt((velk[n] ** 2).sum(axis=1))
lines = ax.plot(x, zk[n].T, c="#0066cc") lines = ax.plot(x, zk[n].T, c="#0066cc")
quiv = [] quiv = []
for i in range(10): for i in range(len(lines)-1):
quiv.append( quiv.append(
ax.quiver( ax.quiver(
x[::50], x[::50],
@ -87,4 +87,4 @@ ani = animation.FuncAnimation(
fig, animate, frames=wl[:, 0].size, interval=20, blit=True fig, animate, frames=wl[:, 0].size, interval=20, blit=True
) )
ani.save(out.joinpath("layers.mp4", codec="h265")) ani.save(out.joinpath("layers.mp4"))

View file

@ -38,7 +38,7 @@ x0 = config.getint("post", "x0")
arg_x0 = np.abs(x - x0).argmin() arg_x0 = np.abs(x - x0).argmin()
t0 = config.getfloat("post", "t0") t0 = config.getfloat("post", "t0")
arg_t0 = np.abs(t - t0).argmin() arg_t0 = np.abs(t - t0).argmin()
dt = config.getfloat("post", "dt") dt = np.diff(t).mean()
f = 1 / dt f = 1 / dt
nperseg = config.getint("post", "nperseg", fallback=None) nperseg = config.getint("post", "nperseg", fallback=None)
log.info(f"Computing reflection coefficient at x={x0}") log.info(f"Computing reflection coefficient at x={x0}")
@ -59,10 +59,10 @@ R = 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]))
/ (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( #R = np.sqrt(
# (1 + G**2 - 2 * G * np.cos(th_eta_u)) # (1 + G**2 - 2 * G * np.cos(th_eta_u))
# / (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"): if config.has_option("post", "compare"):
inp_comp = pathlib.Path(config.get("post", "compare")) inp_comp = pathlib.Path(config.get("post", "compare"))
@ -92,10 +92,10 @@ if config.has_option("post", "compare"):
(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]))
/ (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( #R_ = np.sqrt(
# (1 + G_**2 - 2 * G_ * np.cos(th_eta_u_)) # (1 + G_**2 - 2 * G_ * np.cos(th_eta_u_))
# / (1 + G_**2 + 2 * G_ * np.cos(th_eta_u_)) # / (1 + G_**2 + 2 * G_ * np.cos(th_eta_u_))
# ) #)
# Plotting # Plotting
@ -166,8 +166,8 @@ out = pathlib.Path(config.get("post", "out")).joinpath(f"t{t0}x{x0}")
log.info(f"Saving plots in '{out}'") log.info(f"Saving plots in '{out}'")
out.mkdir(parents=True, exist_ok=True) out.mkdir(parents=True, exist_ok=True)
fig.savefig(out.joinpath("t.png")) fig.savefig(out.joinpath("t.pdf"))
fig_r.savefig(out.joinpath("R.png")) fig_r.savefig(out.joinpath("R.pdf"))
fig_x.savefig(out.joinpath("x.png")) fig_x.savefig(out.joinpath("x.pdf"))
log.info("Finished post-processing") log.info("Finished post-processing")

View file

@ -1,9 +1,12 @@
import logging
import subprocess import subprocess
import tempfile import tempfile
import numpy as np import numpy as np
log = logging.getLogger("read_swash")
class ReadSwash: class ReadSwash:
def __init__(self): def __init__(self):
self._n_x = None self._n_x = None
@ -13,7 +16,9 @@ class ReadSwash:
@classmethod @classmethod
def read_nohead(cls, path): def read_nohead(cls, path):
log.info(f"Replacing \\s with \\n in '{path}'")
subprocess.run(("sed", "-i", r"s/\s\+/\n/g", path)) subprocess.run(("sed", "-i", r"s/\s\+/\n/g", path))
log.info(f"Loading '{path}'")
return np.loadtxt(path) return np.loadtxt(path)
def read_time(self, path): def read_time(self, path):

View file

@ -30,7 +30,7 @@ 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"))) np.save(inp.joinpath("xp"), rsws.read_x(sws_out.joinpath("xp.dat")))
var = { var = {
"dep": rsws.read_scalar, #"dep": rsws.read_scalar,
"botl": rsws.read_const, "botl": rsws.read_const,
"watl": rsws.read_scalar, "watl": rsws.read_scalar,
"pressk": rsws.read_scalar_lay, "pressk": rsws.read_scalar_lay,
@ -38,6 +38,7 @@ var = {
"zk": rsws.read_scalar_lay, "zk": rsws.read_scalar_lay,
"velk": rsws.read_vector_lay, "velk": rsws.read_vector_lay,
"vz": rsws.read_scalar_lay, "vz": rsws.read_scalar_lay,
"vel": rsws.read_vector,
} }
with ThreadPool(len(var)) as pool: with ThreadPool(len(var)) as pool:

View file

@ -1,19 +1,19 @@
$************************* GENERAL *************************************** $************************* GENERAL ***************************************
PROJ 'GW' 'T1' PROJ 'GW' 'T1'
SET NAUT SET NAUT
SET LEVEL 0.5 SET LEVEL 4.5
SET MAXERR 1 SET MAXERR 1
SET DEPMIN 0.001 $SET DEPMIN 0.001
MODE DYN ONED MODE NONST ONED
$************************ GRIDS *************************************** $************************ GRIDS ***************************************
CGRID REG -1450 0 0 1450 0 1450 0 CGRID REG -1450 0 0 1750 0 1750 0
INPGRID BOT REG -1450 0 0 1450 0 1 0 $x0 y0 theta nx-1 ny-1 dx dy INPGRID BOT REG -1450 0 0 1750 0 1 0 $x0 y0 theta nx-1 ny-1 dx dy
VERT 10 $nb couches VERT 2
READ BOTTOM -1 'bathy.dat' 3 0 FREE READ BOTTOM -1 'bathy.dat' 3 0 FREE
INPGRID PORO REG -1450 0 0 1450 0 1 0 INPGRID PORO REG -1450 0 0 1750 0 1 0
INPGRID PSIZ REG -1450 0 0 1450 0 1 0 INPGRID PSIZ REG -1450 0 0 1750 0 1 0
INPGRID HSTRUC REG -1450 0 0 1450 0 1 0 INPGRID HSTRUC REG -1450 0 0 1750 0 1 0
READINP PORO 1 'poro.dat' 3 0 FREE READINP PORO 1 'poro.dat' 3 0 FREE
READINP PSIZ 1 'psize.dat' 3 0 FREE READINP PSIZ 1 'psize.dat' 3 0 FREE
@ -24,7 +24,7 @@ $*********************** BOUNDARIES ****************************************
INIT ZERO INIT ZERO
BOUN SHAP JON SIG PEAK DSPR DEGR BOUN SHAP JON SIG PEAK DSPR DEGR
BOUN SIDE W BTYPE WEAK SMOOT 10 SEC ADDBOUNDWAVE CON SPECT 7.31 16.9 0 0 CYCLE 20 MIN BOUN SIDE W BTYPE WEAK SMOOT 10 SEC ADDBOUNDWAVE CON SPECT 7.31 16.9 0 0 CYCLE 20 MIN
BOUN SIDE E BTYPE SOMMERFELD SPON E 250
$*********************** PHYSICS ******************************************* $*********************** PHYSICS *******************************************