2 possible bugs in Shc boundary conditions

[dcherian]

I think I found 2 bugs in the code for the new Shchepetkin boundary conditions. In u2dbc_im.F for the eastern boundary, the current code has

Code: Select all

              cff=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend  ,j)+                &                                      
     &                            GRID(ng)%h(Iend+1,j)))
              cff1=SQRT(g*cff)
              Cx=dt2d*cff1*0.5_r8*(GRID(ng)%pm(Iend  ,j)+              &                                      
     &                             GRID(ng)%pm(Iend+1,j))

cff1 is sqrt(g/h) but Cx (which is c_tilde in Mason et al. (2010) notation) needs sqrt(g*h). However, in the last line we do require sqrt(g/h) when calculating ubar - corresponding line is

Code: Select all

              ubar(Iend+1,j,kout)=0.5_r8*                              &                                      
     &                            ((1.0_r8-Cx)*ubar(Iend+1,j,know)+    &                                      
     &                             Cx*ubar(Iend,j,know)+               &                                      
     &                             bry_val+                            &                                      
     &                             cff1*(Zx-BOUNDARY(ng)%zeta_east(j)))

So the definition of Cx needs to change to

Code: Select all

              Cx=dt2d*SQRT(g*1/cff)*0.5_r8*(GRID(ng)%pm(Iend  ,j)+        &                                      
     &                             GRID(ng)%pm(Iend+1,j))

I'm not sure if correcting cff to use h+zeta instead of just h is required for good accuracy.

The second bug is for the case Cx > Co. In the code, Zx is reassigned the value of Zx + cff2*cff3 when this condition is true.

Code: Select all

             Zx=(0.5_r8+Cx)*zeta(Iend  ,j,know)+                       &                                      
     &           (0.5_r8-Cx)*zeta(Iend+1,j,know)
              IF (Cx.gt.Co) THEN
                cff2=(1.0_r8-Co/Cx)**2
                cff3=zeta(Iend,j,kout)+                                &                                      
     &               Cx*zeta(Iend+1,j,know)-                           &                                      
     &               (1.0_r8-Cx)*zeta(Iend,j,know)
                Zx=Zx+cff2*cff3
              END IF

But, Zx+cff2*cff3 gives (0.5 + Cx - (1-Cx)*(1-Co/Cx)^2) * zeta(Iend,j,know) + two other terms. The other two terms have the correct coefficient when comparing against eqn (7) in Mason et al. (2010). But for zeta(Iend,j,know) we need (0.5 + Cx - (1+Cx)*(1-Co/Cx)^2) * zeta(Iend,j,know) for the code to agree with eqn (7). So, the required correction is

Code: Select all

                cff3=zeta(Iend,j,kout)+                                &                                      
     &               Cx*zeta(Iend+1,j,know)-                           &                                      
     &               (1.0_r8+Cx)*zeta(Iend,j,know)

Similar corrections are required for all 4 boundaries.

Deepak

Mason, E., Molemaker, J., Shchepetkin, A. F., Colas, F., McWilliams, J. C., & Sangrà, P. (2010). Procedures for offline grid nesting in regional ocean models. Ocean Modelling, 35(1-2), 1–15. doi:10.1016/j.ocemod.2010.05.007

[dcherian]

I can now confirm that my application behaves a lot better when the two changes are made. Earlier, I had a bunch of waves running around my domain. Now, they exit the open boundaries very nicely.

Deepak

[shchepet]

The original code in u2dbc_im.F is as follows, and it is correct:

Code: Select all

          cff=0.5*(h(iend,j)+h(iend+1,j))
          hx=sqrt(g/cff)
          cx=dtfast*cff*hx*0.5*(pm(iend,j)+pm(iend+1,j))

          zx=(0.5+cx)*zeta(iend,j,kstp)+(0.5-cx)*zeta(iend+1,j,kstp)
          if (cx > 0.292893218813452) then
            zx=zx + ( zeta(iend,j,knew) +cx*zeta(iend+1,j,kstp)
     &                               -(1.+cx)*zeta(iend,j,kstp)
     &                           )*(1.-0.292893218813452/cx)**2
          endif

          ubar(iend+1,j,knew)=0.5*( (1.-cx)*ubar(iend+1,j,kstp)
     &                                     +cx*ubar(iend,j,kstp)
#    ifdef M2_FRC_BRY
     &                       +ubar_east(j) +hx*(zx-zeta_east(j))
#    else
     &                   +ubclm(iend+1,j) +hx*(zx-ssh(iend+1,j))
#    endif
     &                                                        )

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.

Notes:

1. hx=sqrt(g/cff), that is division by cff, not multiplication, however

2. cx=dtfast*cff*hx*0.5*(pm(iend,j)+pm(iend+1,j)) 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
surface).

(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
fundamentally limited);

(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.

[dcherian]

Thanks!

The changes I made agree with the code you posted. I like the "radiation cooling" analogy.

Deepak

[shchepet]

OK. Let's just summarize it. I just downloaded the latest version, ROMS 3.7 svn rev. 726, and the relevant piece of
code looks like

Code: Select all

               cff=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend  ,j)+                 &   ! <-- A
     &                            GRID(ng)%h(Iend+1,j)))
              cff1=SQRT(g*cff)                                                ! <-- B
              Cx=dt2d*cff1*0.5_r8*(GRID(ng)%pm(Iend  ,j)+               &    ! <-- C
     &                             GRID(ng)%pm(Iend+1,j))
              Zx=(0.5_r8+Cx)*zeta(Iend  ,j,know)+                       &
     &           (0.5_r8-Cx)*zeta(Iend+1,j,know)
              IF (Cx.gt.Co) THEN
                cff2=(1.0_r8-Co/Cx)**2
                cff3=zeta(Iend,j,kout)+                                 &
     &               Cx*zeta(Iend+1,j,know)-                            &
     &               (1.0_r8-Cx)*zeta(Iend,j,know)                         ! <-- D
                Zx=Zx+cff2*cff3
              END IF
              ubar(Iend+1,j,kout)=0.5_r8*                               &
     &                            ((1.0_r8-Cx)*ubar(Iend+1,j,know)+     &
     &                             Cx*ubar(Iend,j,know)+                &
     &                             bry_val+                             &
     &                             cff1*(Zx-BOUNDARY(ng)%zeta_east(j)))   ! <-- E

where by A,B,C,D,E I highlighted the lines which need to be paid attention to.

There are two errors in the code above: (1) Cx computed by line C is incorrect, and also dimensionally incorrect.
It must be non-dimensional, while it is 1/Length if computed as above; and (2) the three weighting coefficients
for the three free-surface values in computation of cff3 (ending at line D) should add up to 0. Currently they
are not.

The corrected version is

Code: Select all

              cff=0.5_r8*(GRID(ng)%h(Iend  ,j)+                        &   ! <-- A
     &                    GRID(ng)%h(Iend+1,j))
              cff1=SQRT(g/cff)                                                  ! <-- B
              Cx=dt2d*cff1*cff* 0.5_r8*(GRID(ng)%pm(Iend  ,j)+         &      ! <-- C
     &                                   GRID(ng)%pm(Iend+1,j))
              Zx=(0.5_r8+Cx)*zeta(Iend  ,j,know)+                       &
     &           (0.5_r8-Cx)*zeta(Iend+1,j,know)
              IF (Cx.gt.Co) THEN
                cff2=(1.0_r8-Co/Cx)**2
                cff3=zeta(Iend,j,kout)+                                 &
     &               Cx*zeta(Iend+1,j,know)-                            &
     &               (1.0_r8+Cx)*zeta(Iend,j,know)                         ! <-- D
                Zx=Zx+cff2*cff3
              END IF
              ubar(Iend+1,j,kout)=0.5_r8*                               &
     &                            ((1.0_r8-Cx)*ubar(Iend+1,j,know)+     &
     &                             Cx*ubar(Iend,j,know)+                &
     &                             bry_val+                             &
     &                             cff1*(Zx-BOUNDARY(ng)%zeta_east(j)))  ! <-- E

A. changed 1.0_r8/(0.5_r8*(...)) into just 0.5_r8*(...)

B. changed SQRT(g*cff) into SQRT(g/cff)

C. introduced an extra multiplication by cff just after cff1. Note: changes A,B,C applied together
do not actually change the computed value of cff1, however Cx is different.

D. changed (1.0_r8-Cx) into (1.0_r8+Cx) Note: this is unrelated to A,B,C,E. The three weighting
coefficients for the three free-surface values should add up to 0.

E. there is no need for any change, however note that the meaning of cff1 must be preserved in changes A,B.


Similar corrections need to be applied to western edge boundary conditions in u2dbc_im.F and in southern
and northern edges in v2dbc_im.F

AGRIF and UCLA codes are not affected.

[arango]

Indeed, I fixed those bugs. Please update. Thank you for bringing this to my attention. I also corrected the tangent linear, representer, and adjoint versions of u2dbc_im.F and v2dbc_im.F.

[kate]

The WET_DRY case still needs cff=(h+zeta) for cases in which h goes negative. The following sqrt operation behaves badly if cff is negative.