SourceMethod_ESAtmo2D_t Subroutine

  1. SELF_ESAtmo2D_t.f90
  2. SELF_ESAtmo2D_t
  3. SourceMethod_ESAtmo2D_t

public subroutine SourceMethod_ESAtmo2D_t(this)

Souza et al. (2023) non-conservative gravity flux differencing.

The rhov equation carries the body force -rho * partial_y Phi where Phi = gy is the geopotential (state variable index 5; no flux, no source of its own). On a curvilinear mesh, the SBP-EC strong-form flux differencing for the non-conservative product rho * partial_y Phi reads:

[rho * d_y Phi]_i = (1/J_i) sum_r sum_j D_split^r[i_r, j] * 0.5*(Ja^r_y(i) + Ja^r_y(j)) * _log(s_i, s_j) * (Phi_j - Phi_i)

and we set source(rhov) = - that result. The other variables (rho, rhou, rho*theta, Phi) all have zero source.

Arguments

TypeIntentOptionalAttributesName
class(ESAtmo2D_t), intent(inout) :: this

Calls

proc~~sourcemethod_esatmo2d_t~~CallsGraph proc~sourcemethod_esatmo2d_t SourceMethod_ESAtmo2D_t proc~log_mean log_mean proc~sourcemethod_esatmo2d_t->proc~log_mean

Contents

Source Code

SourceMethod_ESAtmo2D_t

Source Code

  subroutine SourceMethod_ESAtmo2D_t(this)
    !! Souza et al. (2023) non-conservative gravity flux differencing.
    !!
    !! The rho*v equation carries the body force -rho * partial_y Phi where
    !! Phi = g*y is the geopotential (state variable index 5; no flux, no
    !! source of its own). On a curvilinear mesh, the SBP-EC strong-form
    !! flux differencing for the non-conservative product rho * partial_y Phi
    !! reads:
    !!
    !!   [rho * d_y Phi]_i = (1/J_i) sum_r sum_j D_split^r[i_r, j]
    !!                       * 0.5*(Ja^r_y(i) + Ja^r_y(j))
    !!                       * <rho>_log(s_i, s_j) * (Phi_j - Phi_i)
    !!
    !! and we set source(rho*v) = - that result. The other variables
    !! (rho, rho*u, rho*theta, Phi) all have zero source.
    implicit none
    class(ESAtmo2D_t),intent(inout) :: this
    ! Local
    integer :: i,j,iEl,nn
    real(prec) :: rho_ijk,phi_ijk,rho_p,phi_p
    real(prec) :: rho_log,Ja_avg,acc,jac

    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                  iEl=1:this%mesh%nElem)

      rho_ijk = this%solution%interior(i,j,iEl,1)
      phi_ijk = this%solution%interior(i,j,iEl,5)
      acc = 0.0_prec

      do nn = 1,this%solution%N+1
        ! xi^1 partner (nn, j) — uses Ja^1_y (d=2, r=1)
        rho_p = this%solution%interior(nn,j,iEl,1)
        phi_p = this%solution%interior(nn,j,iEl,5)
        rho_log = log_mean(rho_ijk,rho_p)
        Ja_avg = 0.5_prec*(this%geometry%dsdx%interior(i,j,iEl,1,2,1)+ &
                           this%geometry%dsdx%interior(nn,j,iEl,1,2,1))
        acc = acc+this%solution%interp%dSplitMatrix(nn,i)* &
              Ja_avg*rho_log*(phi_p-phi_ijk)

        ! xi^2 partner (i, nn) — uses Ja^2_y (d=2, r=2)
        rho_p = this%solution%interior(i,nn,iEl,1)
        phi_p = this%solution%interior(i,nn,iEl,5)
        rho_log = log_mean(rho_ijk,rho_p)
        Ja_avg = 0.5_prec*(this%geometry%dsdx%interior(i,j,iEl,1,2,2)+ &
                           this%geometry%dsdx%interior(i,nn,iEl,1,2,2))
        acc = acc+this%solution%interp%dSplitMatrix(nn,j)* &
              Ja_avg*rho_log*(phi_p-phi_ijk)
      enddo

      jac = this%geometry%J%interior(i,j,iEl,1)
      this%source%interior(i,j,iEl,1) = 0.0_prec
      this%source%interior(i,j,iEl,2) = 0.0_prec
      this%source%interior(i,j,iEl,3) = -acc/jac
      this%source%interior(i,j,iEl,4) = 0.0_prec
      this%source%interior(i,j,iEl,5) = 0.0_prec

    enddo

  endsubroutine SourceMethod_ESAtmo2D_t