Ocean Modeling Discussion

ROMS/TOMS

Search for:
It is currently Thu Sep 21, 2017 3:21 am




Post new topic Reply to topic  [ 7 posts ] 

All times are UTC

Author Message
PostPosted: Thu Mar 13, 2014 3:41 am 
Offline

Joined: Mon Oct 17, 2011 1:41 am
Posts: 12
Location: WHOI
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:
              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:
              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:
              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:
             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:
                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


Top
 Profile  
Reply with quote  
PostPosted: Thu Mar 20, 2014 1:08 am 
Offline

Joined: Mon Oct 17, 2011 1:41 am
Posts: 12
Location: WHOI
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


Top
 Profile  
Reply with quote  
PostPosted: Fri Mar 21, 2014 3:05 am 
Offline

Joined: Fri Nov 14, 2003 4:57 pm
Posts: 170
Location: UCLA, USA
The original code in u2dbc_im.F is as follows, and it is correct:
Code:
          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.


Quote:
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.


Top
 Profile  
Reply with quote  
PostPosted: Wed Mar 26, 2014 3:44 am 
Offline

Joined: Mon Oct 17, 2011 1:41 am
Posts: 12
Location: WHOI
Thanks!

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

Deepak


Top
 Profile  
Reply with quote  
PostPosted: Wed Mar 26, 2014 8:08 pm 
Offline

Joined: Fri Nov 14, 2003 4:57 pm
Posts: 170
Location: UCLA, USA
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:
               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:
              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.


Top
 Profile  
Reply with quote  
PostPosted: Fri Mar 28, 2014 10:03 pm 
Offline
Site Admin
User avatar

Joined: Wed Feb 26, 2003 4:41 pm
Posts: 1011
Location: IMCS, Rutgers University
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.


Top
 Profile  
Reply with quote  
PostPosted: Mon Mar 31, 2014 6:16 pm 
Offline
User avatar

Joined: Wed Jul 02, 2003 5:29 pm
Posts: 3216
Location: IMS/UAF, USA
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.


Top
 Profile  
Reply with quote  
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 7 posts ] 

All times are UTC


Who is online

Users browsing this forum: xiaozhu557 and 4 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
Powered by phpBB® Forum Software © phpBB Group