update ffs ke

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
 
 

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 <phi>/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,28 +134,19 @@ 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
            call center_zed(iv, coeff(ia, :) ,  -nzgrid)
         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)

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 <g>, 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,14 +1470,14 @@ 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;
             !> the sign is consistent with phi appearing on the RHS of the eqn and int g appearing on the LHS.
             !> 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 <g> 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 <g>, so deallocate
          deallocate (gyro_g)
 

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

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)