Skip to content

Commit

Permalink
latest changes
Browse files Browse the repository at this point in the history
  • Loading branch information
@GeorgiaActon
GeorgiaActon committed Jan 16, 2025
1 parent 1db4acf commit e42c6b0
Show file tree
Hide file tree
Showing 12 changed files with 364 additions and 1,383 deletions.
1,156 changes: 0 additions & 1,156 deletions #gyro_averages.f90#
View file

This file was deleted.

1 change: 0 additions & 1 deletion .#gyro_averages.f90
View file

This file was deleted.

2 changes: 2 additions & 0 deletions Makefile.depend
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,6 +177,7 @@ euterpe_interface.o: \
extended_zgrid.o: \
constants.o \
kt_grids.o \
mp.o \
zgrid.o
ffs_solve.o: \
constants.o \
Expand Down Expand Up @@ -332,6 +333,7 @@ init_g.o: \
stella_transforms.o \
system_fortran.o \
text_options.o \
time_advance.o \
vpamu_grids.o \
zgrid.o
inputprofiles_interface.o: \
Expand Down
44 changes: 29 additions & 15 deletions dist_fn.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,17 +134,21 @@ subroutine init_kperp2
use kt_grids, only: akx, aky
use kt_grids, only: zonal_mode
use kt_grids, only: nalpha


use extended_zgrid, only: periodic

implicit none

integer :: iky, ikx

integer :: iz, ia

if (kp2init) return
kp2init = .true.

!> allocate the kperp2 array to contain |k_perp|^2
allocate (kperp2(naky, nakx, nalpha, -nzgrid:nzgrid))

do iky = 1, naky
if (zonal_mode(iky)) then
do ikx = 1, nakx
Expand All @@ -162,16 +166,20 @@ subroutine init_kperp2
end do
end if
end do

! NB: should really avoid this by using higher resolution when reading in VMEC geometry and then
! NB: course-graining if necessary to map onto lower-resolution stella grid
! ensure kperp2 is positive everywhere (only might go negative if using full-flux-surface due to interpolation)
where (kperp2 < 0.0)
kperp2 = 0.0
end where
end where

call enforce_single_valued_kperp2

! do iky = 1, naky
! if (periodic(iky)) kperp2(iky, :, :, nzgrid) = kperp2(iky, :, :, -nzgrid)
! end do
call enforce_single_valued_kperp2

end subroutine init_kperp2

!> init_dkperp2dr allocates and initialises the dkperp2dr array, needed for radial variation
Expand Down Expand Up @@ -229,7 +237,7 @@ subroutine enforce_single_valued_kperp2
use dist_fn_arrays, only: kperp2
use kt_grids, only: naky, nalpha
use zgrid, only: nzgrid
use extended_zgrid, only: neigen, nsegments, ikxmod
use extended_zgrid, only: neigen, nsegments, ikxmod, periodic

implicit none

Expand All @@ -239,15 +247,21 @@ subroutine enforce_single_valued_kperp2
allocate (tmp(nalpha)); tmp = 0.0

do iky = 1, naky
do ie = 1, neigen(iky)
if (nsegments(ie, iky) > 1) then
do iseg = 2, nsegments(ie, iky)
tmp = 0.5 * (kperp2(iky, ikxmod(iseg - 1, ie, iky), :, nzgrid) + kperp2(iky, ikxmod(iseg, ie, iky), :, -nzgrid))
kperp2(iky, ikxmod(iseg, ie, iky), :, -nzgrid) = tmp
kperp2(iky, ikxmod(iseg - 1, ie, iky), :, nzgrid) = tmp
end do
end if
end do
if(.not. periodic(iky)) then
do ie = 1, neigen(iky)
if (nsegments(ie, iky) > 1) then
do iseg = 2, nsegments(ie, iky)
!tmp = 0.5 * (kperp2(iky, ikxmod(iseg - 1, ie, iky), :, nzgrid) + kperp2(iky, ikxmod(iseg, ie, iky), :, -nzgrid))
!kperp2(iky, ikxmod(iseg, ie, iky), :, -nzgrid) = tmp
!kperp2(iky, ikxmod(iseg - 1, ie, iky), :, nzgrid) = tmp

tmp = kperp2(iky, ikxmod(iseg, ie, iky), :, -nzgrid)
kperp2(iky, ikxmod(iseg - 1, ie, iky), :, nzgrid) = tmp

end do
end if
end do
end if
end do

deallocate (tmp)
Expand Down
138 changes: 119 additions & 19 deletions extended_zgrid.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -400,15 +400,63 @@ subroutine map_to_extended_zgrid(it, ie, iky, g, gext, ulim)
curr_shift = curr_shift / phase_shift(iky)
ikx = ikxmod(iseg, ie, iky)
itmod = it_right(itmod)
llim = ulim + 1
ulim = llim + nzed_segment - 1
gext(llim:ulim) = g(ikx, iz_low(iseg) + 1:iz_up(iseg), itmod) * curr_shift
! llim = ulim + 1
! ulim = llim + nzed_segment - 1
! gext(llim:ulim) = g(ikx, iz_low(iseg) + 1:iz_up(iseg), itmod) * curr_shift
llim = ulim
ulim = llim + nzed_segment
!gext(llim) = 0.5 * (g(ikx, iz_low(iseg), itmod) * curr_shift + gext(llim))
!gext(llim + 1: ulim) = g(ikx, iz_low(iseg) + 1:iz_up(iseg), itmod) * curr_shift

gext(llim:ulim) = g(ikx, iz_low(iseg):iz_up(iseg), itmod) * curr_shift
end do
end if
end if

end subroutine map_to_extended_zgrid

subroutine map_to_extended_zgrid2(it, ie, iky, g, gext, ulim)

use zgrid, only: nzgrid
use mp, only: proc0
implicit none

integer, intent(in) :: it, ie, iky
complex, dimension(:, -nzgrid:, :), intent(in) :: g
complex, dimension(:), intent(out) :: gext
integer, intent(out) :: ulim

integer :: iseg, ikx, itmod
integer :: llim
complex :: curr_shift

! avoid double-counting at boundaries between 2pi segments
iseg = 1
curr_shift = 1.
ikx = ikxmod(iseg, ie, iky)
llim = 1; ulim = nzed_segment + 1

if (periodic(iky)) then
gext(llim:ulim - 1) = g(ikx, iz_low(iseg):iz_up(iseg) - 1, it)
gext(ulim) = g(ikx, iz_low(iseg),it) / phase_shift(iky)
else
gext(llim:ulim) = g(ikx, iz_low(iseg):iz_up(iseg), it) * curr_shift
if (nsegments(ie, iky) > 1) then
itmod = it
do iseg = 2, nsegments(ie, iky)
curr_shift = curr_shift / phase_shift(iky)
ikx = ikxmod(iseg, ie, iky)
itmod = it_right(itmod)
llim = ulim
ulim = llim + nzed_segment
gext(llim) = 0.5 * (g(ikx, iz_low(iseg), itmod) * curr_shift + gext(llim))
gext(llim + 1: ulim) = g(ikx, iz_low(iseg) + 1:iz_up(iseg), itmod) * curr_shift
end do
end if
end if

end subroutine map_to_extended_zgrid2

subroutine map_from_extended_zgrid(it, ie, iky, gext, g)

use zgrid, only: nzgrid
Expand All @@ -430,54 +478,106 @@ subroutine map_from_extended_zgrid(it, ie, iky, gext, g)

if (periodic(iky)) then
g(ikx, iz_low(iseg):iz_up(iseg) - 1, it) = gext(llim:ulim - 1)
g(ikx, iz_up(iseg), it) = g(ikx, iz_low(iseg), it) * phase_shift(iky)
g(ikx, iz_up(iseg), it) = g(ikx, iz_low(iseg), it) / phase_shift(iky)
! g(ikx, iz_up(iseg), it) = g(ikx, iz_low(iseg), it) * phase_shift(iky)
else
g(ikx, iz_low(iseg):iz_up(iseg), it) = gext(llim:ulim)
if (nsegments(ie, iky) > 1) then
itmod = it
do iseg = 2, nsegments(ie, iky)
curr_shift = curr_shift * phase_shift(iky)
llim = ulim + 1
ulim = llim + nzed_segment - 1
! llim = ulim + 1
! ulim = llim + nzed_segment - 1
ikx = ikxmod(iseg, ie, iky)
itmod = it_right(itmod)
g(ikx, iz_low(iseg):iz_up(iseg), itmod) = gext(llim - 1:ulim) * curr_shift
! g(ikx, iz_low(iseg):iz_up(iseg), itmod) = gext(llim - 1:ulim) * curr_shift

llim = ulim
ulim = llim + nzed_segment
g(ikx, iz_low(iseg):iz_up(iseg), itmod) = gext(llim:ulim) * curr_shift
g(ikxmod(iseg - 1, ie, iky), iz_up(iseg - 1), itmod) = g(ikx, iz_low(iseg), itmod)

end do
end if
end if

end subroutine map_from_extended_zgrid

subroutine map_from_extended_zgrid2(it, ie, iky, gext, g)

use zgrid, only: nzgrid
use mp, only: proc0
implicit none

integer, intent(in) :: it, ie, iky
complex, dimension(:), intent(in) :: gext
complex, dimension(:, -nzgrid:, :), intent(in out) :: g

integer :: iseg, ikx, itmod
integer :: llim, ulim
complex :: curr_shift

iseg = 1
curr_shift = 1.
ikx = ikxmod(iseg, ie, iky)
llim = 1; ulim = nzed_segment + 1

if (periodic(iky)) then
g(ikx, iz_low(iseg):iz_up(iseg) - 1, it) = gext(llim:ulim - 1)
g(ikx, iz_up(iseg), it) = g(ikx, iz_low(iseg), it) / phase_shift(iky)
else
g(ikx, iz_low(iseg):iz_up(iseg), it) = gext(llim:ulim)
if (nsegments(ie, iky) > 1) then
itmod = it
do iseg = 2, nsegments(ie, iky)
curr_shift = curr_shift * phase_shift(iky)

ikx = ikxmod(iseg, ie, iky)
itmod = it_right(itmod)

llim = ulim
ulim = llim + nzed_segment

!!!! Can add this plus one and it doesnt do anything???? strange
! g(ikx, iz_low(iseg):iz_up(iseg), itmod) = gext(llim:ulim) * curr_shift
g(ikx, iz_low(iseg):iz_up(iseg), itmod) = gext(llim:ulim + 1) * curr_shift

!g(ikxmod(iseg - 1, ie, iky), iz_up(iseg - 1), itmod) = g(ikx, iz_low(iseg), itmod)
!g(ikx, iz_low(iseg):iz_up(iseg), itmod) = gext(llim:ulim) * curr_shift

end do
end if
end if

end subroutine map_from_extended_zgrid2

subroutine enforce_reality (field)

use zgrid, only: nzgrid, ntubes
use kt_grids, only: naky, nakx

use mp, only: proc0 , broadcast
implicit none

complex, dimension(:, :, -nzgrid:, :), intent(inout) :: field
complex, dimension(:), allocatable :: field_ext
integer :: ulim, nz_ext
integer :: ulim, llim, nz_ext
integer :: it, ie, iky
integer :: iseg, ikx, itmod
complex :: curr_shift
integer :: ikx

field(1, 1, :, :) = real(field(1, 1, :, :))
do ikx = 2, nakx / 2 + 1
field(1, ikx, :, :) = 0.5 * (field(1, ikx, :, :) + conjg(field(1, nakx - ikx + 2, :, :)))
field(1, nakx - ikx + 2, :, :) = conjg(field(1, ikx, :, :))
end do

it = 1
! do ikx = 2, nakx / 2 + 1
! field(1, ikx, :, :) = 0.5 * (field(1, ikx, :, :) + conjg(field(1, nakx - ikx + 2, :, :)))
! field(1, nakx - ikx + 2, :, :) = conjg(field(1, ikx, :, :))
! end do
do iky = 1, naky
do ie = 1, neigen(iky)
nz_ext = nsegments(ie, iky) * nzed_segment + 1
allocate (field_ext(nz_ext)); field_ext = 0.0
call map_to_extended_zgrid (it, ie, iky, field(iky, :, :, :), field_ext, ulim)
call map_to_extended_zgrid2 (it, ie, iky, field(iky, :, :, :), field_ext, ulim)
call map_from_extended_zgrid (it, ie, iky, field_ext, field(iky,:, :, :))
deallocate(field_ext)
end do
if (periodic(iky)) field(iky,:,-nzgrid,:) = field(iky,:,nzgrid,:) * phase_shift(iky)
! if (periodic(iky)) field(iky,:,-nzgrid,:) = field(iky,:,nzgrid,:) * phase_shift(iky)
end do

end subroutine enforce_reality
Expand Down
8 changes: 4 additions & 4 deletions ffs_solve.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ subroutine add_correction_ffs (phiin, gin, source_out)
source1 = 0.0
source_out = 0.0
call get_drifts_ffs_itteration (phiin, gin, source1)
! call get_source_ffs_itteration (phiin, gin, phi_bar, source_out)
! call get_source_ffs_itteration (phiin, gin, source_out)

source_out = source_out + source1

Expand All @@ -37,7 +37,7 @@ end subroutine add_correction_ffs

! contains

subroutine get_source_ffs_itteration (phi, g, phi_bar, source)
subroutine get_source_ffs_itteration (phi, g, source)

use mp, only: proc0
use stella_layouts, only: vmu_lo
Expand All @@ -63,7 +63,7 @@ subroutine get_source_ffs_itteration (phi, g, phi_bar, source)
use extended_zgrid, only: enforce_reality
implicit none

complex, dimension(:, :, -nzgrid:, :), intent(in) :: phi, phi_bar
complex, dimension(:, :, -nzgrid:, :), intent(in) :: phi
complex, dimension(:, :, -nzgrid:, :, vmu_lo%llim_proc:), intent(in) :: g
complex, dimension(:, :, -nzgrid:, :, vmu_lo%llim_proc:), intent (out) :: source

Expand Down Expand Up @@ -105,7 +105,7 @@ subroutine get_source_ffs_itteration (phi, g, phi_bar, source)
call get_dzed(iv, g0, dgphi_dz)
!> d phi/dz
! call get_dzed(iv, phi, dphi_dz)
call get_dzed(iv, phi_bar, dphi_dz)
call get_dzed(iv, phi, dphi_dz)
!> dg/dz
call get_dzed(iv, g(:, :, :, :, ivmu), dgdz)
!> get these quantities in real space
Expand Down
Loading

0 comments on commit e42c6b0

Please sign in to comment.