Computing Advective CFL in diag.F

General scientific issues regarding ROMS

Moderators: arango, robertson

Post Reply
Message
Author
User avatar
shchepet
Posts: 185
Joined: Fri Nov 14, 2003 4:57 pm

Computing Advective CFL in diag.F

#1 Post by shchepet » Wed Jul 06, 2011 9:46 pm

Hernan, and all whoever may find it useful,

Some time ago I found that Computation of Courant number for the purpose of
diagnostics in done by averaging velocities toward the tracer point (see around
line 220 of diag.F),

Code: Select all

              my_Cu=0.5_r8*ABS(u(i,j,k,nstp)+u(i+1,j,k,nstp))*          &
     &              dt(ng)*pm(i,j)
              my_Cv=0.5_r8*ABS(v(i,j,k,nstp)+v(i,j+1,k,nstp))*          &
     &              dt(ng)*pn(i,j)
              my_Cw=0.5_r8*ABS(wvel(i,j,k-1)+wvel(i,j,k))*              &
     &              dt(ng)/Hz(i,j,k)
              my_C=my_Cu+my_Cv+my_Cw
              IF (my_C.gt.my_max_C) THEN
                my_max_C =my_C
                my_max_Cu=my_Cu
                my_max_Cv=my_Cv
                my_max_Cw=my_Cw
                my_max_Ci=i
                my_max_Cj=j
                my_max_Ck=k
              END IF
this is approximately correct, but there is a better way to define advective
Courant number consistently with finite-volume interpretation of grid-point
discrete data.

In the code segment below it advective Courant number is defined as the sum
of fluxes directed OUTWARD from grid box Hz divided by its control volume.
Essentially it is the fraction of water replaced within the grid box "Hz" during
one time step.


Why OUTWARD? Because under this definition, Cu=1 is the condition when
flux-split first-order upstream advection scheme looses its positive definiteness
property (hence stability). Of course, other combinations of spatial discrete
advection schemes and time stepping algorithm should become unstable when
Cu exceeds a different critical value, but it is also comparable to 1.

Note that in the case of "Hz"s fixed in time it does not matter whether it is
incoming or outgoing fluxes, because of the sum of incoming is equal to the
sum of outgoing. If "Hz"s are allowed to change, then it matters and it should
be the sum of outgoing. In principle, this definition is applicable to a code
which allows "drying" cells (i.e., one cannot take out more volume from an
about-to-dry cell that this cell has at the beginning of the time step), or,
in fact, any weird/unstructured grids with not necessarily six facets for
each grid box -- it can be any.

This also solves the dilemma which numbers to report: separate in each
direction, Cx,Cy,Cz or some combination of them, Cx+Cy+Cz. As the practical
utility of computing advective Courant number is to diagnose "hot spots" of
ROMS grids, which cause model blow up, and identify the primary reason -- either
exceeding horizontal or vertical velocity (say due to combination of a large
topographic slope and very fine local vertical resolution) , so it makes sense
to separate the two, but reporting two horizontal components is hardly useful:
ROMS grid resolutions are more or less horizontally isotropic most of the time.

Reported values are: "Cu_Adv" is full (tri-dimensional) Courant number;
"i,j,k" are coordinates where its maximum occurs; "Cu_W" is contribution
into "Cu_Adv" due to vertical velocity.

Code: Select all

# ifdef MAX_ADV_CFL
#  ifdef MASKING
              ciV=dt*rmask(i,j)*pm(i,j)*pn(i,j)/Hz(i,j,k)
#  else
              ciV=dt*pm(i,j)*pn(i,j)/Hz(i,j,k)
#  endif
              cw=ciV*( max(W(i,j,k), 0.) - min(W(i,j,k-1), 0.))

              cx=cw+ciV*( max(FlxU(i+1,j,k), 0.) -min(FlxU(i,j,k), 0.)
     &                   +max(FlxV(i,j+1,k), 0.) -min(FlxV(i,j,k), 0.))

              if (cx .gt. my_Cu_Adv) then
                my_Cu_Adv=cx
                my_Cu_W=cw
                my_i_cmax=i
                my_j_cmax=j
                my_k_cmax=k
              endif
# else
...........
Where FlxU, FlxV are the same as Huon, Hvon in Rutgers code.


In theory, the most dramatic difference between the two diagnostics is when
velocity field is a tri-dimensional checkerboard pattern. Then averaging of
velocity components leads to zero, hence zero Courant number(?), while
computing it via sum of outgoing fluxes may lead to a large number.


Hernan, adapt it from http://www.atmos.ucla.edu/~alex/ROMS/tmp/diag.F.

User avatar
arango
Site Admin
Posts: 1117
Joined: Wed Feb 26, 2003 4:41 pm
Location: IMCS, Rutgers University
Contact:

Re: Computing Advective CFL in diag.F

#2 Post by arango » Thu Jul 07, 2011 3:28 pm

Thank you, Sasha. I will take a look latter. I have a lot of stuff accumulated.

Post Reply