The original code in u2dbc_im.F is as follows, and it is correct:
Code: Select all
if (cx > 0.292893218813452) then
zx=zx + ( zeta(iend,j,knew) +cx*zeta(iend+1,j,kstp)
# ifdef M2_FRC_BRY
& +ubar_east(j) +hx*(zx-zeta_east(j))
& +ubclm(iend+1,j) +hx*(zx-ssh(iend+1,j))
you may use it as the reference.
It also comes with a special test problem WAVE_RAD wich can be used to verify
that all 4 boundaries work correctly and result in very little reflection.
, that is division by cff, not multiplication, however
contains an extra multiplication
by cff, hence cx=dtfast*cff*sqrt(g/cff)*.... = dtfast*sqrt(g*cff)*...
, so it contains
sqrt(g*cff)=sqrt(g*h) which is phase speed of barotropic gravity wave. The resultant
overall value of cx is nondimensional, as it should be, because it is Courant number.
3. The above 1. and 2. are constructed in such a way that there is only one division and
one sqrt operation (both division and sqrt are approximately of the theme computational
cost, and are substantially more expensive than multiplicationand addition. This is the
rationale of doing it this way and not the other way around). Note that BOTH hx and cx
are needed by the computations below, and any other variant of their computation
[besides computing sqrt(g/h)] will result in an extra sqrt or an extra division.
4. if (cx > 0.292893218813452) then
maps onto cx > Co in Rutgers code, and zx is
modified inside this logical branch in such a way that the additional term is quadratic
with thespect to cx-Co.
5. It can easily be worked with paper and pencil, that after modification 4. to zx
is applied, zx is a weighted sum of zeta(iend,j,kstp), zeta(iend+1,j,kstp), and
zeta(iend,j,knew) with POSITIVE coefficients in front of all three such that the
three coefficients ADD UP TO 1.
I'm not sure if correcting cff to use h+zeta instead of just h is required for good accuracy.
Generally speaking the answer is NO, but this judgment is close to neutral than categoric.
These kind boundary conditions (I hate to call them "Flather" because calling them this
way is historically incorrect) are derived under the assumptions of
(i) linearization (left- and right-travelling waves are moving independently from each
other); thought the linearization does not necessarily have to be around zeta=0, it can
be viewed as a small perturbation relative to a final background state (e.g., consider
small short waves on top of large-amplitude tidal elevations and depressions of free
(ii) The waves are assumed to propagate in the direction normal to the boundary
(which is never true in practice, hence the physical accuracy of these B.C. is
(iii) The Coriolis term is neglected altogether, so the actual wave speed of surface
gravity wave is slightly larger than sqrt(gh), however there is no evident way to
fix this, because in the presense of Coriolis force waves become dispersive and
these is no self-evident way to decompose them into left and right travelling waves;
The limitations (ii) anf (iii) are offset by the fact that the role of these B.C.
condition is to merely provide "radiative cooling" of the solution, or, simply put,
avoid trapping of surface gravity waves inside the domain. The easiest way to
understand this is to compare this with "clamped" boundares. If normal velocity is
specified at the boundary, then waves will be 100% reflected back into the domain,
and the entire solution will have seiche-mode oscillation in free surface which will
never go away. In contrast, a radiation condition (not necessarily perfect
non-reflection) makes seiche mode go away rather quickly -- it cannot sustain itself
as lon as a significant portion of wave is let out of the domain every thime.