Windwave module

\( \documentclass[12pt]{book} \usepackage{amsmath} \usepackage{epsfig} %\usepackage{graphicx} \topmargin0truemm \oddsidemargin0truemm \textheight=210truemm \flushbottom \footskip15truemm \textwidth=150truemm \renewcommand{\baselinestretch}{1.5} \newcommand{\pd}{\partial} \newcommand{\beq}{\begin{equation}} \newcommand{\beqa}{\begin{eqnarray}} \newcommand{\eeq}{\end{equation}} \newcommand{\eeqa}{\end{eqnarray}} \newcommand{\mi}{\mbox{i}} \newcommand{\kp}{\kappa} \newcommand{\md}{\mbox{d}} \newcommand{\bs}{\boldsymbol} \newcommand{\ba}{\begin{array}} \newcommand{\ea}{\end{array}} \providecommand{\D}{\Delta} \providecommand{\nn}{\nonumber} \title{\sc SELFE Theory Manual} \author{Yinglong Joseph Zhang} % \\Center for Coastal Margin Observation \and Prediction (CMOP)\\ % \\Oregon Health \and Science University} \date{\today} \begin{document} \maketitle \chapter{Transport} \section{Notation on symbols used in this manual} We typically use bold characters to denote vectors and matrices, and unbold characters to denote scalars. \section{Transport equation} \beqa \frac{D T}{D t}&=&\frac{\pd }{\pd z} (\kappa \frac{\pd T}{\pd z}) +\hat{Q},\label{tr1}\\ \kappa \frac{\pd T}{\pd z} &=& \hat{T}, \mbox{ at } z=\eta, \label{tr2} \\ \kappa \frac{\pd T}{\pd z} &=& \hat{T_b}, \mbox{ at } z=-h, \label{tr3} \eeqa where $T$ is a generic scalar variable, and the total derivative is written in a conservative form: \beq \frac{D }{D t}=\frac{\pd }{\pd t}+\nabla \cdot (\bs{u}T) \label{tr4} \eeq and the 3D velocity $\bs{u}$ is divergence free: \beq \nabla \cdot (\bs{u})=0 \label{tr5} \eeq Eq. (\ref{tr5}) guarentees the constancy condition for the transport equation i.e. $T=const$ initially will remain so in the absence of sinks/sources. \section{Finite volume for 3D continuity equation} We use F.V. method to solve the 3D continuity eq. (\ref{tr5}), not only to enforce volume conservation, but also serves as a foundation for mass conservation due to the reason stated above. The continuity eq. is solved prism-wise from bottom to surface; note that the top and bottom faces of the prism are not necessarily horizontal, and so the normal velocity there has horizontal components in it (cf. Fig. \ref{tr6}). \beqa & &\sum_{m=1}^3 \hat{P}_{js(i,m),l} s_{i,m} \frac{q_{js(i,m),l}+q_{js(i,m),l+1}}{2} + \hat{S}_{i,l+1} (\bar{u}_{i,l+1} n_{l+1}^x +\bar{v}_{i,l+1} n_{l+1}^y+\bar{w}_{i,l+1} n_{l+1}^z)- \nn \\ & & \hat{S}_{i,l}(\bar{u}_{i,l} n_{l}^x +\bar{v}_{i,l} n_{l}^y+\bar{w}_{i,l} n_{l}^z)=0 \nn \\ & & (l=k_b^e,...,N_z-1) \label{tr7} \eeqa where $\hat{P}$ and $\hat{S}$ are areas of horizontal and vertical faces of the prism, $i$ is the element number, $js(i,m)$ is the side number, $s_{i,m}$ is a sign related to the outer normal, $k_b^e$ is the bottom index of the element, $(n_l^x,n_l^y,n_l^z)$ is the unit upward normal on the 2 vertical faces, $q$ is the normal velocity at side (with positive direction following internal convention of SELFE, and thus the sign $s_{i,m}$), and $(\bar{u},\bar{v})$ is the horizontal velocity, averaged from 3 sidecenter velocities (which are the 'original' velocities defined in SELFE), on 2 vertical faces. The vertical velocity is then solved from bottom to surface starting from the bottom no-flux b.c. The eq. (\ref{tr7}) can be written in the following compact statement for volume conservation: \beq \sum_{j\in S^+} |Q_j|=\sum_{j\in S^-} |Q_j| \label{tr8} \eeq where $Q=u_n A$ is the flux, $A$ is the area of the face, and $S^\pm$ are the out/inflow faces. This eq. simply states that the incoming and outging fluxes must be balanced out against each other; this will be used below. \section{Upwind scheme} A F.V. discretization of eq. (\ref{tr1}) for a prism shown in Fig. \ref{tr6} is: \beqa T_i^{m+1} &=& T_i^{m}-\frac{\D t'}{V_i}\sum_{j\in S}Q_j T_{j*} +\hat{Q}_{i,k}^{n+1}\D t'+ \nn \\ && \frac{A_i \D t'}{V_i}[\kappa_{i,k}\frac{T_{i,k+1}^{m+1}-T_{i,k}^{m+1}}{\D z_{i,k+1/2}}- \kappa_{i,k-1}\frac{T_{i,k}^{m+1}-T_{i,k-1}^{m+1}}{\D z_{i,k-1/2}}] \nn \\ && (k=k_b^e+1,...,N_z)\label{tr9} \eeqa where $\D t' \neq \D t$ is the trasnport step (subject to Courant condition below), $V_i$ is the volume of the prism, $T_i$ is a shorthand for $T_{i,k}$ (i.e., concentration at prism $(i,k)$). Note that we have treated the diffusion term implicitly. For the sake of brevity we'll drop the source and diffusion terms from now on and focus on the advection term, where different methods diverge. The upwind concentration is defined as: \beq T_{j*}\equiv \tilde{T}_{j*}= \left\{ \ba{ll} T_i, & j\in S^+ \\ T_j, & j\in S^- \ea \right. \eeq where again we have used shorthand $T_j$ for concentration at prism $j$ (i.e. the prism adjacent to $(i,k)$ from face $j$). Eq. (\ref{tr9}) then becomes (only retaining the advective term): \beqa T_i^{m+1} &=& T_i^{m}-\frac{\D t'}{V_i}[\sum_{j\in S^+}|Q_j|T_i-\sum_{j\in S^-}|Q_j|T_j] \nn \\ &=& \left( 1-\frac{\D t'}{V_i}\sum_{j\in S^+}|Q_j|\right) T_i + \frac{\D t'}{V_i} \sum_{j\in S^+}|Q_j|T_j \label{tr10} \eeqa ($T$ without superscript are evaluated at step $m$, i.e. explicitly). We have used the volume conservation statement (\ref{tr8}) here. Therefore the Courant condition is: \beq 1-\frac{\D t'}{V_i}\sum_{j\in S^+}|Q_j| \geq 0 \label{tr11} \eeq SELFE uses this eq. to dynamically adjust time step for transport for each step. Moreover, to improve efficiency, the vertical flux terms in (\ref{tr10}) are treated implicitly, and the corresponding terms are then removed in the Courant condition (\ref{tr11}). This is allowable because upwind is a linear method. \begin{figure} [100] \centering \includegraphics{transport-prism.eps} \caption{Prism and variable locations in SELFE} \label{tr6} \end{figure} \section{TVD scheme} The only difference between TVD and upwind schemes lies in the evaluation of the interfacial concentation: \beq T_{j*}=\tilde{T}_{j*}+\frac{\phi_j}{2}(T_{jD}-\tilde{T}_{j*}) \label{tr12} \eeq where $T_{jD}$ is the downstream concentration, and $0\leq \phi_j \leq 2$ is a limiter function. TVD scheme nominally has 2nd order accuracy due to the anti-diffusion term in (\ref{tr12}). The final eq. for TVD is: \beqa T_i^{m+1} &=& T_i^{m}+\frac{\D t'}{V_i}\sum_{j\in S^-}|Q_j|(T_j-T_i) + \nn \\ &&\frac{\D t'}{V_i}\sum_j |Q_j|\frac{\phi_j}{2}(T_i-T_j)+\mbox{souce}+\mbox{diffusion} \label{tr13} \eeqa and the Courant condition is (option 'AA' in param.in): \beq \D t' \leq \frac{V_i}{\sum_{j\in S^-} |Q_j| (1-\frac{\phi_j}{2}+\delta_i)} \label{tr14} \eeq where \beqa \delta_i &=& \sum_{p\in S^+} \frac{\phi (r_p)}{2r_p} \nn \\ r_p &=& \frac{\sum_{m\in S^-} |Q_m| (T_m-T_i)}{|Q_p|(T_i-T_p)}, p\in S_i^+ \label{tr15} \eeqa The 2nd eq. involves upwind of upwind neighboring prism. Again SELFE automatically calculates the time step accroding to the Courant condition above. The choices for the limiter function include: MINMOD, OSHER, van Leer, Super Bee etc. Since TVD is a nonlinear method, we cannot treat the vertical fluxes implicitly, and so all fluxes have to be treated explicitly. TVD method is therefore more expensive than upwind. A hybrid upwind/TVD, with TVD being used in the deeper depths and upwind in the shallow depths, has been implemented in SELFE to improve efficiency. The user can also manually specify upwind/TVD zones in the domain. \section{Boundary conditions} In either upwiund or TVD schemes, the concentration at the neighboring prism $T_j$ at the open boundary is known. For outflow, $T_j=T_i$ and the signal is advected out of the domain without hindrance. For incoming flow, $T_j$ is specified by the b.c., and SELFE nudges to this value with a relaxation constant, in order to prevent sharp gradient there. \section{Interface for the generic tracer transport module} In SELFE v3.2.x, it's easy to use SELFE's generic tracer transport module by using flag\_model=-1 in param.in. The 3 arrays in selfe\_step.F90, $bdy\_frc,flx\_sf,flx\_bt$ correspond to the source, surface/bottom b.c. in Eqs. (\ref{tr1}-\ref{tr3}). You may use SELFE's conventional way of specifying i.c. (flag\_ic in param.in) or choose to specify your own by modifying sections in selfe\_init.F90. \end{document} \)