FIXME : Add documentation MAB MAB
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