Skip to content

PPM advection #419

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
adcroft merged 4 commits into from
Feb 13, 2017
Merged

PPM advection #419

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d458ec9
Clean up tracer advection PPM implementation
angus-g Feb 6, 2017
8ffb474
Add Colella & Woodward PPM advection
angus-g Feb 6, 2017
dc64bd2
Increase advection stencil with PPM method
angus-g Feb 7, 2017
eb04b32
Begin documentation of tracer advection
angus-g Feb 7, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
149 changes: 79 additions & 70 deletions src/tracer/MOM_tracer_advect.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ module MOM_tracer_advect
!< timing of diagnostic output.
logical :: debug !< If true, write verbose checksums for debugging purposes.
logical :: usePPM !< If true, use PPM instead of PLM
logical :: useHuynh !< If true, use the Huynh scheme for PPM interface values
type(group_pass_type) :: pass_uhr_vhr_t_hprev ! For group pass
end type tracer_advect_CS

Expand Down Expand Up @@ -92,7 +93,10 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg, &
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
landvolfill = 1.0e-20 ! This is arbitrary, but must be positive.
stencil = 2 ! The scheme's stencil; 2 for PLM.
stencil = 2 ! The scheme's stencil; 2 for PLM and PPM:H3

! increase stencil size for Colella & Woodward PPM
if (CS%usePPM .and. .not. CS%useHuynh) stencil = 3

if (.not. associated(CS)) call MOM_error(FATAL, "MOM_tracer_advect: "// &
"tracer_advect_init must be called before advect_tracer.")
Expand Down Expand Up @@ -249,11 +253,11 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg, &

! First, advect zonally.
call advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
isv, iev, jsv-stencil, jev+stencil, k, G, GV, CS%usePPM)
isv, iev, jsv-stencil, jev+stencil, k, G, GV, CS%usePPM, CS%useHuynh)

! Next, advect meridionally.
call advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
isv, iev, jsv, jev, k, G, GV, CS%usePPM)
isv, iev, jsv, jev, k, G, GV, CS%usePPM, CS%useHuynh)

domore_k(k) = 0
do j=jsv-stencil,jev+stencil ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
Expand All @@ -263,11 +267,11 @@ subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, CS, Reg, &

! First, advect meridionally.
call advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
isv-stencil, iev+stencil, jsv, jev, k, G, GV, CS%usePPM)
isv-stencil, iev+stencil, jsv, jev, k, G, GV, CS%usePPM, CS%useHuynh)

! Next, advect zonally.
call advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
isv, iev, jsv, jev, k, G, GV, CS%usePPM)
isv, iev, jsv, jev, k, G, GV, CS%usePPM, CS%useHuynh)

domore_k(k) = 0
do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
Expand Down Expand Up @@ -310,7 +314,7 @@ end subroutine advect_tracer
!> This subroutine does 1-d flux-form advection in the zonal direction using
!! a monotonic piecewise linear scheme.
subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
is, ie, js, je, k, G, GV, usePPM)
is, ie, js, je, k, G, GV, usePPM, useHuynh)
type(ocean_grid_type), intent(inout) :: G
type(verticalGrid_type), intent(in) :: GV
type(tracer_type), dimension(ntr), intent(inout) :: Tr
Expand All @@ -321,7 +325,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
logical, dimension(SZJ_(G),SZK_(G)), intent(inout) :: domore_u
real, intent(in) :: Idt
integer, intent(in) :: ntr, is, ie, js, je,k
logical, intent(in) :: usePPM
logical, intent(in) :: usePPM, useHuynh

real, dimension(SZIB_(G),ntr) :: &
slope_x, & ! The concentration slope per grid point in units of
Expand All @@ -344,11 +348,11 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
! in roundoff and can be neglected, in m.
logical :: do_i(SZIB_(G)) ! If true, work on given points.
logical :: do_any_i
integer :: i, j, m
integer :: i, j, m, i_up
real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
logical :: usePLMslope

usePLMslope = .not. usePPM
usePLMslope = .not. (usePPM .and. useHuynh)

min_h = 0.1*GV%Angstrom
h_neglect = GV%H_subroundoff
Expand Down Expand Up @@ -423,43 +427,43 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &

if (usePPM) then
do m=1,ntr ; do I=is-1,ie
! centre cell depending on upstream direction
if (uhh(I) >= 0.0) then
! Implementation of PPM-H3
Tp = Tr(m)%t(i+1,j,k) ; Tc = Tr(m)%t(i,j,k) ; Tm = Tr(m)%t(i-1,j,k)
i_up = i
else
i_up = i+1
endif

! Implementation of PPM-H3
Tp = Tr(m)%t(i_up+1,j,k) ; Tc = Tr(m)%t(i_up,j,k) ; Tm = Tr(m)%t(i_up-1,j,k)

if (useHuynh) then
aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate
aL = max( min(Tc,Tm), aL) ; aL = min( max(Tc,Tm), aL) ! Bound
aR = ( 5.*Tc + ( 2.*Tp - Tm ) )/6. ! H3 estimate
aR = max( min(Tc,Tp), aR) ; aR = min( max(Tc,Tp), aR) ! Bound
dA = aR - aL ; mA = 0.5*( aR + aL )
if (G%mask2dCu(I,j)*G%mask2dCu(I-1,j)*(Tp-Tc)*(Tc-Tm) <= 0.) then
aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells
elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then
aL = 3.*Tc - 2.*aR
elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then
aR = 3.*Tc - 2.*aL
endif
a6 = 6.*Tc - 3. * (aR + aL) ! Curvature
else
aL = 0.5 * ((Tm + Tc) + (slope_x(i_up-1,m) - slope_x(i_up,m)))
aR = 0.5 * ((Tc + Tp) + (slope_x(i_up,m) - slope_x(i_up+1,m)))
endif

dA = aR - aL ; mA = 0.5*( aR + aL )
if (G%mask2dCu(I_up,j)*G%mask2dCu(I_up-1,j)*(Tp-Tc)*(Tc-Tm) <= 0.) then
aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells
elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then
aL = 3.*Tc - 2.*aR
elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then
aR = 3.*Tc - 2.*aL
endif

a6 = 6.*Tc - 3. * (aR + aL) ! Curvature

if (uhh(I) >= 0.0) then
flux_x(I,m) = uhh(I)*( aR - 0.5 * CFL(I) * ( &
( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) )
( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) )
else
! Implementation of PPM-H3
Tp = Tr(m)%t(i+2,j,k) ; Tc = Tr(m)%t(i+1,j,k) ; Tm = Tr(m)%t(i,j,k)
aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate
aL = max( min(Tc,Tm), aL) ; aL = min( max(Tc,Tm), aL) ! Bound
aR = ( 5.*Tc + ( 2.*Tp - Tm ) )/6. ! H3 estimate
aR = max( min(Tc,Tp), aR) ; aR = min( max(Tc,Tp), aR) ! Bound
dA = aR - aL ; mA = 0.5*( aR + aL )
dA = aR - aL ; mA = 0.5*( aR + aL )
if (G%mask2dCu(I,j)*G%mask2dCu(I+1,j)*(Tp-Tc)*(Tc-Tm) <= 0.) then
aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells
elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then
aL = 3.*Tc - 2.*aR
elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then
aR = 3.*Tc - 2.*aL
endif
a6 = 6.*Tc - 3. * (aR + aL) ! Curvature
flux_x(I,m) = uhh(I)*( aL + 0.5 * CFL(I) * ( &
( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) )
( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) )
endif
enddo ; enddo
else ! PLM
Expand Down Expand Up @@ -573,7 +577,7 @@ end subroutine advect_x
!> This subroutine does 1-d flux-form advection using a monotonic piecewise
!! linear scheme.
subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
is, ie, js, je, k, G, GV, usePPM)
is, ie, js, je, k, G, GV, usePPM, useHuynh)
type(ocean_grid_type), intent(inout) :: G
type(verticalGrid_type), intent(in) :: GV
type(tracer_type), dimension(ntr), intent(inout) :: Tr
Expand All @@ -584,7 +588,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
logical, dimension(SZJB_(G),SZK_(G)), intent(inout) :: domore_v
real, intent(in) :: Idt
integer, intent(in) :: ntr, is, ie, js, je,k
logical, intent(in) :: usePPM
logical, intent(in) :: usePPM, useHuynh

real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
slope_y, & ! The concentration slope per grid point in units of
Expand All @@ -608,11 +612,11 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles.
logical :: do_i(SZIB_(G)) ! If true, work on given points.
logical :: do_any_i
integer :: i, j, m
integer :: i, j, m, j_up
real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
logical :: usePLMslope

usePLMslope = .not. usePPM
usePLMslope = .not. (usePPM .and. useHuynh)

min_h = 0.1*GV%Angstrom
h_neglect = GV%H_subroundoff
Expand Down Expand Up @@ -687,42 +691,43 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &

if (usePPM) then
do m=1,ntr ; do i=is,ie
! centre cell depending on upstream direction
if (vhh(i,J) >= 0.0) then
! Implementation of PPM-H3
Tp = Tr(m)%t(i,j+1,k) ; Tc = Tr(m)%t(i,j,k) ; Tm = Tr(m)%t(i,j-1,k)
j_up = j
else
j_up = j + 1
endif

! Implementation of PPM-H3
Tp = Tr(m)%t(i,j_up+1,k) ; Tc = Tr(m)%t(i,j_up,k) ; Tm = Tr(m)%t(i,j_up-1,k)

if (useHuynh) then
aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate
aL = max( min(Tc,Tm), aL) ; aL = min( max(Tc,Tm), aL) ! Bound
aR = ( 5.*Tc + ( 2.*Tp - Tm ) )/6. ! H3 estimate
aR = max( min(Tc,Tp), aR) ; aR = min( max(Tc,Tp), aR) ! Bound
dA = aR - aL ; mA = 0.5*( aR + aL )
if (G%mask2dCv(i,J)*G%mask2dCv(i,J-1)*(Tp-Tc)*(Tc-Tm) <= 0.) then
aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells
elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then
aL = 3.*Tc - 2.*aR
elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then
aR = 3.*Tc - 2.*aL
endif
a6 = 6.*Tc - 3. * (aR + aL) ! Curvature
else
aL = 0.5 * ((Tm + Tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)))
aR = 0.5 * ((Tc + Tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)))
endif

dA = aR - aL ; mA = 0.5*( aR + aL )
if (G%mask2dCv(i,J_up)*G%mask2dCv(i,J_up-1)*(Tp-Tc)*(Tc-Tm) <= 0.) then
aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells
elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then
aL = 3.*Tc - 2.*aR
elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then
aR = 3.*Tc - 2.*aL
endif

a6 = 6.*Tc - 3. * (aR + aL) ! Curvature

if (vhh(i,J) >= 0.0) then
flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * CFL(i) * ( &
( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) )
( aR - aL ) - a6 * ( 1. - 2./3. * CFL(I) ) ) )
else
! Implementation of PPM-H3
Tp = Tr(m)%t(i,j+2,k) ; Tc = Tr(m)%t(i,j+1,k) ; Tm = Tr(m)%t(i,j,k)
aL = ( 5.*Tc + ( 2.*Tm - Tp ) )/6. ! H3 estimate
aL = max( min(Tc,Tm), aL) ; aL = min( max(Tc,Tm), aL) ! Bound
aR = ( 5.*Tc + ( 2.*Tp - Tm ) )/6. ! H3 estimate
aR = max( min(Tc,Tp), aR) ; aR = min( max(Tc,Tp), aR) ! Bound
dA = aR - aL ; mA = 0.5*( aR + aL )
if (G%mask2dCv(i,J)*G%mask2dCv(i,J+1)*(Tp-Tc)*(Tc-Tm) <= 0.) then
aL = Tc ; aR = Tc ! PCM for local extremum and bounadry cells
elseif ( dA*(Tc-mA) > (dA*dA)/6. ) then
aL = 3.*Tc - 2.*aR
elseif ( dA*(Tc-mA) < - (dA*dA)/6. ) then
aR = 3.*Tc - 2.*aL
endif
a6 = 6.*Tc - 3. * (aR + aL) ! Curvature
flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * CFL(i) * ( &
( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) )
( aR - aL ) + a6 * ( 1. - 2./3. * CFL(I) ) ) )
endif
enddo ; enddo
else ! PLM
Expand Down Expand Up @@ -873,6 +878,10 @@ subroutine tracer_advect_init(Time, G, param_file, diag, CS)
CS%usePPM = .false.
case ("PPM:H3")
CS%usePPM = .true.
CS%useHuynh = .true.
case ("PPM")
CS%usePPM = .true.
CS%useHuynh = .false.
case default
call MOM_error(FATAL, "MOM_tracer_advect, tracer_advect_init: "//&
"Unknown TRACER_ADVECTION_SCHEME = "//trim(mesg))
Expand Down
24 changes: 24 additions & 0 deletions src/tracer/_Advection.dox
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
/*! \page Advection Tracer Advection

MOM6 implements a generalised tracer advection scheme, which is a combination of the modified flux advection scheme (Easter, 1993) with reconstructed tracer distributions. The tracer distributions may be piecewise linear (PLM) or piecewise parabolic (PPM), which may itself use either the Colella and Woodward (CW84) or Huynh (H3) reconstruction.

\section Flux_advection Flux advection
The modified flux advection scheme preserves the tracer mixing ratio in a cell across directional splitting by accounting for changes in mass changes. Fluxes are applied to alternate directions in turn, restricting the applied flux so as not to evacuate all mass out of a cell. Because of this, we need to know the stencil used during the calculation of the reconstruction. Every iteration of the splitting algorithm, cells at the edge of a processor's data domain are invalidated. When this invalidation region extends below the halo, a group pass is required to refresh the halo. A larger stencil (such as for the CW84 reconstruction) therefore introduces more frequent updates, and may impact performance.

\section Tracer_reconstruction Tracer reconstruction
While MOM6 only carries the mean tracer concentration in a cell, a higher order reconstruction is computed for the purpose of advection. Reconstructions are also modified to ensure that monotonicity is preserved (i.e. spurious minima or maxima cannot be introduced).

The piecewise linear (PLM) reconstruction uses the monotonic modified van Leer scheme (Lin et al., 1994). One might think to use the average of the one-sided differences of mean tracer concentration within a cell to calculate the slope of the linear reconstruction, however this method guarantees neither monotonicity, nor positive definiteness. Instead, the method is locally limited to the minimum of this average slope and each of the one-sided slopes, i.e. \f[\Delta \Phi_i = \min\left\{\left|[\Delta \Phi_i]_\text{avg}\right|, 2\left(\Phi_i - \Phi_i^\text{min}\right), 2\left(\Phi_i^\text{max} - \Phi_i\right)\right\}\f]
(where \f$\Phi_i^\text{min}\f$ is the minimum in the 3-point stencil).

In a PPM scheme, for a cell with mean tracer concentration \f$\Phi_i\f$, the values at the left and right interfaces, \f$\Phi_{L,i}\f$ and \f$\Phi_{R,i}\f$ must be estimated. First, an interpolation is used to calculate \f$\Phi_{i-1/2}\f$ and \f$\Phi_{i+1/2}\f$. These values are then modified to preserve monotonicity in each cell, which introduces discontinuities between cell edges (e.g. \f$\Phi_{R,i}\f$ and \f$\Phi_{L,i+1}\f$).

The reconstruction \f$\Phi_i(\xi)\f$ then satisfies three properties:

- total amount of tracer is conserved, \f$\int_{\xi_{i-1/2}}^{\xi_{i+1/2}} \Phi_i(\xi') \,\mathrm d\xi' = \Phi_i\f$
- left interface value matches, \f$\Phi(\xi_{i-1/2}) = \Phi_{L,i}\f$
- right interface value matches, \f$\Phi(\xi_{i+1/2}) = \Phi_{R,i}\f$

There are two methods of reconstruction for a piecewise parabolic (PPM) profile. They differ in the estimate of interface values \f$\Phi_{i+1/2}\f$ prior to monotonicity limiting. The Colella and Woodward (1984) scheme makes use of the limited slope \f$\Delta\Phi_i\f$ from PLM, above. This has the effect of requiring a larger stencil for each reconstruction. On the other hand, the Huynh (1997) scheme reduces the requirement of this stencil, by only examining the tracer concentrations in adjacent cells, at the same time reducing order of accuracy of the reconstruction.

*/