Numerical Issues

[arango]

Use this thread to discuss any numerical issues associated with the documentation in the ROMS/TOMS manual.

[kate]

Many moons ago, I sent out an image about the ROMS tiling. That image is out of date, so here is an update:


Edit: I'll put this one out here too:

[kate]

Tiling where the numbers work out:

Tiling when there is slop, say Mm=34 and NtileJ=4:

[stef]

Hi!

I was wondering if someone could give me hints about where to find documentation of the 3d momentum time-stepping (Adams-Bashforth 3) used in version 2.2. I tried to understand the Shchepetkin-McWilliams paper, and could extract some info about the barotropic LF-AM3, but I think they don't cover the AB3, do they? Is there such a study about the AB3, directed to ocean modelling?

Thanks,
--Stefan

[shchepet]

There is no emphasis for description of the particular main 3D time
stepping procedure used in ROMS v. 1.8....2.2 in the Shchepetkin
McWilliams,2005 ROMS paper because this algorithm was always considered
as a provisional variant to be replaced with a more refined one.

The main 3D procedure of v. 1.8....2.2 can be classified as an AB3-TR
Generalized forward-backward scheme [AB3 = Adams--Bashforth 3rd order
(open parabolic rule); TR = trapezoidal rule) and is summarized as
follows:

Code: Select all

 u^{n+1} = u^n +dt*[(23/12)*RHSU^n - (4/3)*RHSU^{n-1} + (5/12)*RHSU^{n-2}]  

 T^{n+1} = T^n -dt*div[ (1/2)*(U^{n+1}+U^n)*T^{n+1/2} ]

where T^{n+1/2} = (5/12)*T^{n+1,*} + (2/3)*T^n - (1/12)*T^{n-1} (the
interpolation coefficients are recognized as the ones for Adams--Moulton
3rd-order closed parabolic integration rule). In its turn T^{n+1,*} comes
from LF predictor sub-step,

Code: Select all

       T^{n+1,*}=T^{n-1} -2*dt*div[ U^n * T^n ]

which is implemented in practice as a pseudo-compressible finite-volume
step with artificial continuity equation, and, The LF-step is also
combined with AM3 interpolation (this is done in pre_step.F), so that
you actually do not see a naked LF explicitly in the code, but rather
you see

Code: Select all

 T^{n+1/2} = (5/12)*[ T^{n-1} -2*dt*div(U^n * T^n)]
            + (2/3)*T^n -(1/12)*T^{n-1}

or

Code: Select all

  T^{n+1/2} = (1/3)*T^{n-1} + (2/3)*T^n -(5/6)*dt*div(U^n * T^n)

where (1/3) appear in "pre_step.F" as "1/2-Gamma"; (2/3) as "1/2+Gamma";
(5/6) as "1-Gamma"; and Gamma is set to 1/6.

Above for simplicity I have omitted multiplication and division by Hz
and metric factors. Lowercase "u" means velocity; uppercase "U" fluxes
(which map onto Huon and Hvom in the code.

Setting of fluxes for corrector sub-step of tracer equations via
Trapezoidal Rule (TR) occurs in "step3d_uv.F", see

Code: Select all

   Huon(i,j,k)=0.5_r8*(Huon(i,j,k)+u(i,j,k,nnew)*DC(i,k))

where u(i,j,k,nnew) is u^{n+1}, DC(i,k) is new-time step cross-section
(Hz^{n+1}/pn or pm for v-component), and Huon(i,j,k) is previously
computed flux "U^n" (in the actual code Huon(i,j,k) is initially computed
within "set_massflux.F" (in ROMS 2.2; I believe that older codes may have
different file name), and then modified in "step3d_uv.F".


Assuming that the stability limit is governed primarily by the phase speed
of the fastest internal mode (basically the first baroclinic mode), the
details of LF-AM3 stepping for tracer become unimportant because the rate
of change of T is dominated by vertical gradient of T multiplied by
vertical velocity. In this case the whole algorithm maps onto Eq. (2.49)
for Shchepetkin McWilliams, 2005 with beta=5/12, delta=1/2,
gamma=epsilon=0. Its stability limit alpha_max=1.1441551 and its
placement of characteristic roots relatively to unit circle looks similar
to one of Fig. 13, upper left panel (that panel shows AB3-AM3, rather
than AB3-TR version).

The exact figure showing AB3-TR characteristic roots can be seen on page 22 in

http://marine.rutgers.edu/po/Workshops/ ... petkin.pdf

which is presentation on Venice, 2004 workshop.

[stef]

Hello!
Thanks for the reply!

I am also searching for information on ROMS finite volume formulations on the c-grid. Specifically, I have troubles understanding how the variables on the boundary of the gridcells are reconstructed to do the integration of the fluxes. It seems like there are many ways to do the reconstruction, and after reading the code and seeing how it is done, I would like to understand why it's done that way. In case anybody is aware of papers on this topic (finite-volumes on a c-grid), it would help me a lot if you could just post the titles here.

Thanks very much,
--Stefan

[stef]

I just reread the scrum manual and there it just says "finite diffences", not finite volume. I looked at the 2d-momentum equation and if I did everything right, then for a rectilinear grid the finite-volume and centered finite difference formulation is the same thing. But for a curvilinear grid you can't call it finite difference method any more, can you? Because the four cell boundaries can have a different length.

Sorry, I know this might not be ineresting for most users, but I have to write something into my work, so I would appreciate if somebody could tell me whether the term "finite volume" or "finite differences" ist the more correct one.

Suggestions for literature are welcome!

Thanks,
--Stefan

[shchepet]

The difference between "finite differences" and "finite volume" is saddle, but significant.

The main thing is interpretation of the gridded data: "finite-volume" assumes that u(i,j,k) is AVERAGE of field "u" over control volume dV(i,j,k), while finite difference assumes that it is instantaneous value at the location x(i,j,k). This leads to different formulas for things like computing derivatives, although within the second-order accuracy the formulas are mostly agree.

The method of derivation is very different, however. Formally, to derive a finite-volume discrete equation you must integrate you equations of motion over a control volume and, whenever possibly, transform the integral into flux form. Then, technically, it boils down to a three stage procedure: (1) given set of grid-box averages find instantaneous values at grid-box interfaces; (2) compute fluxes (perhaps using nonlinear formulas); and (3) add all fluxes coming in/out each control volume to find change of the amount of quantity there.

Methodologically, operation (3) is always "exact".

Operation (2) is "exact" if flux formulas are linear, which is true only in simplest cases, but generally it requires some approximations to be applied. For example, formally speaking if you have velocity and tracer concentration fields at the interface, then you must multiply them and integrate the product over the interface to get flux. But instead you approximate it by a product of "mean" values of each field averaged over the interface individually, and multiply it by the area of the interface, if you are on a curvilinear grid and areas are changing from one grid element to another. These are not the same, but approximately close to.

Operation (1) is never exact, and is, if fact, root of most errors.


Consider two formulas

a(i+1/2) = -1/16 *a(i-1) + 9/16 *a(i) + 9/16 *a(i+1) -1/16 *a(i+1)

and

a(i+1/2) = -1/12 *a(i-1) + 7/12 *a(i) + 7/12 *a(i+1) -1/12 *a(i+1)

the first one interpolates the value of quantity "a" given at discrete locations x(i), x(i+1), ... to the midpoint x(i+1/2) half-way between x(i) and x(i+1). That is, given instantaneous values, compute instantaneous value.

The second formula assumes that a(i) is the average of the field within the interval x(i-1/2) < x < x(i+1/2), and similarly interprets, a(i+1), a(i+2), and a(i-1). Then it not just interpolates to mid-location, but it performs TRANSLATION from averaged to instantaneous values. As the result, coefficients are different.

What happens if you subtract [a(i+1/2) -a(i-1/2)]/dx ?

If a(i+1/2) is computed using bottom formula, you get the fourth-order accurate formula for approximation for the first derivative. If, on the other hand, you use upper formula, you would not get fourth order accuracy.

If one wants to guarantee things like integral conservation of something, and the grid is curvilinear, then "finite volume" is basically the only way to go.

A good reference to read is

Arakawa and Lamb, 1977 Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model, Meth. Comput. Phys., vol17, 174-267. Academic Press.


Common usage of terms "finite volume" and "finite difference" is loose and often interchangeable, especially if people are talking about structured grids. SCRUM is a finite volume model as well, and so is POM, but nobody emphasizes that.

In contrast, unstructured grid people, use the term "finite volume" as opposite to "finite element" and tend to make very big deal out of it.

Term "control volume method" is interchangeable with "finite volume".