weno7_face_swept Function

public pure function weno7_face_swept(qm3, qm2, qm1, q0, qp1, qp2, qp3, sigma) result(face)

WENO7-Z swept-average face value (u > 0, upwind cell = q0).

Four cubic candidates with optimal weights d=(4,18,12,1)/35. Z-weights: tau7 = |beta0 - beta3|.

Coefficient matrices from Balsara & Shu (2000), Table 1. Smoothness indicators from Balsara & Shu (2000), eq. 2.17.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: qm3

Cell averages at i-3 through i+3.
real(kind=wp), intent(in) :: qm2

Cell averages at i-3 through i+3.
real(kind=wp), intent(in) :: qm1

Cell averages at i-3 through i+3.
real(kind=wp), intent(in) :: q0

Cell averages at i-3 through i+3.
real(kind=wp), intent(in) :: qp1

Cell averages at i-3 through i+3.
real(kind=wp), intent(in) :: qp2

Cell averages at i-3 through i+3.
real(kind=wp), intent(in) :: qp3

Cell averages at i-3 through i+3.
real(kind=wp), intent(in) :: sigma

CFL of upwind cell: |u_eff|*dt/dx.

Return Value real(kind=wp)


Called by

proc~~weno7_face_swept~~CalledByGraph proc~weno7_face_swept weno7_face_swept proc~weno_face_conc_x weno_face_conc_x proc~weno_face_conc_x->proc~weno7_face_swept proc~weno_face_conc_y weno_face_conc_y proc~weno_face_conc_y->proc~weno7_face_swept proc~drain_swept_flux_x_weno drain_swept_flux_x_weno proc~drain_swept_flux_x_weno->proc~weno_face_conc_x proc~drain_swept_flux_y_weno drain_swept_flux_y_weno proc~drain_swept_flux_y_weno->proc~weno_face_conc_y proc~continuity_tracer_drain continuity_tracer_drain proc~continuity_tracer_drain->proc~drain_swept_flux_x_weno proc~continuity_tracer_drain->proc~drain_swept_flux_y_weno proc~ocean_dyn_flush_tracer_window ocean_dyn_flush_tracer_window proc~ocean_dyn_flush_tracer_window->proc~continuity_tracer_drain proc~ocean_dyn_step ocean_dyn_step proc~ocean_dyn_step->proc~continuity_tracer_drain proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~continuity_tracer_drain proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~ocean_dyn_flush_tracer_window proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step proc~engine_step->proc~ocean_dyn_step_split proc~rdb_ocean_set_tracer rdb_ocean_set_tracer proc~rdb_ocean_set_tracer->proc~ocean_dyn_flush_tracer_window

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: A0
real(kind=wp), private :: A1
real(kind=wp), private :: A2
real(kind=wp), private :: A3
real(kind=wp), private :: B0
real(kind=wp), private :: B1
real(kind=wp), private :: B2
real(kind=wp), private :: B3
real(kind=wp), private :: C0
real(kind=wp), private :: C1
real(kind=wp), private :: C2
real(kind=wp), private :: C3
real(kind=wp), private :: D0
real(kind=wp), private :: D1
real(kind=wp), private :: D2
real(kind=wp), private :: D3
real(kind=wp), private :: alpha0
real(kind=wp), private :: alpha1
real(kind=wp), private :: alpha2
real(kind=wp), private :: alpha3
real(kind=wp), private :: alpha_sum
real(kind=wp), private :: beta0
real(kind=wp), private :: beta1
real(kind=wp), private :: beta2
real(kind=wp), private :: beta3
real(kind=wp), private, parameter :: d0_w = 4.0_wp/35.0_wp
real(kind=wp), private, parameter :: d1_w = 18.0_wp/35.0_wp
real(kind=wp), private, parameter :: d2_w = 12.0_wp/35.0_wp
real(kind=wp), private, parameter :: d3_w = 1.0_wp/35.0_wp
real(kind=wp), private, parameter :: eps = 1.0e-36_wp
real(kind=wp), private :: face0
real(kind=wp), private :: face1
real(kind=wp), private :: face2
real(kind=wp), private :: face3
real(kind=wp), private :: omega0
real(kind=wp), private :: omega1
real(kind=wp), private :: omega2
real(kind=wp), private :: omega3
real(kind=wp), private :: p_half
real(kind=wp), private :: pp_half
real(kind=wp), private :: ppp_half
real(kind=wp), private :: pppp_half
real(kind=wp), private :: tau7

Source Code

   pure function weno7_face_swept(qm3, qm2, qm1, q0, qp1, qp2, qp3, sigma) result(face)
      !! WENO7-Z swept-average face value (u > 0, upwind cell = q0).
      !!
      !! Four cubic candidates with optimal weights d=(4,18,12,1)/35.
      !! Z-weights: tau7 = |beta0 - beta3|.
      !!
      !! Coefficient matrices from Balsara & Shu (2000), Table 1.
      !! Smoothness indicators from Balsara & Shu (2000), eq. 2.17.
      !$acc routine seq
      real(wp), intent(in) :: qm3, qm2, qm1, q0, qp1, qp2, qp3
         !! Cell averages at i-3 through i+3.
      real(wp), intent(in) :: sigma
         !! CFL of upwind cell: |u_eff|*dt/dx.
      real(wp) :: face

      ! Optimal weights for WENO7 (BS00)
      real(wp), parameter :: d0_w = 4.0_wp/35.0_wp
      real(wp), parameter :: d1_w = 18.0_wp/35.0_wp
      real(wp), parameter :: d2_w = 12.0_wp/35.0_wp
      real(wp), parameter :: d3_w = 1.0_wp/35.0_wp
      real(wp), parameter :: eps = 1.0e-36_wp

      ! Cubic polynomial coefficients A + B*xi + C*xi^2 + D*xi^3
      ! for each of the 4 stencils (via exact rational matrices from BS00)
      real(wp) :: A0, B0, C0, D0
      real(wp) :: A1, B1, C1, D1
      real(wp) :: A2, B2, C2, D2
      real(wp) :: A3, B3, C3, D3

      ! Smoothness indicators (BS00 eq 2.17 exact)
      real(wp) :: beta0, beta1, beta2, beta3, tau7

      ! Z-weights
      real(wp) :: alpha0, alpha1, alpha2, alpha3, alpha_sum
      real(wp) :: omega0, omega1, omega2, omega3

      ! Swept-average face values per candidate
      real(wp) :: face0, face1, face2, face3

      ! Cubic swept-average helpers: p(1/2), p'(1/2), p''(1/2), p'''(1/2)
      real(wp) :: p_half, pp_half, ppp_half, pppp_half

      ! Stencil 0: offsets (-3,-2,-1,0) — q values: qm3,qm2,qm1,q0
      A0 = (1.0_wp/24.0_wp)*qm3 + (-1.0_wp/6.0_wp)*qm2 &
           + (5.0_wp/24.0_wp)*qm1 + (11.0_wp/12.0_wp)*q0
      B0 = (-7.0_wp/24.0_wp)*qm3 + (11.0_wp/8.0_wp)*qm2 &
           + (-23.0_wp/8.0_wp)*qm1 + (43.0_wp/24.0_wp)*q0
      C0 = (-1.0_wp/2.0_wp)*qm3 + 2.0_wp*qm2 &
           + (-5.0_wp/2.0_wp)*qm1 + 1.0_wp*q0
      D0 = (-1.0_wp/6.0_wp)*qm3 + (1.0_wp/2.0_wp)*qm2 &
           + (-1.0_wp/2.0_wp)*qm1 + (1.0_wp/6.0_wp)*q0

      ! Stencil 1: offsets (-2,-1,0,+1) — q values: qm2,qm1,q0,qp1
      A1 = 0.0_wp*qm2 + (-1.0_wp/24.0_wp)*qm1 &
           + (13.0_wp/12.0_wp)*q0 + (-1.0_wp/24.0_wp)*qp1
      B1 = (5.0_wp/24.0_wp)*qm2 + (-9.0_wp/8.0_wp)*qm1 &
           + (5.0_wp/8.0_wp)*q0 + (7.0_wp/24.0_wp)*qp1
      C1 = 0.0_wp*qm2 + (1.0_wp/2.0_wp)*qm1 &
           + (-1.0_wp)*q0 + (1.0_wp/2.0_wp)*qp1
      D1 = (-1.0_wp/6.0_wp)*qm2 + (1.0_wp/2.0_wp)*qm1 &
           + (-1.0_wp/2.0_wp)*q0 + (1.0_wp/6.0_wp)*qp1

      ! Stencil 2: offsets (-1,0,+1,+2) — q values: qm1,q0,qp1,qp2
      A2 = (-1.0_wp/24.0_wp)*qm1 + (13.0_wp/12.0_wp)*q0 &
           + (-1.0_wp/24.0_wp)*qp1 + 0.0_wp*qp2
      B2 = (-7.0_wp/24.0_wp)*qm1 + (-5.0_wp/8.0_wp)*q0 &
           + (9.0_wp/8.0_wp)*qp1 + (-5.0_wp/24.0_wp)*qp2
      C2 = (1.0_wp/2.0_wp)*qm1 + (-1.0_wp)*q0 &
           + (1.0_wp/2.0_wp)*qp1 + 0.0_wp*qp2
      D2 = (-1.0_wp/6.0_wp)*qm1 + (1.0_wp/2.0_wp)*q0 &
           + (-1.0_wp/2.0_wp)*qp1 + (1.0_wp/6.0_wp)*qp2

      ! Stencil 3: offsets (0,+1,+2,+3) — q values: q0,qp1,qp2,qp3
      A3 = (11.0_wp/12.0_wp)*q0 + (5.0_wp/24.0_wp)*qp1 &
           + (-1.0_wp/6.0_wp)*qp2 + (1.0_wp/24.0_wp)*qp3
      B3 = (-43.0_wp/24.0_wp)*q0 + (23.0_wp/8.0_wp)*qp1 &
           + (-11.0_wp/8.0_wp)*qp2 + (7.0_wp/24.0_wp)*qp3
      C3 = 1.0_wp*q0 + (-5.0_wp/2.0_wp)*qp1 &
           + 2.0_wp*qp2 + (-1.0_wp/2.0_wp)*qp3
      D3 = (-1.0_wp/6.0_wp)*q0 + (1.0_wp/2.0_wp)*qp1 &
           + (-1.0_wp/2.0_wp)*qp2 + (1.0_wp/6.0_wp)*qp3

      ! Smoothness indicators: the published integer-coefficient quadratic
      ! forms (Balsara & Shu 2000, eq. 2.17), one per candidate stencil.
      ! These are the forms validated in the numpy prototype
      ! (local_archive/prototypes/tracer_weno/) — do not substitute re-derivations.
      beta0 = qm3*(547.0_wp*qm3 - 3882.0_wp*qm2 + 4642.0_wp*qm1 - 1854.0_wp*q0) &
              + qm2*(7043.0_wp*qm2 - 17246.0_wp*qm1 + 7042.0_wp*q0) &
              + qm1*(11003.0_wp*qm1 - 9402.0_wp*q0) &
              + 2107.0_wp*q0**2
      beta1 = qm2*(267.0_wp*qm2 - 1642.0_wp*qm1 + 1602.0_wp*q0 - 494.0_wp*qp1) &
              + qm1*(2843.0_wp*qm1 - 5966.0_wp*q0 + 1922.0_wp*qp1) &
              + q0*(3443.0_wp*q0 - 2522.0_wp*qp1) &
              + 547.0_wp*qp1**2
      beta2 = qm1*(547.0_wp*qm1 - 2522.0_wp*q0 + 1922.0_wp*qp1 - 494.0_wp*qp2) &
              + q0*(3443.0_wp*q0 - 5966.0_wp*qp1 + 1602.0_wp*qp2) &
              + qp1*(2843.0_wp*qp1 - 1642.0_wp*qp2) &
              + 267.0_wp*qp2**2
      beta3 = q0*(2107.0_wp*q0 - 9402.0_wp*qp1 + 7042.0_wp*qp2 - 1854.0_wp*qp3) &
              + qp1*(11003.0_wp*qp1 - 17246.0_wp*qp2 + 4642.0_wp*qp3) &
              + qp2*(7043.0_wp*qp2 - 3882.0_wp*qp3) &
              + 547.0_wp*qp3**2

      tau7 = abs(beta0 - beta3)
      alpha0 = d0_w*(1.0_wp + (tau7/(eps + beta0))**2)
      alpha1 = d1_w*(1.0_wp + (tau7/(eps + beta1))**2)
      alpha2 = d2_w*(1.0_wp + (tau7/(eps + beta2))**2)
      alpha3 = d3_w*(1.0_wp + (tau7/(eps + beta3))**2)
      alpha_sum = alpha0 + alpha1 + alpha2 + alpha3
      omega0 = alpha0/alpha_sum
      omega1 = alpha1/alpha_sum
      omega2 = alpha2/alpha_sum
      omega3 = alpha3/alpha_sum

      ! Cubic swept-average over [1/2 - sigma, 1/2]:
      !   q_face = p(1/2) - sigma/2 * p'(1/2) + sigma^2/6 * p''(1/2) - sigma^3/24 * p'''(1/2)
      ! where:
      !   p(1/2)   = A + B/2 + C/4 + D/8
      !   p'(1/2)  = B + C + 3D/4
      !   p''(1/2) = 2C + 3D
      !   p'''(1/2)= 6D

      ! Stencil 0
      p_half = A0 + 0.5_wp*B0 + 0.25_wp*C0 + 0.125_wp*D0
      pp_half = B0 + C0 + 0.75_wp*D0
      ppp_half = 2.0_wp*C0 + 3.0_wp*D0
      pppp_half = 6.0_wp*D0
      face0 = p_half - 0.5_wp*sigma*pp_half &
              + sigma**2/6.0_wp*ppp_half - sigma**3/24.0_wp*pppp_half

      ! Stencil 1
      p_half = A1 + 0.5_wp*B1 + 0.25_wp*C1 + 0.125_wp*D1
      pp_half = B1 + C1 + 0.75_wp*D1
      ppp_half = 2.0_wp*C1 + 3.0_wp*D1
      pppp_half = 6.0_wp*D1
      face1 = p_half - 0.5_wp*sigma*pp_half &
              + sigma**2/6.0_wp*ppp_half - sigma**3/24.0_wp*pppp_half

      ! Stencil 2
      p_half = A2 + 0.5_wp*B2 + 0.25_wp*C2 + 0.125_wp*D2
      pp_half = B2 + C2 + 0.75_wp*D2
      ppp_half = 2.0_wp*C2 + 3.0_wp*D2
      pppp_half = 6.0_wp*D2
      face2 = p_half - 0.5_wp*sigma*pp_half &
              + sigma**2/6.0_wp*ppp_half - sigma**3/24.0_wp*pppp_half

      ! Stencil 3
      p_half = A3 + 0.5_wp*B3 + 0.25_wp*C3 + 0.125_wp*D3
      pp_half = B3 + C3 + 0.75_wp*D3
      ppp_half = 2.0_wp*C3 + 3.0_wp*D3
      pppp_half = 6.0_wp*D3
      face3 = p_half - 0.5_wp*sigma*pp_half &
              + sigma**2/6.0_wp*ppp_half - sigma**3/24.0_wp*pppp_half

      face = omega0*face0 + omega1*face1 + omega2*face2 + omega3*face3
   end function weno7_face_swept