Creating a laplacian sponge layer mechanism

Initial plan:

• Track down where rayleigh friction code is.
• Track down where second-order operators are computed.

Understanding the rayleigh friction code:

      if( .not.flagstruct%RF_fast .and. flagstruct%tau > 0. ) then 
        if ( gridstruct%grid_type<4 .or. gridstruct%bounded_domain .or. is_ideal_case ) then 
!         if ( flagstruct%RF_fast ) then
!            call Ray_fast(abs(dt), npx, npy, npz, pfull, flagstruct%tau, u, v, w,  &
!                          dp_ref, ptop, hydrostatic, flagstruct%rf_cutoff, bd)
!         else
             call Rayleigh_Super(abs(bdt), npx, npy, npz, ks, pfull, phis, flagstruct%tau, u, v, w, pt,  &
                  ua, va, delz, gridstruct%agrid, cp_air, rdgas, ptop, hydrostatic,    &    
                 .not. (gridstruct%bounded_domain .or. is_ideal_case), flagstruct%rf_cutoff, gridstruct, domain, bd)
!         endif
        else
             call Rayleigh_Friction(abs(bdt), npx, npy, npz, ks, pfull, flagstruct%tau, u, v, w, pt,  &
                  ua, va, delz, cp_air, rdgas, ptop, hydrostatic, .true., flagstruct%rf_cutoff, gridstruct, domain, bd)
        endif
      endif

! *** Inline Rayleigh friction here?
   if( flagstruct%RF_fast .and. flagstruct%tau > 0. )  & 
   call Ray_fast(abs(dt), npx, npy, npz, pfull, flagstruct%tau, u, v, w,  &
                      ks, dp_ref, ptop, hydrostatic, flagstruct%rf_cutoff, bd)

Understanding the second-order diffusion operators in FV3

The model I'll be using is divergence damping. Function calls can be found around line 700 of atmos_cubed_sphere/model/dyn_core.F90.

  ! Sponge layers with del-2 damping on divergence, vorticity, w, z, and air mass (delp).
! no special damping of potential temperature in sponge layers
              if ( k==1 ) then 
! Divergence damping:
                 nord_k=0;
                 if (is_ideal_case) then 
                    d2_divg = max(flagstruct%d2_bg, flagstruct%d2_bg_k1)
                 else
                    d2_divg = max(0.01, flagstruct%d2_bg, flagstruct%d2_bg_k1)
                 endif
! Vertical velocity:
                   nord_w=0; damp_w = d2_divg
                   if ( flagstruct%do_vort_damp ) then 
! damping on delp and vorticity:
                        nord_v(k)=0;
#ifndef HIWPP
                        damp_vt(k) = 0.5*d2_divg
#endif
                   endif
                   d_con_k = 0. 
              elseif ( k==2 .and. flagstruct%d2_bg_k2>0.01 ) then 
                   nord_k=0; d2_divg = max(flagstruct%d2_bg, flagstruct%d2_bg_k2)
                   nord_w=0; damp_w = d2_divg
                   if ( flagstruct%do_vort_damp ) then 
                        nord_v(k)=0;
#ifndef HIWPP
                        damp_vt(k) = 0.5*d2_divg
#endif
                   endif
                   d_con_k = 0. 
              elseif ( k==3 .and. flagstruct%d2_bg_k2>0.05 ) then 
                   nord_k=0;  d2_divg = max(flagstruct%d2_bg, 0.2*flagstruct%d2_bg_k2)
                   nord_w=0;  damp_w = d2_divg
                   d_con_k = 0. 
              endif
       endif

It appears that horizontal second-order operators are handled within the call to the explicit d_sw shallow water subroutine. Therefore there are flags so that the damping coefficients are set to appropriate values in the sponge levels.

designing an ideal test case to ensure that this implementation is behaving reasonably.

Use the hydrostatic staniforth white test case? Bring that over, at least.

Use the files I modified for the class in the following way:

Note from 08/06/22

Now that I've looked in the source code in JT's new FV3 dycore branch, it appears that he's already added namelist variables that enable $\nabla^2$ sponge-layer diffusion in the top two layers. The default namelist configuration disables the conversion of dissipated energy from Rayleigh friction to thermal energy, This is controlled by the namelist variable fv3_d_con.

My task now is to benchmark the performance of these two diffusion methods .

Comparing Rayleigh friction when fv3_d_con = 1

Actual vector laplacian

The vector laplacian is somewhat inconveniently defined (this also works for tensors) as $$(\nabla \cdot \nabla) \mathbf{v}.$$

Note that that formulation uses operators that can easily be formulated in coordinate-free notation. The end result we wish to get is, under the assumption that $\mathbf{v}_{k} = 0$, we can calculate the laplacian in terms of 2D $\zeta$ relative vorticity and divergence $\delta$ which are available in the D-grid routines.

If we use some vector calculus identities $\nabla (\nabla \cdot \mathbf{v}) - \nabla \times (\nabla \times \mathbf{v})$

Assume that $\mathbf{v} = \begin{bmatrix} v_1 \\ v_2 \\ 0 \end{bmatrix}$ in which case $\nabla \cdot \mathbf{v} = \delta$ and (using cartesian coordiantes for confidence). We assume we're working at a fixed $r$, and so vertical variation of horizontal wind can be disregarded. We then find that $$\nabla \times \mathbf{v} = \zeta \mathbf{k}$$

then we find that $$\nabla \times \zeta \mathbf{k} = \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ \partial_x & \partial_y & \partial_z \\ 0 & 0 & \zeta \end{vmatrix} = \begin{bmatrix} \partial_y \zeta \\ -\partial_x \zeta \\ 0 \end{bmatrix}$$

Therefore we need to know how to take the gradient of a scalar.

But! It turns out that adding a $(\nabla \cdot \nabla) \mathbf{v}$ term to the RHS of the momentum equation, i.e. $$\begin{align*} \der{u}{t} &= \ldots + (\nu (\nabla \cdot \nabla) \mathbf{v})_1\\ \der{v}{t} &= \ldots + (\nu (\nabla \cdot \nabla) \mathbf{v})_2 \end{align*}$$ just gives us $$\begin{align*} \der{u}{t} &= \ldots + \partial_{x_1} \delta + \partial_{x_2} \zeta\\ \der{v}{t} &= \ldots + \partial_{x_2} \delta - \partial_{x_1} \zeta \end{align*}$$

From page 15 of the fv3 specification we know that we can calculate "cell integrated" quantity on the dual grid: $$D = \frac{1}{A_c} \left[\delta_x (u_c \Delta y_c \sin \alpha) + \delta_y (v_c \Delta x_c \sin \alpha) \right]$$

divergence_corner in sw_core shows how to deal with index hell. It appears that $\delta_{x, y}$ is a finite differencing operation. See line 1788 in the sw_core module file. For volume mean vorticity see line 1249, i.e. wk(i,j) = rarea(i,j)*(vt(i,j)-vt(i,j+1)-ut(i,j)+ut(i+1,j))

When applied to vorticity, the analogue of equations (2.3) and (2.4) in the FV3 document can be found at lines 1676 and 1689. i.e.

fx2(i,j) = gridstruct%del6_v(i,j)*(d2(i,j)-d2(i-1,j))
fy2(i,j) = gridstruct%del6_u(i,j)*(d2(i,j)-d2(i,j-1))

where f(xy)2 seems to represent a flux quantity. I.e. ut, vt in the code are flux quantities.

It doesn't appear that the preprocessor name USE_SG is defined anywhere? I've checked the FV3 github as well as the entire CESM codebase to see if it appears in the build infrastructure.

Synthesizing these breadcrumbs:

Application of 0th order divergence damping involves adding a $\nabla \delta$ term to the RHS of the momentum equation. Therefore the treatment of 0th order vorticity damping should show us how to add what we need.

Use the subroutine del6_vt_flux as a guide.

Using the two (simplified) lines as a first guess,

vort(i,j) = damp*wk(i,j)
fx2(i,j) = gridstruct%del6_v(i,j)*(vort(i-1,j)-vort(i,j))
fy2(i,j) = gridstruct%del6_u(i,j)*(vort(i,j-1)-vort(i,j))
! fx2(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1)) * dy(i,j) * (vort(i-1,j)-vort(i,j)) * rdxc(i,j)
! fy2(i,j) = 0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))* dx(i,j) * (vort(i,j-1)-vort(i,j)) * rdyc(i,j)

then our quantity should look something like

! uses C grid convention, seemingly
vort(i,j) = damp*wk_vort(i,j)
divg(i,j) = damp*wk_divg(i,j)
! lap_nu_v = lap_nu_u almost certainly. 
fx2(i,j) = gridstruct%lap_nu_v(i,j)*(divg(i-1,j)-divg(i,j) + vort(i,j-1)-vort(i,j))
fy2(i,j) = gridstruct%lap_nu_u(i,j)*(divg(i,j-1)-divg(i,j) - (vort(i-1,j)-vort(i,j)))
! this is probably wrong:
! fx2(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1)) * dy(i,j) * (divg(i-1,j)-divg(i,j) + vort(i,j-1)-vort(i,j)) * rdxc(i,j)
! fy2(i,j) = 0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))* dx(i,j) * (divg(i,j-1)-divg(i,j) - (vort(i-1,j)-vort(i,j))) * rdyc(i,j)

But we need to use the spherical grid version. Let's do that first thing tomorrow.