init_vnew Subroutine

private subroutine init_vnew()

FIXME : Add documentation MAB MAB

Arguments

None

Contents

Source Code


Source Code

  subroutine init_vnew
    use species, only: nspec, spec, is_electron_species
    use le_grids, only: negrid, energy => energy_maxw, w => w_maxw
    use kt_grids, only: naky, ntheta0, kperp2
    use theta_grid, only: ntgrid
    use run_parameters, only: zeff, delt_option_switch, delt_option_auto
    use gs2_time, only: tunits
    use constants, only: pi, sqrt_pi
    use gs2_save, only: init_vnm
    real, dimension(negrid):: hee, hsg, local_energy
    integer :: ik, ie, is, it, ig
    integer :: istatus
    real :: k4max
    real :: vl, vr, dv2l, dv2r
    real, dimension (2) :: vnm_saved

    if (delt_option_switch == delt_option_auto) then
       call init_vnm (vnm_saved, istatus)
       if (istatus == 0) vnm_init = vnm_saved
    endif

    ! Initialise vnmult from the values in run_parameters. This is
    ! either what has been restored from the restart file  if
    ! `delt_option == 'check_restart'` or 1 otherwise.
    if (all(vnmult == -1)) then
       vnmult = vnm_init
    end if

    if (const_v) then
       local_energy = 1.0
    else
       local_energy = energy
    end if

    do ie = 1, negrid
       hee(ie) = exp(-local_energy(ie))/sqrt(pi*local_energy(ie)) &
            + (1.0 - 0.5/local_energy(ie))*erf(sqrt(local_energy(ie)))
       !>MAB
       ! hsg is the G of Hirshman and Sigmar
       ! added to allow for momentum conservation with energy diffusion
       hsg(ie) = hsg_func(sqrt(local_energy(ie)))
       !<MAB
    end do

    if (.not.allocated(vnew)) then
       ! Here the ky dependence is just to allow us to scale by tunits.
       ! This saves some operations in later use at expense of more memory
       allocate (vnew(naky,negrid,nspec))
       allocate (vnew_s(naky,negrid,nspec))
       allocate (vnew_D(naky,negrid,nspec))
       allocate (vnew_E(naky,negrid,nspec))
       allocate (delvnew(naky,negrid,nspec))
    end if

    if (ei_coll_only) then
       vnew_s = 0 ; vnew_D = 0 ; vnew_E = 0 ; delvnew = 0.
       do is = 1, nspec
          if (is_electron_species(spec(is))) then
             do ie = 1, negrid
                vnew(:, ie, is) = spec(is)%vnewk / local_energy(ie)**1.5 &
                        * zeff * 0.5 * tunits(:)
             end do
          else
             vnew(:, :, is) = 0.0
          end if
       end do
    else
       do is = 1, nspec
          ! Don't include following lines, as haven't yet accounted for possibility of timestep changing.
          ! Correct collision frequency if applying collisions only every nth timestep:
          !spec(is)%vnewk = spec(is)%vnewk * float(timesteps_between_collisions)

          do ie = 1, negrid
             vnew_s(:, ie, is) = spec(is)%vnewk / sqrt(local_energy(ie)) &
                  * hsg(ie) * 4.0 * tunits(:)
          end do

          if (is_electron_species(spec(is))) then
             do ie = 1, negrid
                vnew(:, ie, is) = spec(is)%vnewk / local_energy(ie)**1.5 &
                     * (zeff + hee(ie)) * 0.5 * tunits(:)
                vnew_D(:, ie, is) = spec(is)%vnewk / local_energy(ie)**1.5 &
                     * hee(ie) * tunits(:)
             end do
          else
             do ie = 1, negrid
                vnew(:, ie, is) = spec(is)%vnewk / local_energy(ie)**1.5 &
                     * hee(ie) * 0.5 * tunits(:)
                vnew_D(:, ie, is) = 2.0 * vnew(:, ie, is)
             end do
          end if
       end do
    end if

    if (conservative) then !Presumably not compatible with const_v so we use energy
       do is = 1, nspec
          do ie = 2, negrid-1
             vr = 0.5*(sqrt(energy(ie+1)) + sqrt(energy(ie)))
             vl = 0.5*(sqrt(energy(ie  )) + sqrt(energy(ie-1)))
             dv2r = (energy(ie+1) - energy(ie)) / (sqrt(energy(ie+1)) - sqrt(energy(ie)))
             dv2l = (energy(ie) - energy(ie-1)) / (sqrt(energy(ie)) - sqrt(energy(ie-1)))

             vnew_E(:, ie, is) =  spec(is)%vnewk * tunits * &
                  vnew_E_conservative(vr, vl, dv2r, dv2l) / (sqrt_pi * w(ie))
             delvnew(:, ie, is) = spec(is)%vnewk * tunits * &
                  delvnew_conservative(vr, vl) / (sqrt_pi * sqrt(energy(ie)) * w(ie))
          end do

          ! boundary at v = 0
          ie = 1
          vr = 0.5 * (sqrt(energy(ie+1)) + sqrt(energy(ie)))
          vl = 0.0
          dv2r = (energy(ie+1) - energy(ie)) / (sqrt(energy(ie+1)) - sqrt(energy(ie)))
          dv2l = 0.0
          vnew_E(:, ie, is) =  spec(is)%vnewk * tunits * &
               vnew_E_conservative(vr, vl, dv2r, dv2l) / (sqrt_pi * w(ie))
          delvnew(:, ie, is) = spec(is)%vnewk * tunits * &
               delvnew_conservative(vr, vl) / (sqrt_pi * sqrt(energy(ie)) * w(ie))

          ! boundary at v -> infinity
          ie = negrid
          vr = 0.0
          vl = 0.5 * (sqrt(energy(ie)) + sqrt(energy(ie-1)))
          dv2r = 0.0
          dv2l = (energy(ie) - energy(ie-1)) / (sqrt(energy(ie)) - sqrt(energy(ie-1)))
          vnew_E(:, ie, is) =  spec(is)%vnewk * tunits * &
               vnew_E_conservative(vr, vl, dv2r, dv2l) / (sqrt_pi * w(ie))
          delvnew(:, ie, is) = spec(is)%vnewk * tunits * &
               delvnew_conservative(vr, vl) / (sqrt_pi * sqrt(energy(ie)) * w(ie))
       end do
    else
       do is = 1, nspec
          do ie = 2, negrid-1
             vnew_E(:, ie, is) = local_energy(ie) * (vnew_s(:, ie, is) * &
                  (2.0 - 0.5 / local_energy(ie)) - 2.0 * vnew_D(:, ie, is))
             delvnew(:, ie, is) = vnew_s(:, ie, is) - vnew_D(:, ie, is)
          end do
       end do
    end if

       ! add hyper-terms inside collision operator
!BD: Warning!
!BD: For finite magnetic shear, this is different in form from what appears in hyper.f90
!BD: because kperp2 /= akx**2 + aky**2;  there are cross terms that are dropped in hyper.f90
!BD: Warning!
!BD: Also: there is no "grid_norm" option here and the exponent is fixed to 4 for now
    if (hyper_colls) then
       if (.not. allocated(vnewh)) allocate (vnewh(-ntgrid:ntgrid,ntheta0,naky,nspec))
       k4max = (maxval(kperp2))**2
       do is = 1, nspec
          do ik = 1, naky
             do it = 1, ntheta0
                vnewh(:,it,ik,is) = spec(is)%nu_h * kperp2(:,it,ik)**2/k4max
             end do
          end do
       end do
    end if

  contains
    elemental real function vnew_E_conservative(vr, vl, dv2r, dv2l) result(vnE)
      real, intent(in) :: vr, vl, dv2l, dv2r
      vnE = (vl * exp(-vl**2) * dv2l * hsg_func(vl) - vr * exp(-vr**2) * dv2r * hsg_func(vr))
    end function vnew_E_conservative

    elemental real function delvnew_conservative(vr, vl) result(delVn)
      real, intent(in) :: vr, vl
      delVn = (vl * exp(-vl**2) * hsg_func(vl) - vr * exp(-vr**2) * hsg_func(vr))
    end function delvnew_conservative

  end subroutine init_vnew