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

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

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

revision 1.7 by adcroft, Wed Sep 26 21:59:33 2001 UTC revision 1.22 by heimbach, Tue Jan 22 20:50:41 2008 UTC
# Line 2  Line 2 
2  % $Name$  % $Name$
3    
4  \section{Spatial discretization of the dynamical equations}  \section{Spatial discretization of the dynamical equations}
5    \begin{rawhtml}
6    <!-- CMIREDIR:spatial_discretization_of_dyn_eq: -->
7    \end{rawhtml}
8    
9  Spatial discretization is carried out using the finite volume  Spatial discretization is carried out using the finite volume
10  method. This amounts to a grid-point method (namely second-order  method. This amounts to a grid-point method (namely second-order
11  centered finite difference) in the fluid interior but allows  centered finite difference) in the fluid interior but allows
12  boundaries to intersect a regular grid allowing a more accurate  boundaries to intersect a regular grid allowing a more accurate
13  representation of the position of the boundary. We treat the  representation of the position of the boundary. We treat the
14  horizontal and veritical directions as seperable and differently.  horizontal and vertical directions as separable and differently.
   
 \input{part2/notation}  
15    
16    
17  \subsection{The finite volume method: finite volumes versus finite difference}  \subsection{The finite volume method: finite volumes versus finite difference}
18    \begin{rawhtml}
19    <!-- CMIREDIR:finite_volume: -->
20    \end{rawhtml}
21    
22    
23    
24  The finite volume method is used to discretize the equations in  The finite volume method is used to discretize the equations in
25  space. The expression ``finite volume'' actually has two meanings; one  space. The expression ``finite volume'' actually has two meanings; one
26  is the method of cut or instecting boundaries (shaved or lopped cells  is the method of embedded or intersecting boundaries (shaved or lopped
27  in our terminology) and the other is non-linear interpolation methods  cells in our terminology) and the other is non-linear interpolation
28  that can deal with non-smooth solutions such as shocks (i.e. flux  methods that can deal with non-smooth solutions such as shocks
29  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
30  conservation laws to which the {\it weak solution} is a solution on  of the conservation laws to which the {\it weak solution} is a
31  each finite volume of (sub-domain). The weak solution can be  solution on each finite volume of (sub-domain). The weak solution can
32  constructed outof piece-wise constant elements or be  be constructed out of piece-wise constant elements or be
33  differentiable. The differentiable equations can not be satisfied by  differentiable. The differentiable equations can not be satisfied by
34  piece-wise constant functions.  piece-wise constant functions.
35    
# Line 37  the lengths $\Delta x_i$: Line 43  the lengths $\Delta x_i$:
43  \begin{displaymath}  \begin{displaymath}
44  \Delta x \partial_t \theta + \delta_i ( F ) = 0  \Delta x \partial_t \theta + \delta_i ( F ) = 0
45  \end{displaymath}  \end{displaymath}
46  is exact if $\theta(x)$ is peice-wise constant over the interval  is exact if $\theta(x)$ is piece-wise constant over the interval
47  $\Delta x_i$ or more generally if $\theta_i$ is defined as the average  $\Delta x_i$ or more generally if $\theta_i$ is defined as the average
48  over the interval $\Delta x_i$.  over the interval $\Delta x_i$.
49    
# Line 57  recovers the same ODE's resulting from f Line 63  recovers the same ODE's resulting from f
63  interior of a fluid. Differences arise at boundaries where a boundary  interior of a fluid. Differences arise at boundaries where a boundary
64  is not positioned on a regular or smoothly varying grid. This method  is not positioned on a regular or smoothly varying grid. This method
65  is used to represent the topography using lopped cell, see  is used to represent the topography using lopped cell, see
66  \cite{Adcroft98}. Subtle difference also appear in more than one  \cite{adcroft:97}. Subtle difference also appear in more than one
67  dimension away from boundaries. This happens because the each  dimension away from boundaries. This happens because the each
68  direction is discretized independantly in the finite difference method  direction is discretized independently in the finite difference method
69  while the integrating over finite volume implicitly treats all  while the integrating over finite volume implicitly treats all
70  directions simultaneously. Illustration of this is given in  directions simultaneously. Illustration of this is given in
71  \cite{Adcroft02}.  \cite{ac:02}.
72    
73  \subsection{C grid staggering of variables}  \subsection{C grid staggering of variables}
74    
# Line 79  equations. } Line 85  equations. }
85  The basic algorithm employed for stepping forward the momentum  The basic algorithm employed for stepping forward the momentum
86  equations is based on retaining non-divergence of the flow at all  equations is based on retaining non-divergence of the flow at all
87  times. This is most naturally done if the components of flow are  times. This is most naturally done if the components of flow are
88  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}.
89    
90  Fig. \ref{fig:cgrid3d} shows the components of flow ($u$,$v$,$w$)  Fig. \ref{fig:cgrid3d} shows the components of flow ($u$,$v$,$w$)
91  staggered in space such that the zonal component falls on the  staggered in space such that the zonal component falls on the
92  interface between continiuty cells in the zonal direction. Similarly  interface between continuity cells in the zonal direction. Similarly
93  for the meridional and vertical directions.  The continiuty cell is  for the meridional and vertical directions.  The continuity cell is
94  synonymous with tracer cells (they are one and the same).  synonymous with tracer cells (they are one and the same).
95    
96    
# Line 115  grid data: ({\em model/inc/GRID.h}) Line 121  grid data: ({\em model/inc/GRID.h})
121    
122    
123  \subsection{Horizontal grid}  \subsection{Horizontal grid}
124    \label{sec:spatial_discrete_horizontal_grid}
125    
126  \begin{figure}  \begin{figure}
127  \begin{center}  \begin{center}
# Line 132  grid lines indicate the tracer cell boun Line 139  grid lines indicate the tracer cell boun
139  grid for all panels. a) The area of a tracer cell, $A_c$, is bordered  grid for all panels. a) The area of a tracer cell, $A_c$, is bordered
140  by the lengths $\Delta x_g$ and $\Delta y_g$. b) The area of a  by the lengths $\Delta x_g$ and $\Delta y_g$. b) The area of a
141  vorticity cell, $A_\zeta$, is bordered by the lengths $\Delta x_c$ and  vorticity cell, $A_\zeta$, is bordered by the lengths $\Delta x_c$ and
142  $\Delta y_c$. c) The area of a u cell, $A_c$, is bordered by the  $\Delta y_c$. c) The area of a u cell, $A_w$, is bordered by the
143  lengths $\Delta x_v$ and $\Delta y_f$. d) The area of a v cell, $A_c$,  lengths $\Delta x_v$ and $\Delta y_f$. d) The area of a v cell, $A_s$,
144  is bordered by the lengths $\Delta x_f$ and $\Delta y_u$.}  is bordered by the lengths $\Delta x_f$ and $\Delta y_u$.}
145  \label{fig:hgrid}  \label{fig:hgrid}
146  \end{figure}  \end{figure}
147    
148  The model domain is decomposed into tiles and within each tile a  The model domain is decomposed into tiles and within each tile a
149  quasi-regular grid is used. A tile is the basic unit of domain  quasi-regular grid is used. A tile is the basic unit of domain
150  decomposition for parallelization but may be used whether parallized  decomposition for parallelization but may be used whether parallelized
151  or not; see section \ref{sect:tiles} for more details. Although the  or not; see section \ref{sect:domain_decomposition} for more details.
152  tiles may be patched together in an unstructured manner  Although the tiles may be patched together in an unstructured manner
153  (i.e. irregular or non-tessilating pattern), the interior of tiles is  (i.e. irregular or non-tessilating pattern), the interior of tiles is
154  a structered grid of quadrilateral cells. The horizontal coordinate  a structured grid of quadrilateral cells. The horizontal coordinate
155  system is orthogonal curvilinear meaning we can not necessarily treat  system is orthogonal curvilinear meaning we can not necessarily treat
156  the two horizontal directions as seperable. Instead, each cell in the  the two horizontal directions as separable. Instead, each cell in the
157  horizontal grid is described by the length of it's sides and it's  horizontal grid is described by the length of it's sides and it's
158  area.  area.
159    
160  The grid information is quite general and describes any of the  The grid information is quite general and describes any of the
161  available coordinates systems, cartesian, spherical-polar or  available coordinates systems, cartesian, spherical-polar or
162  curvilinear. All that is necessary to distinguish between the  curvilinear. All that is necessary to distinguish between the
163  coordinate systems is to initialize the grid data (discriptors)  coordinate systems is to initialize the grid data (descriptors)
164  appropriately.  appropriately.
165    
166  In the following, we refer to the orientation of quantities on the  In the following, we refer to the orientation of quantities on the
# Line 179  rAc(i,j)} and {\bf DYg(i,j)} positioned Line 186  rAc(i,j)} and {\bf DYg(i,j)} positioned
186  Fig.~\ref{fig:hgrid}b shows the vorticity cell. The length of the  Fig.~\ref{fig:hgrid}b shows the vorticity cell. The length of the
187  southern edge, $\Delta x_c$, western edge, $\Delta y_c$ and surface  southern edge, $\Delta x_c$, western edge, $\Delta y_c$ and surface
188  area, $A_\zeta$, presented in the vertical are stored in arrays {\bf  area, $A_\zeta$, presented in the vertical are stored in arrays {\bf
189  DXg}, {\bf DYg} and {\bf rAz}.  DXc}, {\bf DYc} and {\bf rAz}.
190  \marginpar{$A_\zeta$: {\bf rAz}}  \marginpar{$A_\zeta$: {\bf rAz}}
191  \marginpar{$\Delta x_c$: {\bf DXc}}  \marginpar{$\Delta x_c$: {\bf DXc}}
192  \marginpar{$\Delta y_c$: {\bf DYc}}  \marginpar{$\Delta y_c$: {\bf DYc}}
# Line 292  using\-Cartes\-ianGrid} in namelist {\em Line 299  using\-Cartes\-ianGrid} in namelist {\em
299  spacing can be set to uniform via scalars {\bf dXspacing} and {\bf  spacing can be set to uniform via scalars {\bf dXspacing} and {\bf
300  dYspacing} in namelist {\em PARM04} or to variable resolution by the  dYspacing} in namelist {\em PARM04} or to variable resolution by the
301  vectors {\bf DELX} and {\bf DELY}. Units are normally  vectors {\bf DELX} and {\bf DELY}. Units are normally
302  meters. Non-dimensional coordinates can be used by interpretting the  meters. Non-dimensional coordinates can be used by interpreting the
303  gravitational constant as the Rayleigh number.  gravitational constant as the Rayleigh number.
304    
305  \subsubsection{Spherical-polar coordinates}  \subsubsection{Spherical-polar coordinates}
# Line 346  The vertical grid is calculated in subro Line 353  The vertical grid is calculated in subro
353  INI\_VERTICAL\_GRID} and specified via the vector {\bf DELR} in  INI\_VERTICAL\_GRID} and specified via the vector {\bf DELR} in
354  namelist {\em PARM04}. The units of ``r'' are either meters or Pascals  namelist {\em PARM04}. The units of ``r'' are either meters or Pascals
355  depending on the isomorphism being used which in turn is dependent  depending on the isomorphism being used which in turn is dependent
356  only on the choise of equation of state.  only on the choice of equation of state.
357    
358  There are alternative namelist vectors {\bf DELZ} and {\bf DELP} which  There are alternative namelist vectors {\bf DELZ} and {\bf DELP} which
359  dictate whether z- or  dictate whether z- or
# Line 360  vertical grid descriptors are stored in Line 367  vertical grid descriptors are stored in
367    
368  The above grid (Fig.~\ref{fig:vgrid}a) is known as the cell centered  The above grid (Fig.~\ref{fig:vgrid}a) is known as the cell centered
369  approach because the tracer points are at cell centers; the cell  approach because the tracer points are at cell centers; the cell
370  centers are mid-way between the cell interfaces. An alternative, the  centers are mid-way between the cell interfaces.
371  vertex or interface centered approach, is shown in  This discretization is selected when the thickness of the
372    levels are provided ({\bf delR}, parameter file {\em data},
373    namelist {\em PARM04})
374    An alternative, the vertex or interface centered approach, is shown in
375  Fig.~\ref{fig:vgrid}b. Here, the interior interfaces are positioned  Fig.~\ref{fig:vgrid}b. Here, the interior interfaces are positioned
376  mid-way between the tracer nodes (no longer cell centers). This  mid-way between the tracer nodes (no longer cell centers). This
377  approach is formally more accurate for evaluation of hydrostatic  approach is formally more accurate for evaluation of hydrostatic
378  pressure and vertical advection but historically the cell centered  pressure and vertical advection but historically the cell centered
379  approach has been used. An alternative form of subroutine {\em  approach has been used. An alternative form of subroutine {\em
380  INI\_VERTICAL\_GRID} is used to select the interface centered approach  INI\_VERTICAL\_GRID} is used to select the interface centered approach
381  but no run time option is currently available.  This form requires to specify $Nr+1$ vertical distances {\bf delRc}
382    (parameter file {\em data}, namelist {\em PARM04}, e.g.
383    {\em verification/ideal\_2D\_oce/input/data})
384    corresponding to surface to center, $Nr-1$ center to center, and center to
385    bottom distances.
386    %but no run time option is currently available.
387    
388  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
389  {\em S/R INI\_VERTICAL\_GRID} ({\em  {\em S/R INI\_VERTICAL\_GRID} ({\em
# Line 386  $\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\e Line 401  $\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\e
401    
402    
403  \subsection{Topography: partially filled cells}  \subsection{Topography: partially filled cells}
404    \begin{rawhtml}
405    <!-- CMIREDIR:topo_partial_cells: -->
406    \end{rawhtml}
407    
408  \begin{figure}  \begin{figure}
409  \begin{center}  \begin{center}
# Line 399  thickness of the open side is given by $ Line 417  thickness of the open side is given by $
417  \label{fig:hfacs}  \label{fig:hfacs}
418  \end{figure}  \end{figure}
419    
420  \cite{Adcroft97} presented two alternatives to the step-wise finite  \cite{adcroft:97} presented two alternatives to the step-wise finite
421  difference representation of topography. The method is known to the  difference representation of topography. The method is known to the
422  engineering community as {\em intersecting boundary method}. It  engineering community as {\em intersecting boundary method}. It
423  involves allowing the boundary to intersect a grid of cells thereby  involves allowing the boundary to intersect a grid of cells thereby
424  modifying the shape of those cells intersected. We suggested allowing  modifying the shape of those cells intersected. We suggested allowing
425  the topgoraphy to take on a peice-wise linear representation (shaved  the topography to take on a piece-wise linear representation (shaved
426  cells) or a simpler piecewise constant representaion (partial step).  cells) or a simpler piecewise constant representation (partial step).
427  Both show dramatic improvements in solution compared to the  Both show dramatic improvements in solution compared to the
428  traditional full step representation, the piece-wise linear being the  traditional full step representation, the piece-wise linear being the
429  best. However, the storage requirements are excessive so the simpler  best. However, the storage requirements are excessive so the simpler
# Line 419  the thickness of a level is determined a Line 437  the thickness of a level is determined a
437  \marginpar{$h_s$: {\bf hFacS}}  \marginpar{$h_s$: {\bf hFacS}}
438  The physical thickness of a tracer cell is given by $h_c(i,j,k) \Delta  The physical thickness of a tracer cell is given by $h_c(i,j,k) \Delta
439  r_f(k)$ and the physical thickness of the open side is given by  r_f(k)$ and the physical thickness of the open side is given by
440  $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
441  $h_s$ are used to describe the geometry: {\bf hFacC}, {\bf hFacW} and  $h_s$ are used to describe the geometry: {\bf hFacC}, {\bf hFacW} and
442  {\bf hFacS} respectively. These are calculated in subroutine {\em  {\bf hFacS} respectively. These are calculated in subroutine {\em
443  INI\_MASKS\_ETC} along with there reciprocals {\bf RECIP\_hFacC}, {\bf  INI\_MASKS\_ETC} along with there reciprocals {\bf RECIP\_hFacC}, {\bf
# Line 456  $h_s^{-1}$: {\bf RECIP\_hFacS} ({\em GRI Line 474  $h_s^{-1}$: {\bf RECIP\_hFacS} ({\em GRI
474    
475    
476  \section{Continuity and horizontal pressure gradient terms}  \section{Continuity and horizontal pressure gradient terms}
477    \begin{rawhtml}
478    <!-- CMIREDIR:continuity_and_horizontal_pressure: -->
479    \end{rawhtml}
480    
481    
482  The core algorithm is based on the ``C grid'' discretization of the  The core algorithm is based on the ``C grid'' discretization of the
483  continuity equation which can be summarized as:  continuity equation which can be summarized as:
# Line 470  continuity equation which can be summari Line 492  continuity equation which can be summari
492  \end{eqnarray}  \end{eqnarray}
493  where the continuity equation has been most naturally discretized by  where the continuity equation has been most naturally discretized by
494  staggering the three components of velocity as shown in  staggering the three components of velocity as shown in
495  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$
496  are the lengths between tracer points (cell centers). The grid lengths  are the lengths between tracer points (cell centers). The grid lengths
497  $\Delta x_g$, $\Delta y_g$ are the grid lengths between cell  $\Delta x_g$, $\Delta y_g$ are the grid lengths between cell
498  corners. $\Delta r_f$ and $\Delta r_c$ are the distance (in units of  corners. $\Delta r_f$ and $\Delta r_c$ are the distance (in units of
# Line 483  A}_c$.  The factors $h_w$ and $h_s$ are Line 505  A}_c$.  The factors $h_w$ and $h_s$ are
505  \marginpar{$h_s$: {\bf hFacS}}  \marginpar{$h_s$: {\bf hFacS}}
506    
507  The last equation, the discrete continuity equation, can be summed in  The last equation, the discrete continuity equation, can be summed in
508  the vertical to yeild the free-surface equation:  the vertical to yield the free-surface equation:
509  \begin{equation}  \begin{equation}
510  {\cal A}_c \partial_t \eta + \delta_i \sum_k \Delta y_g \Delta r_f h_w  {\cal A}_c \partial_t \eta + \delta_i \sum_k \Delta y_g \Delta r_f h_w
511  u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v = {\cal  u + \delta_j \sum_k \Delta x_g \Delta r_f h_s v = {\cal
# Line 494  addition of volume due to excess precipi Line 516  addition of volume due to excess precipi
516  evaporation and only enters the top-level of the {\em ocean} model.  evaporation and only enters the top-level of the {\em ocean} model.
517    
518  \section{Hydrostatic balance}  \section{Hydrostatic balance}
519    \begin{rawhtml}
520    <!-- CMIREDIR:hydrostatic_balance: -->
521    \end{rawhtml}
522    
523  The vertical momentum equation has the hydrostatic or  The vertical momentum equation has the hydrostatic or
524  quasi-hydrostatic balance on the right hand side. This discretization  quasi-hydrostatic balance on the right hand side. This discretization
# Line 502  derived from the buoyancy equation exact Line 527  derived from the buoyancy equation exact
527  from the pressure gradient terms when forming the kinetic energy  from the pressure gradient terms when forming the kinetic energy
528  equation.  equation.
529    
530  In the ocean, using z-ccordinates, the hydrostatic balance terms are  In the ocean, using z-coordinates, the hydrostatic balance terms are
531  discretized:  discretized:
532  \begin{equation}  \begin{equation}
533  \epsilon_{nh} \partial_t w  \epsilon_{nh} \partial_t w
# Line 522  the atmosphere. Line 547  the atmosphere.
547    
548  The difference in approach between ocean and atmosphere occurs because  The difference in approach between ocean and atmosphere occurs because
549  of the direct use of the ideal gas equation in forming the potential  of the direct use of the ideal gas equation in forming the potential
550  energy conversion term $\alpha \omega$. The form of these consversion  energy conversion term $\alpha \omega$. The form of these conversion
551  terms is discussed at length in \cite{Adcroft01}.  terms is discussed at length in \cite{adcroft:02}.
552    
553  Because of the different representation of hydrostatic balance between  Because of the different representation of hydrostatic balance between
554  ocean and atmosphere there is no elegant way to represent both systems  ocean and atmosphere there is no elegant way to represent both systems
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22