user manual

WIKI/_build/html/_sources/nearcom_formulation.rst.txt

@@ -277,3 +277,251 @@ where :math:`p` is the bed porosity, :math:`f_{mor}` represents a morphology fac
 
 
 
+
+*****************************************
+**Numerical schemes for SHORECIRC**
+*****************************************
+
+Compact form of governing equations
+#########################################
+
+The compact form of governing equations can be written as
+
+.. math:: \frac{\partial {\bf \Psi}}{\partial t} + \nabla \cdot {\bf \Theta}({\bf \Psi}) = {\bf S}
+   :label: tvd
+
+where :math:`{\bf \Psi}` and :math:`{\bf \Theta} ({\bf \Psi})` are  the vector of conserved variables and  the flux vector function, respectively, and are given by
+
+.. math:: {\bf \Psi} = \left( \begin{array}{c} J\eta \\
+          JP \\
+          JQ \end{array} \right) 
+
+.. math:: {\bf \Theta} = \left( \begin{array}{c} JP {\bf i} + JQ {\bf j} \\
+          \left[JPu +\frac{1}{2}g (\eta^2 + 2 \eta h)JL^1_1 \right ] {\bf i} + \left[JQu +\frac{1}{2}g (\eta^2 + 2 \eta h)JL^2_1 \right ] {\bf j} \\
+          \left[JPv +\frac{1}{2}g (\eta^2 + 2 \eta h)JL^1_2 \right ]  {\bf i} + \left[JQv +\frac {1}{2}g (\eta^2 + 2 \eta h)JL^2_2 \right ] {\bf j} \end{array}\right) .
+
+.. math:: {\bf S} = \left( \begin{array}{c} 0\\
+          \frac{1}{\rho} \left(F_{\mbox{wx}}  -\tau^b_x + \tau^s_x -f_x,  - ROT_x \right) \\
+          \frac{1}{\rho} \left(F_{\mbox{wy}}  -\tau^b_y + \tau^s_y  -f_y - ROT_y \right) \end{array} \right)
+
+where :math:`(u,v)` represent Cartesian components of velocity, :math:`(P,Q) = (Hu^1, Hu^2)`, the contravariant components of flow flux, (:math:`\tau^b_x, \tau^b_y`), (:math:`\tau^b_x, \tau^b_y`), (:math:`f_x, f_y`) and (:math:`ROT_x, ROT_y`) represent the bottom friction, wind stress,  Coriolis forcing and the rest of terms in :math:`x` and :math:`y` directions, respectively. 
+(:math:`F_{\mbox{wx}}, F_{\mbox{wy}}`) are the gradient of radiation stresses expressed by
+
+.. math:: F_{\mbox{wx}} = -  \left [ \frac{\partial }{\partial \xi^1} \left (S_{xx} J L^1_1+S_{xy} J L^1_2 \right) +  [ \frac{\partial }{\partial \xi^2} \left (S_{xx} J L^2_1+S_{xy} J L^2_2 \right)
+ \right ]
+
+.. math:: F_{\mbox{wy}} = -  \left [ \frac{\partial }{\partial \xi^1} \left (S_{xy} J L^1_1+S_{yy} J L^1_2 \right) + [ \frac{\partial }{\partial \xi^2} \left (S_{xy} J L^2_1+S_{yy} J L^2_2 \right) \right ]
+
+The reason to put the wave forcing as a source term rather than a flux term is that the forcing  is calculated directly from SWAN in a curvilinear grid. 
+
+
+
+
+
+Spatial discretization
+###############################
+
+A combined finite-volume and finite-difference method was applied to the spatial discretization. 
+For the flux terms and the first-order derivative terms, MUSCL-TVD schemes from the 2nd-order to 4th-order  are implemented in the present model. The multi-order MUSCL-TVD scheme can be written in a compact form  according to Erduran et al. (2005) who modified Yamamoto et al.'s (1995)  fourth-order approach. 
+In :math:`x` -direction, for example, the combined form of the interface construction can be written as follows:
+
+.. math:: \phi^L_{i+1/2} = \phi_i +\frac{1}{4} \left[ (1-\kappa_1)\chi(r) \Delta^* \phi_{i-1/2} + (1+\kappa_1) \chi(1/r) \Delta^* \phi_{i+1/2} \right ]
+   :label: left
+
+.. math:: \phi^R_{i-1/2} = \phi_i -\frac{1}{4} \left[ (1+\kappa_1)\chi(r) \Delta^* \phi_{i-1/2} + (1-\kappa_1) \chi(1/r) \Delta^* \phi_{i+1/2} \right ]
+   :label: right
+
+where :math:`\phi^L_{i+1/2}` is the constructed value at the left-hand side of the interface :math:`i+\frac{1}{2}` and :math:`\phi^R_{i-1/2}` is the value at the right-hand side of the interface :math:`i - \frac{1}{2}`.  
+The values of :math:`\Delta^* \phi` are evaluated as follows:
+
+.. math::
+   \begin{array}{l l l}
+      \Delta^* \phi_{i+1/2} = \Delta \phi_{i+1/2} - \kappa_2 \Delta^3 \bar{\phi}_{i+1/2} /6, \\
+      \Delta \phi_{i+1/2} = \phi_{i+1} - \phi_{i}, \\
+      \Delta^3 \bar{\phi}_{i+1/2} = \Delta \bar{\phi}_{i+3/2} - 2 \Delta \bar{\phi}_{i+1/2} + \Delta \bar{\phi}_{i-1/2}, \\
+      \Delta \bar{\phi}_{i-1/2} =\mbox{minmod} (\Delta \phi_{i-1/2}, \Delta \phi_{i+1/2}, \Delta \phi_{i+3/2}), \\
+      \Delta \bar{\phi}_{i+1/2} =\mbox{minmod} (\Delta \phi_{i+1/2}, \Delta \phi_{i+3/2}, \Delta \phi_{i-1/2}), \\
+      \Delta \bar{\phi}_{i+3/2} =\mbox{minmod} (\Delta \phi_{i+3/2}, \Delta \phi_{i-1/2}, \Delta \phi_{i+1/2}) 
+   \end{array}
+   :label: minmod
+
+In :eq:`minmod`, minmod represents the minmod limiter and is given by
+
+.. math:: \mbox{minmod} (j,k,l) = \mbox{sign} (j) \mbox{max} \{ 0, \mbox{min} [|j|, 2 \mbox{sign} (j) k, 2 \mbox{sign} (j) l ] \} .
+
+:math:`\kappa_1` and :math:`\kappa_2` in :eq:`left` and :eq:`right` are  control parameters for orders of the scheme in the compact form.
+The complete form with :math:`(\kappa_1, \kappa_2) = (1/3, 1)` is the fourth-order scheme given by Yamamoto et al. (1995). 
+:math:`(\kappa_1, \kappa_2) = (1/3, 0)` yields a third-order scheme, while the second-order scheme can be retrieved using :math:`(\kappa_1, \kappa_2) = (-1, 0)`. 
+
+:math:`\chi(r)` in :eq:`left` and :eq:`right` is the limiter function. 
+The original scheme introduced by Yamamoto et al. (1998) uses the Minmod limiter as used in :eq:`minmod`. Erduran et al. (2005) found that the use of the van-Leer limiter for the third-order scheme gives more accurate results. 
+Their finding was confirmed by  using the present model in the benchmark tests for wave runup conducted by Tehranirad et al. (2011).  
+The van-Leer limiter can be expressed as
+
+.. math:: \chi(r) = \frac{r+|r|}{1+r}
+
+where
+
+.. math:: r = \frac{\Delta^* \phi_{i+1/2}}{\Delta ^* \phi_{i-1/2}}.
+
+The numerical fluxes are computed using a HLL approximate Riemann solver
+
+.. math:: {\bf \Theta} ({\bf \Psi}^L, {\bf \Psi}^R)= \left \{ \begin{array}{ll} {\bf \Theta} ({\bf \Psi}^L) & \mbox{if} \ \ \ s_L \ge 0 \\
+          {\bf \Theta}^* ({\bf \Psi}^L, {\bf \Psi}^R) & \mbox{if} \ \ \ s_L < 0 < s_R \\
+          {\bf \Theta}({\bf \Psi}^R) & \mbox{if} \ \ \ s_R \le 0,\end{array}  \right.
+
+where
+
+.. math:: {\bf \Theta^*} ({\bf \Psi}^L, {\bf \Psi}^R) = \frac{s_R {\bf \Theta} ({\bf \Psi}^L) -s_L {\bf \Theta}({\bf \Psi}^R) + s_L s_R ({\bf \Psi}^R - {\bf \Psi}^L)}{s_R - s_L}
+
+The wave speeds of the Riemann solver are given by
+
+.. math:: s_L= \mbox{min} ({\bf V}^L \cdot {\bf n} - \sqrt{g (h+\eta)^L} , u_s -\sqrt{\varphi_s}),  
+   :label: sl
+
+.. math:: s_R= \mbox{max} ({\bf V}^R \cdot {\bf n} + \sqrt{g (h+\eta)^R},  u_s +\sqrt{\varphi_s}),  
+   :label: sr
+
+in which :math:`u_s` and :math:`\varphi_s` are estimated as
+
+.. math:: u_s =\frac{1}{2} ({\bf V}^L + {\bf V}^R)\cdot {\bf n} + \sqrt{g (\eta + h)^L} - \sqrt{g(\eta+h)^R}
+
+.. math:: \sqrt{\varphi_s} = \frac{\sqrt{g (\eta + h)^L}+\sqrt{g (\eta + h)^R}}{2} +\frac{({\bf V}^L - {\bf V}^R)\cdot {\bf n}}{4}
+
+and :math:`{\bf n}` is the normalized side vector for a cell face.
+
+
+
+
+
+Time stepping
+##########################
+
+The third-order Strong Stability-Preserving (SSP) Runge-Kutta scheme for nonlinear spatial discretization (Gottlieb et al., 2001) was adopted for time stepping. The scheme is given by
+
+.. math::
+   \begin{array}{l l l}
+      {\bf \Psi}^{(1)} = {\bf \Psi}^{n}  + \Delta t (- \nabla \cdot {\bf \Theta} ({\bf \Psi}^n) + {\bf S}^{(1)} ) \\
+      {\bf \Psi}^{(2)} = \frac{3}{4}{\bf \Psi}^{n}  + \frac{1}{4} \left[   {\bf \Psi}^{(1)} +  \Delta t \left (- \nabla \cdot {\bf \Theta} ({\bf \Psi}^{(1)} ) + {\bf S}^{(2)} \right) \right] \\
+      {\bf \Psi}^{n+1} = \frac{1}{3}{\bf \Psi}^{n}  + \frac{2}{3} \left[   {\bf \Psi}^{(2)} +  \Delta t \left (- \nabla \cdot {\bf \Theta} ({\bf \Psi}^{(2)} ) + {\bf S}^{n+1} \right) \right]
+   \end{array}
+   :label: runge
+
+in which :math:`{\bf \Psi}^{n}` denotes :math:`{\bf \Psi}` at time level :math:`n`. 
+:math:`{\bf \Psi}^{(1)}` and :math:`{\bf \Psi}^{(2)}` are values at intermediate stages in the Runge-Kutta integration. 
+
+  An adaptive time step is chosen, following the Courant-Friedrichs-Lewy (CFL) criterion:
+
+.. math:: \Delta t = C  \mbox{min} \left ( \mbox{min} \frac{\Delta x}{|u_{i,j}| + \sqrt{g (h_{i,j} +\eta_{i,j})}},  \mbox{min} \frac{\Delta y}{|v_{i,j}| + \sqrt{g (h_{i,j} +\eta_{i,j})}} \right )
+   :label: cfl
+
+where :math:`C` is the Courant number and :math:`C=0.5` was used in the following examples.  
+
+
+
+
+
+Wetting-drying schemes for shallow water
+##############################################
+
+The wetting-drying scheme for modeling a moving boundary is straightforward. 
+The normal flux :math:`{\bf n} \cdot {\bf M}` at the cell interface of a dry cell is set to zero. 
+A mirror boundary condition is applied to the high-order MUSCL-TVD scheme. 
+It may be noted that the wave speeds of the Riemann solver  :eq:`sl` and :eq:`sr` for a dry cell are modified as 
+
+
+.. math:: s_L= {\bf V}^L \cdot {\bf n} - \sqrt{g (h+\eta)^L} ,  \ \ \ \    s_R= {\bf V}^L \cdot {\bf n} +2  \sqrt{g (h+\eta)^L}  \  \  \  \mbox{(right dry cell)}
+
+and
+
+.. math:: s_L= {\bf V}^R \cdot {\bf n} - \sqrt{g (h+\eta)^R} ,  \ \ \ \    s_R= {\bf V}^R \cdot {\bf n} +2  \sqrt{g (h+\eta)^R}  \  \  \  \mbox{(left dry cell)}
+
+
+
+
+
+Boundary conditions
+##############################
+
+To incorporate tide and river inflow into the circulation module, two types of open boundary conditions are implemented. 
+One is the surface clamped boundary condition that reads at tidal open boundaries            
+
+.. math:: \eta= \eta_0    \ \ \ \ \ \ \mbox{at tidal open boundaries}
+
+where :math:`\eta_0` is measured or predicted surface elevations at open boundaries. 
+The other is the specified flux boundary condition that is usually used at river boundaries:
+
+.. math:: P^\alpha = Q^{\mbox{flux}}_\beta L^\alpha_\beta
+
+where :math:`Q^{\mbox{flux}}` is the specified volume flux per unit width. 
+
+
+
+
+
+
+*************************************************
+**SWAN-SHORECIRC coupling and parallelization**
+*************************************************
+
+For the SHORECIRC code, we used a domain decomposition technique to subdivide the problem into multiple regions and assign each subdomain to a separate processor core. 
+Each subdomain region contains an overlapping area of ghost cells, three-row deep, as required by the highest-order ( the fourth order) MUSCL-TVD scheme in the model. 
+The Message Passing Interface (MPI) with non-blocking communication is used to exchange  data in the overlapping region between neighboring processors. 
+
+For the SWAN code, we used the existing domain decomposition scheme which is slightly different from the SHORECIRC code.
+The version 40.51AB of SWAN uses a single-direction domain decomposition and takes into account number of dry points in the grid splitting. Because of the different domain decomposition between SHORECIRC and SWAN, it is difficult to directly pass variables between the decomposed domains. 
+We now used a two-step mapping method which is  to first gather a passing variable into the global domain, and then to distribute it into each sub-domain. 
+Dry points in SWAN are set to be wet with a small water depth (1 mm) in order to make an equal split of a computational domain. 
+
+:numref:`domain` shows an example of domain decomposition in SWAN and SHORECIRC. 
+
+To investigate performance of the parallel program, numerical simulations of an 
+idealized  case are tested with different numbers of processors  on a Linux cluster located at University of Delaware. 
+The test case is set up in a numerical grid  of 1800 :math:`\times` 1800 cells.  
+:numref:`fig1` shows the model speedup versus number of processors. 
+It can be seen that performance scales nearly proportional to the number of processors, with some delay caused by inefficiencies in parallelization, such as inter-processor communication time.
+
+
+.. figure:: figures/parallel.jpg
+   :name: domain
+   :scale: 40%
+   :align: center
+
+   An example of domain decomposition in SWAN and SHORECIRC/SEDIMENT.
+
+
+.. figure:: figures/speedup.jpg
+   :name: fig1
+   :scale: 35%
+   :align: center
+
+   Variation in model performance with number of processors for a 1800 x 1800 domain.  Straight line indicates arithmetic speedup. Actual performance is shown in the curved line **(the figure needs to be updated using the coupled model !!!)**.
+
+
+
+
+
+
+
+*******************
+**Users' Manual**
+*******************
+
+Program outline and flow chart
+##################################
+
+The code was written using Fortran 90 with the c preprocessor (cpp) statements for separation of the source code. Arrays are dynamically allocated at runtime. Precision is selected using the *selected_real_kind* Fortran intrinsic function defined in the makefile.  The default precision is single. 
+
+The present version of NearCoM-TVD includes a number of options including (1) choice of serial or parallel code (2) Cartesian or curvilinear coordinate, (3) samples.
+
+The flow chart is shown in :numref:`chart`. 
+
+.. figure:: figures/chart.jpg
+   :name: chart
+   :scale: 50%
+   :align: center
+
+   Flow chart of the main program.
+
+
+

WIKI/_build/html/index.html

@@ -55,6 +55,19 @@ <h1>NEARCOM-TVD<a class="headerlink" href="#nearcom-tvd" title="Permalink to thi
 <li class="toctree-l2"><a class="reference internal" href="nearcom_formulation.html#seabed-evolution-equation">Seabed evolution equation</a></li>
 </ul>
 </li>
+<li class="toctree-l1"><a class="reference internal" href="nearcom_formulation.html#numerical-schemes-for-shorecirc"><strong>Numerical schemes for SHORECIRC</strong></a><ul>
+<li class="toctree-l2"><a class="reference internal" href="nearcom_formulation.html#compact-form-of-governing-equations">Compact form of governing equations</a></li>
+<li class="toctree-l2"><a class="reference internal" href="nearcom_formulation.html#spatial-discretization">Spatial discretization</a></li>
+<li class="toctree-l2"><a class="reference internal" href="nearcom_formulation.html#time-stepping">Time stepping</a></li>
+<li class="toctree-l2"><a class="reference internal" href="nearcom_formulation.html#wetting-drying-schemes-for-shallow-water">Wetting-drying schemes for shallow water</a></li>
+<li class="toctree-l2"><a class="reference internal" href="nearcom_formulation.html#boundary-conditions">Boundary conditions</a></li>
+</ul>
+</li>
+<li class="toctree-l1"><a class="reference internal" href="nearcom_formulation.html#swan-shorecirc-coupling-and-parallelization"><strong>SWAN-SHORECIRC coupling and parallelization</strong></a></li>
+<li class="toctree-l1"><a class="reference internal" href="nearcom_formulation.html#users-manual"><strong>Users’ Manual</strong></a><ul>
+<li class="toctree-l2"><a class="reference internal" href="nearcom_formulation.html#program-outline-and-flow-chart">Program outline and flow chart</a></li>
+</ul>
+</li>
 </ul>
 </div>
 </div>