integrate_moment_lec Subroutine

private subroutine integrate_moment_lec(lo, g, total)

Perform an integral over velocity space whilst in the LE_LAYOUT in which we have ensured that all of velocity space is local. As such we don't need any calls to MPI reduction routines. Note that this means the processors for different distributed spatial points (x,y) don't know the results at other points.

Arguments

Type IntentOptional Attributes Name
type(le_layout_type), intent(in) :: lo
complex, intent(in), dimension (:,:,lo%llim_proc:) :: g
complex, intent(out), dimension (lo%llim_proc:) :: total

Contents

Source Code


Source Code

  subroutine integrate_moment_lec (lo, g, total)
    use layouts_type, only: le_layout_type
    use gs2_layouts, only: ig_idx, it_idx, ik_idx, is_idx
    use kt_grids, only: kwork_filter
    use array_utils, only: zero_array
    implicit none
    type (le_layout_type), intent (in) :: lo
    complex, dimension (:,:,lo%llim_proc:), intent (in) :: g
    complex, dimension (lo%llim_proc:), intent (out) :: total
    integer :: ixi, ie, ile, ig, it, ik, is, nxup
    call zero_array(total)
    if (nxi > 2 * ng2) then !Could be grid_has_trapped_particles?
       nxup = nxi + 1
    else
       nxup = nxi
    end if

    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP PRIVATE(ile, is, it, ik, ig, ie, ixi) &
    !$OMP SHARED(lo, kwork_filter, negrid, nxup, w, wxi, g) &
    !$OMP REDUCTION(+ : total) &
    !$OMP SCHEDULE(guided)
    do ile = lo%llim_proc, lo%ulim_proc
       it = it_idx (lo,ile)
       ik = ik_idx (lo,ile)
       if (kwork_filter(it, ik)) cycle
       is = is_idx (lo,ile)
       ig = ig_idx (lo,ile)
       do ie = 1, negrid
!CMR, 2/10/2013:
!   nxi+1 limit on do loop below is CRUCIAL, as its stores phase space point
!   corresponding to g_lo (il=nlambda, isign=2).
!   This MUST contribute to the v-space integral, but is NOT
!   needed in collision operator as EQUIVALENT to g_lo(il=nlambda, isign=2).
!   (In collisions at ig=0, both of these points are EXACTLY equivalent, xi=0.)
!MRH actually nxi+1 is needed in the collision operator for consistency 16/08/2018  
          do ixi = 1, nxup
             total(ile) = total(ile) + w(ie, is) * wxi(ixi, ig) * g(ixi, ie, ile)
          end do
       end do
    end do
    !$OMP END PARALLEL DO
    ! No need for communication since all velocity grid points are together
    ! and each prcessor does not touch the unset place
    ! They actually don't need to keep all 4D array
    ! Do we stay in le_layout for total?
    ! --- ile contains necessary and sufficient information for (ig,it,ik,is)

  end subroutine integrate_moment_lec