Dear All,
I was well aware of this feature of SPEM/SCRUM/ROMS codes for years (at least 5 years) and do not consider this feature as a bug. Neither I consider the latest fix as an urgent measure to be done immediately, and I do not see any reason for panic or statements like
I cannot believe that we missed it for so many years.
Furthermore, implementing this fix alone (just replacing the
metrics.F, as proposed by Hernan)
without changing several places thoughout the rest of the code is dangerous because it may potentially destroy the conservation properties in the momentum equation and utimately threatens stability of the code (see below).
SPEM/SCRUM/ROMS codes (all of them) were indeed very sloppy in handling curvilinear grids. There are numerous examples where consecutive averaging of already averaged quantities takes place, resulting in potential loss of accuracy. For example, inside the grid generation software (this is
seagrid for those who remember) the metric coefficients
pm=1/dx, and
pn=1/dy are computed first at their naturally on the C-grid (at
U- and
V-points, respectively) and then averaged onto RHO-points: algebraic averaging of
pm,
pn, e.g.,
Code: Select all
pm_(i,j) = [ pm_(i,j) + pm_(i+1,j) ]/2
hence harmonic averaging of associated grid spacings
dx,
dy. These averaged values of
pm,
pn are then stored in the grid file, and whenever necessary are
averaged again, see for example routine
metrics, old or new, does not matter. Effectively, we are computing
dx(i+1/2,j) using 3-grid-interval stencil, rather than the natural 1-interval stencil.
Regarding to the bug fix proposed by Hernan, basically what we talking here is replacing harmonic averaging of
dx_ and
dy_ with algebraic averaging,
Code: Select all
dx_(i,j) + dx_(i+1,j)
dx_(i+1/2,j) = -----------------------
2
vs.
2
dx_(i+1/2,j) = ---------------------------
1/dx_(i,j) + 1/dx_(i+1,j)
These two agree with each other with second order of accuracy and the difference between them is asymptotically small with controlling smallness parameter
Code: Select all
dx_(i,j) - dx_(i+1,j)
epsilon = -----------------------
dx_(i,j) + dx_(i+1,j)
(
second order agreement means that the difference vanishes as second power of epsilon above).
From the accuracy point of view, neither formula has a decisive advantage over the other. Consider, for example, grid with exponential stretching:
where
j is grid index,
X is length scale, while
N controls stretching. Then
Code: Select all
dx_(j) = 2X * sinh(1/(2N)) * exp(j/N)
and, similarly, the
exact value of halfway midpoint
dx_(j+1/2)
Code: Select all
dx_(j+1/2) = 2X * sinh(1/(2N)) * exp((j+1/2)/N)
on the other hand, algebraic averaging of
dx_(j) and
dx_(j+1) yields
Code: Select all
dx_(j)+dx_(j+1)
----------------- = 2X * sinh(1/(2N)) * cosh(1/(2N)) * exp((j+1/2)/N)
2
which has an extra multiplier
cosh(1/(2N)); while harmonic averaging yields
Code: Select all
2 2X * sinh(1/(2N))
----------------------- = ------------------- * exp((j+1/2)/N)
1/dx_(j) + 1/dx_(j+1) cosh(1/(2N))
where the same term appears in the denominator. One cannot say that one approximation is better than the other, however one can verify that their mean produces a more accurate approximation than the other two:
Code: Select all
algebraic + harmonic [cosh(1/(2N))]^2 + 1
---------------------- = [dx_(j+1/2)]_"exact" * ----------------------
2 2 * cosh(1/(2N))
where the extra multiplier on the right behaves like
Code: Select all
[cosh(1/(2N))]^2 + 1 1
---------------------- = 1 + --- * (1/(2N))^4 + ...
2 * cosh(1/(2N)) 8
In this case, truncation errors of algebraic and harmonic averages are comparable in order, but have the opposite sign, so that they cancel each other.
One of the features of harmonic averaging is that it produces values which are systematically less than the algebraic averaging (the results are equal only if the numbers to be averaged are the same). This leads to the fact that the sum of all areas as seen from above of
U- and
V- control volumes is systematically less than the sum of
RHO-areas. In principle, this effect does not jeopardize the conservation properties of ROMS code, because we multiply by area and divide by exactly the same number, for example in
U-equation it is
Code: Select all
a multiplier of om_u(i,j)*om_u(i,j)
followed by multiplier
Code: Select all
0.25*(pm(i,j)+pm(i-1,j))*(pn(i,j)+pn(i-1,j))
See, for example,
step2d and tracer computation of
rhs_ubarand its subsequent final use to compute the updated value of
ubar(i,j,knew).
As long ad it is guaranteed that the product of the two multipliers above is equal to unity, the conservation properties are respected, regardless of the details how the multipliers are computed. The
old way of doing it is defining
Code: Select all
om_u(i,j)=2./(pm(i,j)+pm(i-1,j))
and
on_u(i,j)=2./(pn(i,j)+pn(i-1,j))
However conservation properties are lost, if, for example, one replaces the
metrics.F routine with the one proposed by Hernan and leaving
step2d unchanged.
Code: Select all
recall that 0.25*(pm(i,j)+pm(i-1,j))*(1./pm(i,j)+1./pm(i-1,j) > 1
with 1 achieved only if
pm(i,j)=pm(i-1,j). This kind of multiplier may show up if one implements
metrics.F as suggested by Hernan, while keeping the rest of the code as before. Depending on specific version of the code, it may be equivalent to multiplying u and ubar by this factor at every time step, hence instability. I cheched ROMS 2.1, and in that code actually only the right-hand-side are multiplied and divided by horizontal surfaces, while ubar and u are not; this is in contrast to the algorithm for tracer equations, where
t(:,:,:,:) are also multiplied and divided then during time stepping. So such instability would not occur on ROMS 2.1. The conservation properies, however, will be lost in any code.
Where do we go from here?
It is my vision that the correct way to address the problem is to reconsider the whole technological chain of grid generation and the format of storing of grid data in the files. The current strategy is to store a lot of redundant information, but at the same time allow irreversible loss of information, for example averaging of
pn and
pm to RHO-points inside grid-generation package.
Grid generation:
We should recall that SPEM/SCRUM/ROMS codes employ orthogonal curvilinear coordinates. This means that there is exists a conformal mapping procedure between physical coordinates
X,
Y and the corresponding logical coordinates
xi,
eta (basically these are grid-index coordinate, with, perhaps, only one scaling factor deciding aspect ratio on logical coordinate space.
Since this is conformal mapping,
X(xi,eta) and
Y(xi,eta) are harmonic functions of the logical coordinates,
Code: Select all
d^2 X d^2 X d^2 X d^2 X
-------- + -------- = 0 and -------- + -------- = 0
d xi^2 d eta^2 d xi^2 d eta^2
which are naturally discretized as follows:
Code: Select all
at RHO-points,
Xu_(i+1/2,j) + Xu_(i-1/2,j) - 2 Xr_(i,j)
-----------------------------------------
dxi^2
Xv_(i,j+1/2) + Xv_(i,j-1/2) - 2 Xr_(i,j)
+ ----------------------------------------- = 0
deta^2
where
dxi and
deta are globally defined constants (actually only ratio of
dxi/deta makes sense, because there is no right-hand-side, and the equations can be scaled by an arbitrary number).
Similarly,
Code: Select all
U-points
Xr_(i+1,j) + Xr_(i,j) - 2 Xu_(i+1/2,j)
---------------------------------------
dxi^2
Xp_(i+1/2,j+1/2) + Xp_(i+1/2,j-1/2) - 2 Xu_(i+1/2,j)
+ ------------------------------------------------------ = 0
deta^2
V-points
Xp_(i+1/2,j+1/2) + Xp_(i-1/2,j+1/2) - 2 Xv_(i,j+1/2)
-----------------------------------------------------
dxi^2
Xr_(i,j+1) + Xr_(i,j) - 2 Xv_(i,j+1/2)
+ ---------------------------------------- = 0
deta^2
PSI-points
Xv_(i+1,j+1/2) + Xv_(i,j+1/2) - 2 Xp_(i+1/2,j+1/2)
---------------------------------------------------
dxi^2
Xu_(i+1/2,j+1) + Xu_(i+1/2,j) - 2 Xp_(i+1/2,j+1/2)
+ ---------------------------------------------------- = 0
deta^2
The above notation places the four relevant variables on the grid as follows:
Code: Select all
Xu_(i-1/2,j+1) Xr_(i,j+1) Xu_(i+1/2,j+1)
Xp_(i-1/2,j+1/2) Xv_(i,j+1/2) Xp_(i+1/2,j+1/2)
Xu_(i-1/2,j) Xr_(i,j) Xu_(i+1/2,j)
Xp_(i-1/2,j-1/2) Xv_(i,j-1/2) Xp_(i+1/2,j-1/2)
Where basically one can define an array of twice the size in each direction and just solve for regular 5-point Laplacian. (Y is discretized in a similar way).
This leads to the natural interpolation rules for
U-,
V-, and
PSI-points, as well as subsequent computation of all other metric terms: suppose we store
X and both
RHO- and
PSI-points in the file. Then,
X at
U- and
V-points, can be reconstructed as
Code: Select all
[[Xr_(i,j)+Xr_(i+1,j)]/dxi^2 + [Xp_(i+1/2,j+1/2)+Xp_(i+1/2,j 1/2)]/deta^2
Xu(i+1/2,j) = ---------------------------------------------------------------------------
2/dxi^2 + 2/deta^2
similarly
[[Xp_(i+1/2,j+1/2)+Xp_(i-1/2,j+1/2)]/dxi^2 + [[Xr_(i,j)+Xr-(i,j+1)]/deta^2
Xv(i,j+1/2) = ---------------------------------------------------------------------------
2/dxi^2 + 2/deta^2
which are
exact in sense that it identically respects the harmonic property (zero-Laplacian) defined on a 5-point stencil specified above.
Once all four values,
Xr,
Xu,
Xv,
Xp are known, all relevant grid spacings can be computed by simple differencing without any averaging at all, for example,
Code: Select all
om_u(i+1/2,j)=Xr(i+1,j)-Xr(i,j)
pm_u=1./[Xr(i+1,j)-Xr(i,j)]
Basically it all translates into storing the four arrays, say
Xr,
Xp,
Yr,
Yp (alternatively
Xu,
Xv,
Yu,
Yv) in the grid file and ...that is it. Everything else (
pm,
pn,
Coriolis,
angles,
areas, etc... ) can be computed from them on the fly. In the case of spherical coordinates either
lon-
lat coordinates, or their Mercator projections can be stored.
Doing everything accurately is possible, and does not lead to extra computational burden inside the hydrodynamic code of ROMS, but once again, it is not just code: the whole supporting infrastructure (grid generation, data file formats, pre- post processing matlab scripts, etc) should adopt a different standard.
Alexander Shchepetkin
.
.
