/[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.2 by adcroft, Wed Aug 8 22:19:02 2001 UTC revision 1.3 by adcroft, Thu Aug 9 16:26:43 2001 UTC
# Line 38  grid data: ({\em model/inc/GRID.h}) Line 38  grid data: ({\em model/inc/GRID.h})
38    
39  The finite volume method is used to discretize the equations in  The finite volume method is used to discretize the equations in
40  space. The expression ``finite volume'' actually has two meanings; one  space. The expression ``finite volume'' actually has two meanings; one
41  involves invocation of the weak formulation (e.g. integral  is the method of cut or instecting boundaries (shaved or lopped cells
42  formulation) and the other involves non-linear expressions for  in our terminology) and the other is non-linear interpolation methods
43  interpolation of data (e.g. flux limiters). We use both but they can  that can deal with non-smooth solutions such as shocks (i.e. flux
44  and will be ddiscussed seperately. The finite volume method discretizes by invoking the weak formulation of the equations or integral form. For example, the 1-D advection-diffusion equation:  limiters for advection). Both make use of the integral form of the
45    conservation laws to which the {\it weak solution} is a solution on
46    each finite volume of (sub-domain). The weak solution can be
47    constructed outof piece-wise constant elements or be
48    differentiable. The differentiable equations can not be satisfied by
49    piece-wise constant functions.
50    
51    As an example, the 1-D constant coefficient advection-diffusion
52    equation:
53  \begin{displaymath}  \begin{displaymath}
54  \partial_t \theta + \partial_x ( u \theta - \kappa \partial_x \theta ) = 0  \partial_t \theta + \partial_x ( u \theta - \kappa \partial_x \theta ) = 0
55  \end{displaymath}  \end{displaymath}
56  can be discretized by integrating of finite lengths $\Delta x$:  can be discretized by integrating over finite sub-domains, i.e.
57    the lengths $\Delta x_i$:
58  \begin{displaymath}  \begin{displaymath}
59  \Delta x \partial_t \theta + \delta_i ( F ) = 0  \Delta x \partial_t \theta + \delta_i ( F ) = 0
60  \end{displaymath}  \end{displaymath}
61  is exact but where the flux  is exact if $\theta(x)$ is peice-wise constant over the interval
62    $\Delta x_i$ or more generally if $\theta_i$ is defined as the average
63    over the interval $\Delta x_i$.
64    
65    The flux, $F_{i-1/2}$, must be approximated:
66  \begin{displaymath}  \begin{displaymath}
67  F = u \overline{\theta} - \frac{\kappa}{\Delta x_c} \partial_i \theta  F = u \overline{\theta} - \frac{\kappa}{\Delta x_c} \partial_i \theta
68  \end{displaymath}  \end{displaymath}
69  is approximate. The method for obtained $\overline{\theta}$ is  and this is where truncation errors can enter the solution. The
70  unspecified and non-lienar finite volume methods can be invoked.  method for obtaining $\overline{\theta}$ is unspecified and a wide
71    range of possibilities exist including centered and upwind
72  ....  INCOMPLETE  interpolation, polynomial fits based on the the volume average
73    definitions of quantities and non-linear interpolation such as
74    flux-limiters.
75    
76    Choosing simple centered second-order interpolation and differencing
77    recovers the same ODE's resulting from finite differencing for the
78    interior of a fluid. Differences arise at boundaries where a boundary
79    is not positioned on a regular or smoothly varying grid. This method
80    is used to represent the topography using lopped cell, see
81    \cite{Adcroft98}. Subtle difference also appear in more than one
82    dimension away from boundaries. This happens because the each
83    direction is discretized independantly in the finite difference method
84    while the integrating over finite volume implicitly treats all
85    directions simultaneously. Illustration of this is given in
86    \cite{Adcroft02}.
87    
88  \subsection{C grid staggering of variables}  \subsection{C grid staggering of variables}
89    
# Line 122  curvilinear. All that is necessary to di Line 149  curvilinear. All that is necessary to di
149  coordinate systems is to initialize the grid data (discriptors)  coordinate systems is to initialize the grid data (discriptors)
150  appropriately.  appropriately.
151    
 \marginpar{Caution!}  
152  In the following, we refer to the orientation of quantities on the  In the following, we refer to the orientation of quantities on the
153  computational grid using geographic terminology such as points of the  computational grid using geographic terminology such as points of the
154  compass. This is purely for convenience but should note be confused  compass.
155    \marginpar{Caution!}
156    This is purely for convenience but should note be confused
157  with the actual geographic orientation of model quantities.  with the actual geographic orientation of model quantities.
158    
 \marginpar{$A_c$: {\bf rAc}}  
 \marginpar{$\Delta x_g$: {\bf DXg}}  
 \marginpar{$\Delta y_g$: {\bf DYg}}  
159  Fig.~\ref{fig:hgrid}a shows the tracer cell (synonymous with the  Fig.~\ref{fig:hgrid}a shows the tracer cell (synonymous with the
160  continuity cell). The length of the southern edge, $\Delta x_g$,  continuity cell). The length of the southern edge, $\Delta x_g$,
161  western edge, $\Delta y_g$ and surface area, $A_c$, presented in the  western edge, $\Delta y_g$ and surface area, $A_c$, presented in the
162  vertical are stored in arrays {\bf DXg}, {\bf DYg} and {\bf rAc}. The  vertical are stored in arrays {\bf DXg}, {\bf DYg} and {\bf rAc}.
163  ``g'' suffix indicates that the lengths are along the defining grid  \marginpar{$A_c$: {\bf rAc}}
164  boundaries. The ``c'' suffix associates the quantity with the cell  \marginpar{$\Delta x_g$: {\bf DXg}}
165  centers. The quantities are staggered in space and the indexing is  \marginpar{$\Delta y_g$: {\bf DYg}}
166  such that {\bf DXg(i,j)} is positioned to the south of {\bf rAc(i,j)}  The ``g'' suffix indicates that the lengths are along the defining
167  and {\bf DYg(i,j)} positioned to the west.  grid boundaries. The ``c'' suffix associates the quantity with the
168    cell centers. The quantities are staggered in space and the indexing
169    is such that {\bf DXg(i,j)} is positioned to the south of {\bf
170    rAc(i,j)} and {\bf DYg(i,j)} positioned to the west.
171    
 \marginpar{$A_\zeta$: {\bf rAz}}  
 \marginpar{$\Delta x_c$: {\bf DXc}}  
 \marginpar{$\Delta y_c$: {\bf DYc}}  
172  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
173  southern edge, $\Delta x_c$, western edge, $\Delta y_c$ and surface  southern edge, $\Delta x_c$, western edge, $\Delta y_c$ and surface
174  area, $A_\zeta$, presented in the vertical are stored in arrays {\bf  area, $A_\zeta$, presented in the vertical are stored in arrays {\bf
175  DXg}, {\bf DYg} and {\bf rAz}.  The ``z'' suffix indicates that the  DXg}, {\bf DYg} and {\bf rAz}.
176  lengths are measured between the cell centers and the ``$\zeta$'' suffix  \marginpar{$A_\zeta$: {\bf rAz}}
177  associates points with the vorticity points. The quantities are  \marginpar{$\Delta x_c$: {\bf DXc}}
178  staggered in space and the indexing is such that {\bf DXc(i,j)} is  \marginpar{$\Delta y_c$: {\bf DYc}}
179  positioned to the north of {\bf rAc(i,j)} and {\bf DYc(i,j)} positioned  The ``z'' suffix indicates that the lengths are measured between the
180  to the east.  cell centers and the ``$\zeta$'' suffix associates points with the
181    vorticity points. The quantities are staggered in space and the
182    indexing is such that {\bf DXc(i,j)} is positioned to the north of
183    {\bf rAc(i,j)} and {\bf DYc(i,j)} positioned to the east.
184    
 \marginpar{$A_w$: {\bf rAw}}  
 \marginpar{$\Delta x_v$: {\bf DXv}}  
 \marginpar{$\Delta y_f$: {\bf DYf}}  
185  Fig.~\ref{fig:hgrid}c shows the ``u'' or western (w) cell. The length of  Fig.~\ref{fig:hgrid}c shows the ``u'' or western (w) cell. The length of
186  the southern edge, $\Delta x_v$, eastern edge, $\Delta y_f$ and  the southern edge, $\Delta x_v$, eastern edge, $\Delta y_f$ and
187  surface area, $A_w$, presented in the vertical are stored in arrays  surface area, $A_w$, presented in the vertical are stored in arrays
188  {\bf DXv}, {\bf DYf} and {\bf rAw}.  The ``v'' suffix indicates that  {\bf DXv}, {\bf DYf} and {\bf rAw}.
189  the length is measured between the v-points, the ``f'' suffix  \marginpar{$A_w$: {\bf rAw}}
190  indicates that the length is measured between the (tracer) cell faces  \marginpar{$\Delta x_v$: {\bf DXv}}
191  and the ``w'' suffix associates points with the u-points (w stands for  \marginpar{$\Delta y_f$: {\bf DYf}}
192  west). The quantities are staggered in space and the indexing is such  The ``v'' suffix indicates that the length is measured between the
193  that {\bf DXv(i,j)} is positioned to the south of {\bf rAw(i,j)} and  v-points, the ``f'' suffix indicates that the length is measured
194  {\bf DYf(i,j)} positioned to the east.  between the (tracer) cell faces and the ``w'' suffix associates points
195    with the u-points (w stands for west). The quantities are staggered in
196    space and the indexing is such that {\bf DXv(i,j)} is positioned to
197    the south of {\bf rAw(i,j)} and {\bf DYf(i,j)} positioned to the east.
198    
 \marginpar{$A_s$: {\bf rAs}}  
 \marginpar{$\Delta x_f$: {\bf DXf}}  
 \marginpar{$\Delta y_u$: {\bf DYu}}  
199  Fig.~\ref{fig:hgrid}d shows the ``v'' or southern (s) cell. The length of  Fig.~\ref{fig:hgrid}d shows the ``v'' or southern (s) cell. The length of
200  the northern edge, $\Delta x_f$, western edge, $\Delta y_u$ and  the northern edge, $\Delta x_f$, western edge, $\Delta y_u$ and
201  surface area, $A_s$, presented in the vertical are stored in arrays  surface area, $A_s$, presented in the vertical are stored in arrays
202  {\bf DXf}, {\bf DYu} and {\bf rAs}.  The ``u'' suffix indicates that  {\bf DXf}, {\bf DYu} and {\bf rAs}.
203  the length is measured between the u-points, the ``f'' suffix  \marginpar{$A_s$: {\bf rAs}}
204  indicates that the length is measured between the (tracer) cell faces  \marginpar{$\Delta x_f$: {\bf DXf}}
205  and the ``s'' suffix associates points with the v-points (s stands for  \marginpar{$\Delta y_u$: {\bf DYu}}
206  south). The quantities are staggered in space and the indexing is such  The ``u'' suffix indicates that the length is measured between the
207  that {\bf DXf(i,j)} is positioned to the north of {\bf rAs(i,j)} and  u-points, the ``f'' suffix indicates that the length is measured
208  {\bf DYu(i,j)} positioned to the west.  between the (tracer) cell faces and the ``s'' suffix associates points
209    with the v-points (s stands for south). The quantities are staggered
210    in space and the indexing is such that {\bf DXf(i,j)} is positioned to
211    the north of {\bf rAs(i,j)} and {\bf DYu(i,j)} positioned to the west.
212    
213  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
214  {\em S/R INI\_CARTESIAN\_GRID} ({\em  {\em S/R INI\_CARTESIAN\_GRID} ({\em
# Line 347  $\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\e Line 375  $\Delta r_c^{-1}$: {\bf RECIP\_DRc} ({\e
375    
376  \end{minipage} }  \end{minipage} }
377    
378    
379    \subsection{Topography: partially filled cells}
380    
381    \begin{figure}
382    \centerline{
383    \resizebox{4.5in}{!}{\includegraphics{part2/vgrid-xz.eps}}
384    }
385    \caption{
386    A schematic of the x-r plane showing the location of the
387    non-dimensional fractions $h_c$ and $h_w$. The physical thickness of a
388    tracer cell is given by $h_c(i,j,k) \Delta r_f(k)$ and the physical
389    thickness of the open side is given by $h_w(i,j,k) \Delta r_f(k)$.}
390    \label{fig:hfacs}
391    \end{figure}
392    
393    \cite{Adcroft97} presented two alternatives to the step-wise finite
394    difference representation of topography. The method is known to the
395    engineering community as {\em intersecting boundary method}. It
396    involves allowing the boundary to intersect a grid of cells thereby
397    modifying the shape of those cells intersected. We suggested allowing
398    the topgoraphy to take on a peice-wise linear representation (shaved
399    cells) or a simpler piecewise constant representaion (partial step).
400    Both show dramatic improvements in solution compared to the
401    traditional full step representation, the piece-wise linear being the
402    best. However, the storage requirements are excessive so the simpler
403    piece-wise constant or partial-step method is all that is currently
404    supported.
405    
406    Fig.~\ref{fig:hfacs} shows a schematic of the x-r plane indicating how
407    the thickness of a level is determined at tracer and u points.
408    \marginpar{$h_c$: {\bf hFacC}}
409    \marginpar{$h_w$: {\bf hFacW}}
410    \marginpar{$h_s$: {\bf hFacS}}
411    The physical thickness of a tracer cell is given by $h_c(i,j,k) \Delta
412    r_f(k)$ and the physical thickness of the open side is given by
413    $h_w(i,j,k) \Delta r_f(k)$. Three 3-D discriptors $h_c$, $h_w$ and
414    $h_s$ are used to describe the geometry: {\bf hFacC}, {\bf hFacW} and
415    {\bf hFacS} respectively. These are calculated in subroutine {\em
416    INI\_MASKS\_ETC} along with there reciprocals {\bf RECIP\_hFacC}, {\bf
417    RECIP\_hFacW} and {\bf RECIP\_hFacS}.
418    
419    The non-dimensional fractions (or h-facs as we call them) are
420    calculated from the model depth array and then processed to avoid tiny
421    volumes. The rule is that if a fraction is less than {\bf hFacMin}
422    then it is rounded to the nearer of $0$ or {\bf hFacMin} or if the
423    physical thickness is less than {\bf hFacMinDr} then it is similarly
424    rounded. The larger of the two methods is used when there is a
425    conflict. By setting {\bf hFacMinDr} equal to or larger than the
426    thinnest nominal layers, $\min{(\Delta z_f)}$, but setting {\bf
427    hFacMin} to some small fraction then the model will only lop thick
428    layers but retain stability based on the thinnest unlopped thickness;
429    $\min{(\Delta z_f,\mbox{\bf hFacMinDr})}$.
430    
431    \fbox{ \begin{minipage}{4.75in}
432    {\em S/R INI\_MASKS\_ETC} ({\em model/src/ini\_masks\_etc.F})
433    
434    $h_c$: {\bf hFacC} ({\em GRID.h})
435    
436    $h_w$: {\bf hFacW} ({\em GRID.h})
437    
438    $h_s$: {\bf hFacS} ({\em GRID.h})
439    
440    $h_c^{-1}$: {\bf RECIP\_hFacC} ({\em GRID.h})
441    
442    $h_w^{-1}$: {\bf RECIP\_hFacW} ({\em GRID.h})
443    
444    $h_s^{-1}$: {\bf RECIP\_hFacS} ({\em GRID.h})
445    
446    \end{minipage} }
447    
448    
449  \subsection{Continuity and horizontal pressure gradient terms}  \subsection{Continuity and horizontal pressure gradient terms}
450    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22