/[MITgcm]/manual/s_algorithm/text/spatial-discrete.tex

Diff of /manual/s_algorithm/text/spatial-discrete.tex

Parent Directory | Revision Log | View Revision Graph Revision Graph | View Patch Patch

-revision 1.5 by adcroft,
Thu Aug  9 20:45:27 2001 UTC
+revision 1.19 by jmc,
Wed Apr  5 01:16:27 2006 UTC
 Line 2
  % $Name$
  \section{Spatial discretization of the dynamical equations}
+ \begin{rawhtml}
+ <!-- CMIREDIR:spatial_discretization_of_dyn_eq: -->
+ \end{rawhtml}
  Spatial discretization is carried out using the finite volume
  method. This amounts to a grid-point method (namely second-order
  centered finite difference) in the fluid interior but allows
  boundaries to intersect a regular grid allowing a more accurate
  representation of the position of the boundary. We treat the
- horizontal and veritical directions as seperable and differently.
+ horizontal and vertical directions as separable and differently.
  \input{part2/notation}
  \subsection{The finite volume method: finite volumes versus finite difference}
+ \begin{rawhtml}
+ <!-- CMIREDIR:finite_volume: -->
+ \end{rawhtml}
  The finite volume method is used to discretize the equations in
  space. The expression ``finite volume'' actually has two meanings; one
- is the method of cut or instecting boundaries (shaved or lopped cells
+ is the method of embedded or intersecting boundaries (shaved or lopped
- in our terminology) and the other is non-linear interpolation methods
+ cells in our terminology) and the other is non-linear interpolation
- that can deal with non-smooth solutions such as shocks (i.e. flux
+ methods that can deal with non-smooth solutions such as shocks
- limiters for advection). Both make use of the integral form of the
+ (i.e. flux limiters for advection). Both make use of the integral form
- conservation laws to which the {\it weak solution} is a solution on
+ of the conservation laws to which the {\it weak solution} is a
- each finite volume of (sub-domain). The weak solution can be
+ solution on each finite volume of (sub-domain). The weak solution can
- constructed outof piece-wise constant elements or be
+ be constructed out of piece-wise constant elements or be
  differentiable. The differentiable equations can not be satisfied by
  piece-wise constant functions.
-Line 37 
 the lengths $\Delta x_i$:
+Line 45 
 the lengths $\Delta x_i$:
  \begin{displaymath}
  \Delta x \partial_t \theta + \delta_i ( F ) = 0
  \end{displaymath}
- is exact if $\theta(x)$ is peice-wise constant over the interval
+ is exact if $\theta(x)$ is piece-wise constant over the interval
  $\Delta x_i$ or more generally if $\theta_i$ is defined as the average
  over the interval $\Delta x_i$.
-Line 57 
 recovers the same ODE's resulting from f
+Line 65 
 recovers the same ODE's resulting from f
  interior of a fluid. Differences arise at boundaries where a boundary
  is not positioned on a regular or smoothly varying grid. This method
  is used to represent the topography using lopped cell, see
- \cite{Adcroft98}. Subtle difference also appear in more than one
+ \cite{adcroft:97}. Subtle difference also appear in more than one
  dimension away from boundaries. This happens because the each
- direction is discretized independantly in the finite difference method
+ direction is discretized independently in the finite difference method
  while the integrating over finite volume implicitly treats all
  directions simultaneously. Illustration of this is given in
- \cite{Adcroft02}.
+ \cite{ac:02}.
  \subsection{C grid staggering of variables}
  \begin{figure}
- \centerline{ \resizebox{!}{2in}{ \includegraphics{part2/cgrid3d.eps}} }
+ \begin{center}
+ \resizebox{!}{2in}{ \includegraphics{part2/cgrid3d.eps}}
+ \end{center}
  \caption{Three dimensional staggering of velocity components. This
  facilitates the natural discretization of the continuity and tracer
  equations. }
-Line 77 
 equations. }
+Line 87 
 equations. }
  The basic algorithm employed for stepping forward the momentum
  equations is based on retaining non-divergence of the flow at all
  times. This is most naturally done if the components of flow are
- staggered in space in the form of an Arakawa C grid \cite{Arakawa70}.
+ staggered in space in the form of an Arakawa C grid \cite{arakawa:77}.
  Fig. \ref{fig:cgrid3d} shows the components of flow ($u$,$v$,$w$)
  staggered in space such that the zonal component falls on the
- interface between continiuty cells in the zonal direction. Similarly
+ interface between continuity cells in the zonal direction. Similarly
- for the meridional and vertical directions.  The continiuty cell is
+ for the meridional and vertical directions.  The continuity cell is
  synonymous with tracer cells (they are one and the same).
-Line 113 
 grid data: ({\em model/inc/GRID.h})
+Line 123 
 grid data: ({\em model/inc/GRID.h})
  \subsection{Horizontal grid}
+ \label{sec:spatial_discrete_horizontal_grid}
  \begin{figure}
- \centerline{ \begin{tabular}{cc}
+ \begin{center}
+ \begin{tabular}{cc}
    \raisebox{1.5in}{a)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Ac.eps}}
  & \raisebox{1.5in}{b)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Az.eps}}
  \\
    \raisebox{1.5in}{c)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Au.eps}}
  & \raisebox{1.5in}{d)}\resizebox{!}{2in}{ \includegraphics{part2/hgrid-Av.eps}}
- \end{tabular} }
+ \end{tabular}
+ \end{center}
  \caption{
  Staggering of horizontal grid descriptors (lengths and areas). The
  grid lines indicate the tracer cell boundaries and are the reference
-Line 136 
 is bordered by the lengths $\Delta x_f$
+Line 149 
 is bordered by the lengths $\Delta x_f$
  The model domain is decomposed into tiles and within each tile a
  quasi-regular grid is used. A tile is the basic unit of domain
- decomposition for parallelization but may be used whether parallized
+ decomposition for parallelization but may be used whether parallelized
- or not; see section \ref{sect:tiles} for more details. Although the
+ or not; see section \ref{sect:domain_decomposition} for more details.
- tiles may be patched together in an unstructured manner
+ Although the tiles may be patched together in an unstructured manner
  (i.e. irregular or non-tessilating pattern), the interior of tiles is
- a structered grid of quadrilateral cells. The horizontal coordinate
+ a structured grid of quadrilateral cells. The horizontal coordinate
  system is orthogonal curvilinear meaning we can not necessarily treat
- the two horizontal directions as seperable. Instead, each cell in the
+ the two horizontal directions as separable. Instead, each cell in the
  horizontal grid is described by the length of it's sides and it's
  area.
  The grid information is quite general and describes any of the
  available coordinates systems, cartesian, spherical-polar or
  curvilinear. All that is necessary to distinguish between the
- coordinate systems is to initialize the grid data (discriptors)
+ coordinate systems is to initialize the grid data (descriptors)
  appropriately.
  In the following, we refer to the orientation of quantities on the
-Line 288 
 using\-Cartes\-ianGrid} in namelist {\em
+Line 301 
 using\-Cartes\-ianGrid} in namelist {\em
  spacing can be set to uniform via scalars {\bf dXspacing} and {\bf
  dYspacing} in namelist {\em PARM04} or to variable resolution by the
  vectors {\bf DELX} and {\bf DELY}. Units are normally
- meters. Non-dimensional coordinates can be used by interpretting the
+ meters. Non-dimensional coordinates can be used by interpreting the
  gravitational constant as the Rayleigh number.
  \subsubsection{Spherical-polar coordinates}
-Line 313 
 other grids, the horizontal grid descrip
+Line 326 
 other grids, the horizontal grid descrip
  \subsection{Vertical grid}
  \begin{figure}
- \centerline{ \begin{tabular}{cc}
+ \begin{center}
+   \begin{tabular}{cc}
    \raisebox{4in}{a)} \resizebox{!}{4in}{
    \includegraphics{part2/vgrid-cellcentered.eps}} & \raisebox{4in}{b)}
    \resizebox{!}{4in}{ \includegraphics{part2/vgrid-accurate.eps}}
- \end{tabular} }
+ \end{tabular}
+ \end{center}
  \caption{Two versions of the vertical grid. a) The cell centered
  approach where the interface depths are specified and the tracer
  points centered in between the interfaces. b) The interface centered
-Line 340 
 The vertical grid is calculated in subro
+Line 355 
 The vertical grid is calculated in subro
  INI\_VERTICAL\_GRID} and specified via the vector {\bf DELR} in
  namelist {\em PARM04}. The units of ``r'' are either meters or Pascals
  depending on the isomorphism being used which in turn is dependent
- only on the choise of equation of state.
+ only on the choice of equation of state.
  There are alternative namelist vectors {\bf DELZ} and {\bf DELP} which
  dictate whether z- or
-Line 354 
 vertical grid descriptors are stored in
+Line 369 
 vertical grid descriptors are stored in
  The above grid (Fig.~\ref{fig:vgrid}a) is known as the cell centered
  approach because the tracer points are at cell centers; the cell
- centers are mid-way between the cell interfaces. An alternative, the
+ centers are mid-way between the cell interfaces.
- vertex or interface centered approach, is shown in
+ This discretization is selected when the thickness of the
+ levels are provided ({\bf delR}, parameter file {\em data},
+ namelist {\em PARM04})
+ An alternative, the vertex or interface centered approach, is shown in
  Fig.~\ref{fig:vgrid}b. Here, the interior interfaces are positioned
  mid-way between the tracer nodes (no longer cell centers). This
  approach is formally more accurate for evaluation of hydrostatic
  pressure and vertical advection but historically the cell centered
  approach has been used. An alternative form of subroutine {\em
  INI\_VERTICAL\_GRID} is used to select the interface centered approach
- but no run time option is currently available.
+ This form requires to specify $Nr+1$ vertical distances {\bf delRc}
+ (parameter file {\em data}, namelist {\em PARM04}, e.g.
+ {\em verification/ideal\_2D\_oce/input/data})
+ corresponding to surface to center, $Nr-1$ center to center, and center to
+ bottom distances.
+ %but no run time option is currently available.
  \fbox{ \begin{minipage}{4.75in}
  {\em S/R INI\_VERTICAL\_GRID} ({\em
-Line 380 
 $\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\e
+Line 403 
 $\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\e
  \subsection{Topography: partially filled cells}
+ \begin{rawhtml}
+ <!-- CMIREDIR:topo_partial_cells: -->
+ \end{rawhtml}
  \begin{figure}
- \centerline{
+ \begin{center}
  \resizebox{4.5in}{!}{\includegraphics{part2/vgrid-xz.eps}}
- }
+ \end{center}
  \caption{
  A schematic of the x-r plane showing the location of the
  non-dimensional fractions $h_c$ and $h_w$. The physical thickness of a
-Line 393 
 thickness of the open side is given by $
+Line 419 
 thickness of the open side is given by $
  \label{fig:hfacs}
  \end{figure}
- \cite{Adcroft97} presented two alternatives to the step-wise finite
+ \cite{adcroft:97} presented two alternatives to the step-wise finite
  difference representation of topography. The method is known to the
  engineering community as {\em intersecting boundary method}. It
  involves allowing the boundary to intersect a grid of cells thereby
  modifying the shape of those cells intersected. We suggested allowing
- the topgoraphy to take on a peice-wise linear representation (shaved
+ the topography to take on a piece-wise linear representation (shaved
- cells) or a simpler piecewise constant representaion (partial step).
+ cells) or a simpler piecewise constant representation (partial step).
  Both show dramatic improvements in solution compared to the
  traditional full step representation, the piece-wise linear being the
  best. However, the storage requirements are excessive so the simpler
-Line 413 
 the thickness of a level is determined a
+Line 439 
 the thickness of a level is determined a
  \marginpar{$h_s$: {\bf hFacS}}
  The physical thickness of a tracer cell is given by $h_c(i,j,k) \Delta
  r_f(k)$ and the physical thickness of the open side is given by
- $h_w(i,j,k) \Delta r_f(k)$. Three 3-D discriptors $h_c$, $h_w$ and
+ $h_w(i,j,k) \Delta r_f(k)$. Three 3-D descriptors $h_c$, $h_w$ and
  $h_s$ are used to describe the geometry: {\bf hFacC}, {\bf hFacW} and
  {\bf hFacS} respectively. These are calculated in subroutine {\em
  INI\_MASKS\_ETC} along with there reciprocals {\bf RECIP\_hFacC}, {\bf
-Line 450 
 $h_s^{-1}$: {\bf RECIP\_hFacS} ({\em GRI
+Line 476 
 $h_s^{-1}$: {\bf RECIP\_hFacS} ({\em GRI
  \section{Continuity and horizontal pressure gradient terms}
+ \begin{rawhtml}
+ <!-- CMIREDIR:continuity_and_horizontal_pressure: -->
+ \end{rawhtml}
  The core algorithm is based on the ``C grid'' discretization of the
  continuity equation which can be summarized as:
  \begin{eqnarray}
- \partial_t u + \frac{1}{\Delta x_c} \delta_i \left. \frac{ \partial \Phi}{\partial r}\right|_{s} \eta + \frac{\epsilon_{nh}}{\Delta x_c} \delta_i \Phi_{nh}' & = & G_u - \frac{1}{\Delta x_c} \delta_i \Phi_h' \\
+ \partial_t u + \frac{1}{\Delta x_c} \delta_i \left. \frac{ \partial \Phi}{\partial r}\right|_{s} \eta + \frac{\epsilon_{nh}}{\Delta x_c} \delta_i \Phi_{nh}' & = & G_u - \frac{1}{\Delta x_c} \delta_i \Phi_h' \label{eq:discrete-momu} \\
- \partial_t v + \frac{1}{\Delta y_c} \delta_j \left. \frac{ \partial \Phi}{\partial r}\right|_{s} \eta + \frac{\epsilon_{nh}}{\Delta y_c} \delta_j \Phi_{nh}' & = & G_v - \frac{1}{\Delta y_c} \delta_j \Phi_h' \\
+ \partial_t v + \frac{1}{\Delta y_c} \delta_j \left. \frac{ \partial \Phi}{\partial r}\right|_{s} \eta + \frac{\epsilon_{nh}}{\Delta y_c} \delta_j \Phi_{nh}' & = & G_v - \frac{1}{\Delta y_c} \delta_j \Phi_h' \label{eq:discrete-momv} \\
- \epsilon_{nh} \left( \partial_t w + \frac{1}{\Delta r_c} \delta_k \Phi_{nh}' \right) & = & \epsilon_{nh} G_w + \overline{b}^k - \frac{1}{\Delta r_c} \delta_k \Phi_{h}' \\
+ \epsilon_{nh} \left( \partial_t w + \frac{1}{\Delta r_c} \delta_k \Phi_{nh}' \right) & = & \epsilon_{nh} G_w + \overline{b}^k - \frac{1}{\Delta r_c} \delta_k \Phi_{h}' \label{eq:discrete-momw} \\
  \delta_i \Delta y_g \Delta r_f h_w u +
  \delta_j \Delta x_g \Delta r_f h_s v +
  \delta_k {\cal A}_c w & = & {\cal A}_c \delta_k (P-E)_{r=0}
-Line 464 
 continuity equation which can be summari
+Line 494 
 continuity equation which can be summari
  \end{eqnarray}
  where the continuity equation has been most naturally discretized by
  staggering the three components of velocity as shown in
- Fig.~\ref{fig-cgrid3d}. The grid lengths $\Delta x_c$ and $\Delta y_c$
+ Fig.~\ref{fig:cgrid3d}. The grid lengths $\Delta x_c$ and $\Delta y_c$
  are the lengths between tracer points (cell centers). The grid lengths
  $\Delta x_g$, $\Delta y_g$ are the grid lengths between cell
  corners. $\Delta r_f$ and $\Delta r_c$ are the distance (in units of
-Line 477 
 A}_c$.  The factors $h_w$ and $h_s$ are
+Line 507 
 A}_c$.  The factors $h_w$ and $h_s$ are
  \marginpar{$h_s$: {\bf hFacS}}
  The last equation, the discrete continuity equation, can be summed in
- the vertical to yeild the free-surface equation:
+ the vertical to yield the free-surface equation:
  \begin{equation}
- {\cal A}_c \partial_t \eta + \delta_i \sum_k \Delta y_g \Delta r_f h_w u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v =
+ {\cal A}_c \partial_t \eta + \delta_i \sum_k \Delta y_g \Delta r_f h_w
- {\cal A}_c(P-E)_{r=0}
+ u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v = {\cal
+ A}_c(P-E)_{r=0} \label{eq:discrete-freesurface}
  \end{equation}
  The source term $P-E$ on the rhs of continuity accounts for the local
  addition of volume due to excess precipitation and run-off over
  evaporation and only enters the top-level of the {\em ocean} model.
  \section{Hydrostatic balance}
+ \begin{rawhtml}
+ <!-- CMIREDIR:hydrostatic_balance: -->
+ \end{rawhtml}
  The vertical momentum equation has the hydrostatic or
  quasi-hydrostatic balance on the right hand side. This discretization
-Line 495 
 derived from the buoyancy equation exact
+Line 529 
 derived from the buoyancy equation exact
  from the pressure gradient terms when forming the kinetic energy
  equation.
- In the ocean, using z-ccordinates, the hydrostatic balance terms are
+ In the ocean, using z-coordinates, the hydrostatic balance terms are
  discretized:
  \begin{equation}
  \epsilon_{nh} \partial_t w
  + g \overline{\rho'}^k + \frac{1}{\Delta z} \delta_k \Phi_h' = \ldots
+ \label{eq:discrete_hydro_ocean}
  \end{equation}
  In the atmosphere, using p-coordinates, hydrostatic balance is
  discretized:
  \begin{equation}
  \overline{\theta'}^k + \frac{1}{\Delta \Pi} \delta_k \Phi_h' = 0
+ \label{eq:discrete_hydro_atmos}
  \end{equation}
  where $\Delta \Pi$ is the difference in Exner function between the
  pressure points. The non-hydrostatic equations are not available in
-Line 513 
 the atmosphere.
+Line 549 
 the atmosphere.
  The difference in approach between ocean and atmosphere occurs because
  of the direct use of the ideal gas equation in forming the potential
- energy conversion term $\alpha \omega$. The form of these consversion
+ energy conversion term $\alpha \omega$. The form of these conversion
- terms is discussed at length in \cite{Adcroft01}.
+ terms is discussed at length in \cite{adcroft:02}.
  Because of the different representation of hydrostatic balance between
  ocean and atmosphere there is no elegant way to represent both systems

 Legend:



Removed from v.1.5
 


changed lines


 
Added in v.1.19
 Legend:



Removed from v.1.5
 


changed lines


 
Added in v.1.19
-Removed from v.1.5
+Added in v.1.19

	ViewVC Help
Powered by ViewVC 1.1.22