get_verr Subroutine

public subroutine get_verr(errest, erridx, phi, bpar)

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.

Arguments

Type IntentOptional 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


Contents

Source Code


Source Code

  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