calculate_potentials_from_nonadiabatic_dfn Subroutine

  1. dist_fn.fpp
  2. dist_fn
  3. calculate_potentials_from_nonadiabatic_dfn

public subroutine calculate_potentials_from_nonadiabatic_dfn(f_in, phi, apar, bpar, gf_lo, local_only)

Uses

Calculates the potentials consistent with the passed nonadiabatic distribution function, f_in. Note this is not what GS2 usually evolves as g/gnew but the adjusted version.

This is closely related to get_init_field which does the same job but works with the g/gnew modified distribution function instead.

Currently we force potentials to zero if they are not included in the simulation. We should provide a method to calculate these fields even if they are not included in the evolution (e.g. for diagnostics, collisions etc.). Perhaps an optional flag here would be enough.

Arguments

Type IntentOptional Attributes Name
complex, intent(in), dimension(-ntgrid:, :, g_lo%llim_proc:) :: f_in
complex, intent(out), dimension(-ntgrid:, :, :) :: phi
complex, intent(out), dimension(-ntgrid:, :, :) :: apar
complex, intent(out), dimension(-ntgrid:, :, :) :: bpar
logical, optional :: gf_lo
logical, optional :: local_only

Contents

Source Code

calculate_potentials_from_nonadiabatic_dfn

Source Code

  subroutine calculate_potentials_from_nonadiabatic_dfn(f_in, phi, apar, bpar, &
       gf_lo, local_only)
    use gs2_layouts, only: g_lo
    use theta_grid, only: ntgrid, bmag
    use run_parameters, only: has_phi, has_apar, has_bpar, beta, tite
    use kt_grids, only: naky, ntheta0, inv_kperp2
    use species, only: spec, has_electron_species, has_ion_species
    use optionals, only: get_option_with_default
    use dist_fn_arrays, only: antot, antota, antotp
    implicit none
    complex, dimension(-ntgrid:, :, g_lo%llim_proc:), intent(in) :: f_in
    complex, dimension(-ntgrid:, :, :), intent(out) :: phi, apar, bpar
    logical, optional :: gf_lo, local_only
    real :: phi_weight
    real, dimension(-ntgrid:ntgrid) :: bpar_weight
    logical :: local_gf_lo, local_local_only
    integer :: it, ik, ig
    integer :: it_llim, it_ulim, ik_llim, ik_ulim

    local_gf_lo = get_option_with_default(gf_lo, .false.)
    local_local_only = get_option_with_default(local_only, .false.)

    ! Get the weighted velocity space integrals of the passed
    ! distribution function.
    if(local_gf_lo) then
       call getan_nogath_from_dfn (f_in, antot, antota, antotp, local_only)
    else
       call getan_from_dfn (f_in, antot, antota, antotp, local_only)
    end if

    if (local_local_only) then
       it_llim = g_lo%it_min ; it_ulim = g_lo%it_max
       ik_llim = g_lo%ik_min ; ik_ulim = g_lo%ik_max
    else
       it_llim = 1 ; it_ulim = ntheta0
       ik_llim = 1 ; ik_ulim = naky
    end if

    if (has_phi) then
       ! We could add in kperp2*poisfac to the weight to retain this feature
       if (.not. has_electron_species(spec) .or. .not. has_ion_species(spec) ) then
          phi_weight = 1.0 / (sum(spec%dens*spec%z*spec%z / spec%temp) + tite)
       else
          phi_weight = 1.0 / sum(spec%dens*spec%z*spec%z / spec%temp)
       end if

       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(it, ik) &
       !$OMP SHARED(ik_llim, ik_ulim, it_llim, it_ulim, phi, antot, phi_weight) &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = ik_llim, ik_ulim
          do it = it_llim, it_ulim
             phi(:, it, ik) = antot(:, it, ik) * phi_weight
          end do
       end do
       !$OMP END PARALLEL DO
    else
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(it, ik) &
       !$OMP SHARED(ik_llim, ik_ulim, it_llim, it_ulim, phi) &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = ik_llim, ik_ulim
          do it = it_llim, it_ulim
             phi(:, it, ik) = 0.0
          end do
       end do
       !$OMP END PARALLEL DO
    end if

    if (has_apar) then
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(it, ik, ig) &
       !$OMP SHARED(ik_llim, ik_ulim, it_llim, it_ulim, ntgrid, inv_kperp2, apar, antota) &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = ik_llim, ik_ulim
          do it = it_llim, it_ulim
             do ig = -ntgrid, ntgrid
                apar(ig, it, ik) = antota(ig, it, ik) * inv_kperp2(ig, it, ik)
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    else
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(it, ik) &
       !$OMP SHARED(ik_llim, ik_ulim, it_llim, it_ulim, apar) &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = ik_llim, ik_ulim
          do it = it_llim, it_ulim
             apar(:, it, ik) = 0.0
          end do
       end do
       !$OMP END PARALLEL DO
    end if

    if (has_bpar) then
       bpar_weight = -beta / bmag**2
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(it, ik) &
       !$OMP SHARED(ik_llim, ik_ulim, it_llim, it_ulim, bpar, bpar_weight, antotp) &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = ik_llim, ik_ulim
          do it = it_llim, it_ulim
             bpar(:, it, ik) = antotp(:, it, ik) * bpar_weight
          end do
       end do
       !$OMP END PARALLEL DO
    else
       !$OMP PARALLEL DO DEFAULT(none) &
       !$OMP PRIVATE(it, ik) &
       !$OMP SHARED(ik_llim, ik_ulim, it_llim, it_ulim, bpar) &
       !$OMP COLLAPSE(2) &
       !$OMP SCHEDULE(static)
       do ik = ik_llim, ik_ulim
          do it = it_llim, it_ulim
             bpar(:, it, ik) = 0.0
          end do
       end do
       !$OMP END PARALLEL DO
    end if

  end subroutine calculate_potentials_from_nonadiabatic_dfn