sourcemethod_LinearEuler2D_PML_t Subroutine

  subroutine sourcemethod_LinearEuler2D_PML_t(this)
    !! Hu (2001) unsplit PML source term, evaluated per node. In the
    !! interior (sigma_x = sigma_y = 0) this leaves the acoustic
    !! variables untouched and the auxiliaries integrate q in time but
    !! never re-enter the dynamics (since the coupling coefficient
    !! sigma_x*sigma_y is zero). Inside the PML, the (sigma_x + sigma_y)
    !! damping plus the sigma_x*sigma_y*phi term produce the correct
    !! perfectly-matched behaviour for the linear Euler 2D system.
    !!
    !! Note: we override sourcemethod rather than source2d because the
    !! per-node sigma_x(i,j,iel), sigma_y(i,j,iel) lookups are not
    !! reachable from the pure source2d(s,dsdx) signature.
    implicit none
    class(LinearEuler2D_PML_t),intent(inout) :: this
    ! Local
    integer :: i,j,iel
    real(prec) :: sx,sy
    real(prec) :: s(1:this%nvar)

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

      sx = this%sigma_x%interior(i,j,iel,1)
      sy = this%sigma_y%interior(i,j,iel,1)
      s = this%solution%interior(i,j,iel,1:this%nvar)

      this%source%interior(i,j,iel,1) = -(sx+sy)*s(1)-sx*sy*s(6)
      this%source%interior(i,j,iel,2) = -(sx+sy)*s(2)-sx*sy*s(7)
      this%source%interior(i,j,iel,3) = -(sx+sy)*s(3)-sx*sy*s(8)
      this%source%interior(i,j,iel,4) = -(sx+sy)*s(4)-sx*sy*s(9)
      this%source%interior(i,j,iel,5) = 0.0_prec
      this%source%interior(i,j,iel,6) = s(1)
      this%source%interior(i,j,iel,7) = s(2)
      this%source%interior(i,j,iel,8) = s(3)
      this%source%interior(i,j,iel,9) = s(4)

    enddo

  endsubroutine sourcemethod_LinearEuler2D_PML_t