invert_rhs_1 Subroutine

private subroutine invert_rhs_1(phinew, bparnew, source)

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

Arguments

Type IntentOptional Attributes Name
complex, intent(in), dimension (-ntgrid:,:,:) :: phinew
complex, intent(in), dimension (-ntgrid:,:,:) :: bparnew
complex, intent(inout), dimension (-ntgrid:,:,g_lo%llim_proc:) :: source

Contents

Source Code


Source Code

  subroutine invert_rhs_1 (phinew, bparnew, source)
    use dist_fn_arrays, only: gnew, vperp2, aj1, aj0, set_h_zero, g
    use theta_grid, only: ntgrid
    use le_grids, only: ng2, 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: fbpar, fphi, ieqzip
    use kt_grids, only: kwork_filter
    use species, only: spec, nonmaxw_corr
    use mp, only: mp_abort
    implicit none
    complex, dimension (-ntgrid:,:,:), intent (in) :: phinew, bparnew
    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (in out) :: source
    integer :: iglo
    integer :: ig, ik, it, il, ie, is
    integer :: ilmin
    complex :: beta1
    complex :: adjleft, adjright
    logical :: use_pass_homog, speriod_flag
    logical :: is_wfb, is_passing, is_trapped
    logical :: is_self_periodic, is_passing_like

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(iglo, ik, it, il, ie, is, use_pass_homog, speriod_flag, &
    !$OMP ilmin, is_wfb, is_passing, is_trapped, is_self_periodic, &
    !$OMP adjright, adjleft, ig, beta1, is_passing_like) &
    !$OMP SHARED(g_lo, kwork_filter, ieqzip, is_ttp, start_from_previous_solution, &
    !$OMP gnew, phinew, bparnew, fphi, fbpar, boundary_option_switch, g, &
    !$OMP aky, ng2, forbid, g_h, zero_forbid, mixed_wfb, &
    !$OMP connections, trapped_wfb, passing_wfb, wfb, ntgrid, nonad_zero, vperp2, &
    !$OMP aj1, aj0, spec, nonmaxw_corr, l_links, r_links, 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

       if (use_pass_homog) then
          ilmin = 1
       else
          ilmin = ng2 + 1
       end if

       ! 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
       !/////////////////////////

       !CMR, 18/4/2012:
       ! What follows is a selectable improved parallel bc for passing particles.
       !                                            (prompted by Greg Hammett)
       ! Original bc is: g_gs2 = gnew = 0 at ends of domain:
       !   ONLY implies zero incoming particles in nonadiabatic delta f if phi=bpar=0
       ! Here ensure ZERO incoming particles in nonadiabatic delta f at domain ends
       !  (exploits code used in subroutine g_adjust to transform g_wesson to g_gs2)
       if ( nonad_zero .and. is_passing_like ) then
          !Only apply the new boundary condition to the leftmost
          !cell for sign going from left to right
          if (l_links(it, ik) .eq. 0) then
             adjleft = 2.0*vperp2(-ntgrid,iglo)*aj1(-ntgrid,iglo) &
                  *bparnew(-ntgrid,it,ik)*fbpar &
                  + spec(is)%zt*phinew(-ntgrid,it,ik)*aj0(-ntgrid,iglo) &
                  *fphi*nonmaxw_corr(ie,is)
             gnew(-ntgrid,1,iglo) = -adjleft
          end if
          !Only apply the new boundary condition to the rightmost
          !cell for sign going from right to left
          if (r_links(it, ik) .eq. 0) then
             adjright = 2.0*vperp2(ntgrid,iglo)*aj1(ntgrid,iglo) &
                  *bparnew(ntgrid,it,ik)*fbpar &
                  + spec(is)%zt*phinew(ntgrid,it,ik)*aj0(ntgrid,iglo) &
                  *fphi*nonmaxw_corr(ie,is)
             gnew(ntgrid,2,iglo) = -adjright
          end if
       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