diff --git a/extended_zgrid.f90 b/extended_zgrid.f90 index 0869cbada5..71ab1b4eb1 100644 --- a/extended_zgrid.f90 +++ b/extended_zgrid.f90 @@ -75,7 +75,7 @@ subroutine init_extended_zgrid !> Usually set to zero for standard local simulation, but can !> have an effect for global simulations and simulations with low !> magnetic shear that use periodic boundary conditions everywhere - phase_shift = exp(zi * aky * phase_shift_angle) + phase_shift = exp( zi * aky * phase_shift_angle) if (boundary_option_switch == boundary_option_linked .or. boundary_option_switch == boundary_option_linked_stellarator) then @@ -418,7 +418,8 @@ subroutine map_from_extended_zgrid(it, ie, iky, gext, g) curr_shift = 1. ikx = ikxmod(iseg, ie, iky) llim = 1; ulim = nzed_segment + 1 - g(ikx, iz_low(iseg):iz_up(iseg), it) = gext(llim:ulim) + 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) @@ -462,7 +463,8 @@ subroutine enforce_reality (field) 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 diff --git a/ffs_solve.f90 b/ffs_solve.f90 index 4bbf448008..8a0c502045 100644 --- a/ffs_solve.f90 +++ b/ffs_solve.f90 @@ -102,13 +102,10 @@ subroutine get_source_ffs_itteration (phi, g, source) call gyro_average(phi, g0, j0_ffs(:, :, :, ivmu)) !> Full d /dz -!! call get_dgdz_centered(g0, ivmu, dgphi_dz) call get_dzed(iv, g0, dgphi_dz) !> d phi/dz -!! call get_dgdz_centered(phi, ivmu, dphi_dz) call get_dzed(iv, phi, dphi_dz) !> dg/dz -!! call get_dgdz(g(:, :, :, :, ivmu), ivmu, dgdz) call get_dzed(iv, g(:, :, :, :, ivmu), dgdz) !> get these quantities in real space do it = 1, ntubes @@ -137,8 +134,6 @@ subroutine get_source_ffs_itteration (phi, g, source) call center_zed(iv, coeff2(ia, :) , -nzgrid) end do g2y = spread(spread(coeff, 2, ikx_max), 4, ntubes) * g2y + spread(spread(coeff2, 2, ikx_max), 4, ntubes) * g0y -! g2y = spread(spread(stream_correction(:,:,iv,is), 2, ikx_max), 4, ntubes) & -! * (g2y + g0y * spread(spread(maxwell_mu_avg(:, :, imu, is), 2, ikx_max),4, ntubes ) * scratch) coeff = stream_store_full (:,:,iv,is) * (maxwell_mu(:, :, imu, is) - maxwell_mu_avg(:, :, imu, is)) * scratch do ia = 1, ny @@ -146,19 +141,12 @@ subroutine get_source_ffs_itteration (phi, g, source) end do g3y = spread(spread(coeff, 2, ikx_max), 4, ntubes) * g1y -! g3y = spread(spread(stream_store_full (:,:,iv,is) * (maxwell_mu(:, :, imu, is) & -!! - maxwell_mu_avg(:, :, imu, is)), 2, ikx_max),4, ntubes) & -! * g1y * scratch - coeff = stream_store_full (:,:,iv,is) * maxwell_mu(:, :, imu, is) do ia = 1, ny call center_zed(iv, coeff(ia, :) , -nzgrid) end do g0y = spread(spread(coeff, 2, ikx_max), 4, ntubes) * (g0y - g1y) * scratch -! g0y = spread(spread(stream_store_full (:,:,iv,is) * maxwell_mu(:, :, imu, is) , 2, ikx_max),4, ntubes) & -! * (g0y - g1y) * scratch -! g0y = 0.0 g0y = g0y + g2y + g3y do it = 1, ntubes @@ -169,7 +157,6 @@ subroutine get_source_ffs_itteration (phi, g, source) end do call enforce_reality (source(:,:,:,:,ivmu)) -!! call center_zed(iv, source(:,:,:,:,ivmu)) end do deallocate(g0) diff --git a/fields.fpp b/fields.fpp index 23d2439bc4..3758baf21b 100644 --- a/fields.fpp +++ b/fields.fpp @@ -781,9 +781,9 @@ contains if (.not. allocated(gamtot)) allocate (gamtot(naky, nakx, -nzgrid:nzgrid)); gamtot = 0. 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 - end if +! if (zonal_mode(1) .and. akx(1) < epsilon(0.) .and. has_electron_species(spec)) then + gamtot(1, 1, :) = 0.0 + !end if if (.not. has_electron_species(spec)) then ia = 1 @@ -812,7 +812,7 @@ contains ! call test_band_lu_factorisation (gam0_ffs, lu_gam0_ffs) call band_lu_factorisation_ffs(gam0_ffs, lu_gam0_ffs) end if - + deallocate (wgts) deallocate (kperp2_swap) deallocate (aj0_alpha, gam0_alpha) @@ -1316,7 +1316,6 @@ contains call gyro_average(phiold, source2, gam0_ffs) source2 = source2 - gamtot_t * phiold - source = source - source2 where (gamtot_t < epsilon(0.0)) @@ -1330,7 +1329,7 @@ contains source(1, 1, :, :) = 0.0 end if -! call enforce_reality(source) + call enforce_reality(source) deallocate(source2, gamtot_t) @@ -1385,10 +1384,10 @@ contains !> store result in phi, which will be further modified below to account for polarization term call sum_allreduce(source) - allocate(source_copy(naky,nakx, -nzgrid:nzgrid, ntubes)) ; source_copy = 0.0 - source_copy = spread(source, 4, ntubes) - call enforce_reality (source_copy) - deallocate(source_copy) + ! allocate(source_copy(naky,nakx, -nzgrid:nzgrid, ntubes)) ; source_copy = 0.0 + ! source_copy = spread(source, 4, ntubes) + ! call enforce_reality (source_copy) + ! deallocate(source_copy) !> no longer need , so deallocate deallocate (gyro_g, gyro_g2) @@ -1437,6 +1436,7 @@ contains integer :: it, ia complex :: tmp + complex, dimension(:, :, :, :), allocatable :: source_copy allocate (source(naky, nakx, -nzgrid:nzgrid)); source = 0.0 if (fphi > epsilon(0.0)) then @@ -1470,7 +1470,6 @@ contains phi(1, 1, :, :) = 0.0 end if -! call enforce_reality (phi) else !> calculate the contribution to quasineutrality coming from the velocity space !> integration of the guiding centre distribution function g; @@ -1478,6 +1477,7 @@ contains !> this is returned in source if (debug) write (*, *) 'fields::advance_fields::get_fields_ffs::get_g_integral_contribution' call get_g_integral_contribution(g, source) + !> use sum_s int d3v and QN to solve for phi !> NB: assuming here that ntubes = 1 for FFS sim if (debug) write (*, *) 'fields::advance_fields::get_phi_ffs' @@ -1608,13 +1608,13 @@ contains call sum_allreduce(source) !!> Better fix when only on the implicit solve - if (present(implicit_solve)) then - allocate(source_copy(naky,nakx, -nzgrid:nzgrid, ntubes)) ; source_copy = 0.0 - source_copy = spread(source, 4, ntubes) - call enforce_reality (source_copy) - source = source_copy (:,:,:,1) - deallocate(source_copy) - end if + ! if (present(implicit_solve)) then + ! allocate(source_copy(naky,nakx, -nzgrid:nzgrid, ntubes)) ; source_copy = 0.0 + ! source_copy = spread(source, 4, ntubes) + ! call enforce_reality (source_copy) + ! source = source_copy (:,:,:,1) + ! deallocate(source_copy) + ! end if !> no longer need , so deallocate deallocate (gyro_g) diff --git a/gyro_averages.f90 b/gyro_averages.f90 index f98e4cab86..d35e5cbb73 100644 --- a/gyro_averages.f90 +++ b/gyro_averages.f90 @@ -11,7 +11,6 @@ module gyro_averages public :: find_max_required_kalpha_index public :: j1_ffs public :: j0_const, j0_B_const - public :: j0max_const private @@ -52,7 +51,7 @@ module gyro_averages logical :: debug = .false. type(coupled_alpha_type), dimension(:, :, :, :), allocatable :: j1_ffs - real, dimension(:, :, :, :), allocatable :: j0_const, j0_B_const, j0max_const + real, dimension(:, :, :, :), allocatable :: j0_const, j0_B_const contains @@ -293,9 +292,6 @@ subroutine init_bessel_ffs complex, dimension(:, :), allocatable :: j0_const_in_kalpha, j0_B_const_in_kalpha complex, dimension(:, :), allocatable :: j0_const_c, j0_B_const_c - real, dimension(:), allocatable :: j0max - complex, dimension(:, :), allocatable :: j0max_const_in_kalpha, j0max_const_c - ! call open_output_file (j0_ffs_unit, '.j0_ffs') ! call open_output_file (j0_B_ffs_unit, '.j0_over_B_ffs') @@ -331,11 +327,6 @@ subroutine init_bessel_ffs allocate (j0_const(naky, nakx, -nzgrid:nzgrid, vmu_lo%llim_proc:vmu_lo%ulim_alloc)); j0_const = 0.0 allocate (j0_B_const(naky, nakx, -nzgrid:nzgrid, vmu_lo%llim_proc:vmu_lo%ulim_alloc)); j0_B_const = 0.0 - allocate (j0max(nalpha)) - allocate (j0max_const_in_kalpha(naky_all, ikx_max)); j0max_const_in_kalpha = 0.0 - allocate (j0max_const_c(naky, nakx)); j0max_const_c = 0.0 - allocate (j0max_const(naky, nakx, -nzgrid:nzgrid, vmu_lo%llim_proc:vmu_lo%ulim_alloc)); j0max_const = 0.0 - ia_max_j0_count = 0; ia_max_j0_B_count = 0 do iz = -nzgrid, nzgrid ! if (proc0) write (*, *) 'calculating Fourier coefficients needed for gyro-averaging with alpha variation; zed index: ', iz @@ -361,7 +352,6 @@ subroutine init_bessel_ffs !> compute J_0*B*exp(-v^2), needed when integrating g over v-space in Maxwell's equations, !> due to B in v-space Jacobian j0_B(ia) = aj0_alpha(ia) * bmag(ia, iz) - j0max = aj0_alpha(ia) * maxwell_vpa(iv, is) * maxwell_mu(ia, iz, imu, is) aj1_alpha(ia) = j1(arg) end do @@ -371,10 +361,9 @@ subroutine init_bessel_ffs call transform_alpha2kalpha(aj0_alpha, aj0_kalpha) call transform_alpha2kalpha(j0_B, j0_B_kalpha) call transform_alpha2kalpha(aj1_alpha, aj1_kalpha) - + j0_const_in_kalpha(iky, ikx) = aj0_kalpha(1) j0_B_const_in_kalpha(iky, ikx) = j0_B_kalpha(1) - j0max_const_in_kalpha(iky, ikx) = j0max(1) !> given the Fourier coefficients aj0_kalpha, calculate the minimum number of coefficients needed, !> called j0_ffs%max_idx, to ensure that the relative error in the total spectral energy is below a specified tolerance @@ -412,8 +401,6 @@ subroutine init_bessel_ffs j0_const(:, :, iz, ivmu) = real(j0_const_c) call swap_kxky_back_ordered(j0_B_const_in_kalpha, j0_B_const_c) j0_B_const(:, :, iz, ivmu) = real(j0_B_const_c) - call swap_kxky_back_ordered(j0max_const_in_kalpha, j0max_const_c) - j0max_const(:, :, iz, ivmu) = real(j0max_const_c) end do end do @@ -423,8 +410,6 @@ subroutine init_bessel_ffs deallocate (j0_const_in_kalpha, j0_const_c) deallocate (j0_B_const_in_kalpha, j0_B_const_c) - deallocate (j0max_const_in_kalpha, j0max_const_c, j0max) - !> calculate the reduction factor of Fourier modes !> used to represent J0 !> avoid overflow by converting integers to reals before multiplying @@ -629,7 +614,6 @@ subroutine finish_bessel if (allocated(j0_B_const)) deallocate (j0_B_const) if (allocated(j0_const)) deallocate (j0_const) - if (allocated(j0max_const)) deallocate (j0max_const) bessinit = .false. @@ -1053,7 +1037,6 @@ subroutine band_lu_factorisation_ffs(gam0, lu_gam0) use common_types, only: coupled_alpha_type, gam0_ffs_type use zgrid, only: nzgrid use kt_grids, only: ikx_max, naky_all, naky - implicit none type(coupled_alpha_type), dimension(:, :, -nzgrid:), intent(in) :: gam0 diff --git a/implicit_solve.f90 b/implicit_solve.f90 index 8ce49c6979..16f9fb4474 100644 --- a/implicit_solve.f90 +++ b/implicit_solve.f90 @@ -116,6 +116,7 @@ subroutine advance_implicit_terms(g, phi, apar, bpar) !> 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 if (driftkinetic_implicit .and. (itt .NE. 1)) then + fields_updated = .false. call advance_fields(g2, phi, apar, bpar, dist=trim(dist_choice)) phi_old = phi end if @@ -197,7 +198,6 @@ subroutine advance_implicit_terms(g, phi, apar, bpar) else call update_pdf end if - !! change !!error = 0.0 itt = itt + 1 @@ -476,8 +476,8 @@ subroutine add_streaming_contribution_phi ! and store in dummy variable z_scratch z_scratch = maxwell_mu(ia, :, imu, is) call center_zed(iv, z_scratch, -nzgrid) - else - z_scratch = 1.0 + else + z_scratch = 1.0 end if end if ! multiply by Maxwellian factor @@ -527,7 +527,7 @@ subroutine add_drifts_contribution_phi ikx = ikx_from_izext(izext) iz = iz_from_izext(izext) scratch(izext) = zi * scratch(izext) * (akx(ikx) * wdriftx_phi(ia, iz, ivmu) & - + aky(iky) * (wdrifty_phi(ia, iz, ivmu) + wstar(ia, iz, ivmu))) + + aky(iky) * (wdrifty_phi(ia, iz, ivmu) + wstar(ia, iz, ivmu))) end do call center_zed(iv, scratch, 1, periodic(iky)) @@ -1250,7 +1250,7 @@ subroutine invert_parstream_response(phi, apar, bpar) ! unpack phi (and apar if evolving) from the field vector; ! also apply phase shift at periodic point phi(iky, ikx, -nzgrid:nzgrid - 1, it) = fld(:nzp) - phi(iky, ikx, nzgrid, it) = phi(iky, ikx, -nzgrid, it) / phase_shift(iky) + phi(iky, ikx, nzgrid, it) = phi(iky, ikx, -nzgrid, it) * phase_shift(iky) if (include_apar) then apar(iky, ikx, -nzgrid:nzgrid - 1, it) = fld(offset_apar + 1:nzp + offset_apar) apar(iky, ikx, nzgrid, it) = apar(iky, ikx, -nzgrid, it) / phase_shift(iky) diff --git a/init_g.f90 b/init_g.f90 index 77e76a2c93..fd1406c260 100644 --- a/init_g.f90 +++ b/init_g.f90 @@ -338,43 +338,44 @@ 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 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 - end do + ! 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 + + ! call map_from_extended_zgrid(it, ie, iky, phiext, phi(iky, :, :, :)) + ! deallocate(phiext, zed_ext) + ! end do + ! end do + ! 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)) / width0)**2) * cmplx(1.0, 1.0) - ! end do + do iz = -nzgrid, nzgrid + phi(:, :, iz,1) = exp(-((zed(iz)) / width0)**2) * cmplx(1.0, 1.0) + end do phi(1, 1, :, :) = 0.0 - - ! if (zonal_mode(1)) then - ! if (abs(akx(1)) < epsilon(0.0)) then - ! phi(1, 1, :, :) = 0.0 - ! end if - - ! if (reality) then - ! do ikx = 1, nakx - ikx_max - ! phi(1, nakx - ikx + 1, :, :) = conjg(phi(1, ikx + 1, :, :)) - ! end do - ! end if - ! end if + + if (zonal_mode(1)) then + if (abs(akx(1)) < epsilon(0.0)) then + phi(1, 1, :, :) = 0.0 + end if + + if (reality) then + do ikx = 1, nakx - ikx_max + phi(1, nakx - ikx + 1, :, :) = conjg(phi(1, ikx + 1, :, :)) + end do + end if + end if allocate (g_swap(naky_all, ikx_max)) allocate (phiy(ny, ikx_max, -nzgrid:nzgrid)) ; phiy = 0.0 @@ -382,7 +383,7 @@ subroutine ginit_default_ffs call swap_kxky(phi(:, :, iz, 1), g_swap) call transform_ky2y(g_swap, phiy(:, :, iz)) end do - + allocate (g0x(ny, ikx_max, -nzgrid:nzgrid, vmu_lo%llim_proc:vmu_lo%ulim_alloc)); g0x = 0.0 do ivmu = vmu_lo%llim_proc, vmu_lo%ulim_proc iv = iv_idx(vmu_lo, ivmu) @@ -390,16 +391,17 @@ subroutine ginit_default_ffs is = is_idx(vmu_lo, ivmu) do iy = 1, ny do ikx = 1, ikx_max - do iz = -nzgrid, nzgrid + do iz = -nzgrid, nzgrid g0x(iy, ikx, iz, ivmu) = phiinit * 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 end do end do - + + gnew = 0.0 do ivmu = vmu_lo%llim_proc, vmu_lo%ulim_proc do iz = -nzgrid, nzgrid call transform_y2ky(g0x(:, :, iz, ivmu), g_swap) @@ -408,6 +410,7 @@ subroutine ginit_default_ffs end do gnew = spread(gnew(:,:,:,1,:), 4, ntubes) gnew (1,1,:,:,:) = 0.0 + do ivmu = vmu_lo%llim_proc, vmu_lo%ulim_proc call enforce_reality (gnew(:, :, :, :, ivmu)) end do