The GKE equation and Poissons Equation should be correct to Delta^1 delta^1 in the local limit, assuming rare collisions in the equilibrium solution. This hasn't been thoroughly tested yet. Further corrections are needed when the rare collisions assumption is lifted. Should be able to account for this by passing the g_neo distribution to all operators on the LHS of the GKE, rather than just the time derivative which is currently the case.

Author: OJ Robertson

STELLA_CODE/arrays/arrays.f90

@@ -28,6 +28,7 @@ module arrays
    public :: initialised_wstar1
    public :: initialised_wpol
    public :: initialised_neo_curv_drift
+   public :: initialised_neo_mag_drift
 
    !----------------------------------------------------------------------------
    ! For the Gyrokinetic Equation
@@ -75,6 +76,8 @@ module arrays
    public :: wpol
    public :: neocurvx
    public :: neocurvy
+   public :: neomagx
+   public :: neomagy
 
    !----------------------------------------------------------------------------
    ! For the Field Equations
@@ -83,8 +86,6 @@ module arrays
    ! Arrays used to calculate the fields for electrostatic simulations
    public :: denominator_fields
    public :: denominator_fields_h
-   ! For NEO's higher order corrections. 
-   public :: denominator_fields_neo_adiab
 
    ! Arrays used to calculate the fields for electrostatic simulations
    ! considering a Modified Boltzmann Response for the electrons
@@ -120,6 +121,7 @@ module arrays
    logical :: initialised_wstar1
    logical :: initialised_wpol
    logical :: initialised_neo_curv_drift
+   logical :: initialised_neo_mag_drift
 
    !----------------------------------------------------------------------------
    ! For the Gyrokinetic Equation
@@ -159,6 +161,7 @@ module arrays
    real, dimension(:, :, :), allocatable :: neo_dchidz_coeff
    real, dimension(:, :, :), allocatable :: wstar1, wpol
    real, dimension(:, :, :), allocatable :: neocurvx, neocurvy
+   real, dimension(:, :, :), allocatable :: neomagx, neomagy
 
    !----------------------------------------------------------------------------
    ! For the Field Equations

STELLA_CODE/calculations/calculations_timestep.f90

@@ -64,6 +64,7 @@ subroutine init_cfl
       use arrays, only: neo_chi_coeff, neo_apar_coeff, neo_dchidz_coeff
       use arrays, only: wstar1, wpol
       use arrays, only: neocurvx, neocurvy
+      use arrays, only: neomagx, neomagy
 
       ! Collisions
       use dissipation_and_collisions, only: include_collisions, collisions_implicit
@@ -90,10 +91,11 @@ subroutine init_cfl
       real :: cfl_dt_neo_chi_coeff, cfl_dt_neo_apar_coeff, cfl_dt_neo_dchidz_coeff
       real :: cfl_dt_wstar1, cfl_dt_wpol
       real :: cfl_dt_neocurvx, cfl_dt_neocurvy
+      real :: cfl_dt_neomagx, cfl_dt_neomagy
       real :: neo_chi_coeff_max, neo_apar_coeff_max, neo_dchidz_coeff_max
       real :: wstar1_max, wpol_max
       real :: neocurvx_max, neocurvy_max      
-
+      real :: neomagx_max, neomagy_max
 
       !-------------------------------------------------------------------------
 
@@ -236,7 +238,7 @@ subroutine init_cfl
 
       ! Check that the introduction of the neoclassical wpol drive doesn't break the CFL condition.
       if (neoclassical_is_enabled()) then
-         !  Only calculate the CFL constaint if there are non-zero akx present. 
+         ! Only calculate the CFL constaint if there are non-zero akx present. 
           if (maxval(abs(akx)) > epsilon(0.0)) then
               wpol_max = maxval(abs(wpol))
               if (nproc > 1) then
@@ -273,6 +275,33 @@ subroutine init_cfl
           end if
       end if
 
+
+      ! Check that the introduction of the neoclassical neomagy doesn't break the CFL condition.
+      if (neoclassical_is_enabled()) then
+          neomagy_max = maxval(abs(neomagy))
+          if (nproc > 1) then
+              call max_allreduce(neomagy_max)
+          end if
+
+          cfl_dt_neomagy = abs(code_dt) / max(maxval(abs(aky)) * neomagy_max, zero)
+          cfl_dt_linear = min(cfl_dt_linear, cfl_dt_neomagy)
+      end if
+
+      ! Check that the introduction of the neoclassical neocurvx drive doesn't break the CFL condition.
+      if (neoclassical_is_enabled()) then
+          ! Only calculate the CFL constaint if there are non-zero akx present. 
+          if (maxval(abs(akx)) > epsilon(0.0)) then
+              neomagx_max = maxval(abs(neomagx))
+              if (nproc > 1) then
+                  call max_allreduce(neomagx_max)
+              end if
+
+              cfl_dt_neomagx = abs(code_dt) / max(maxval(abs(akx)) * neomagx_max, zero)
+              cfl_dt_linear = min(cfl_dt_linear, cfl_dt_neomagx)
+          end if
+      end if
+
+
       ! Get the minimum value of <cfl_dt_linear> on all processors
       if (runtype_option_switch == runtype_multibox) call scope(allprocs)
       call min_allreduce(cfl_dt_linear)
@@ -295,7 +324,9 @@ subroutine init_cfl
          if (neoclassical_is_enabled()) write (*, '(A12,ES12.4)') 'wstar1: ', cfl_dt_wstar1
          if (neoclassical_is_enabled() .and. maxval(abs(akx)) > epsilon(0.0)) write (*, '(A12,ES12.4)') 'wpol: ', cfl_dt_wpol
          if (neoclassical_is_enabled()) write (*, '(A12,ES12.4)') 'neocurvy: ', cfl_dt_neocurvy
-         if (neoclassical_is_enabled() .and. maxval(abs(akx)) > epsilon(0.0)) write (*, '(A12,ES12.4)') 'neocurvx: ', cfl_dt_neocurvx 
+         if (neoclassical_is_enabled() .and. maxval(abs(akx)) > epsilon(0.0)) write (*, '(A12,ES12.4)') 'neocurvx: ', cfl_dt_neocurvx
+         if (neoclassical_is_enabled()) write (*, '(A12,ES12.4)') 'neomagy: ', cfl_dt_neomagy
+         if (neoclassical_is_enabled() .and. maxval(abs(akx)) > epsilon(0.0)) write (*, '(A12,ES12.4)') 'neomagx: ', cfl_dt_neomagx 
          write (*, '(A12,ES12.4)') '   total: ', cfl_dt_linear
          write (*, *)
       end if

STELLA_CODE/calculations/calculations_tofrom_ghf.f90

@@ -846,14 +846,17 @@ subroutine g_or_gbar_to_gbarneo_kxkyz(g0, phi, apar, bpar, facchi)
             ! Calculate the phi contribution to gyoraveraged fields.
 
             field = fphi * facchi * spec(is)%zt * phi(iky, ikx, iz, it) * spread(maxwell_vpa(:, is), 2, nmu) &
-               * spread(maxwell_mu(ia, iz, :, is), 1, nvpa) * maxwell_fac(is) 
+               * spread(maxwell_mu(ia, iz, :, is), 1, nvpa)
 
             ! If apar is present, we must account for this.
             if (include_apar) then 
                 field = field - 2.0 * facchi * spec(is)%zt * apar(iky, ikx, iz, it) * spread(vpa, 2, nmu) * spread(maxwell_vpa(:, is), 2, nmu) &
-                * spread(maxwell_mu(ia, iz, :, is), 1, nvpa) * maxwell_fac(is) * spec(is)%stm_psi0 
+                * spread(maxwell_mu(ia, iz, :, is), 1, nvpa) * spec(is)%stm_psi0 
             end if
 
+            ! Multiply by the neoclassical distribution factor. 
+            ! field = field * ( ( 0.5/bmag(ia, iz) ) * dneo_h_dmu(iz, ivmu, 1) - ( neo_h(iz, ivmu, 1) - spec(is)%z * neo_phi(iz) ) )
+
             ! Gyroaverage.
             call gyro_average(field, ikxkyz, gyro_averaged_field)
 
@@ -872,8 +875,11 @@ subroutine g_or_gbar_to_gbarneo_kxkyz(g0, phi, apar, bpar, facchi)
                 is = is_idx(kxkyz_lo, ikxkyz)
 
                 field = 4.0 * facchi * spread(mu, 1, nvpa) * bpar(iky, ikx, iz, it) * spread(maxwell_vpa(:, is), 2, nmu) &
-                * spread(maxwell_mu(ia, iz, :, is), 1, nvpa) * maxwell_fac(is)
+                * spread(maxwell_mu(ia, iz, :, is), 1, nvpa)
 
+                ! Multiply by the neoclassical distribution factor. 
+                ! field = field * ( ( 0.5/bmag(ia, iz) ) * dneo_h_dmu(iz, ivmu, 1) - ( neo_h(iz, ivmu, 1) - spec(is)%z * neo_phi(iz) ) )
+
                 ! Gyroaverage the J_1 contribution.
                 call gyro_average_j1(field, ikxkyz, gyro_averaged_field)
 
@@ -975,7 +981,8 @@ subroutine g_or_gbar_to_gbarneo_vmu_single(ivmu, g0, phi, apar, bpar, facchi)
                ! First calculate the gyroaveraged field being used, <>. 
                ! Calculate the phi contribution to gyoraveraged fields.
 
-               field = spec(is)%zt * facchi * phi (:, :, iz, it) * maxwell_vpa(iv, is) * maxwell_mu(ia, iz, imu, is) * maxwell_fac(is)
+               field = spec(is)%zt * facchi * phi (:, :, iz, it) * maxwell_vpa(iv, is) * maxwell_mu(ia, iz, imu, is) &
+               * maxwell_fac(is)
 
                ! If apar is present, we must account for this.
                if (include_apar) then 
@@ -985,8 +992,7 @@ subroutine g_or_gbar_to_gbarneo_vmu_single(ivmu, g0, phi, apar, bpar, facchi)
 
                ! Mutliply by the neoclassical distribution factor. 
 
-               field = field * ( ( 0.5/bmag(ia, iz) ) * dneo_h_dmu(iz, ivmu, 1) - ( neo_h(iz, ivmu, 1) - spec(is)%z * neo_phi(iz, 1) ) )
-
+               field = field * ( ( 0.5/bmag(ia, iz) ) * dneo_h_dmu(iz, ivmu, 1) - ( neo_h(iz, ivmu, 1) - spec(is)%z * neo_phi(iz) ) )
                ! Gyroaverage the J₀ contribution.
                call gyro_average(field, iz, ivmu, gyro_averaged_field)
 
@@ -1002,8 +1008,11 @@ subroutine g_or_gbar_to_gbarneo_vmu_single(ivmu, g0, phi, apar, bpar, facchi)
            do it = 1, ntubes
                do iz = -nzgrid, nzgrid
 
-                   field = 4.0 * facchi * mu(imu) * bpar(:, :, iz, it) * maxwell_vpa(iv, is) * maxwell_mu(ia, iz, imu, is) * maxwell_fac(is)
+                   field = 4.0 * facchi * mu(imu) * bpar(:, :, iz, it) * maxwell_vpa(iv, is) * maxwell_mu(ia, iz, imu, is) &
+                   * maxwell_fac(is)
 
+                    field = field * ( ( 0.5/bmag(ia, iz) ) * dneo_h_dmu(iz, ivmu, 1) - ( neo_h(iz, ivmu, 1) - spec(is)%z * neo_phi(iz) ) )
+
                    ! Gyroaverage the J_1 contribution.
                    call gyro_average_j1(field, iz, ivmu, gyro_averaged_field)
 

STELLA_CODE/field_equations/field_equations_fluxtube.f90

@@ -628,7 +628,7 @@ subroutine init_field_equations_fluxtube
       use parallelisation_layouts, only: iz_idx, it_idx, ikx_idx, iky_idx, is_idx
       use parallelisation_layouts, only: iv_idx, imu_idx
       ! Arrays
-      use arrays, only: denominator_fields, denominator_fields_neo_adiab, denominator_fields_MBR
+      use arrays, only: denominator_fields, denominator_fields_MBR
       use arrays, only: denominator_fields_h, denominator_fields_MBR_h, efac, efacp
       use arrays_gyro_averages, only: aj0v
 
@@ -659,15 +659,15 @@ subroutine init_field_equations_fluxtube
       ! NEO data for higher order corrections. 
       use neoclassical_terms_neo, only: neoclassical_is_enabled
       use neoclassical_terms_neo, only: neo_phi, neo_h, neo_phi, dneo_h_dmu
+      use neoclassical_terms_neo, only: neo_dens 
 
       implicit none
 
       ! Local variables
       integer :: ikxkyz, iz, it, ikx, iky, ia, iv, imu, ivmu
       integer :: is, is_inner, is_outer
-      real :: tmp, tmp_neo_adiab, wgt
+      real :: tmp, wgt
       real, dimension(:, :), allocatable :: g0
-      real, dimension(:, :), allocatable :: g0_neo_adiab
 
       !-------------------------------------------------------------------------
 
@@ -696,9 +696,9 @@ subroutine init_field_equations_fluxtube
          ! ============================================================================================================= !
          !
          ! If NEO's corrections are enabled, the denominator for phi calculations picks up a correction. 
-         ! This is given by: 
+         ! This is given by: ... 
          !
-         ! sum_s Z_s n_s [ (2B/sqrt(pi)) int dvpa int dmu J_0 * g ] 
+         ! 
          !
          !
          ! ============================================================================================================= !         
@@ -724,7 +724,7 @@ subroutine init_field_equations_fluxtube
                     imu = imu_idx(vmu_lo, ivmu)
                     iv = iv_idx(vmu_lo, ivmu)
 
-                    g0(iv, imu) = g0(iv, imu) * ( 1 + neo_h(iz, ivmu, 1) - spec(is_inner)%z * neo_phi(iz, 1) - ( 0.5/bmag(ia, iz) ) * dneo_h_dmu(iz, ivmu, 1) )
+                    g0(iv, imu) = g0(iv, imu) * ( 1 + neo_h(iz, ivmu, 1) - spec(is_inner)%z * neo_phi(iz) - ( 0.5/bmag(ia, iz) ) * dneo_h_dmu(iz, ivmu, 1) )
                 end do
             end if
 
@@ -892,56 +892,26 @@ subroutine init_field_equations_fluxtube
 
             ! If NEO's corrections are enabled, compute the higher order correction to the denominator. 
             if(neoclassical_is_enabled()) then 
-                allocate(g0_neo_adiab(nvpa, nmu))
-                allocate(denominator_fields_neo_adiab(nakx, naky, -nzgrid:nzgrid))
-                denominator_fields_neo_adiab = 0.0
-
-                 do ikxkyz = kxkyz_lo%llim_proc, kxkyz_lo%ulim_proc
-                     it = it_idx(kxkyz_lo, ikxkyz)
-                     if (it /= 1) cycle
-                     iky = iky_idx(kxkyz_lo, ikxkyz)
-                     ikx = ikx_idx(kxkyz_lo, ikxkyz)
-                     iz = iz_idx(kxkyz_lo, ikxkyz)
-                     is_outer = is_idx(kxkyz_lo, ikxkyz)
-                     if (is_outer /= 2) cycle
-                 
-                     ! Calculate exp(v²) for each (kx,ky,z).
-                     g0_neo_adiab = spread(maxwell_vpa(:, is_outer), 2, nmu) * spread(maxwell_mu(ia, iz, :, is_outer), 1, nvpa) * maxwell_fac(is_outer)
-
-                     ! Iterate over velocity space.
-                     do ivmu = vmu_lo%llim_proc, vmu_lo%ulim_proc
-                         is_inner = is_idx(vmu_lo, ivmu)
-                         if (is_inner /= is_outer) cycle
-                         imu = imu_idx(vmu_lo, ivmu)
-                         iv = iv_idx(vmu_lo, ivmu)
-
-                         g0_neo_adiab(iv, imu) = g0_neo_adiab(iv, imu) * ( neo_h(iz, ivmu, 1) - spec(is_inner)%z * neo_phi(iz, 1) )
-                     end do
-
-                     ! Calculate denominator_fields_neo_adiab[iky,ikz,iz].
-                     call integrate_vmu(g0_neo_adiab, iz, tmp_neo_adiab)
-                     denominator_fields_neo_adiab(ikx, iky, iz) = denominator_fields_neo_adiab(ikx, iky, iz) + tmp_neo_adiab
-                 end do
-
-                 ! Sum the values on all processors and send them to <proc0>.
-                 call sum_allreduce(denominator_fields_neo_adiab)
-
-                 ! Calculate the approproate denominators for g and h in the presence of the higher order equilibriu. 
-                 denominator_fields = denominator_fields + efac * ( 1 + denominator_fields_neo_adiab ) 
-                 ! denominator_fields_h = denominator_fields_h + efac * ( 1 + denominator_fields_neo_adiab )
-             else 
-                 ! Otherwise, calculate the approporiate denominators for g and h in the presence of the Maxwellian equilbirium. 
-                 efacp = efac * (spec(ion_species)%tprim - spec(ion_species)%fprim)
+                ! Calculate the approproate denominators for g and h in the presence of the higher order equilibrium.
+                do iz = -nzgrid, nzgrid 
+                    efacp = efac * (spec(ion_species)%tprim - spec(ion_species)%fprim)
+
+                    denominator_fields(:, :, iz) = denominator_fields(:, :, iz) + efac * ( 1 + neo_dens(iz, 2) )
+                    ! denominator_fields_h(:, :, iz) = denominator_fields_h(:, :, iz) + efac * ( 1 + neo_dens(iz, 2) )
+                end do
+            else 
+                ! Otherwise, calculate the approporiate denominators for g and h in the presence of the Maxwellian equilbirium. 
+                efacp = efac * (spec(ion_species)%tprim - spec(ion_species)%fprim)
 
-                 ! Add the contribution of adiabatic electrons to <denominator_fields>
-                 ! Calculate denominator_fields = sum_(s not e) (Z_s² n_s/T_s) (1 - Gamma0) + (n_e/T_e)
-                 ! This is the factor that multiplies phi in the final expression. 
-                 denominator_fields = denominator_fields + efac
+                ! Add the contribution of adiabatic electrons to <denominator_fields>
+                ! Calculate denominator_fields = sum_(s not e) (Z_s² n_s/T_s) (1 - Gamma0) + (n_e/T_e)
+                ! This is the factor that multiplies phi in the final expression. 
+                denominator_fields = denominator_fields + efac
 
-                 ! Add the contribution of adiabatic electrons to <denominator_fields_h>
-                 ! Calculate denominator_fields_h = sum_(s not e) (Z_s² n_s/T_s) + (n_e/T_e)
-                 denominator_fields_h = denominator_fields_h + efac
-             end if
+                ! Add the contribution of adiabatic electrons to <denominator_fields_h>
+                ! Calculate denominator_fields_h = sum_(s not e) (Z_s² n_s/T_s) + (n_e/T_e)
+                denominator_fields_h = denominator_fields_h + efac
+            end if
 
             !---------------------------------------------------------------------
             !                 Modified Boltzmann Response (MBR)
@@ -978,8 +948,6 @@ subroutine init_field_equations_fluxtube
 
          ! Deallocate temporary arrays.
          if (allocated(g0)) deallocate (g0)
-         if (allocated(g0_neo_adiab)) deallocate (g0_neo_adiab)
-
       end if
 
    end subroutine init_field_equations_fluxtube

STELLA_CODE/geometry/geometry.f90

@@ -86,6 +86,9 @@ module geometry
    public :: aref, bref, twist_and_shift_geo_fac
    public :: q_as_x, get_x_to_rho, gfac
    public :: dVolume, grad_x_grad_y_end
+
+   ! For NEO's neoclassical calculations, we need b_dot_gradtheta. 
+   public :: b_dot_gradtheta
 
    ! Flux tube only needs b_dot_gradz_avg(z)
    public :: b_dot_gradz_avg
@@ -142,6 +145,9 @@ module geometry
    ! Geometric quantities for full flux surface
    real, dimension(:, :), allocatable :: b_dot_gradz
 
+   ! Geometric quantites for NEO's higher order corrections. 
+   real, dimension(:, :), allocatable :: b_dot_gradtheta
+
    ! Geometric quantities for the momentum flux
    real, dimension(:, :), allocatable :: gradzeta_gradx_R2overB2
    real, dimension(:, :), allocatable :: gradzeta_grady_R2overB2
@@ -402,6 +408,9 @@ subroutine get_geometry_arrays_from_VMEC(nalpha, naky)
       real, dimension(:, :), allocatable :: grad_x_grad_x, grad_y_grad_y, grad_y_grad_x
       real, dimension(:, :), allocatable :: gradzeta_gradpsit_R2overB2, gradzeta_gradalpha_R2overB2
 
+      ! For NEO's neoclassical corrections. 
+      real, dimension(:, :), allocatable :: b_dot_gradtheta_arr
+
       ! Local variables
       real :: rho, shat, iota, field_period_ratio
       real :: dpsidpsit, dpsidx, dpsitdx, dxdpsit 
@@ -1449,6 +1458,9 @@ subroutine finish_geometry
       if (allocated(gradzeta_grady_R2overB2)) deallocate (gradzeta_grady_R2overB2)
       if (allocated(b_dot_gradzeta_RR)) deallocate (b_dot_gradzeta_RR)
 
+      ! Needed for NEO's higher order corrections.
+      ! if (allocated(b_dot_gradtheta_arr)) deallocate (b_dot_gradtheta_arr)
+
       ! Only initialise once
       initialised_geometry = .false.
 

STELLA_CODE/geometry/geometry_miller.f90

@@ -620,7 +620,8 @@ subroutine get_miller_geometry(nzgrid, zed_in, zed_equal_arc, &
 
       ! Calculate <b_dot_gradtheta> = b . ∇θ (formerly <gradpar>)
       b_dot_gradtheta = dpsip_drho / (bmag * jacrho)
-      
+     
+
       ! Calculate <b_dot_gradB> = b . ∇B (formerly <gradparB>)
       b_dot_gradB = b_dot_gradtheta * dB_dtheta
 
@@ -798,7 +799,7 @@ subroutine interpolate_functions
             call geo_spline(theta, d_gradxdotgrady_drho, zed_arc, d_gradxdotgrady_drho_out)
             call geo_spline(theta, d_gradxdotgradx_drho, zed_arc, d_gradxdotgradx_drho_out)
             call geo_spline(theta, djac_drho / dpsip_drho, zed_arc, djac_drho_out)
-
+            
             deallocate (zed_arc)
          else
             if (debug) write (*, *) 'geometry_miller::zed_equal_arc=.false.'
@@ -945,7 +946,7 @@ subroutine allocate_arrays(nr, nz)
       ! For the neoclassical terms
       allocate (gds23(-nz:nz)); gds23 = 0.0
       allocate (gds24(-nz:nz)); gds24 = 0.0
-      
+
       ! The R(r, θ) and Z(r, θ) positions of a Miller equilibrium
       ! Here we have nr = 3 for the radial derivatives
       allocate (Rr(nr, -nz:nz)); Rr = 0.0
@@ -1135,7 +1136,7 @@ subroutine deallocate_arrays
          if (allocated(delta_theta)) deallocate (delta_theta)
          if (allocated(bmag_psi0)) deallocate (bmag_psi0)
          if (allocated(gradrho_psi0)) deallocate (gradrho_psi0)
-         
+
       end subroutine deallocate_arrays
 
    end subroutine finish_miller_geometry