SetPMLProfile_LinearEuler2D_PML_t Subroutine

  subroutine SetPMLProfile_LinearEuler2D_PML_t(this,x_interior_min,x_interior_max, &
                                               y_interior_min,y_interior_max, &
                                               pml_width,sigma_max,ramp_exponent)
    !! Populate the per-node sigma_x, sigma_y fields. Only nodes inside
    !! elements whose material name starts with this%pml_material_prefix
    !! receive non-zero damping; all other nodes are forced to zero so
    !! the interior solution remains unmodified.
    implicit none
    class(LinearEuler2D_PML_t),intent(inout) :: this
    real(prec),intent(in) :: x_interior_min,x_interior_max
    real(prec),intent(in) :: y_interior_min,y_interior_max
    real(prec),intent(in) :: pml_width
    real(prec),intent(in) :: sigma_max
    integer,intent(in),optional :: ramp_exponent
    ! Local
    integer :: i,j,iel,matid,p,prefix_len
    real(prec) :: x,y,dx,dy,sx,sy
    character(LEN=SELF_MESH_MATNAME_LENGTH) :: matname
    logical :: is_pml

    this%pml_x_min = x_interior_min
    this%pml_x_max = x_interior_max
    this%pml_y_min = y_interior_min
    this%pml_y_max = y_interior_max
    this%pml_width = pml_width
    this%pml_sigma_max = sigma_max
    if(present(ramp_exponent)) then
      this%pml_ramp_exponent = ramp_exponent
    endif
    p = this%pml_ramp_exponent
    prefix_len = len_trim(this%pml_material_prefix)

    ! Default everything to zero, then fill the PML elements.
    this%sigma_x%interior(:,:,:,1) = 0.0_prec
    this%sigma_y%interior(:,:,:,1) = 0.0_prec

    do iel = 1,this%mesh%nElem
      matid = this%mesh%elemMaterial(iel)
      matname = this%mesh%materialNames(matid)
      is_pml = (len_trim(matname) >= prefix_len) .and. &
               (matname(1:prefix_len) == this%pml_material_prefix(1:prefix_len))
      if(.not. is_pml) cycle

      do j = 1,this%solution%N+1
        do i = 1,this%solution%N+1
          x = this%geometry%x%interior(i,j,iel,1,1)
          y = this%geometry%x%interior(i,j,iel,1,2)

          ! Distance into PML measured from the interior box edge.
          dx = max(0.0_prec,x_interior_min-x,x-x_interior_max)
          dy = max(0.0_prec,y_interior_min-y,y-y_interior_max)

          ! Polynomial ramp normalised by pml_width so sigma -> sigma_max
          ! at the far edge of a layer of thickness pml_width.
          sx = sigma_max*(dx/pml_width)**p
          sy = sigma_max*(dy/pml_width)**p

          this%sigma_x%interior(i,j,iel,1) = sx
          this%sigma_y%interior(i,j,iel,1) = sy
        enddo
      enddo
    enddo

    call this%sigma_x%UpdateDevice()
    call this%sigma_y%UpdateDevice()

  endsubroutine SetPMLProfile_LinearEuler2D_PML_t