Determine trapped pitch angle weights using "new" polynomial interpolation.
The new method uses Lagrange polynomial interpolation to provide
an accurate approximation of the integration over the trapped
region. The old method effectively uses the trapezium rule to
integrate. In both cases the integration variable is the parallel
velocity spanning the range [-v_||wfb, v||_wfb]. Both methods
give zero weight to the wfb at the locations where B == Bmax.
Generally the old method appears to be accurate to at most single
precision, whilst the new method is more accurate, reaching the
double precision round off level even with just a few trapped
pitch angles. It can be noted that the the old method can reach
the same level of accuracy at some theta grid points but not
all. Further work may be able to optimise the grid to avoid these
locations. There are a number of reasons why we may not want to
make new_trap_int default to true, including :
1. It is not compatible with radau_gauss_grid = T (which is the
default).
2. One might expect high order Lagrange interpolation to be
unstable so may anticipate that the new method will not work
well at higher ntheta. Empirically it seems this does not
become an issue here. Further, the input nmax can be used
to limit the maximum Lagrange order used by splitting the
domain.
3. For large ntheta the new method leads to slow
initialisation. This is due to the calculation of the error
coefficients for use with get_verr. This could perhaps be
optmised.
Both points 2 and 3 could be avoided to some extent if the
trapped pitch angle grid resolution is decoupled from the theta
grid resolution. It may be possible to make the new approach
compatible with radau_gauss_grid or we may decide
radau_gauss_grid should default to false.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real, | intent(in), | dimension(:) | :: | lambda_grid |
The lambda grid points. |
|
| real, | intent(inout), | dimension(-ntgrid:, :) | :: | weights |
The integration weights for the lambda grid. Note this is
intent |
subroutine setup_trapped_lambda_grids_new_trap_int(lambda_grid, weights)
use theta_grid, only: ntgrid, bmag, bmax
use file_utils, only: open_output_file, close_output_file, error_unit
use mp, only: proc0
implicit none
!> The lambda grid points.
real, dimension(:), intent(in) :: lambda_grid
!> The integration weights for the lambda grid. Note this is
!> intent `in out` as we only set a portion of the grid so may
!> want to keep the rest of the input values.
real, dimension(-ntgrid:, :), intent(in out) :: weights
real, dimension (:), allocatable :: ytmp, yb, wb
real, dimension (:,:), allocatable :: wberr
integer :: npts, ntrap
real :: llim, ulim, wgt_tmp
logical :: eflag
integer :: ig, il, ndiv, divmax
! Initialise the error flag
eflag = .false.
! max number of trapped particles (occurs at outboard midplane)
ntrap = nlambda - ng2
! wlterr contains weights for less accurate trapped integrals (for error estimation)
allocate(wlterr(-ntgrid:ntgrid,nlambda,ntrap))
wlterr = 0.0
! Find the integration weights for each theta grid point
do ig = -ntgrid, ntgrid
! First we count how many trapped lambda bounce points we
! have to consider at this point. In other words how many
! bounce points are not forbidden for the current theta
! grid point.
! npts is the number of lambda_grid values in the trapped integral (varies with theta)
npts = 0
do il = ng2+1, nlambda
if (1.0 - lambda_grid(il)*bmag(ig) > -bouncefuzz) then
npts = npts + 1
end if
end do
! If there are any valid bounce points then we need to
! calculate the weights. If npts == 1 then we probably also expect
! the integration weights to be 0 so we could probably skip all of
! this work for npts = 1 as well.
if (npts > 0) then
! ytmp is an array containing pitch angle grid points (for vpa >= 0)
! These are the v||/v points for each lambda between wfb (sets maximum
! v||/v here) and the lambda which bounces at this point (v||/v = 0)
! yb is an array containing the full set of [-v||/v, v||/v] grid points.
! Constructed directly from ytmp by copying and flipping sign+order
! wb is an array containing the integration weights calculated for the
! given v||/v grid.
allocate(ytmp(npts), yb(2*npts-1), wb(2*npts-1))
! wberr is an array used to hold the weights used in the error estimation
! in [[trap_error]] called by [[get_verr]]
allocate(wberr(2*npts-1,npts))
ytmp = 0.0; yb = 0.0; wb = 0.0; wberr = 0.0
! loop computes transformed variable of integration (v||/v)
do il = ng2+1, ng2+npts
ytmp(il-ng2) = sqrt(max(1 - lambda_grid(il)*bmag(ig), 0.0))
end do
! define array (yb) with pitch-angle gridpts
! corresponding to both positive and negative vpa
if (npts > 1) yb(:npts-1) = -ytmp(:npts-1)
yb(npts:) = ytmp(npts:1:-1)
! get grid point weights for trapped particle integral
! Note : Here ulim and llim are the upper and lower vpar
! values for the valid trapped pitch angles at this theta
! location. In other words ulim == v_||_wfb(theta(ig)) =
! sqrt(max(1.0-lambda_grid(ng2+1)*bmag(ig),0.0)) = yb(2*npts-1)
! We could probably replace ulim and llim with yb(2*npts-1)
! and yb(1) respectively.
ulim = sqrt(max(1.0-bmag(ig)/bmax,0.0))
llim = -ulim
! Call get_weights to find the integration weights for the v||/v grid
! Note we don't call this if ulim == 0 as the grid then has no extent.
! This situation can arise when bmag(ig) == bmax where we would find
! npts = 1, corresponding to just the wfb being a valid lambda. This
! means that at this location the wfb has zero integration weight.
! Note we currently ignore the error flag and other returned values other
! than the weights. We should probably check this to at least warn the
! user if something looks badly behaved.
if (ulim > 0.0) then
if (new_trap_int_split) then
! We split the weights calculation into two symmetric domains from
! -v||_wfb to 0 and from 0 to v||_wfb. This has been shown to have
! favourable properties when compared to treating the full -v||_wfb
! to v||_wfb domain in a single call. In particular significantly more
! accurate results when integrating certain functions (see mod vpa test
! in the le_grids_integrate unit tests).
call get_weights (nmax, llim, 0., yb(:npts), wb(:npts), ndiv, divmax, eflag)
if (eflag .and. proc0) then
write(error_unit(), '("Error flag set by first call to get_weights in setup_trapped_lambda_grids")')
end if
! Store the v|| = 0 weight for the lower domain as this will be
! clobbered by the subsequent get_weights call so we need to add
! this on afterwards
wgt_tmp = wb(npts)
call get_weights (nmax, 0., ulim, yb(npts:), wb(npts:), ndiv, divmax, eflag)
if (eflag .and. proc0) then
write(error_unit(), '("Error flag set by second call to get_weights in setup_trapped_lambda_grids")')
end if
wb(npts) = wb(npts) + wgt_tmp
else
call get_weights (nmax, llim, ulim, yb, wb, ndiv, divmax, eflag)
if (eflag .and. proc0) then
write(error_unit(), '("Error flag set by first call to get_weights in setup_trapped_lambda_grids")')
end if
end if
end if
! Calculate the weights for use in the trapped
! integration error estimation code. This can be quite
! expensive so we would like to to only call this if the
! error estimation code is active (i.e. if we're going to
! call [[get_verr]]).
if (npts > 1) call get_trapped_lambda_grid_error_estimate_weights(yb, llim, ulim, npts, wberr)
! Convert from v||/v weights to lambda
! weights. Essentially just sum up the postive and
! negtive v|| points corresponding to the same lambda
! point. Note we skip the v|| == 0 point currently to ensure we
! don't double count.
if (npts > 1) then
do il = ng2+1, ng2+npts-1
! take into account possible asymmetry of weights about xi = 0
! due to unequal # of grid points per integration interval
! Whilst the v||/v grid should always be symmetric
! the algorithm in [[get_weights]] could
! potentially lead to asymmetry in the resulting
! weights. This may be an indication that things
! aren't behaving well.
weights(ig,il) = wb(il-ng2) + wb(2*npts-il+ng2)
wlterr(ig,il,:npts) = wberr(il-ng2,:) + wberr(2*npts-il+ng2,:)
end do
end if
! avoid double counting of gridpoint at vpa=0
weights(ig,ng2+npts) = wb(npts)
wlterr(ig,ng2+npts,:npts) = wberr(npts,:)
deallocate(ytmp, yb, wb, wberr)
end if
end do
end subroutine setup_trapped_lambda_grids_new_trap_int