2042563664

Committed 16 Sep 2025 01:04PM UTC coverage: 10.677% (+0.07%) from 10.607%

Build # 2042563664

Build Type

push

gitlab-ci

Committed by David Dickinson

Commit Message

Merged in minor/tidy_up_some_gs2_io_routines (pull request #1173)

Run Details

4712 of 44131 relevant lines covered (10.68%)

126484.23 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/src/diagnostics/diagnostics_zonal_transfer.f90

!> calculate tau(kx,ky), the transfer of free energy by the non-linearity in fourier space
module diagnostics_zonal_transfer
  implicit none

  private

  public :: init_diagnostics_transfer, finish_diagnostics_transfer
  public :: calculate_zonal_transfer, write_zonal_transfer

  ! free energy transfer due to nonlinearity
  complex, dimension(:,:), allocatable :: zonal_transfer
  logical :: initialised = .false.

contains
  !> FIXME : Add documentation
  subroutine init_diagnostics_transfer
    use kt_grids, only: ntheta0, naky
    implicit none
    if(initialised) return
    initialised = .true.
    allocate(zonal_transfer(ntheta0, naky))
  end subroutine init_diagnostics_transfer

  subroutine finish_diagnostics_transfer
    implicit none
    if (.not. initialised) return
    initialised = .false.
    deallocate(zonal_transfer)
  end subroutine finish_diagnostics_transfer
  
  !> calculate nonlinear term appearing in GKE for now assume
  !> electrostatic. Why do we need this rather than just using
  !> nonlinear_terms:calculate_current_nl_source_and_cfl_limit?
  subroutine non_linear_term (g1)
    use theta_grid, only: ntgrid, kxfac
    use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx
    use gs2_layouts, only: yxf_lo
    use dist_fn_arrays, only: g, g_adjust, to_g_gs2, from_g_gs2
    use species, only: spec
    use gs2_transforms, only: transform2, inverse2
    use kt_grids, only: aky, akx
    use fields_arrays, only: phi, bpar
    use constants, only: zi
    implicit none
    complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent(out) :: g1
    integer :: i, j, iglo, ik, it, ig, is
    real, dimension(:,:), allocatable :: banew, gbnew, bracketnew

    allocate (banew(yxf_lo%ny,yxf_lo%llim_proc:yxf_lo%ulim_alloc)) ; banew = 0.
    allocate (gbnew(yxf_lo%ny,yxf_lo%llim_proc:yxf_lo%ulim_alloc)) ; gbnew = 0.
    allocate (bracketnew(yxf_lo%ny,yxf_lo%llim_proc:yxf_lo%ulim_alloc)) ; bracketnew = 0.

    ! Form g1=i*kx*phi
    call load_kx_phi

    ! Transform to real space
    ! 2D fft of g1 is returned as banew
    call transform2 (g1, banew)

    ! Form g1=i*ky*g
    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       ik = ik_idx(g_lo,iglo)
       g1(:,:,iglo)=g(:,:,iglo)*zi*aky(ik)
    end do

    ! Transform to real space
    call transform2 (g1, gbnew)

    ! Calculate (d phi /dx).(dg/dy)
    do j = yxf_lo%llim_proc, yxf_lo%ulim_proc
       do i = 1, yxf_lo%ny
          bracketnew(i,j) = banew(i,j)*gbnew(i,j)*kxfac
       end do
    end do

    ! Form g1=i*ky*chi
    call load_ky_phi

    !Transform to real space
    call transform2 (g1, banew)

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       it = it_idx(g_lo,iglo)
       g1(:,:,iglo)=g(:,:,iglo)*zi*akx(it)
    enddo

    ! Transform to real space
    call transform2 (g1, gbnew)

    ! Calculate (d phi /dy).(dg/dx) and subtract from (d phi /dx).(dg/dy)
    do j = yxf_lo%llim_proc, yxf_lo%ulim_proc
       do i = 1, yxf_lo%ny
          bracketnew(i,j) = bracketnew(i,j) - banew(i,j)*gbnew(i,j)*kxfac
       end do
    end do

    ! Transform nonlinearity back to spectral space
    ! g1 contains Fourier coefficients associated with nonlinearity
    call inverse2 (bracketnew, g1)

    call g_adjust(g,phi,bpar, direction = from_g_gs2)

    do iglo = g_lo%llim_proc, g_lo%ulim_proc
       is = is_idx(g_lo,iglo)
       g1(:,:,iglo)=g1(:,:,iglo)*conjg(g(:,:,iglo))*spec(is)%temp
    end do

    call g_adjust(g,phi,bpar, direction = to_g_gs2)

    deallocate (banew, gbnew, bracketnew)

  contains

    !> FIXME : Add documentation
    subroutine load_kx_phi
      use dist_fn_arrays, only: aj0  !< bessel function from the fourier coeff of average
      implicit none
      complex :: fac

      do iglo = g_lo%llim_proc, g_lo%ulim_proc
         it = it_idx(g_lo,iglo)
         ik = ik_idx(g_lo,iglo)
         do ig = -ntgrid, ntgrid
            fac = zi*akx(it)*aj0(ig,iglo)*phi(ig,it,ik)
            g1(ig,1,iglo) = fac
            g1(ig,2,iglo) = fac
         end do
      end do

    end subroutine load_kx_phi

    !> FIXME : Add documentation    
    subroutine load_ky_phi
      use dist_fn_arrays, only: aj0
      implicit none
      complex :: fac

      do iglo = g_lo%llim_proc, g_lo%ulim_proc
         it = it_idx(g_lo,iglo)
         ik = ik_idx(g_lo,iglo)
         do ig = -ntgrid, ntgrid
            fac = zi*aky(ik)*aj0(ig,iglo)*phi(ig,it,ik)
            g1(ig,1,iglo) = fac
            g1(ig,2,iglo) = fac
         end do
      end do

    end subroutine load_ky_phi

  end subroutine non_linear_term

  !> Calculates tau, the transfer of free energy as a function of (kx,ky) and
  !> stores in module level zonal_transfer
  subroutine calculate_zonal_transfer
    use mp, only: proc0, broadcast
    use theta_grid, only: ntgrid, field_line_average
    use gs2_layouts, only: g_lo, ik_idx, it_idx
    use species, only: nspec
    use kt_grids, only: ntheta0, naky
    use le_grids, only: integrate_moment
    implicit none
    complex, dimension (:,:,:), allocatable :: g0
    complex, dimension (:,:,:), allocatable :: inter
    complex, dimension (:,:,:,:), allocatable :: inter_s
    integer :: is
    
    allocate (g0(-ntgrid:ntgrid,2,g_lo%llim_proc:g_lo%ulim_alloc))
    allocate (inter(-ntgrid:ntgrid,ntheta0,naky)) ; inter = 0.
    allocate (inter_s(-ntgrid:ntgrid,ntheta0,naky,nspec)) ; inter_s = 0.

    ! obtain the nonlinear term
    call non_linear_term(g0)

    ! integrate over velocities
    ! can only be called after w and wl are assigned in init_le_grids
    call integrate_moment(g0, inter_s)

    if (proc0) then
       ! sum over species
       do is=1,nspec
          inter = inter + inter_s(:,:,:,is)
       end do

       ! integrate over theta, must be done after velocity space integration
       zonal_transfer = field_line_average(inter)
    end if
    
    call broadcast (zonal_transfer)
    deallocate (g0, inter, inter_s)
  end subroutine calculate_zonal_transfer

  !> FIXME : Add documentation  
  subroutine write_zonal_transfer(file_id, nout)
    use gs2_io, only: nc_write_zonal_transfer
    use mp, only: proc0
    implicit none
    integer, intent(in) :: file_id, nout
    if (.not. proc0) return
    call nc_write_zonal_transfer(file_id, nout, zonal_transfer)
  end subroutine write_zonal_transfer
end module diagnostics_zonal_transfer

1	!> calculate tau(kx,ky), the transfer of free energy by the non-linearity in fourier space
2	module diagnostics_zonal_transfer
3	implicit none
4
5	private
6
7	public :: init_diagnostics_transfer, finish_diagnostics_transfer
8	public :: calculate_zonal_transfer, write_zonal_transfer
9
10	! free energy transfer due to nonlinearity
11	complex, dimension(:,:), allocatable :: zonal_transfer
12	logical :: initialised = .false.
13
14	contains
15	!> FIXME : Add documentation
16	subroutine init_diagnostics_transfer	×
17	use kt_grids, only: ntheta0, naky
18	implicit none
19	if(initialised) return	×
20	initialised = .true.	×
21	allocate(zonal_transfer(ntheta0, naky))	×
22	end subroutine init_diagnostics_transfer
23
24	subroutine finish_diagnostics_transfer	×
25	implicit none
26	if (.not. initialised) return	×
27	initialised = .false.	×
28	deallocate(zonal_transfer)	×
29	end subroutine finish_diagnostics_transfer
30
31	!> calculate nonlinear term appearing in GKE for now assume
32	!> electrostatic. Why do we need this rather than just using
33	!> nonlinear_terms:calculate_current_nl_source_and_cfl_limit?
34	subroutine non_linear_term (g1)	×
35	use theta_grid, only: ntgrid, kxfac
36	use gs2_layouts, only: g_lo, ik_idx, it_idx, is_idx
37	use gs2_layouts, only: yxf_lo
38	use dist_fn_arrays, only: g, g_adjust, to_g_gs2, from_g_gs2
39	use species, only: spec
40	use gs2_transforms, only: transform2, inverse2
41	use kt_grids, only: aky, akx
42	use fields_arrays, only: phi, bpar
43	use constants, only: zi
44	implicit none
45	complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent(out) :: g1
46	integer :: i, j, iglo, ik, it, ig, is
47	real, dimension(:,:), allocatable :: banew, gbnew, bracketnew	×
48
49	allocate (banew(yxf_lo%ny,yxf_lo%llim_proc:yxf_lo%ulim_alloc)) ; banew = 0.	×
50	allocate (gbnew(yxf_lo%ny,yxf_lo%llim_proc:yxf_lo%ulim_alloc)) ; gbnew = 0.	×
51	allocate (bracketnew(yxf_lo%ny,yxf_lo%llim_proc:yxf_lo%ulim_alloc)) ; bracketnew = 0.	×
52
53	! Form g1=ikxphi
54	call load_kx_phi	×
55
56	! Transform to real space
57	! 2D fft of g1 is returned as banew
58	call transform2 (g1, banew)	×
59
60	! Form g1=ikyg
61	do iglo = g_lo%llim_proc, g_lo%ulim_proc	×
62	ik = ik_idx(g_lo,iglo)	×
63	g1(:,:,iglo)=g(:,:,iglo)ziaky(ik)	×
64	end do
65
66	! Transform to real space
67	call transform2 (g1, gbnew)	×
68
69	! Calculate (d phi /dx).(dg/dy)
70	do j = yxf_lo%llim_proc, yxf_lo%ulim_proc	×
71	do i = 1, yxf_lo%ny	×
72	bracketnew(i,j) = banew(i,j)gbnew(i,j)kxfac	×
73	end do
74	end do
75
76	! Form g1=ikychi
77	call load_ky_phi	×
78
79	!Transform to real space
80	call transform2 (g1, banew)	×
81
82	do iglo = g_lo%llim_proc, g_lo%ulim_proc	×
83	it = it_idx(g_lo,iglo)	×
84	g1(:,:,iglo)=g(:,:,iglo)ziakx(it)	×
85	enddo
86
87	! Transform to real space
88	call transform2 (g1, gbnew)	×
89
90	! Calculate (d phi /dy).(dg/dx) and subtract from (d phi /dx).(dg/dy)
91	do j = yxf_lo%llim_proc, yxf_lo%ulim_proc	×
92	do i = 1, yxf_lo%ny	×
93	bracketnew(i,j) = bracketnew(i,j) - banew(i,j)gbnew(i,j)kxfac	×
94	end do
95	end do
96
97	! Transform nonlinearity back to spectral space
98	! g1 contains Fourier coefficients associated with nonlinearity
99	call inverse2 (bracketnew, g1)	×
100
101	call g_adjust(g,phi,bpar, direction = from_g_gs2)	×
102
103	do iglo = g_lo%llim_proc, g_lo%ulim_proc	×
104	is = is_idx(g_lo,iglo)	×
105	g1(:,:,iglo)=g1(:,:,iglo)conjg(g(:,:,iglo))spec(is)%temp	×
106	end do
107
108	call g_adjust(g,phi,bpar, direction = to_g_gs2)	×
109
110	deallocate (banew, gbnew, bracketnew)	×
111
112	contains
113
114	!> FIXME : Add documentation
115	subroutine load_kx_phi	×
116	use dist_fn_arrays, only: aj0 !< bessel function from the fourier coeff of average
117	implicit none
118	complex :: fac
119
120	do iglo = g_lo%llim_proc, g_lo%ulim_proc	×
121	it = it_idx(g_lo,iglo)	×
122	ik = ik_idx(g_lo,iglo)	×
123	do ig = -ntgrid, ntgrid	×
124	fac = ziakx(it)aj0(ig,iglo)*phi(ig,it,ik)	×
125	g1(ig,1,iglo) = fac	×
126	g1(ig,2,iglo) = fac	×
127	end do
128	end do
129
130	end subroutine load_kx_phi	×
131
132	!> FIXME : Add documentation
133	subroutine load_ky_phi	×
134	use dist_fn_arrays, only: aj0
135	implicit none
136	complex :: fac
137
138	do iglo = g_lo%llim_proc, g_lo%ulim_proc	×
139	it = it_idx(g_lo,iglo)	×
140	ik = ik_idx(g_lo,iglo)	×
141	do ig = -ntgrid, ntgrid	×
142	fac = ziaky(ik)aj0(ig,iglo)*phi(ig,it,ik)	×
143	g1(ig,1,iglo) = fac	×
144	g1(ig,2,iglo) = fac	×
145	end do
146	end do
147
148	end subroutine load_ky_phi	×
149
150	end subroutine non_linear_term
151
152	!> Calculates tau, the transfer of free energy as a function of (kx,ky) and
153	!> stores in module level zonal_transfer
154	subroutine calculate_zonal_transfer	×
155	use mp, only: proc0, broadcast
156	use theta_grid, only: ntgrid, field_line_average	×
157	use gs2_layouts, only: g_lo, ik_idx, it_idx
158	use species, only: nspec
159	use kt_grids, only: ntheta0, naky
160	use le_grids, only: integrate_moment
161	implicit none
162	complex, dimension (:,:,:), allocatable :: g0	×
163	complex, dimension (:,:,:), allocatable :: inter	×
164	complex, dimension (:,:,:,:), allocatable :: inter_s	×
165	integer :: is
166
167	allocate (g0(-ntgrid:ntgrid,2,g_lo%llim_proc:g_lo%ulim_alloc))	×
168	allocate (inter(-ntgrid:ntgrid,ntheta0,naky)) ; inter = 0.	×
169	allocate (inter_s(-ntgrid:ntgrid,ntheta0,naky,nspec)) ; inter_s = 0.	×
170
171	! obtain the nonlinear term
172	call non_linear_term(g0)	×
173
174	! integrate over velocities
175	! can only be called after w and wl are assigned in init_le_grids
176	call integrate_moment(g0, inter_s)	×
177
178	if (proc0) then	×
179	! sum over species
180	do is=1,nspec	×
181	inter = inter + inter_s(:,:,:,is)	×
182	end do
183
184	! integrate over theta, must be done after velocity space integration
185	zonal_transfer = field_line_average(inter)	×
186	end if
187
188	call broadcast (zonal_transfer)	×
189	deallocate (g0, inter, inter_s)	×
190	end subroutine calculate_zonal_transfer	×
191
192	!> FIXME : Add documentation
193	subroutine write_zonal_transfer(file_id, nout)	×
194	use gs2_io, only: nc_write_zonal_transfer
195	use mp, only: proc0
196	implicit none
197	integer, intent(in) :: file_id, nout
198	if (.not. proc0) return	×
199	call nc_write_zonal_transfer(file_id, nout, zonal_transfer)	×
200	end subroutine write_zonal_transfer
201	end module diagnostics_zonal_transfer

gyrokinetics / gs2 / 2042563664

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous