Error estimate obtained by comparing standard integral with less-accurate integral
Estimate of the (1) absolute and (2) relative errors resulting from velocity space integrals in the calculation of the following quantities in the given dimensions: (1) k phi, energy (2) k phi, untrapped pitch angles (3) k phi, trapped pitch angles, (4) k apar, energy, (5) k apar, untrapped angles. Relative errors should be < 0.1.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real, | intent(out), | dimension (:,:) | :: | errest |
Estimated error. FIXME: Make this fixed size. |
|
integer, | intent(out), | dimension (:,:) | :: | erridx |
Indices of maximum error. FIXME: Make this fixed size |
|
complex, | intent(in), | dimension (-ntgrid:,:,:) | :: | phi |
Electrostatic potential and parallel magnetic field |
|
complex, | intent(in), | dimension (-ntgrid:,:,:) | :: | bpar |
Electrostatic potential and parallel magnetic field |
subroutine get_verr (errest, erridx, phi, bpar)
use mp, only: broadcast
use le_grids, only: integrate_species
use le_grids, only: eint_error, trap_error, lint_error, wdim
use le_grids, only: ng2, nlambda, new_trap_int, grid_has_trapped_particles
use theta_grid, only: ntgrid
use kt_grids, only: ntheta0, naky, aky, akx
use species, only: nspec, spec
use dist_fn_arrays, only: gnew, aj0, vpa, g_adjust, g_work, to_g_gs2, from_g_gs2
use run_parameters, only: has_phi, has_apar, beta
use gs2_layouts, only: g_lo
use collisions, only: init_lorentz, init_ediffuse, init_lorentz_conserve
use collisions, only: etol, ewindow, etola, ewindowa, init_diffuse_conserve
use collisions, only: vnmult, vary_vnew
use array_utils, only: zero_array
implicit none
integer :: ig, it, ik, iglo, isgn, ntrap
!> Indices of maximum error.
!> FIXME: Make this fixed size
integer, dimension (:,:), intent (out) :: erridx
!> Estimated error.
!> FIXME: Make this fixed size.
real, dimension (:,:), intent (out) :: errest
!> Electrostatic potential and parallel magnetic field
complex, dimension (-ntgrid:,:,:), intent (in) :: phi, bpar
complex, dimension (:,:,:), allocatable :: phi_app, apar_app
complex, dimension (:,:,:,:), allocatable :: phi_e, phi_l, phi_t
complex, dimension (:,:,:,:), allocatable :: apar_e, apar_l, apar_t
real, dimension (:), allocatable :: wgt, errtmp
integer, dimension (:), allocatable :: idxtmp
real :: vnmult_target
logical :: trap_flag=.true. !Value doesn't really matter as just used in present() checks
allocate(wgt(nspec))
allocate(errtmp(2))
allocate(idxtmp(3))
!This should probably be something like if(first) then ; call broadcast ; first=.false. ; endif
!as we only deallocate kmax during finish_dist_fn.
!Really wdim should just be broadcast when we first calculate it -- currently slightly complicated
!as calculation triggered in gs2_diagnostics and only called by proc0.
if (.not. allocated(kmax)) call broadcast (wdim) ! Odd-looking statement. BD
if (has_phi) then
allocate(phi_app(-ntgrid:ntgrid,ntheta0,naky))
allocate(phi_e(-ntgrid:ntgrid,ntheta0,naky,wdim))
allocate(phi_l(-ntgrid:ntgrid,ntheta0,naky,ng2))
end if
if (has_apar) then
allocate(apar_app(-ntgrid:ntgrid,ntheta0,naky))
allocate(apar_e(-ntgrid:ntgrid,ntheta0,naky,wdim))
allocate(apar_l(-ntgrid:ntgrid,ntheta0,naky,ng2))
end if
! first call to g_adjust converts gyro-averaged dist. fn. (g)
! into nonadiabatic part of dist. fn. (h)
call g_adjust (gnew, phi, bpar, direction = from_g_gs2)
! take gyro-average of h at fixed total position (not g.c. position)
if (has_phi) then
do iglo = g_lo%llim_proc, g_lo%ulim_proc
! TEMP FOR TESTING -- MAB
! il = il_idx(g_lo,iglo)
do isgn = 1, 2
do ig=-ntgrid, ntgrid
! TEMP FOR TESTING -- MAB
g_work(ig,isgn,iglo) = aj0(ig,iglo)*gnew(ig,isgn,iglo)
! g_work(ig,isgn,iglo) = sqrt(max(0.0,1.0-bmag(ig)*al(il)))
! g_work = 1.
end do
end do
end do
wgt = spec%z*spec%dens
call integrate_species (g_work, wgt, phi_app)
! integrates dist fn of each species over v-space
! after dropping an energy grid point and returns
! phi_e, which contains the integral approximations
! to phi for each point dropped
call eint_error (g_work, wgt, phi_e)
! integrates dist fn of each species over v-space
! after dropping an untrapped lambda grid point and returns phi_l.
! phi_l contains ng2 approximations for the integral over lambda that
! come from dropping different pts from the gaussian quadrature grid
call lint_error (g_work, wgt, phi_l)
! next loop gets error estimate for trapped particles, if there are any
if (grid_has_trapped_particles() .and. new_trap_int) then
ntrap = nlambda - ng2
allocate(phi_t(-ntgrid:ntgrid,ntheta0,naky,ntrap))
call zero_array(phi_t)
call trap_error (g_work, wgt, phi_t)
end if
end if
if (has_apar) then
do iglo = g_lo%llim_proc, g_lo%ulim_proc
do isgn = 1, 2
do ig=-ntgrid, ntgrid
g_work(ig,isgn,iglo) = aj0(ig,iglo)*vpa(ig,isgn,iglo)*gnew(ig,isgn,iglo)
end do
end do
end do
wgt = 2.0*beta*spec%z*spec%dens*sqrt(spec%temp/spec%mass)
call integrate_species (g_work, wgt, apar_app)
call eint_error (g_work, wgt, apar_e)
call lint_error (g_work, wgt, apar_l)
if (grid_has_trapped_particles() .and. new_trap_int) then
ntrap = nlambda - ng2
allocate(apar_t(-ntgrid:ntgrid,ntheta0,naky,ntrap))
call zero_array(apar_t)
call trap_error (g_work, wgt, apar_t)
end if
end if
! second call to g_adjust converts from h back to g
call g_adjust (gnew, phi, bpar, direction = to_g_gs2)
if (.not. allocated(kmax)) then
allocate (kmax(ntheta0, naky))
do ik = 1, naky
do it = 1, ntheta0
kmax(it,ik) = max(akx(it),aky(ik))
end do
end do
end if
errest = 0.0
erridx = 0
if (has_phi) then
call estimate_error (phi_app, phi_e, kmax, errtmp, idxtmp)
errest(1,:) = errtmp
erridx(1,:) = idxtmp
call estimate_error (phi_app, phi_l, kmax, errtmp, idxtmp)
errest(2,:) = errtmp
erridx(2,:) = idxtmp
if (grid_has_trapped_particles() .and. new_trap_int) then
call estimate_error (phi_app, phi_t, kmax, errtmp, idxtmp, trap_flag)
errest(3,:) = errtmp
erridx(3,:) = idxtmp
deallocate (phi_t)
end if
end if
if (has_apar) then
call estimate_error (apar_app, apar_e, kmax, errtmp, idxtmp)
errest(4,:) = errtmp
erridx(4,:) = idxtmp
call estimate_error (apar_app, apar_l, kmax, errtmp, idxtmp)
errest(5,:) = errtmp
erridx(5,:) = idxtmp
end if
if (vary_vnew) then
vnmult_target = vnmult(2)
if (errest(1,2) > etol + ewindow .or. errest(4,2) > etola + ewindowa) then
call get_vnewk (vnmult(2), vnmult_target, increase)
else if (errest(1,2) < etol - ewindow .and. errest(4,2) < etola - ewindowa) then
call get_vnewk (vnmult(2), vnmult_target, decrease)
end if
call init_ediffuse (vnmult_target)
call init_diffuse_conserve
end if
if (vary_vnew) then
vnmult_target = vnmult(1)
if (errest(2,2) > etol + ewindow .or. errest(3,2) > etol + ewindow &
.or. errest(5,2) > etola + ewindowa) then
call get_vnewk (vnmult(1), vnmult_target, increase)
else if (errest(2,2) < etol - ewindow .and. errest(5,2) < etola - ewindowa .and. &
(errest(3,2) < etol - ewindow .or. .not.trap_flag)) then
call get_vnewk (vnmult(1), vnmult_target, decrease)
end if
call init_lorentz (vnmult_target)
call init_lorentz_conserve
end if
deallocate (wgt, errtmp, idxtmp)
if (has_phi) deallocate(phi_app, phi_e, phi_l)
if (has_apar) deallocate(apar_app, apar_e, apar_l)
end subroutine get_verr