Skip to content

Commit

Permalink
Latest update
Browse files Browse the repository at this point in the history
  • Loading branch information
@GeorgiaActon
GeorgiaActon committed Dec 17, 2024
1 parent 25246c4 commit f4fecfc
Show file tree
Hide file tree
Showing 4 changed files with 54 additions and 58 deletions.
4 changes: 2 additions & 2 deletions extended_zgrid.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -459,11 +459,11 @@ subroutine enforce_reality (field)
deallocate(field_ext)
end do
end do
! if (periodic(iky)) field(iky,:,nzgrid,:) = field(iky,:,-nzgrid,:)
if (periodic(iky)) field(iky,:,nzgrid,:) = field(iky,:,-nzgrid,:)
end do

! better way around!
field(1, :, nzgrid, :) = field(1, :, -nzgrid, :)
! field(1, :, nzgrid, :) = field(1, :, -nzgrid, :)
! field(1, :, -nzgrid, :) = field(1, :, nzgrid, :)
end subroutine enforce_reality

Expand Down
55 changes: 23 additions & 32 deletions fields.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -780,15 +780,20 @@ contains
if (.not. associated(gam0_ffs_corr(iky, ikx, iz)%fourier)) &
allocate (gam0_ffs_corr(iky, ikx, iz)%fourier(gam0_ffs_corr(iky, ikx, iz)%max_idx))

if(ikx == 1 .and. iky == naky) then
gam0_const(iky, ikx, iz) = 0.0
gam0_ffs_corr(iky, ikx, iz)%fourier = 0.0
else
gam0_const(iky, ikx, iz) = gam0_ffs(iky, ikx, iz)%fourier(1)
gam0_ffs_corr(iky, ikx, iz)%fourier(1) = 0.0
gam0_ffs_corr(iky, ikx, iz)%fourier(2:) = gam0_ffs(iky, ikx, iz)%fourier(2:)
end if

! if(ikx == 1 .and. iky == naky) then
! gam0_const(iky, ikx, iz) = 0.0
! gam0_ffs_corr(iky, ikx, iz)%fourier = gam0_ffs(iky, ikx, iz)%fourier
! else
! gam0_const(iky, ikx, iz) = gam0_ffs(iky, ikx, iz)%fourier(1)
! gam0_ffs_corr(iky, ikx, iz)%fourier(1) = 0.0
! gam0_ffs_corr(iky, ikx, iz)%fourier(2:) = gam0_ffs(iky, ikx, iz)%fourier(2:)
! end if

gam0_const(iky, ikx, iz) = gam0_ffs(iky, ikx, iz)%fourier(1)
! gam0_ffs_corr(iky, ikx, iz)%fourier(1) = gam0_const(iky, ikx, iz)
gam0_ffs_corr(iky, ikx, iz)%fourier(1) = 0.0
gam0_ffs_corr(iky, ikx, iz)%fourier(2:) = gam0_ffs(iky, ikx, iz)%fourier(2:)

end do
end do
end do
Expand All @@ -805,7 +810,7 @@ contains
gamtot = real(gamtot_con)
!> TODO-GA: move this to adiabatic response factor
! if (zonal_mode(1) .and. akx(1) < epsilon(0.) .and. has_electron_species(spec)) then
!gamtot(1, 1, :) = 0.0
gamtot(1, 1, :) = 0.0
!end if

if (.not. has_electron_species(spec)) then
Expand Down Expand Up @@ -1012,7 +1017,7 @@ contains
call get_fields_ffs(g, phi, apar, implicit_solve=.true.)
else
if (debug) write (*, *) 'fields::advance_fields::get_fields_ffs'
call get_fields_ffs(g, phi, apar, implicit_solve=.true.)
call get_fields_ffs(g, phi, apar)
end if
else
call get_fields_vmulo(g, phi, apar, bpar, dist)
Expand Down Expand Up @@ -1336,41 +1341,27 @@ contains
source = 0.0
call get_g_integral_contribution_source(gold, source(:,:,:,1) )
do iz = -nzgrid, nzgrid
call gyro_average(phiold(:,:,iz,1), source2(:,:,iz,1), gam0_ffs_corr(:,:,iz))
end do
call gyro_average(phiold, source2, gam0_ffs_corr)
source = spread(source(:,:,:,1), 4, ntubes)
source2 = spread(source2(:,:,:,1), 4, ntubes)
write(*,*) 'source', source(1,2,-nzgrid,1)
write(*,*) 'phiold', phiold(1,2,-nzgrid,1)
write(*,*) 'gam0_ffs_corr', gam0_ffs_corr(naky,2,-nzgrid)%fourier
write(*,*) 'source2' , source2(1,2,-nzgrid,1)
source = source - source2
write(*,*) 'source before gamtot', source(1,2,-nzgrid,1)
write(*,*) 'gamtot', gamtot(1,2,-nzgrid)
! source = source - source2
where (gamtot_t < epsilon(0.0))
source= 0.0
elsewhere
source = source / gamtot_t
end where
write(*,*) 'source after gamtot', source(1,2,-nzgrid,1)
! do iz = -nzgrid, nzgrid
! do iky = 1, naky
! do ikx = 1, nakx
! if(abs(gamtot(iky, ikx, iz)) .LT. 1E-005) source = 0.0
! end do
! end do
! end do
source(1,:,:,:) = 0.0
if (any(gamtot(1, 1, :) < epsilon(0.))) source(1, 1, :, :) = 0.0
if (akx(1) < epsilon(0.)) then
source(1, 1, :, :) = 0.0
end if
source(1, 1, :, :) = 0.0
!!!! call enforce_reality(source)
! call enforce_reality(source)
deallocate(source2, gamtot_t)
Expand Down
8 changes: 4 additions & 4 deletions implicit_solve.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -147,8 +147,8 @@ subroutine advance_implicit_terms(g, phi, apar, bpar)
!> NB the 'g' here is g_inh^{n+1, i+1}
call advance_fields(g, phi, apar, bpar, dist=trim(dist_choice), implicit_solve=.true.)
!> g2 = g^{n+1, i}
!> phi_old = phi^{n+1, i}
call get_fields_source(g2, phi_old, fields_source_ffs)
!> phi_old = phi^{n+1, i}
call get_fields_source(g2, phi_old, fields_source_ffs)
phi = phi + fields_source_ffs
else
call advance_fields(g, phi, apar, bpar, dist=trim(dist_choice))
Expand Down Expand Up @@ -186,11 +186,11 @@ subroutine advance_implicit_terms(g, phi, apar, bpar)
!> solving for full g = g_{inh} + g_{hom}
modify = .true.
call update_pdf(modify)
g2 = g
!> save the incoming pdf and phi, as they will be needed later
!> this will become g^{n+1, i} -- the g from the previous iteration
fields_updated = .false.
call advance_fields(g2, phi, apar, bpar, dist=trim(dist_choice))
call advance_fields(g, phi, apar, bpar, dist=trim(dist_choice))
g2 = g
phi_old = phi
else
call update_pdf
Expand Down
45 changes: 25 additions & 20 deletions init_g.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -336,33 +336,33 @@ subroutine ginit_default_ffs
right = .not. left

it = 1
do iky = 1, naky
do it = 1, ntubes
do ie = 1, neigen(iky)
nz_ext = nsegments(ie, iky) * nzed_segment + 1
allocate (phiext(nz_ext))
allocate (zed_ext(nz_ext))
zed_ext = [((2 * pi * j) / (nz_ext+1), j=-nz_ext/2 , nz_ext/2)]
! do iky = 1, naky
! do it = 1, ntubes
! do ie = 1, neigen(iky)
! nz_ext = nsegments(ie, iky) * nzed_segment + 1
! allocate (phiext(nz_ext))
! allocate (zed_ext(nz_ext))
! zed_ext = [((2 * pi * j) / (nz_ext+1), j=-nz_ext/2 , nz_ext/2)]

do iz = 1, nz_ext
phiext(iz) = exp(-(zed_ext(iz) /width0)**2) * cmplx(1.0, 1.0)
end do
! do iz = 1, nz_ext
! phiext(iz) = exp(-(zed_ext(iz) /width0)**2) * cmplx(1.0, 1.0)
! end do

call map_from_extended_zgrid(it, ie, iky, phiext, phi(iky, :, :, :))
deallocate(phiext, zed_ext)
end do
end do
! call map_from_extended_zgrid(it, ie, iky, phiext, phi(iky, :, :, :))
! deallocate(phiext, zed_ext)
! end do
! end do
! end do

do iz = -nzgrid, nzgrid
phi(:, :, iz,1) = exp(-((zed(iz) - zed0) / width0)**2) * cmplx(1.0, 1.0)
end do

if (chop_side) then
if (left) phi(:, :, :-1, :) = 0.0
if (right) phi(:, :, 1:, :) = 0.0
end if

! do iz = -nzgrid, nzgrid
! phi(:, :, iz,1) = exp(-((zed(iz) - zed0) / width0)**2) * cmplx(1.0, 1.0)
! end do

if (zonal_mode(1)) then
if (abs(akx(1)) < epsilon(0.0)) then
phi(1, 1, :, :) = 0.0
Expand Down Expand Up @@ -390,9 +390,9 @@ subroutine ginit_default_ffs
do iy = 1, ny
do ikx = 1, ikx_max
do iz = -nzgrid, nzgrid
g0x(iy, ikx, iz, ivmu) = phiinit * phiy(iy, ikx, iz) / abs(spec(is)%z) &
g0x(iy, ikx, iz, ivmu) = phiy(iy, ikx, iz) / abs(spec(is)%z) &
* (den0 + 2.0 * zi * vpa(iv) * upar0) &
* maxwell_mu(iy, iz, imu, is) &
! * maxwell_mu(iy, iz, imu, is) &
* maxwell_vpa(iv, is) * maxwell_fac(is)
end do
end do
Expand All @@ -413,6 +413,11 @@ subroutine ginit_default_ffs
do ivmu = vmu_lo%llim_proc, vmu_lo%ulim_proc
call enforce_reality (gnew(:, :, :, :, ivmu))
end do

! Normalise and scale
gnew = gnew / sum(real(gnew)**2 + aimag(gnew)**2)
gnew = gnew * phiinit

deallocate (g_swap)

gvmu = 0.
Expand Down

0 comments on commit f4fecfc

Please sign in to comment.