FIXME : Add documentation
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real, | intent(out), | dimension (:) | :: | aa | ||
real, | intent(out), | dimension (:) | :: | bb | ||
real, | intent(out), | dimension (:) | :: | cc | ||
real, | intent(out), | dimension (:) | :: | dd | ||
real, | intent(out), | dimension (:) | :: | hh | ||
integer, | intent(in) | :: | ig | |||
integer, | intent(in) | :: | ik | |||
integer, | intent(in) | :: | it | |||
integer, | intent(in) | :: | ie | |||
integer, | intent(in) | :: | is |
subroutine get_lorentz_matrix (aa, bb, cc, dd, hh, ig, ik, it, ie, is)
use species, only: spec
use le_grids, only: al, energy => energy_maxw, xi, ng2
use le_grids, only: jend, al, il_is_passing, il_is_wfb
use gs2_time, only: code_dt, tunits
use kt_grids, only: kperp2
use theta_grid, only: bmag
use warning_helpers, only: is_zero
implicit none
real, dimension (:), intent (out) :: aa, bb, cc, dd, hh
integer, intent (in) :: ig, ik, it, ie, is
integer :: il, je, te, te2, teh
real :: slb0, slb1, slb2, slbl, slbr, vhyp, vn, vnh, vnc, ee, deltaxi
je = jend(ig)
!
!CMR, 17/2/2014:
! te, te2, teh: indices in xi, which runs from +1 -> -1.
! te : index of minimum xi value >= 0.
! te2 : total #xi values = index of minimum xi value (= -1)
! teh : index of minimum xi value > 0.
! teh = te if no bouncing particle at this location
! OR teh = te-1 if there is a bouncing particle
!
if (il_is_passing(je) .or. (il_is_wfb(je) .and. special_wfb_lorentz)) then
!CMRDDGC, 17/2/2014:
! This clause is appropriate for Lorentz collisons with
! SPECIAL (unphysical) treatment of wfb at its bounce point
te = ng2
te2 = 2*ng2
teh = ng2
else
!CMRDDGC, 17/2/2014:
! This clause is appropriate for Lorentz collisons with
! STANDARD treatment of wfb at its bounce point
te = je
te2 = 2*je-1
teh = je-1
end if
if (collision_model_switch == collision_model_lorentz_test) then
vn = vnmult(1) * abs(spec(is)%vnewk) * tunits(ik)
vnc = 0.
vnh = 0.
else
if (hyper_colls) then
vhyp = vnewh(ig, it, ik, is)
else
vhyp = 0.0
end if
! vnc and vnh only needed when heating is true
vnc = vnmult(1) * vnew(ik, ie, is)
if (hypermult) then
vn = vnmult(1) * vnew(ik, ie, is) * (1 + vhyp)
vnh = vhyp * vnc
else
vn = vnmult(1) * vnew(ik, ie, is) + vhyp
vnh = vhyp
end if
end if
aa = 0.0 ; bb = 0.0 ; cc = 0.0 ; dd = 0.0 ; hh = 0.0
select case (lorentz_switch)
case (lorentz_scheme_default)
do il = 2, te-1
slb0 = safe_sqrt(1.0 - bmag(ig) * al(il-1))
slb1 = safe_sqrt(1.0 - bmag(ig) * al(il))
slb2 = safe_sqrt(1.0 - bmag(ig) * al(il+1))
slbl = (slb1 + slb0) * 0.5 ! xi(j-1/2)
slbr = (slb1 + slb2) * 0.5 ! xi(j+1/2)
ee = 0.5 * energy(ie)*(1 + slb1**2) * kperp2(ig,it,ik) * cfac &
/ (bmag(ig) * spec(is)%zstm)**2
! coefficients for tridiagonal matrix:
cc(il) = 2.0 * vn * code_dt * (1 - slbr**2) / &
(pitch_weights(il, ig) * (slb2 - slb1))
aa(il) = 2.0 * vn * code_dt * (1 - slbl**2) / &
(pitch_weights(il, ig) * (slb1 - slb0))
bb(il) = 1.0 - (aa(il) + cc(il)) + ee * vn * code_dt
! coefficients for entropy heating calculation
if (heating) then
dd(il) = vnc * (-2.0 * (1 - slbr**2) / &
(pitch_weights(il, ig) * (slb2 - slb1)) + ee)
hh(il) = vnh * (-2.0 * (1 - slbr**2) / &
(pitch_weights(il, ig) * (slb2 - slb1)) + ee)
end if
end do
! boundary at xi = 1
slb1 = safe_sqrt(1.0 - bmag(ig) * al(1))
slb2 = safe_sqrt(1.0 - bmag(ig) * al(2))
slbr = (slb1 + slb2) * 0.5
ee = 0.5 * energy(ie) * (1 + slb1**2) * kperp2(ig,it,ik) * cfac &
/ (bmag(ig) * spec(is)%zstm)**2
cc(1) = 2.0 * vn * code_dt * (1.0 - slbr**2) &
/ (pitch_weights(1, ig) * (slb2 - slb1))
aa(1) = 0.0
bb(1) = 1.0 - cc(1) + ee * vn * code_dt
if (heating) then
dd(1) = vnc * (-2.0 * (1 - slbr**2) &
/ (pitch_weights(1, ig) * (slb2 - slb1)) + ee)
hh(1) = vnh * (-2.0 * (1 - slbr**2) &
/ (pitch_weights(1, ig) * (slb2 - slb1)) + ee)
end if
! boundary at xi = 0
il = te
slb0 = safe_sqrt(1.0 - bmag(ig) * al(il-1))
if (te == ng2) then
slb1 = safe_sqrt(1.0 - bmag(ig) * al(il))
slb2 = -slb1
else
slb1 = 0.0
slb2 = -slb0
end if
slbl = (slb1 + slb0) * 0.5
slbr = (slb1 + slb2) * 0.5
ee = 0.5 * energy(ie) * (1 + slb1**2) * kperp2(ig,it,ik) * cfac &
/ (bmag(ig) * spec(is)%zstm)**2
!CMR, 6/3/2014:
! STANDARD treatment of pitch angle scattering must resolve T-P boundary.
! NEED special_wfb= .false. to resolve T-P boundary at wfb bounce point
! (special_wfb= .true. AVOIDS TP boundary at wfb bounce point)
!
! Original code (pre-r2766) used eq.(42) Barnes et al, Phys Plasmas 16, 072107
! (2009), with pitch angle weights to enhance accuracy in derivatives.
! NB THIS FAILS at wfb bounce point, giving aa=cc=infinite,
! because weight wl(ig,il)=0 at wfb bounce point.
! UNPHYSICAL as d/dxi ((1-xi^2)g) IS NOT resolved numerically for wl=0.
! MUST accept limitations of the grid resolution and USE FINITE coefficients!
! FIX here by setting a FINITE width of the trapped region at wfb B-P
! deltaxi=xi(ig,ng2)-xi(ig,ng2+2)
! ASIDE: NB deltaxi=wl is actually 2*spacing in xi !!!
! which explains upfront factor 2 in definition of aa, cc
deltaxi = pitch_weights(il, ig)
if (te == je) deltaxi = 2 * deltaxi ! MRH vpar = 0 fix
! MRH appropriate when endpoint (te) is a vpar = 0 point (je)
! factor of 2 required above because xi=0 is a repeated point
! on vpar > 0, vpar < 0 grids, and hence has half the weight associated
! to it at each appearance compared to what the weight should be
! as calculated in the continuum of points xi = (-1,1)
if ((.not. special_wfb_lorentz) .and. (is_zero(deltaxi)) .and. il_is_wfb(il)) then
deltaxi = xi(ng2, ig) - xi(ng2 + 2, ig)
endif
cc(il) = 2.0 * vn * code_dt * (1 - slbr**2) / (deltaxi * (slb2 - slb1))
aa(il) = 2.0 * vn * code_dt * (1 - slbl**2) / (deltaxi * (slb1 - slb0))
bb(il) = 1.0 - (aa(il) + cc(il)) + ee * vn * code_dt
if (heating) then
dd(il) = vnc * (-2.0 * (1.0 - slbr**2) / (deltaxi * (slb2 - slb1)) + ee)
hh(il) = vnh * (-2.0 * (1.0 - slbr**2) / (deltaxi * (slb2 - slb1)) + ee)
end if
!CMRend
case (lorentz_scheme_old)
do il = 2, te-1
slb0 = safe_sqrt(1.0 - bmag(ig) * al(il-1))
slb1 = safe_sqrt(1.0 - bmag(ig) * al(il))
slb2 = safe_sqrt(1.0 - bmag(ig) * al(il+1))
slbl = (slb1 + slb0) * 0.5 ! xi(j-1/2)
slbr = (slb1 + slb2) * 0.5 ! xi(j+1/2)
ee = 0.5 * energy(ie) * (1 + slb1**2) * kperp2(ig,it,ik) * cfac &
/ (bmag(ig) * spec(is)%zstm)**2
! coefficients for tridiagonal matrix:
cc(il) = -vn * code_dt * (1 - slbr**2) / ((slbr - slbl) * (slb2 - slb1))
aa(il) = -vn * code_dt * (1 - slbl**2) / ((slbr - slbl) * (slb1 - slb0))
bb(il) = 1.0 - (aa(il) + cc(il)) + ee * vn * code_dt
! coefficients for entropy heating calculation
if (heating) then
dd(il) = vnc * ((1 - slbr**2) / ((slbr - slbl) * (slb2 - slb1)) + ee)
hh(il) = vnh * ((1 - slbr**2) / ((slbr - slbl) * (slb2 - slb1)) + ee)
end if
end do
! boundary at xi = 1
slb0 = 1.0
slb1 = safe_sqrt(1.0 - bmag(ig) * al(1))
slb2 = safe_sqrt(1.0 - bmag(ig) * al(2))
slbl = (slb1 + slb0) * 0.5
slbr = (slb1 + slb2) * 0.5
ee = 0.5 * energy(ie) * (1 + slb1**2) * kperp2(ig,it,ik) * cfac &
/ (bmag(ig) * spec(is)%zstm)**2
cc(1) = -vn * code_dt * (-1.0 - slbr) / (slb2 - slb1)
aa(1) = 0.0
bb(1) = 1.0 - (aa(1) + cc(1)) + ee * vn * code_dt
if (heating) then
dd(1) = vnc * ((1 - slbr**2) / ((slbr - slbl) * (slb2 - slb1)) + ee)
hh(1) = vnh * ((1 - slbr**2) / ((slbr - slbl) * (slb2 - slb1)) + ee)
end if
! boundary at xi = 0
il = te
slb0 = safe_sqrt(1.0 - bmag(ig) * al(il-1))
if (te == ng2) then
slb1 = safe_sqrt(1.0 - bmag(ig) * al(il))
slb2 = -slb1
else
slb1 = 0.0
slb2 = -slb0
end if
slbl = (slb1 + slb0) * 0.5
slbr = (slb1 + slb2) * 0.5
ee = 0.5 * energy(ie) * (1 + slb1**2) * kperp2(ig,it,ik) * cfac &
/ (bmag(ig) * spec(is)%zstm)**2
cc(il) = -vn * code_dt * (1 - slbr**2) / ((slbr - slbl) * (slb2 - slb1))
aa(il) = -vn * code_dt * (1 - slbl**2) / ((slbr - slbl) * (slb1 - slb0))
bb(il) = 1.0 - (aa(il) + cc(il)) + ee * vn * code_dt
if (heating) then
dd(il) = vnc * ((1 - slbr**2) / ((slbr - slbl) * (slb2 - slb1)) + ee)
hh(il) = vnh * ((1 - slbr**2) / ((slbr - slbl) * (slb2 - slb1)) + ee)
end if
end select
! assuming symmetry in xi, fill in the rest of the arrays.
aa(te+1:te2) = cc(teh:1:-1)
bb(te+1:te2) = bb(teh:1:-1)
cc(te+1:te2) = aa(teh:1:-1)
if (heating) then
dd(te+1:te2) = dd(teh:1:-1)
hh(te+1:te2) = hh(teh:1:-1)
end if
end subroutine get_lorentz_matrix