Skip to content
GitLab
Commit a4dc26f5 authored by Jérémy Baudry's avatar Jérémy Baudry
Browse files

ajout documentation

parent a7863271
Hide whitespace changes
Inline Side-by-side
documentation/doc.tex 0 → 100644
View file @ a4dc26f5
\documentclass[10pt,a4paper]{report}
\usepackage[utf8x]{inputenc}
\usepackage{ucs}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{float}
\usepackage{graphicx}
\usepackage{listings}
\usepackage{color} %May be necessary if you want to color links
\usepackage{hyperref}
\hypersetup{
colorlinks=true, %set true if you want colored links
linktoc=all, %set to all if you want both sections and subsections linked
linkcolor=blue, %choose some color if you want links to stand out
}
%\usepackage[obeyspaces]{url}
\newcommand{\inlinecode}{\texttt}
%\author{Dany Dumont, Eliott Bismuth, Jérémy Baudry}
\title{User Manual and Documentation of WIM}
\begin{document}
\maketitle
\tableofcontents
\chapter{Introduction}
WIM (Wave in-ice Model) is a simple conceptual 1D spectral model designed to assess some processes evolving in the interactions between waves and sea ice.
\chapter{Physical aspects}
The model capture two main features of the interactions between waves and sea ice:
The attenuation of wave energy by ice and the fracturation of the ice by waves.
This section review the physical processes implemented in WIM.
\section{The basics: The wave spectrum}
A record of the sea surface elevation $\eta(t)$ can be considered to be the sum of a large number of harmonic waves with random amplitude and phase.
\begin{equation}
\eta(t)=\sum_{i=1}^{N}\underline{a_i}cos(2\pi f_it+\underline{\alpha_i})
\end{equation}
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5]{harmonic}
\caption{A wave record is a sum of many harmonic waves with different amplitudes and phases }
\end{center}
\end{figure}
The phases and amplitudes being random variables, they can be characterised with their probability density function. For each frequency, the phase $\underline{\alpha_i}$ is uniformly distributed between 0 and 2$\pi$ and the amplitude $\underline{a_i}$ is rayleigh distributed.
\begin{equation}
P(a_i)=\frac{\pi}{2}\frac{a_i}{\mu_i^2}\text{exp}\left( -\frac{\pi a_i^2}{4\mu_i^2}\right)
\end{equation}
where $\mu_i$ is the expected value of the amplitude $\mu_i=E\lbrace\underline{a_i}\rbrace$.
The function that shows this expected value over each frequency is called the amplitude spectrum.
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5]{spectre1}
\caption{amplitude spectrum}
\end{center}
\end{figure}
But it is more relevant to consider the variance density spectrum
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5]{spectre2}
\caption{variance density spectrum}
\end{center}
\end{figure}
By multiplying the variance density spectrum by $\rho g$, we obtain the \emph{energy spectrum}
\begin{equation}
S_{energy}=\rho gS_{variance}
\end{equation}
In WIM, the energy spectrum is described using an analytical function. There are many differents equation used to model the wave spectrum, in WIM it is possible to choose between a JONSWAP or a Bretschneider type spectrum. (see section 4.2.3)
\section{Advection}
The wave spectrum is advected by solving
\begin{equation}
\frac{\partial S}{\partial t}+c_g\frac{\partial S}{\partial x}=0
\end{equation}
\section{Attenuation}
The attenuation source term is defined as:
\begin{equation}
R_{ice}=\widehat{\alpha} S
\end{equation}
where $\widehat{\alpha}$ is a dimensional attenuation coefficient defined such as:
\begin{equation}
\widehat{\alpha}=\frac{\alpha c}{\left\langle D\right\rangle }
\end{equation}
where $\alpha$ is the dimensionless attenuation coefficient, i.e. the
(average) amount of attenuation per individual floe, which is a function of ice thickness and wave period.
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5]{alpha}
\caption{Dimensionless energy attenuation coefficient $\alpha$ given bythe model of kohout et meylan (2008) }
\end{center}
\end{figure}
\section{floe breaking}
The model of Williams et al.(2013) is used to parametrize the wave induced floe breaking. \\
The flexural strain imposed by passing waves on ice is defined as:
\begin{equation}
\varepsilon=\frac{h}{2}\frac{\partial^2\eta}{\partial x^2}
\end{equation}
where $\eta$ is the sea surface elevation,$h$ the ice thickness and $x$ the horizontal distance.
We can also define the significant strain
\begin{equation}
E_s=2\sqrt{\left\langle \varepsilon^2\right\rangle}
\end{equation}
The mean square value is $\left\langle \varepsilon^2\right\rangle=m_0\left[ \varepsilon\right] $
where $m_0\left[ \varepsilon\right] $ is the 0th moment of the strain spectrum.
For a sinusoidal wave of the form $A_{ice}sin\left(k_{ice}x-\omega t\right) $ the maximum strain is \begin{equation}
E=\frac{h}{2}k_{ice}^2A_{ice}(\omega)
\end{equation}
where $k_{ice}$ is the wavelength in ice and $A_{ice}$ is the amplitude in the ice.
We can then define the 0th moment of the strain spectrum
\begin{equation}
m_0\left[ \varepsilon\right]=\int S(\omega)E^2(\omega)d\omega
\end{equation}
If the significant strain $E_s$ exceed a certain threshold $E_s>\varepsilon_c$ then the ice breaks and the maximumfloesiwe is set to
\begin{equation}
D_{max}=max\left[ \lambda_w/2,D_{min}\right]
\end{equation}
where $D_{min}$ is the minimum floe size and $\lambda_w$ the dominant wavelength defined as $ \lambda_w=2\pi/k_w$ and $k_w=k_{ice}(2\pi/T_w)$.
$T_w$ is the dominant period of the spectrum defined as:
\begin{equation}
T_w=2\pi\sqrt{\frac{m_0\left[\eta \right] }{m_2\left[\eta \right]}}
\end{equation}
where $m_n\left[\eta \right]$ is the nth moment of the wave spectrum defined as:
\begin{equation}
m_n\left[\eta \right]=\int\omega^nS(\omega)d\omega
\end{equation}
\chapter{Getting started}
\section{Get WIM source code}
The WIM source code is available on the plateform \url{https://gitlasso.uqar.ca/}.
To download WIM from this plateform, please execute the following steps:
\begin{enumerate}
\item create the destination folder where the source code will be download by typing \inlinecode{mkdir <directory name>}
\item go to the destination folder \inlinecode{cd <directory name>} and execute \inlinecode{ git clone https://gitlasso.uqar.ca/bauj0001/WIM2.git .}
\end{enumerate}
\subsection{Contents}
There are 3 subdirectories:
\begin{enumerate}
\item \path{src}: This directory contains all the source files of the routines and subroutines of WIM and the makefile to link and compile all source files
\item \path{nml}: Contains the namelist file which allow to change model parameters without recompiling the code
\item \path{output}: Contains the outputs in NETCDF format.
\end{enumerate}
\section{Compile and execute WIM}
Compilation in WIM is relatively easy. The following steps should however be executed correctly to compile the code without errors.
\begin{enumerate}
\item Open the \textbf{makefile} in \path{WIM/src/}
\item Choose the compiler that will be used to produce the executable by changing the variable COMPILER. By default, the compiler used is GFORTRAN but it should work with other compilers like INTEL FORTRAN
\item WIM uses the NETCDF library to save outputs. Ensure first that the fortran NETCDF library is currently installed on your machine. If it is not the case, you can download it at \url{http://www.unidata.ucar.edu/downloads/netcdf/index.jsp}. Once the NETCDF library is installed on your machine, you should link it to WIM by indicating the path of the includes files, modules and libraries in the variables NETCDFINC, NETCDFMOD, NETCDFLIB.
\item You can now compile the code by typing \inlinecode{make} in \path{src/}. The executable will be created in the parent directory \path{WIM/}.You can also type \inlinecode{make clean} or \inlinecode{make mrproper} to remove the executable, the objects .o and modules .mod if needed.
\item Now execute WIM by typing \inlinecode{./<executable name>}
\end{enumerate}
\section{Beginning with WIM: using namelists}
Designing a simulation with WIM is quite easy. The WIM model uses \textit{namelists} allowing to change some parameter values without recompiling the source code. afterchanging parameters in the namelists, the values will be read by the routine \textbf{parameters.f90} (see \textit{routines details} for more informations). All namelists are written in a single file \path{nml/parameter.nml}. The file is composed by 4 namelists:
\begin{itemize}
\item \textbf{model paramete}r: This namelist allow to change some global parameters of the simulations such as the siwe of the domain, the name of the output directory, the name of the simulation etc..
\item \textbf{waves parameters}: Specifies the waves parameters such as the significant height, the peak period and some features like allowing dispersion or not.
\item \textbf{spectrum parameters}: Specifies the parameters to construct the wave spectrum such as the number of frequency bin, the minimal and maximal frequency.
\item \textbf{ice parameters}: Specifies some ice features like the ice concentration, thickness, initial mean diameters of the floes etc..
\end{itemize}
After having changed the parameters in the namelists, you can now run WIM.
\chapter{Implementation}
\section{General}
\subsection{General algorithm}
\begin{verbatim}
read namelists
allocate memory
initialize
do time loop
advection
do space loop
attenuation
floe breaking
compute new FSD
end space loop
end time loop
write outputs
deallocate memory
\end{verbatim}
\subsection{Solving ODEs}
Time ordinary differential equations are solved in WIM using the first-order Euler method.
For a ODE of the general form
\begin{equation}
\frac{dX}{dt}=F(X,t)
\end{equation}
The euler method give:
\begin{equation}
X_{n+1}=X_n + F(X_n,t_n)\Delta t
\end{equation}
\section{Routine details}
\subsection{main.f90}
This is the main program of WIM. It calls all the other subroutines and does the time Loop. It also contains the subroutine \inlinecode{progress} which displays a progress bar in the consol while the model is running.
\subsection{parameters.f90}
The module \path{parameters.f90} contains all variable and parameters definitions used in the model. All others routines call this module with the Fortran function \inlinecode{USE PARAMETERS}. This module also contains some other subroutines:
\begin{enumerate}
\item \textbf{read namelists}: Read parameters value specified in the differents namelists contained in the file \inlinecode{parameters.nml} (see Beginning with WIM: using namelists)
\item \textbf{array allocation}: Allocate memory for the different arrays and matrix.
\end{enumerate}
\subsection{initialization.f90}
In this routine, the initial spectrum and ice transect are built and initial values for arrays are set.
\subsubsection{initial spectrum}
It is possible de choose between different methods to build the initial spectrum by setting the parameter \inlinecode{init spec} in the namelist \inlinecode{spectrum parameters}. If \inlinecode{init spec}=1 then a JONSWAP spectrum is built according to
\begin{equation}
S_{init}=\alpha H_s^2\frac{f_p^4}{f^5}\text{exp}\left[\frac{-5}{4}\left(\frac{f_p}{f} \right)^4 \right] \gamma^{\text{exp}\left[ \frac{-(f-fp)^2}{2\sigma^2f_p^2}\right] }
\end{equation}
where
\[\alpha=\frac{0.0624}{0.23+0.0336\gamma-\frac{0.185}{1.9+\gamma}}\]
and
\[
\sigma=
\begin{cases}
0.07, & \text{if } f\leq f_p \\
0.09 , & \text{if } f> f_p
\end{cases}\]
(see \url{http://hmf.enseeiht.fr/travaux/CD0910/bei/beiere/groupe4/node/59}).\\
If \inlinecode{init spec}=2 then a Bretschneider spectrum is built according to
\begin{equation}
S_{init}=\frac{1.25H_s^2T^5}{8\pi T_p^2}e^{-1.25(T/T_p)^4}
\end{equation}
where $T_p$ is the peak period and $H_s$ the significant height
\subsubsection{Ice transect}
Parameters for ice transect can be set in the namelist \inlinecode{ice parameters}. Ice concentration is uniform through the domain
\begin{equation}
C(x)=
\begin{cases}
0, & \text{if } x<X_1 \\
\inlinecode{cice} , & \text{if } X_1<x<X_2
\end{cases}
\end{equation}
If \inlinecode{ice thick}=0 then the ice thickness is constant along the transect. if \inlinecode{ice thick}=1 the ice thickness depends on the distance from the ice edge according to
\begin{equation}
h(x)=
\begin{cases}
0, & \text{if } x<X_1 \\
h_{\infty}\left(0.1+0.9\left(1-e^{-\left(x-X_1/X_h \right) }\right) \right) , & \text{if } X_1<x<X_2
\end{cases}
\end{equation}
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5]{icetransect}
\caption{Schematical visualization of the ice transect}
\end{center}
\end{figure}
\subsection{advection.f90}
Here the wave spectrum is advected through the domain by solving
\begin{equation}
\frac{\partial S}{\partial t}+c_g\frac{\partial S}{\partial x}=0
\end{equation}.
This equation is discretized using the second order Lax-Wendroff sheme with a superbee flux limiter.
\subsubsection{Lax-Wendroff without Flux limiter}
The Lax-Wendroff sheme is a \textit{multistep} method where the quantity $S_i^{t+1}$ is calculated using the quantities $S_{i-\frac{1}{2}}^{t+\frac{1}{2}}$
and $S_{i+\frac{1}{2}}^{t+\frac{1}{2}}$.\\
Using the central difference sheme, the quantity $S_i^{t+1}$ can be approximated by
\begin{equation}
S_i^{t+1}=S_i^{t}-\frac{c_g\Delta t}{\Delta x}\left(S_{i+\frac{1}{2}}^{t+\frac{1}{2}}-S_{i-\frac{1}{2}}^{t+\frac{1}{2}} \right)
\end{equation}
The both quantities $S_{i\pm\frac{1}{2}}^{t+\frac{1}{2}}$ are obtained by the Lax method:
\begin{equation}
S_{i-\frac{1}{2}}^{t+\frac{1}{2}}=\frac{1}{2}\left(S_i^{t}+S_{i-1}^{t}\right)-\frac{c_g\Delta t}{2\Delta x}\left(S_i^{t}-S_{i-1}^{t} \right)
\end{equation}
\begin{equation}
S_{i+\frac{1}{2}}^{t+\frac{1}{2}}=\frac{1}{2}\left(S_i^{t}+S_{i+1}^{t}\right)-\frac{c_g\Delta t}{2\Delta x}\left(S_{i+1}^{t}-S_{i}^{t} \right)
\end{equation}
\begin{figure}[H]
\includegraphics[scale=0.6]{laxW}
\caption{Schematical visualization of the Lax-Wendroff method}
\end{figure}
The Lax-Wendroff sheme is then:
\begin{equation}
S^{t+1}_i=S^{t}_i-\frac{c_g\Delta t}{2\Delta x}\left(S^{t}_{i+1}-S^{t}_{i-1} \right)+\frac{c_g^2\Delta t^2}{2\Delta x^2}\left(S^{t}_{i+1}-2S^{t}_{i}+S^{t}_{i-1} \right)
\end{equation}
\subsubsection{The Superbee Flux limiter}
It can be seen above that the Lax-Wendroff sheme is a combination of a 1-order term and a 2-order term. To avoid oscillation near steps due to this second term, we add a flux limiter which will "weight" the two terms.
\begin{equation}
\phi(\theta)=\text{max}\left[0,\text{min}\left(2\theta,1 \right),\text{min}\left(\theta,2 \right) \right]
\end{equation}
where $\theta$ is given by
\begin{equation}
\theta=\frac{S^t_{i+1}-S^t_{i}}{S^t_{i}-S^t_{i-1}}
\end{equation}
The Lax-Wendroff with superbee flux limiter is then
\begin{equation}
S^{t+1}_i=S^{t}_i-\phi(\theta)\frac{c_g\Delta t}{2\Delta x}\left(S^{t}_{i+1}-S^{t}_{i-1} \right)+\phi(\theta)\frac{c_g^2\Delta t^2}{2\Delta x^2}\left(S^{t}_{i+1}-2S^{t}_{i}+S^{t}_{i-1} \right)
\end{equation}
\subsubsection{Boundary conditions}
Boundary conditions are set according to $\Delta S=0$
\subsubsection{Wave dispersion}
In WIM it is possible to allow or not wave dispersion by changing the parameter \inlinecode{disp} in the namelist \inlinecode{waves parameters}.
If \inlinecode{disp} take the value 1, dispersion is allowed and $c_g$ depend on the wave period.
If \inlinecode{disp} take the value 0, the spectrum is advected at all frequencies at the same speed $c_g=max[c_g(\omega)]$
\subsection{attenuation.f90}
In this routine, the dimensionless attenuation coefficient $\alpha$ from Kohout and Meylan (2008) is approximated and the attenuation source term is computed
\subsection{fsd build.f90}
In this routine the average floe size $\left\langle D\right\rangle $ for wave attenuation calculation. It is possible to choose between two methods by setting the parameter \inlinecode{FSD sheme} in the namelist \inlinecode{model parameters}.
If \inlinecode{FSD sheme}$=0$ then the method proposed by Dumont et al (2011) is used.
\begin{equation}
\left\langle D\right\rangle=\frac{\sum_{m=0}^{M}(\xi^2f)^m\xi^{-m}D_{max}}{\sum_{m=0}^{M}(\xi^2f)^m}
\end{equation}
where $f$ is the breaking probability (ie. the fragility of the floes), $D_{min}$ the minimum floe size, $\xi$ a parameter which determine the number of pieces each floe will be fragmented into and $M$ the number of fragmentation steps defined as:
\begin{equation}
M=\text{log}_\xi(D_{max}/D_{min})
\end{equation}
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.4]{break}
\caption{Schematical visualization of the ice transect}
\end{center}
\end{figure}
if \inlinecode{FSD sheme}$=1$, a power law is used to compute $\left\langle D\right\rangle $
\begin{equation}
\left\langle D\right\rangle=A^{-1}\int_{D_{min}}^{D_{max}}D^{-\gamma +1}
\end{equation}
where $\gamma$ is the exponent of the power law and $A$ a coefficient defined as
\begin{equation}
A=\frac{1}{1-\gamma}\left( D_{max}^{1-\gamma}- D_{min}^{1-\gamma}\right)
\end{equation}
\subsection{write output.f90}
In this routine, the output file is created in NETCDF format. Here is an example of an output file:
\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5]{output}
\caption{Example of an output file}
\end{center}
\end{figure}
\chapter{Graphical User Interface}
\begin{center}
This section is not yet available
\end{center}
\end{document}
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment