FIXME : Add documentation
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real, | intent(in), | optional | :: | vnmult_target |
subroutine init_lorentz (vnmult_target)
use le_grids, only: negrid, jend, ng2, nlambda, al, il_is_passing, il_is_wfb
use le_grids, only: setup_trapped_lambda_grids_old_finite_difference, setup_passing_lambda_grids
use gs2_layouts, only: ig_idx, ik_idx, ie_idx, is_idx, it_idx, lz_lo, le_lo
use theta_grid, only: ntgrid
use species, only: has_hybrid_electron_species, spec
use array_utils, only: zero_array
implicit none
real, intent (in), optional :: vnmult_target
integer :: ig, ixi, it, ik, ie, is, je, te2, ile, ilz
real, dimension (:), allocatable :: aa, bb, cc, dd, hh
real, dimension (:, :), allocatable :: local_weights
allocate (aa(nxi_lim), bb(nxi_lim), cc(nxi_lim), dd(nxi_lim), hh(nxi_lim))
if (.not. allocated(pitch_weights)) then
! Start by just copying the existing weights
allocate(local_weights(-ntgrid:ntgrid, nlambda))
local_weights = 0.0
! We have to recaculate the passing weights here as when Radau-Gauss
! grids are used the WFB (il=ng2+1) is treated as both passing and
! trapped. Otherwise we could just copy the passing weights from wl.
call setup_passing_lambda_grids(al, local_weights)
call setup_trapped_lambda_grids_old_finite_difference(al, local_weights)
allocate(pitch_weights(nlambda, -ntgrid:ntgrid))
pitch_weights = transpose(local_weights)
end if
call init_vpdiff
if(use_le_layout) then
if (.not.allocated(c1le)) then
allocate (c1le (nxi_lim, negrid, le_lo%llim_proc:le_lo%ulim_alloc))
allocate (betaale (nxi_lim, negrid, le_lo%llim_proc:le_lo%ulim_alloc))
allocate (qle (nxi_lim, negrid, le_lo%llim_proc:le_lo%ulim_alloc))
if (heating) then
allocate (d1le (nxi_lim, negrid, le_lo%llim_proc:le_lo%ulim_alloc))
allocate (h1le (nxi_lim, negrid, le_lo%llim_proc:le_lo%ulim_alloc))
call zero_array(d1le)
call zero_array(h1le)
end if
vnmult(1) = max(1.0, vnmult(1))
endif
call zero_array(c1le)
call zero_array(betaale)
call zero_array(qle)
if (present(vnmult_target)) then
vnmult(1) = max (vnmult_target, 1.0)
end if
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(ile, ik, it, is, ig, je, te2, ie, aa, bb, cc, dd, hh, ixi) &
!$OMP SHARED(le_lo, jend, special_wfb_lorentz, ng2, negrid, spec, c1le, &
!$OMP d1le, h1le, qle, betaale) &
!$OMP SCHEDULE(static)
do ile = le_lo%llim_proc, le_lo%ulim_proc
ig = ig_idx(le_lo,ile)
ik = ik_idx(le_lo,ile)
it = it_idx(le_lo,ile)
is = is_idx(le_lo,ile)
je = jend(ig)
!if (je <= ng2+1) then
! MRH the above line is likely the original cause of the wfb bug in collisions
! by treating collisions arrays here as if there are only passing particles
! when there are only passing particles plus wfb (and no other trapped particles)
! introduced an inconsistency in the tri-diagonal solve for theta =+/- pi
! Fixed below.
! Here te2 is the total number of unique valid xi at this point.
! For passing particles this is 2*ng2 whilst for trapped we subtract one
! from the number of valid pitch angles due to the "degenerate" duplicate
! vpar = 0 point.
if (il_is_passing(je) .or. (il_is_wfb(je) .and. special_wfb_lorentz) ) then ! MRH
te2 = 2 * ng2
else
te2 = 2 * je - 1
end if
do ie = 1, negrid
call get_lorentz_matrix (aa, bb, cc, dd, hh, ig, ik, it, ie, is)
if (has_hybrid_electron_species(spec)) &
call set_hzero_lorentz_collisions_matrix(aa, bb, cc, je, ik, is)
c1le(:, ie, ile) = cc
if (allocated(d1le)) then
d1le(:, ie, ile) = dd
h1le(:, ie, ile) = hh
end if
qle(1, ie, ile) = 0.0
betaale(1, ie, ile) = 1.0 / bb(1)
do ixi = 1, te2-1
qle(ixi+1, ie, ile) = aa(ixi+1) * betaale(ixi, ie, ile)
betaale(ixi+1, ie, ile) = 1.0 / ( bb(ixi+1) - qle(ixi+1, ie, ile) &
* c1le(ixi, ie, ile) )
end do
qle(te2+1:, ie, ile) = 0.0
betaale(te2+1:, ie, ile) = 0.0
end do
end do
!$OMP END PARALLEL DO
else
if (.not.allocated(c1)) then
allocate (c1(nxi_lim,lz_lo%llim_proc:lz_lo%ulim_alloc))
allocate (betaa(nxi_lim,lz_lo%llim_proc:lz_lo%ulim_alloc))
allocate (ql(nxi_lim,lz_lo%llim_proc:lz_lo%ulim_alloc))
if (heating) then
allocate (d1 (nxi_lim,lz_lo%llim_proc:lz_lo%ulim_alloc))
allocate (h1 (nxi_lim,lz_lo%llim_proc:lz_lo%ulim_alloc))
call zero_array(d1)
call zero_array(h1)
end if
vnmult(1) = max(1.0, vnmult(1))
end if
call zero_array(c1)
call zero_array(betaa)
call zero_array(ql)
if (present(vnmult_target)) then
vnmult(1) = max (vnmult_target, 1.0)
end if
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(ilz, ik, it, is, ig, je, te2, ie, aa, bb, cc, dd, hh, ixi) &
!$OMP SHARED(lz_lo, jend, special_wfb_lorentz, spec, c1, d1, h1, betaa, ql, ng2) &
!$OMP SCHEDULE(static)
do ilz = lz_lo%llim_proc, lz_lo%ulim_proc
is = is_idx(lz_lo,ilz)
ik = ik_idx(lz_lo,ilz)
it = it_idx(lz_lo,ilz)
ie = ie_idx(lz_lo,ilz)
ig = ig_idx(lz_lo,ilz)
je = jend(ig)
!if (je <= ng2+1) then
! MRH the above line is likely the original cause of the wfb bug in collisions
! by treating collisions arrays here as if there are only passing particles
! when there are only passing particles plus wfb (and no other trapped particles)
! introduced an inconsistency in the tri-diagonal solve for theta =+/- pi
! Fixed below
if (il_is_passing(je) .or. (il_is_wfb(je) .and. special_wfb_lorentz) ) then ! MRH
te2 = 2 * ng2
else
te2 = 2 * je - 1
end if
call get_lorentz_matrix (aa, bb, cc, dd, hh, ig, ik, it, ie, is)
if (has_hybrid_electron_species(spec)) &
call set_hzero_lorentz_collisions_matrix(aa, bb, cc, je, ik, is)
c1(:, ilz) = cc
if (allocated(d1)) then
d1(:, ilz) = dd
h1(:, ilz) = hh
end if
ql(1, ilz) = 0.0
betaa(1, ilz) = 1.0 / bb(1)
do ixi = 1, te2-1
ql(ixi+1, ilz) = aa(ixi+1) * betaa(ixi, ilz)
betaa(ixi+1, ilz) = 1.0 / (bb(ixi+1) - ql(ixi+1, ilz) * c1(ixi, ilz))
end do
ql(te2+1:, ilz) = 0.0
c1(te2+1:, ilz) = 0.0
betaa(te2+1:, ilz) = 0.0
end do
!$OMP END PARALLEL DO
end if
deallocate (aa, bb, cc, dd, hh)
end subroutine init_lorentz