Calculates the "new" distribution function for the passed source calculated using the passed fields and current distribution in g.
Also calculates the homogeneous solution in some cases.
Note that further work may be required to produce the final distribution function. In particular with link boundary conditions we may need to add in the homogeneous solution to ensure continuity across linked boundaries.
The routine is roughly organised as: 1. Calculate constants 2. Loop over domain - Calculate flags - Set initial conditions - Set parallel boundary conditions - Solve for homogeneous and inhomogeneous solution - Ensure continuity (e.g. periodicity in theta and in sigma at bounce points) 3. Enforce parity if required. !!!!!!!!!!!!!!!!!!!!!!!!! time advance vpar < 0 ! !!!!!!!!!!!!!!!!!!!!!!!!! inhomogeneous part: gnew r=ainv=0 if forbid(ig,il) or forbid(ig+1,il), so gnew=0 in forbidden region and at upper bounce point !!!!!!!!!!!!!!!!!!!!!!!!! time advance vpar > 0 ! !!!!!!!!!!!!!!!!!!!!!!!!! First set BCs for trapped particles at lower bounce point
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
complex, | intent(in), | dimension (-ntgrid:,:,:) | :: | phinew | ||
complex, | intent(in), | dimension (-ntgrid:,:,:) | :: | bparnew | ||
complex, | intent(inout), | dimension (-ntgrid:,:,g_lo%llim_proc:) | :: | source |
subroutine invert_rhs_1 (phinew, bparnew, source)
use dist_fn_arrays, only: gnew, set_h_zero, g
use theta_grid, only: ntgrid
use le_grids, only: forbid, is_passing_hybrid_electron, passing_wfb, trapped_wfb, &
il_is_wfb, is_ttp, il_is_passing, il_is_trapped, &
is_upper_bounce_point, is_lower_bounce_point, mixed_wfb
use kt_grids, only: aky
use gs2_layouts, only: g_lo, ik_idx, it_idx, il_idx, ie_idx, is_idx
use run_parameters, only: ieqzip
use kt_grids, only: kwork_filter
implicit none
complex, dimension (-ntgrid:,:,:), intent (in) :: phinew, bparnew
complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in out) :: source
integer :: iglo, ig, ik, it, il, ie, is
complex :: beta1
logical :: use_pass_homog, speriod_flag, is_self_periodic
logical :: is_wfb, is_passing, is_trapped, is_passing_like
!$OMP PARALLEL DO DEFAULT(none) &
!$OMP PRIVATE(iglo, ik, it, il, ie, is, ig, use_pass_homog, speriod_flag, &
!$OMP is_wfb, is_passing, is_trapped, is_self_periodic, is_passing_like, beta1) &
!$OMP SHARED(g_lo, kwork_filter, ieqzip, is_ttp, start_from_previous_solution, &
!$OMP gnew, phinew, bparnew, boundary_option_switch, g, &
!$OMP aky, forbid, g_h, zero_forbid, mixed_wfb, l_links, r_links, &
!$OMP connections, trapped_wfb, passing_wfb, wfb, ntgrid, r, ainv, source) &
!$OMP SCHEDULE(static)
do iglo = g_lo%llim_proc, g_lo%ulim_proc
!/////////////////////////
!// Constants, flags etc.
!/////////////////////////
ik = ik_idx(g_lo,iglo)
it = it_idx(g_lo,iglo)
!Skip work if we're not interested in this ik and it
if(kwork_filter(it,ik)) cycle
if(ieqzip(it,ik)) cycle
il = il_idx(g_lo,iglo)
ie = ie_idx(g_lo,iglo)
is = is_idx(g_lo,iglo)
!/////////////////////////
!// Handle hybrid electrons when passing
!/////////////////////////
if (is_passing_hybrid_electron(is, ik, il)) then
call set_h_zero (gnew, phinew, bparnew, iglo)
cycle
end if
!CMR, 17/4/2012:
! use_pass_homog = T iff one of following applies:
! boundary_option_self_periodic
! OR boundary_option_linked
! OR aky=0
! if use_pass_homog = T, compute homogeneous solution (g1) for passing.
! otherwise ONLY compute inhomogenous solution for passing particles.
!
! speriod_flag = T iff boundary_option_linked AND aky=0
! if speriod_flag = T, apply self-periodic bcs for passing and wfb.
select case (boundary_option_switch)
case (boundary_option_self_periodic)
use_pass_homog = .true.
case (boundary_option_linked)
use_pass_homog = .true.
speriod_flag = aky(ik) == 0.0
case default
use_pass_homog = .false.
end select
use_pass_homog = use_pass_homog .or. aky(ik) == 0.0
! ng2+1 is WFB.
is_wfb = il_is_wfb(il)
! Is this particle passing. Note we exclude the wfb here as we can choose
! how we treat this as given by passing_wfb, trapped_wfb etc.
is_passing = il_is_passing(il)
! Should this particle be treated as a passing particle
is_passing_like = is_passing .or. (is_wfb .and. passing_wfb)
! Is this particle trapped? Note we exclude the wfb here as we can choose
! how we treat this as given by passing_wfb, trapped_wfb etc.
is_trapped = il_is_trapped(il)
! Decide if this point is self periodic
!CMR, 17/4/2012:
! self_periodic bc is applied as follows:
! if boundary_option_linked = .T.
! isolated wfb (ie no linked domains)
! passing + wfb if aky=0
! else (ie boundary_option_linked = .F.)
! wfb
! passing and wfb particles if boundary_option_self_periodic = T
!
!CMR, 6/8/2014:
! Parallel BC for wfb is as follows:
! ballooning space
! wfb ALWAYS self-periodic (-ntgrid:ntgrid)
! linked fluxtube
! wfb self-periodic (-ntgrid:ntgrid) ONLY if ISOLATED,
! OR if iky=0
! OTHERWISE wfb treated as passing for now in (-ntgrid:ntgrid)
! store homogeneous and inhomog solutions
! AND set amplitude of homog later in invert_rhs_linked
! to satisfy self-periodic bc in fully linked SUPER-CELL
if (boundary_option_switch == boundary_option_linked) then
is_self_periodic = (speriod_flag .and. is_passing) & ! Passing particle with ky=0
.or. (speriod_flag .and. is_wfb .and. .not. trapped_wfb) & ! Non-trapped wfb treatment of wfb with ky = 0
.or. (is_wfb .and. .not. connections(iglo)%neighbor & ! Isolated (i.e. ky/=0 cell with no connections) wfb
.and. mixed_wfb ) ! treated as neither passing or trapped particle (i.e. unique)
else
is_self_periodic = (use_pass_homog .and. is_passing_like) & ! ky=0 or with self-periodic boundaries and passing like particles
.or. (is_wfb .and. mixed_wfb) ! Mixed wfb is always self periodic when not linked
end if
!/////////////////////////
!// Initialisation
!/////////////////////////
if (start_from_previous_solution) then
! Initialise gnew to g. Note we probably only need to do this when
! istep == istep_last (i.e. on the second half of a timestep). We
! may also be able to get a slight improvement by copying over the
! non-adiabatic part of g and then forming gnew using the new fields.
! For example we could do the equivalent of
! call g_adjust(gnew, phi, bpar, from_g_gs2)
! call g_adjust(gnew, phinew, bparnew, to_g_gs2)
! after copying g into gnew.
gnew(:, :, iglo) = g(:, :, iglo)
else
gnew(:, :, iglo) = 0.0
end if
!/////////////////////////
!// Boundaries
!/////////////////////////
if (is_passing_like .and. .not. is_self_periodic) then
call apply_parallel_boundary_conditions(gnew, source, it, ik, ie, is, iglo, &
phinew, bparnew, l_links(it, ik) == 0, r_links(it, ik) == 0)
end if
!/////////////////////////
!// Time advance
!/////////////////////////
!!!!!!!!!!!!!!!!!!!!!!!!!!!
! time advance vpar < 0 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!
! inhomogeneous part: gnew
! r=ainv=0 if forbid(ig,il) or forbid(ig+1,il), so gnew=0 in forbidden
! region and at upper bounce point
do ig = ntgrid-1, -ntgrid, -1
gnew(ig,2,iglo) = -gnew(ig+1,2,iglo)*r(ig,2,iglo) + ainv(ig,2,iglo)*source(ig,2,iglo)
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!
! time advance vpar > 0 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!
! First set BCs for trapped particles at lower bounce point
if (is_trapped) then
! match boundary conditions at lower bounce point
! Note we can almost change this range to the full theta domain
! now we use is_lower_bounce_point. We would require a slight
! change in how we handle copying the bounce point value into the source
! as currently we set it at one point before the bounce point. We might
! be able to set gnew at the bounce point directly with some small
! changes to ainv etc.
do ig = -ntgrid, ntgrid - 1
! Skip for ttp
if (is_ttp(ig, il)) cycle
if (is_lower_bounce_point(ig+1, il)) then
!CMR: init_invert_rhs sets ainv=1 at lower bounce point for trapped
! source below ensures gnew(lb,1,iglo)=gnew(lb,2,iglo)
! where lb is the lower bounce point.
source(ig,1,iglo) = gnew(ig+1,2,iglo)
end if
end do
end if
! If trapped wfb enforce the bounce condition at the lower bounce point
! Note this treats the wfb as bouncing at the end of the parallel domain.
! For nperiod > 1 this means the wfb does not bounce at the interior points
! where bmag = bmax
if (is_wfb .and. trapped_wfb) then
ig = -ntgrid
gnew(ig,1,iglo) = gnew(ig,2,iglo)
endif
! time advance vpar > 0 inhomogeneous part
do ig = -ntgrid, ntgrid - 1
if (.not. is_ttp(ig, il)) then
! Use two statements here to *prevent* erroneous vectorization of this loop by gfortran 7.2
gnew(ig + 1, 1, iglo) = ainv(ig, 1, iglo) * source(ig, 1, iglo)
gnew(ig + 1, 1, iglo) = gnew(ig + 1, 1, iglo) - gnew(ig, 1, iglo) * r(ig, 1, iglo)
else
! If we're totally trapped then just enforce continuity in sigma.
! Note we don't allow for totally trapped particles at ig = +/- ntgrid
! so it doesn't matter that the outer loop does not cover the full range.
if (.not. forbid(ig, il)) gnew(ig, 1, iglo) = gnew(ig, 2, iglo)
end if
end do
!/////////////////////////
!// Periodicity at bounce point and parallel boundary
!/////////////////////////
if (is_self_periodic) call self_periodic(gnew, g_h, iglo, is, ie, it, ik, is_wfb, phinew, bparnew)
! add correct amount of homogeneous solution for trapped particles to satisfy boundary conditions
!CMR, 24/7/2014:
!Revert to looping from il>= ng2+2, i.e. exclude wfb as:
! (1) wfb bc is already handled above
! (2) forbid never true for wfb, so including ng2+1 in loop is pointless.
if (is_trapped) then
beta1 = 0.0
do ig = ntgrid-1, -ntgrid, -1
! If totally trapped or forbidden then cycle
if (forbid(ig, il) .or. is_ttp(ig, il)) then
beta1 = 0.0
cycle !CMR: to avoid pointless arithmetic later in loop
else if (forbid(ig+1,il)) then !This could be is_upper_bounce_point(ig, il)?
beta1 = (gnew(ig,1,iglo) - gnew(ig,2,iglo))/(1.0 - g_h(ig, 1, iglo))
end if
gnew(ig,1,iglo) = gnew(ig,1,iglo) + beta1*g_h(ig, 1, iglo)
gnew(ig,2,iglo) = gnew(ig,2,iglo) + beta1*g_h(ig, 2, iglo)
end do
end if
! If wfb is trapped, enforce the bounce condition at the upper bounce point
! by combining homogeneous and inhomogeneous solutions
! Note this treats the wfb as bouncing at the end of the parallel grid.
! When nperiod > 1 this means we do not consider the wfb to bounce at the
! interior points with bmag=bmax. This treatment works out slightly similar
! to the self_periodic treatment but mixes the different sigma.
if (is_wfb .and. trapped_wfb) then
ig = ntgrid
! Note here we find the difference between gnew at different sigma
! but the same theta, whilst in self_periodic we find the difference
! between different theta at the same sigma.
! Note if trapped_wfb then g1(ntgrid, 2) == 1.0 so the denominator
! is equivalent to g1(ig,2) - g1(ig,1) here.
beta1 = (gnew(ig,1,iglo) - gnew(ig,2,iglo))/(1.0 - g_h(ig, 1, iglo))
gnew(:, :, iglo) = gnew(:, :, iglo) + beta1*g_h(:, :, iglo)
end if
!CMR,DD, 25/7/2014:
! Not keen on following kludge zeroing forbidden region
! Introduced new flag: zero_forbid defaults to .false.
! Tested and default is fine linearly, expect should work nonlinearly,
! (Can easily restore old behaviour by setting: zero_forbid=.true.)
! zero out spurious gnew outside trapped boundary
if(zero_forbid)then
where (forbid(:,il))
gnew(:,1,iglo) = 0.0
gnew(:,2,iglo) = 0.0
end where
endif
end do
!$OMP END PARALLEL DO
if (def_parity) call enforce_parity(parity_redist)
end subroutine invert_rhs_1