From d458ec9474709a0ce840889cea0d5ebdfd705aee Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Mon, 6 Feb 2017 13:26:05 -0500 Subject: [PATCH 1/4] Clean up tracer advection PPM implementation There was some code duplication in the implementation of PPM tracer advection, which was eliminated in favour of adding an additional conditional check. --- src/tracer/MOM_tracer_advect.F90 | 125 ++++++++++++++----------------- 1 file changed, 58 insertions(+), 67 deletions(-) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index d2be24dcc3..a8e1b40cfa 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -344,7 +344,7 @@ 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 @@ -423,43 +423,38 @@ 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 + 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) + + 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_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 - ! 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) - 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 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 @@ -608,7 +603,7 @@ 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 @@ -687,42 +682,38 @@ 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 + 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) + + 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_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 - ! 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) - 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)*( 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 From 8ffb474c2a6eb760a1435684d840de425f909f04 Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Mon, 6 Feb 2017 13:51:49 -0500 Subject: [PATCH 2/4] Add Colella & Woodward PPM advection This scheme uses slope data from neighbouring cells to calculate the edge values that are used for constructing the interpolating parabola within a cell. This differs from the Huynh (PPM:H3) scheme in that the stencil is larger. --- src/tracer/MOM_tracer_advect.F90 | 51 +++++++++++++++++++++----------- 1 file changed, 33 insertions(+), 18 deletions(-) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index a8e1b40cfa..21335b8452 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -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 @@ -249,11 +250,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 @@ -263,11 +264,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 @@ -310,7 +311,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 @@ -321,7 +322,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 @@ -348,7 +349,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & 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 @@ -433,10 +434,15 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! 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) - 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 + 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 + 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 @@ -568,7 +574,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 @@ -579,7 +585,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 @@ -607,7 +613,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & 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 @@ -692,10 +698,15 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & ! 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) - 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 + 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 + 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 @@ -864,6 +875,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)) From dc64bd2e052d1dd0e3999c3392976a51a9918ebd Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Tue, 7 Feb 2017 16:14:20 -0500 Subject: [PATCH 3/4] Increase advection stencil with PPM method The non-Huyhn advection reconstruction requires a larger stencil to compute interface values. To ensure the group pass mechanism updates when the stencil is invalid, we change the stencil value to 3 for this advection method. --- src/tracer/MOM_tracer_advect.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 21335b8452..55a8aa9247 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -93,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.") From eb04b3231b73c96c63d77f1dcee99cfef5b9caa9 Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Tue, 7 Feb 2017 16:32:03 -0500 Subject: [PATCH 4/4] Begin documentation of tracer advection This is an incomplete introduction to the tracer advection algorithm in MOM6. It tries to cover aspects of the flux advection algorithm, including numerical considerations, as well as the subcell reconstruction methods that are used for calculating the actual fluxes. --- src/tracer/_Advection.dox | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 src/tracer/_Advection.dox diff --git a/src/tracer/_Advection.dox b/src/tracer/_Advection.dox new file mode 100644 index 0000000000..7f68afb8b1 --- /dev/null +++ b/src/tracer/_Advection.dox @@ -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. + +*/ \ No newline at end of file