Commit 46467502 authored by Matthieu Flamand's avatar Matthieu Flamand

Manual added in doc/Manual and a note on how to implement the...

Manual added in doc/Manual and a note on how to implement the Saint-Venant-Exner approach in a coupled way in SLIM is also added
parent 4362afd1
Pipeline #2386 failed with stage
in 2 minutes and 22 seconds
start = 0;
end = 10;
//Create the points to draw the contour
Point(1) = {start, 0.605, 0, 1.0};
Point(2) = {start, -0.605, 0, 1.0};
Point(3) = {start, 0.2, 0, 1.0};
Point(4) = {start, 0.1492, 0, 1.0};
Point(5) = {start, -0.1492, 0, 1.0};
Point(6) = {start, -0.2, 0, 1.0};
Point(7) = {end, 0.605, 0, 1.0};
Point(8) = {end, -0.605, 0, 1.0};
Point(9) = {end, 0.2, 0, 1.0};
Point(10) = {end, -0.1492, 0, 1.0};
Point(11) = {end, 0.1492, 0, 1.0};
Point(12) = {end, -0.2, 0, 1.0};
//Create the lines between the points
Line(1) = {7, 1};
Line(2) = {2, 8};
Line(3) = {1, 3};
Line(4) = {3, 4};
Line(5) = {4, 5};
Line(6) = {5, 6};
Line(7) = {6, 2};
Line(8) = {9, 7};
Line(9) = {11, 9};
Line(10) = {10, 11};
Line(11) = {12, 10};
Line(12) = {8, 12};
Line(13) = {3, 9};
Line(14) = {4, 11};
Line(15) = {5, 10};
Line(16) = {6, 12};
//Defines the boundary of every surface. The elements will be generated separately on each surface
Line Loop(17) = {13, 8, 1, 3};
Plane Surface(18) = {17};
Line Loop(19) = {14, 9, -13, 4};
Plane Surface(20) = {19};
Line Loop(21) = {15, 10, -14, 5};
Plane Surface(22) = {21};
Line Loop(23) = {16, 11, -15, 6};
Plane Surface(24) = {23};
Line Loop(25) = {2, 12, -16, 7};
Plane Surface(26) = {25};
//Defines the physical tags that will be used to impose the boundary conditions
Physical Line("Wall") = {1, 2};
Physical Line("In") = {3, 4, 5, 6, 7};
Physical Line("Out") = {8, 9, 10, 11, 12};
Physical Surface("Domain") = {18, 20, 22, 24, 26};
//Other parameters
Mesh.CharacteristicLengthFactor = 0.1;
Mesh.Algorithm = 6;
import slimPre
import mySlim as slim
from dgpy.scripts import slim_private
import numpy
##################
## Parameters ##
##################
#File parameters
meshfile = "CompoundChannel.msh"
output_bath = "output/bath" #Output folder for the bathymetry
output_nc = "netcdf/" #Output folder for NETcdf files
mesh_start = 0
mesh_end = 10
#Time parameters
dt = 1 #Time step
tstart = 0 #Starting time for the hydrodynamics
tend = 160 #End time of the simulation
tinit = 80 #Starting time for the tracer equation
tstab = 50 #Time at which the boundary conditions reach their final value
export = dt
#Bathymetry parameters
b = 0.605
bfp = 0.405
h = 0.0752
hfp = 0.0508
s0 = 0.0019
level = 0.0119
#Parameters boundary conditions
q = 0.015
slip = 0.996
#Flow parameters
smagorinsky_coefficient = 0.4
manning_mc = 0.010
manning_w = 0.001
manning_fp = 0.016
#Parameters tracer
concentrationIn = 0.0
concentration = 0.0
source = 1.0
#Parameters sediment
initial_bottom_concentration = 0.0
windU = 1e-4
windV = 1e-4
rhosediment = 2650
porosity = 0.4
rhosedimentbottom = rhosediment*(1.-porosity)
d50 = 91e-6
d90 = 0.0004
shields = 0.0119
diffusivity = 1e-6
g = 9.81
fluxMc = 1
diffusivity_mode_constant = False
diffusivity_mc = 1.2e-3
diffusivity_w = 1.0e-3
diffusivity_fp = 0.4e-3
smagorinsky_coefficient = 0.3
h = 0.063
manning_mc = 0.010
manning_w = 0.01
manning_fp = 0.012
q = 0.015
level = 0.005
#The functions to calculate various parameters are defined. They will be used later
def bathf(x,y) :
bathi = 0.
if abs(y) >= (b-bfp) :
bathi = h-hfp
elif abs (y) >= (b-bfp-hfp) :
bathi = h-(abs(y)-(b-bfp-hfp))
else :
bathi = h
return bathi+x*s0-level
def calcEta(x,y) :
return -x*s0+level
def calcUInit(x,y) :
return 0.
def calcManning(x,y):
n = manning_mc
if abs(y) >= (b-bfp) :
n = manning_fp
elif abs (y) >= (b-bfp-hfp) :
n = manning_w
else :
n = manning_mc
return n
def calcDiffusivity(x,y):
if abs(y) >= (b-bfp) :
d = diffusivity_fp
elif abs (y) >= (b-bfp-hfp) :
d = diffusivity_mc + (abs(y)-(b-bfp-hfp))/(hfp)*(diffusivity_fp-diffusivity_mc)
else :
d = diffusivity_mc
return d
def calcC(x,y) :
c=concentration
return c
def calcSource(x,y):
c = 0
if abs(y-0)<=(b-bfp-hfp)+0.0001 and abs(x)<0.1:
c=0.8*source/(0.1*(b-bfp-hfp)*2*h*1000)*700*q
return c
def calcSed(x,y) :
sed = initial_bottom_concentration
return sed
def calcQ(t0, i,dt):
t = t0 + i*dt
discharge = q
if t<tstab:
discharge = q/tstab*t
return discharge
def calcEtaIn(t0,i,dt):
return calcEta(mesh_start,0)
def calcEtaOut(t0,i,dt):
return calcEta(mesh_end,0)
#Parameters for the sediment
factorW= 0.01
wMax=0.0001
u0erosion = 0.4
u0deposition = u0erosion
w0=10
a1=0.000001
n=4
a01=0.028
a02=0.144
omega2=2.4538
eros0Fact=0.245
if not slim_private.path.exists(output_nc):
slim_private.makedirs(output_nc)
###################
### pre process ###
###################
#Create the object of type mesh
mesh = slimPre.Mesh(meshfile)
#Create the slimPre region
region = slimPre.Region(mesh)
#array of shape (n,3), n being the number of nodes in the region, with the coordinates of the mesh nodes
xyz = region.coordinates
#The x and y coordinates are stored in a NETcdf file
slimPre.write_file(output_nc + 'coordinates.nc', region=region, time=None, data=[('x',xyz[:,0]),('y',xyz[:,1])])
#######################
### Flow parameters ###
#######################
#Initiates the vectors containing the flow parameters on the whole domain
bath = numpy.zeros(xyz.shape[0])
manning_vector = numpy.zeros(xyz.shape[0])
diffusivity_vector = numpy.zeros(xyz.shape[0])
#Fills the flow parameters as defined in the functions previously
for i in range(xyz.shape[0]) :
bath[i] = bathf(xyz[i,0], xyz[i,1])
manning_vector[i] = calcManning(xyz[i,0], xyz[i,1])
diffusivity_vector[i] = calcDiffusivity(xyz[i,0], xyz[i,1])
#The flow parameters are saved in NETcdf files, the standard input for SLIM simulations
slimPre.write_file(output_nc + 'bath.nc', region=region, time=None, data=[('bath',bath)])
slimPre.write_file(output_nc + 'manning.nc', region=region, time=None, data=[('manning',manning_vector)])
slimPre.write_file(output_nc + 'sediment.nc', region=region, time=None, data=[('diffusivity',diffusivity_vector)])
#Export a field in the netcdf format to msh format and gives the same name as the original mesh file but in the folder "output/bath"
slimPre.netcdf_to_msh(meshfile, output_nc + "bath.nc", "bath", output_bath)
#####################
### open boundary ###
#####################
#The object of type time which stores each time step is created
time_vector = numpy.arange(0,tend+2*dt,dt,numpy.float64)
time = slimPre.Time(time_vector)
#The vectors that will hold the boundary conditions are set to zero
#They have the same size as the time. There is one value for each time step
etaOut = numpy.zeros_like(time_vector)
discharge_vector = numpy.zeros_like(time_vector)
#Fills the boundary conditions as defined previously
for i in range(time_vector.shape[0]) :
etaOut[i] = calcEtaOut(tstart,i,dt)
discharge_vector[i] = calcQ(tstart,i,dt)
concentrationIn_vector=concentrationIn*numpy.ones_like(time_vector)
concentrationOut_vector=concentration*numpy.ones_like(time_vector)
#The boundary conditions are saved in two NETcdf files. One for the hydrodynamics, one for the sediments
#The NETcdf files are applied to the time (time=time)
slimPre.write_file(output_nc + 'boundaryIn.nc', region=None, time=time, data=[('dischargeIn', discharge_vector)])
slimPre.write_file(output_nc + 'boundaryOut.nc', region=None, time=time, data=[('etaOut', etaOut)])
slimPre.write_file(output_nc + 'sedimentBoundaryIn.nc', region=None, time=time, data=[('c', concentrationIn_vector)])
slimPre.write_file(output_nc + 'sedimentBoundaryOut.nc', region=None, time=time, data=[('c', concentrationOut_vector)])
#########################
### initial condition ###
#########################
#Initiates the vectors for the initial conditions to zero
etaInit = numpy.zeros(xyz.shape[0])
uInit = numpy.zeros(xyz.shape[0])
vInit = numpy.zeros(xyz.shape[0])
cInit = numpy.zeros(xyz.shape[0])
bottomConcentration_vector = numpy.zeros(xyz.shape[0])
source_vector = numpy.zeros(xyz.shape[0])
#Fills the vector the way the initial conditions were previously defined
for i in range(xyz.shape[0]) :
etaInit[i] = calcEta(xyz[i,0], xyz[i,1])
uInit[i] = calcUInit(xyz[i,0], xyz[i,1])
cInit[i] = calcC(xyz[i,0], xyz[i,1])
bottomConcentration_vector[i] = calcSed(xyz[i,0], xyz[i,1])
source_vector[i] = calcSource(xyz[i,0], xyz[i,1])
windU_vector = windU*numpy.ones(xyz.shape[0])
windV_vector = windV*numpy.ones(xyz.shape[0])
#The initial conditions are saved in two NETcdf files. One for the hydrodynamics, one for the sediments
#The NETcdf file is now applied at the domain (region=region)
slimPre.write_file(output_nc + 'init.nc', region=region, time=None, data=[('eta', etaInit),('u', uInit),('v', vInit)])
slimPre.write_file(output_nc + 'sedimentInit.nc', region=region, time=None, data=[('c', cInit),('initial_bottom_concentration', bottomConcentration_vector),('source', source_vector),('windU', windU_vector),('windV', windV_vector)])
###########
### run ###
###########
#The object of time domain is created by giving it the NETcdf file containing the bathymetry
domain = slim.Domain(meshfile,(output_nc + 'bath.nc','bath'))
#The system of the shallow water equations is created and the optional terms are added
eq = slim.ShallowWater2d(domain, "implicit",initial_time=tstart, final_time=tend, export_every_sub_time_step = False)
eq.set_viscosity("smagorinsky",smagorinsky_coefficient)
eq.set_dissipation("manning", coefficient=(output_nc + 'manning.nc','manning'))
#The initial and boundary conditions are applied
eq.set_initial_condition((output_nc + 'init.nc','eta'),(output_nc + 'init.nc','u'),(output_nc + 'init.nc','v'))
eq.set_boundary_open('In', discharge=(output_nc + 'boundaryIn.nc','dischargeIn'))
eq.set_boundary_open('Out', sse=(output_nc + 'boundaryOut.nc','etaOut'))
eq.set_boundary_coast('Wall', slip_factor=slip)
#The system of the tracer equation is created and the optional terms are added
eq2 = slim.ShallowWaterTracer2d(domain, "implicit",eq, name="tracer", offline = False, initial_time=tinit, final_time=tend)
eq2.set_diffusivity("variable", variable_diffusivity=(output_nc + 'sediment.nc','diffusivity'))
eq2.set_source((output_nc + 'sedimentInit.nc','source'))
eq2.set_sediment((output_nc + 'sedimentInit.nc','initial_bottom_concentration'),(output_nc + 'sedimentInit.nc','windU'),(output_nc + 'sedimentInit.nc','windV'), rhosediment, porosity, d50, d90, shields, False, dt, factorW, wMax, u0erosion, u0deposition, w0, a01, a02, a1, n, omega2, eros0Fact)
#The initial and boundary conditions are applied
eq2.set_initial_condition((output_nc + 'sedimentInit.nc','c'))
eq2.set_boundary_open('In', (output_nc + 'sedimentBoundaryIn.nc','c'))
eq2.set_boundary_open('Out', (output_nc + 'sedimentBoundaryOut.nc','c'))
eq2.set_boundary_coast('Wall')
#Create an object of type Loop = temporal Solver
loop = slim.Loop(maximum_time_step=dt, export_time=export)
#Add the equations to solve
loop.add_equation(eq)
loop.add_equation(eq2)
#The simulation starts
loop.run()
slimPre.exit(0)
This diff is collapsed.
%
% file: localoperator.tex
% author: Victor Brena
% description: Briefly describes properties of the local operator.
%
\chapter{Generating a mesh with GMSH}
\label{app:app01}
The .geo file used to generate the mesh of the compound channel using GMSH is presented hereunder as an example on how to use SLIM.
\lstinputlisting[language=C, frame=none, backgroundcolor=\color{white}, commentstyle=\color{orange},keywordstyle=\color{cyan}, numberstyle=\footnotesize\color{darkgray},stringstyle=\color{purple}, , numbers=left]{chapters/appendices/Compound_Channel.geo}
\chapter{Running a simulation with SLIM}
\label{app:app02}
The python script used for the compound channel is presented hereunder as an example on how to use SLIM.
\lstinputlisting[language=Python, frame=none, backgroundcolor=\color{white}, commentstyle=\color{orange},keywordstyle=\color{cyan}, numberstyle=\footnotesize\color{darkgray},stringstyle=\color{purple},numbers=left ]{chapters/appendices/Compound_Channel.py}
\
\ No newline at end of file
%% Creator: Inkscape 0.91_64bit, www.inkscape.org
%% PDF/EPS/PS + LaTeX output extension by Johan Engelen, 2010
%% Accompanies image file 'Reference.pdf' (pdf, eps, ps)
%%
%% To include the image in your LaTeX document, write
%% \input{<filename>.pdf_tex}
%% instead of
%% \includegraphics{<filename>.pdf}
%% To scale the image, write
%% \def\svgwidth{<desired width>}
%% \input{<filename>.pdf_tex}
%% instead of
%% \includegraphics[width=<desired width>]{<filename>.pdf}
%%
%% Images with a different path to the parent latex file can
%% be accessed with the `import' package (which may need to be
%% installed) using
%% \usepackage{import}
%% in the preamble, and then including the image with
%% \import{<path to file>}{<filename>.pdf_tex}
%% Alternatively, one can specify
%% \graphicspath{{<path to file>/}}
%%
%% For more information, please see info/svg-inkscape on CTAN:
%% http://tug.ctan.org/tex-archive/info/svg-inkscape
%%
\begingroup%
\makeatletter%
\providecommand\color[2][]{%
\errmessage{(Inkscape) Color is used for the text in Inkscape, but the package 'color.sty' is not loaded}%
\renewcommand\color[2][]{}%
}%
\providecommand\transparent[1]{%
\errmessage{(Inkscape) Transparency is used (non-zero) for the text in Inkscape, but the package 'transparent.sty' is not loaded}%
\renewcommand\transparent[1]{}%
}%
\providecommand\rotatebox[2]{#2}%
\ifx\svgwidth\undefined%
\setlength{\unitlength}{400.61204325bp}%
\ifx\svgscale\undefined%
\relax%
\else%
\setlength{\unitlength}{\unitlength * \real{\svgscale}}%
\fi%
\else%
\setlength{\unitlength}{\svgwidth}%
\fi%
\global\let\svgwidth\undefined%
\global\let\svgscale\undefined%
\makeatother%
\begin{picture}(1,0.81300162)%
\put(0,0){\includegraphics[width=\unitlength,page=1]{Reference.pdf}}%
\put(0.6257294,0.44969967){\color[rgb]{0,0,0}\makebox(0,0)[lt]{\begin{minipage}{0.09699441\unitlength}\raggedright $h$\end{minipage}}}%
\put(0.38666662,0.69675022){\color[rgb]{0,0,0}\makebox(0,0)[lt]{\begin{minipage}{0.14549168\unitlength}\raggedleft $\eta$\end{minipage}}}%
\put(0,0){\includegraphics[width=\unitlength,page=2]{Reference.pdf}}%
\put(0.76776947,0.70055246){\color[rgb]{0,0,0}\makebox(0,0)[lt]{\begin{minipage}{0.0655596\unitlength}\raggedright $z$\end{minipage}}}%
\put(0.86657263,0.62381709){\color[rgb]{0,0,0}\makebox(0,0)[lt]{\begin{minipage}{0.06757681\unitlength}\raggedright $x$\end{minipage}}}%
\end{picture}%
\endgroup%
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!-- Created with Inkscape (http://www.inkscape.org/) -->
<svg
xmlns:dc="http://purl.org/dc/elements/1.1/"
xmlns:cc="http://creativecommons.org/ns#"
xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns="http://www.w3.org/2000/svg"
xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
width="141.32703mm"
height="114.8991mm"
viewBox="0 0 500.76505 407.1228"
id="svg2"
version="1.1"
inkscape:version="0.91 r13725"
sodipodi:docname="Reference.svg">
<defs
id="defs4">
<marker
inkscape:isstock="true"
style="overflow:visible"
id="marker5191"
refX="0"
refY="0"
orient="auto"
inkscape:stockid="Arrow2Lend">
<path
transform="matrix(-1.1,0,0,-1.1,-1.1,0)"
d="M 8.7185878,4.0337352 -2.2072895,0.01601326 8.7185884,-4.0017078 c -1.7454984,2.3720609 -1.7354408,5.6174519 -6e-7,8.035443 z"
style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.625;stroke-linejoin:round;stroke-opacity:1"
id="path5193"
inkscape:connector-curvature="0" />
</marker>
<marker
inkscape:isstock="true"
style="overflow:visible"
id="marker5133"
refX="0.0"
refY="0.0"
orient="auto"
inkscape:stockid="Arrow2Lstart">
<path
transform="scale(1.1) translate(1,0)"
d="M 8.7185878,4.0337352 L -2.2072895,0.016013256 L 8.7185884,-4.0017078 C 6.9730900,-1.6296469 6.9831476,1.6157441 8.7185878,4.0337352 z "
style="fill-rule:evenodd;stroke-width:0.625;stroke-linejoin:round;stroke:#000000;stroke-opacity:1;fill:#000000;fill-opacity:1"
id="path5135" />
</marker>
<marker
inkscape:stockid="Arrow2Lstart"
orient="auto"
refY="0.0"
refX="0.0"
id="Arrow2Lstart"
style="overflow:visible"
inkscape:isstock="true"
inkscape:collect="always">
<path
id="path4189"
style="fill-rule:evenodd;stroke-width:0.625;stroke-linejoin:round;stroke:#000000;stroke-opacity:1;fill:#000000;fill-opacity:1"
d="M 8.7185878,4.0337352 L -2.2072895,0.016013256 L 8.7185884,-4.0017078 C 6.9730900,-1.6296469 6.9831476,1.6157441 8.7185878,4.0337352 z "
transform="scale(1.1) translate(1,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lstart"
orient="auto"
refY="0.0"
refX="0.0"
id="Arrow1Lstart"
style="overflow:visible"
inkscape:isstock="true">
<path
id="path4171"
d="M 0.0,0.0 L 5.0,-5.0 L -12.5,0.0 L 5.0,5.0 L 0.0,0.0 z "
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt;stroke-opacity:1;fill:#000000;fill-opacity:1"
transform="scale(0.8) translate(12.5,0)" />
</marker>
<marker
inkscape:stockid="Arrow2Lend"
orient="auto"
refY="0"
refX="0"
id="marker4595"
style="overflow:visible"
inkscape:isstock="true"
inkscape:collect="always">
<path
inkscape:connector-curvature="0"
id="path4597"
style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.625;stroke-linejoin:round;stroke-opacity:1"
d="M 8.7185878,4.0337352 -2.2072895,0.01601326 8.7185884,-4.0017078 c -1.7454984,2.3720609 -1.7354408,5.6174519 -6e-7,8.035443 z"
transform="matrix(-1.1,0,0,-1.1,-1.1,0)" />
</marker>
<marker
inkscape:isstock="true"
style="overflow:visible"
id="marker4497"
refX="0"
refY="0"
orient="auto"
inkscape:stockid="Arrow2Lend"
inkscape:collect="always">
<path
transform="matrix(-1.1,0,0,-1.1,-1.1,0)"
d="M 8.7185878,4.0337352 -2.2072895,0.01601326 8.7185884,-4.0017078 c -1.7454984,2.3720609 -1.7354408,5.6174519 -6e-7,8.035443 z"
style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.625;stroke-linejoin:round;stroke-opacity:1"
id="path4499"
inkscape:connector-curvature="0" />
</marker>
<marker
inkscape:isstock="true"
style="overflow:visible"
id="marker4469"
refX="0"
refY="0"
orient="auto"
inkscape:stockid="Arrow2Lend">
<path
transform="matrix(-1.1,0,0,-1.1,-1.1,0)"
d="M 8.7185878,4.0337352 -2.2072895,0.01601326 8.7185884,-4.0017078 c -1.7454984,2.3720609 -1.7354408,5.6174519 -6e-7,8.035443 z"
style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.625;stroke-linejoin:round;stroke-opacity:1"
id="path4471"
inkscape:connector-curvature="0" />
</marker>
<marker
inkscape:isstock="true"
style="overflow:visible"
id="marker4447"
refX="0"
refY="0"
orient="auto"
inkscape:stockid="Arrow2Lend">
<path
transform="matrix(-1.1,0,0,-1.1,-1.1,0)"
d="M 8.7185878,4.0337352 -2.2072895,0.01601326 8.7185884,-4.0017078 c -1.7454984,2.3720609 -1.7354408,5.6174519 -6e-7,8.035443 z"
style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.625;stroke-linejoin:round;stroke-opacity:1"
id="path4449"
inkscape:connector-curvature="0" />
</marker>
<marker
inkscape:stockid="Arrow2Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow2Lend"
style="overflow:visible"
inkscape:isstock="true"
inkscape:collect="always">
<path
id="path4192"
style="fill:#000000;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:0.625;stroke-linejoin:round;stroke-opacity:1"
d="M 8.7185878,4.0337352 -2.2072895,0.01601326 8.7185884,-4.0017078 c -1.7454984,2.3720609 -1.7354408,5.6174519 -6e-7,8.035443 z"
transform="matrix(-1.1,0,0,-1.1,-1.1,0)"
inkscape:connector-curvature="0" />
</marker>
</defs>
<sodipodi:namedview
id="base"
pagecolor="#ffffff"
bordercolor="#666666"
borderopacity="1.0"
inkscape:pageopacity="0.0"
inkscape:pageshadow="2"
inkscape:zoom="1.979899"
inkscape:cx="402.4924"
inkscape:cy="292.3113"
inkscape:document-units="px"