Routine for getting the weights used in estimating the error on trapped particle integrals as a part of trap_error
This routine is all to calculate wberr/wlterr which is just used in trap_error, which in turn is only used in get_verr. This is used as a part of the vary_vnew adaptive collisionality code triggered through the diagnostics. This routine is rather expensive for high theta resolution - measured in the order of 4-5 minutes for ntheta=256. As such we might want to consider skipping this if the variable vnewk code is not active. As we do this for each npt and each theta, and npt ~ntheta/2 we might expect this to scale quadratically with ntheta. In practice it seems to scale somewhat worse than this. This may in part be due to an increase in the number of iterations required in get_weights such that this may scale closer to ntheta cubed.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real, | intent(in), | dimension(:) | :: | yb | ||
| real, | intent(in) | :: | llim | |||
| real, | intent(in) | :: | ulim | |||
| integer, | intent(in) | :: | npts | |||
| real, | intent(inout), | dimension(:, :) | :: | wberr |
subroutine get_trapped_lambda_grid_error_estimate_weights(yb, llim, ulim, npts, wberr)
implicit none
real, dimension(:), intent(in) :: yb
real, intent(in) :: llim, ulim
integer, intent(in) :: npts
real, dimension(:, :), intent(in out) :: wberr
real, dimension (:), allocatable :: yberr, wberrtmp
integer :: ix
logical :: eflag
integer :: divmaxerr, ndiverr
! Can't get error estimate for npts = 0 or 1 as we don't have any points
! that we can drop.
if (npts > 1) then
do ix=1,npts
if (ix == 1) then
! drop the first and last grid points from the integral
allocate (yberr(2*npts-3),wberrtmp(2*npts-3))
yberr = 0.0; wberrtmp = 0.0
yberr = yb(2:2*npts-2)
else if (ix == npts) then
! drop the vpa=0 grid point from the integral
allocate (yberr(2*npts-2),wberrtmp(2*npts-2))
yberr = 0.0; wberrtmp = 0.0
yberr(:npts-1) = yb(:npts-1)
yberr(npts:) = yb(npts+1:)
else
! drop the grid points corresponding to ix and its negative from the integral
allocate (yberr(2*npts-3),wberrtmp(2*npts-3))
yberr = 0.0; wberrtmp = 0.0
yberr(:ix-1) = yb(:ix-1)
yberr(ix:2*npts-ix-2) = yb(ix+1:2*npts-ix-1)
yberr(2*npts-ix-1:) = yb(2*npts-ix+1:)
end if
call get_weights (nmax, llim, ulim, yberr, wberrtmp, ndiverr, divmaxerr, eflag)
! insert a weight of zero into indices corresponding to ix and its conjugate
if (ix == 1) then
wberr(2:2*npts-2,1) = wberrtmp
else if (ix == npts) then
wberr(:npts-1,npts) = wberrtmp(:npts-1)
wberr(npts+1:,npts) = wberrtmp(npts:)
else
wberr(:ix-1,ix) = wberrtmp(:ix-1)
wberr(ix+1:2*npts-ix-1,ix) = wberrtmp(ix:2*npts-ix-2)
wberr(2*npts-ix+1,ix) = wberrtmp(2*npts-ix-1)
end if
deallocate (yberr,wberrtmp)
end do
end if
end subroutine get_trapped_lambda_grid_error_estimate_weights