Calculate the characteristic fluxes for all elements for a single mode.
Having all of these calculations in a subroutine made it easier to ensure that all local variables are private (for OpenMP parallelization) and to avoid segmentation faults.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(cgf), | intent(in) | :: | rho | A complex grid function containing the time derivative variable, . |
||
| type(cgf), | intent(in) | :: | pi | A complex grid function containing the radial derivative variable, . |
||
| integer(kind=ip), | intent(in) | :: | nmodes | The total number of modes. |
||
| integer(kind=ip), | intent(in) | :: | mode | The mode for which the flux is calculated. |
||
| type(rgf), | intent(in) | :: | r | A real grid function that contains the Scharzschild radial coordinate. |
||
| type(scal_schw_eff), | intent(inout) | :: | effs | The effective source object. Needed to add and subtract the derivatives of the singular field when calculating fluxes at the world-tube boundary. |
||
| type(tdc) | :: | time_dc | The time dependent coordinate object. Needed to transform the derivatives of the field when calculating fluxes at the boundary of the time dependent coordinate region. |
|||
| type(rgfb), | intent(in), | dimension(:) | :: | e_lambda | A 1d array of size (2) of real boundary grid functions that contains the characteristic speeds at the element boundary. |
|
| type(rgfb), | intent(in), | dimension(:,:) | :: | e_s | A 2d array of size (2,2) of real boundary grid functions that contains the matrix used to convert from characteristic to evolved variables. |
|
| type(rgfb), | intent(in), | dimension(:,:) | :: | e_sinv | A 2d array of size (2,2) of real boundary grid functions that contains the matrix used to convert from evolved to characteristic variables. |
|
| type(rgf), | intent(in), | dimension(:) | :: | e_coeffs | A 1d array of real grid functions that contains the coefficients for the wave equation. |
|
| type(ref_element), | intent(in) | :: | r_elem | The reference element. |
||
| type(cgf), | intent(inout) | :: | rho_flux | On output contains the component of the characteristic flux. |
||
| type(cgf), | intent(inout) | :: | pi_flux | On output contains the component of the characteristic flux. |
subroutine get_elem_flux ( rho, pi, nmodes, mode, r, effs, time_dc, &
e_lambda, e_s, e_sinv, e_coeffs, r_elem, &
rho_flux, pi_flux )
!! Calculate the characteristic fluxes for all elements for a single mode.
!!
!! Having all of these calculations in a subroutine made it easier to ensure
!! that all local variables are private (for OpenMP parallelization) and to
!! avoid segmentation faults.
use grid, only : drdrho, Tminus_ind, Tplus_ind
use grid_function
use parameters, only : order, use_world_tube, use_particle, &
use_generic_orbit, q_charge
use world_tube, only : wt
implicit none
type(cgf), intent(in) :: rho
!! A complex grid function containing the time derivative variable,
!! \(\Upsilon\).
type(cgf), intent(in) :: pi
!! A complex grid function containing the radial derivative variable,
!! \(\Pi\).
integer(ip), intent(in) :: nmodes
!! The total number of modes.
integer(ip), intent(in) :: mode
!! The mode for which the flux is calculated.
type(rgf), intent(in) :: r
!! A real grid function that contains the Scharzschild radial coordinate.
type(scal_schw_eff), intent(inout) :: effs
!! The effective source object. Needed to add and subtract the derivatives
!! of the singular field when calculating fluxes at the world-tube
!! boundary.
type(tdc) :: time_dc
!! The time dependent coordinate object. Needed to transform the
!! derivatives of the field when calculating fluxes at the boundary of
!! the time dependent coordinate region.
type(rgfb), dimension(:), intent(in) :: e_lambda
!! A 1d array of size (2) of real boundary grid functions that contains
!! the characteristic speeds at the element boundary.
type(rgfb), dimension(:,:), intent(in) :: e_s
!! A 2d array of size (2,2) of real boundary grid functions that contains
!! the matrix used to convert from characteristic to evolved variables.
type(rgfb), dimension(:,:), intent(in) :: e_sinv
!! A 2d array of size (2,2) of real boundary grid functions that contains
!! the matrix used to convert from evolved to characteristic variables.
type(rgf), dimension(:), intent(in) :: e_coeffs
!! A 1d array of real grid functions that contains the coefficients for the
!! wave equation.
type(ref_element), intent(in) :: r_elem
!! The reference element.
type(cgf), intent(inout) :: rho_flux
!! On output contains the \(\Upsilon\) component of the characteristic
!! flux.
type(cgf), intent(inout) :: pi_flux
!! On output contains the \(\Pi\) component of the characteristic flux.
integer(ip), dimension(2) :: ind
integer(ip) :: i, j, l, ne, np
complex(wp), dimension(2,2) :: uint, uext
real(wp), dimension(2,2) :: lambda
real(wp), dimension(2,2,2) :: s, sinv
complex(wp), dimension(2,2) :: flux_int
complex(wp), dimension(1,2) :: psi_sing
complex(wp) :: psi_t, psi_rstar, psi_lambda, psi_xi
logical :: debug_output = .false.
ne = rho%n
ind(1) = 1
! Loop over all elements.
do i = 1, ne
np = rho%elems(i)%order + 1
ind(2) = np
! Get the interior data.
uint(1,1) = rho%elems(i)%var(1)
uint(2,1) = rho%elems(i)%var(np)
uint(1,2) = pi%elems(i)%var(1)
uint(2,2) = pi%elems(i)%var(np)
! Get the exterior data
if (i>1) then
uext(1,1) = rho%elems(i-1)%var(np)
uext(1,2) = pi%elems(i-1)%var(np)
else
! Zero exterior data at the horizon.
uext(1,:) = 0.0_wp
end if
if (i<ne) then
uext(2,1) = rho%elems(i+1)%var(1)
uext(2,2) = pi%elems(i+1)%var(1)
else
! Zero exterior data at Scri+.`
uext(2,:) = 0.0_wp
endif
! If we have a world tube and use an effective source we need to subtract
! and/or add the singular field appropriately.
if ( use_world_tube .and. use_particle ) then
associate ( r1 => r%elems(i)%var(1), &
r2 => r%elems(i)%var(np), &
b1 => wt%boundary_info%elems(i)%bvar(1), &
b2 => wt%boundary_info%elems(i)%bvar(2) )
! Use the worldtube object's boundary information to determine if we are
! at a worldtube boundary.
if ( b1 /= izero ) then
call effs%get_dsingular_dt ( r1, mode, psi_sing(1,1:) )
call effs%get_dsingular_dr ( r1, mode, psi_sing(1,2:) )
! If we are on a generic orbit transform the derivatives of the singular
! field to time dependent coordinates before subtracting or adding (b1 is
! either -1 or 1 depending on whether the current element is inside or
! outside of the worldtube).
if ( use_generic_orbit ) then
psi_t = psi_sing(1,1)
psi_rstar = psi_sing(1,2)
call time_dc%tortoise_to_tdc ( i, 1, psi_t, &
psi_rstar, &
psi_lambda, psi_xi )
psi_sing(1,1) = psi_lambda
psi_sing(1,2) = psi_xi
end if
uext(1,:) = uext(1,:) + q_charge*b1*psi_sing(1,:)
end if
! Use the worldtube object's boundary information to determine if we are
! at a worldtube boundary.
if ( b2 /= izero ) then
call effs%get_dsingular_dt ( r2, mode, psi_sing(1,1:) )
call effs%get_dsingular_dr ( r2, mode, psi_sing(1,2:) )
! If we are on a generic orbit transform the derivatives of the singular
! field to time dependent coordinates before subtracting or adding (b1 is
! either -1 or 1 depending on whether the current element is inside or
! outside of the worldtube).
if ( use_generic_orbit ) then
psi_t = psi_sing(1,1)
psi_rstar = psi_sing(1,2)
call time_dc%tortoise_to_tdc ( i, 2, psi_t, &
psi_rstar, &
psi_lambda, psi_xi )
psi_sing(1,1) = psi_lambda
psi_sing(1,2) = psi_xi
end if
uext(2,:) = uext(2,:) + q_charge*b2*psi_sing(1,:)
end if
end associate
end if
! If we are on the boundary of the time dependent coordinate region, we
! have to transform the exterior data (either to or from time dependent
! coordinates) to match the interior data.
if ( use_generic_orbit ) then
if ( i == Tminus_ind(1) ) then
call time_dc%tdc_to_tortoise ( i+1, 1, uext(2,1), &
uext(2,2), psi_t, &
psi_rstar )
uext(2,1) = psi_t
uext(2,2) = psi_rstar
end if
if ( i == Tminus_ind(2) ) then
call time_dc%tortoise_to_tdc ( i, 1, uext(1,1), &
uext(1,2), psi_lambda, &
psi_xi )
uext(1,1) = psi_lambda
uext(1,2) = psi_xi
end if
if ( i == Tplus_ind(1) ) then
call time_dc%tortoise_to_tdc ( i, 2, uext(2,1), &
uext(2,2), psi_lambda, &
psi_xi )
uext(2,1) = psi_lambda
uext(2,2) = psi_xi
end if
if ( i == Tplus_ind(2) ) then
call time_dc%tdc_to_tortoise ( i-1, 2, uext(1,1), &
uext(1,2), psi_t, &
psi_rstar )
uext(1,1) = psi_t
uext(1,2) = psi_rstar
end if
end if
! Store the eigen-data in local variables.
do j = 1, 2
lambda(:,j) = e_lambda(j)%elems(i)%bvar(:)
do l = 1, 2
s(:,l,j) = e_s(l,j)%elems(i)%bvar(:)
sinv(:,l,j) = e_sinv(l,j)%elems(i)%bvar(:)
end do
! Calculat the interior flux.
flux_int(j,:) = flux ( uint(j,:), &
e_coeffs(2)%elems(i)%var(ind(j)), &
e_coeffs(1)%elems(i)%var(ind(j)) )
end do
! Pass the information to the characteristic flux routine from the
! reference element object and get the numerical flux back (the lift matrix
! has already been applied.
flux_result(:,:,mode) = r_elem%characteristic_flux ( 2, order, uint, uext, &
flux_int, lambda, s, sinv, &
debug_output )
! The only thing left to do is to store the result after multiplying with the
! Jacobian.
do j = 1, np
rho_flux%elems(i)%var(j) = drdrho%elems(i)%var(j)*flux_result(j,1,mode)
pi_flux%elems(i)%var(j) = drdrho%elems(i)%var(j)*flux_result(j,2,mode)
end do
end do
end subroutine get_elem_flux