SphericalSoundWave_LinearEuler2D_t Subroutine

  subroutine SphericalSoundWave_LinearEuler2D_t(this,rhoprime,Lr,x0,y0,c)
    !! This subroutine sets the initial condition for a weak blast wave
    !! problem. The initial condition is given by
    !!
    !! \begin{equation}
    !! \begin{aligned}
    !! u &= 0 \\
    !! v &= 0 \\
    !! p &= \rho' c^2 \exp\left( -\ln(2) \frac{(x-x_0)^2 + (y-y_0)^2}{L_r^2} \right)
    !! \end{aligned}
    !! \end{equation}
    !!
    !! `rhoprime` is the amplitude of the (diagnostic) density anomaly that the
    !! pressure pulse corresponds to through the acoustic relation p = rho*c^2;
    !! the density anomaly itself is not a solution variable.
    !! The sound speed `c` (passed as an argument since it is not a
    !! scalar model attribute) and the background density `this%rho0` are set
    !! uniformly across the domain.
    implicit none
    class(LinearEuler2D_t),intent(inout) :: this
    real(prec),intent(in) ::  rhoprime,Lr,x0,y0,c
    ! Local
    integer :: i,j,iEl
    real(prec) :: x,y,rho,r

    print*,__FILE__," : Configuring weak blast wave initial condition. "
    print*,__FILE__," : rhoprime = ",rhoprime
    print*,__FILE__," : Lr = ",Lr
    print*,__FILE__," : x0 = ",x0
    print*,__FILE__," : y0 = ",y0
    print*,__FILE__," : c = ",c

    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
                  iel=1:this%mesh%nElem)
      x = this%geometry%x%interior(i,j,iEl,1,1)-x0
      y = this%geometry%x%interior(i,j,iEl,1,2)-y0
      r = sqrt(x**2+y**2)

      rho = (rhoprime)*exp(-log(2.0_prec)*r**2/Lr**2)

      this%solution%interior(i,j,iEl,1) = 0.0_prec
      this%solution%interior(i,j,iEl,2) = 0.0_prec
      this%solution%interior(i,j,iEl,3) = rho*c*c
      this%solution%interior(i,j,iEl,4) = c
      this%solution%interior(i,j,iEl,5) = this%rho0 ! uniform background density

    enddo

    call this%ReportMetrics()
    call this%solution%UpdateDevice()

  endsubroutine SphericalSoundWave_LinearEuler2D_t