Merge pull request #187 from mabarnes/feature-gyroaverages

Feature gyroaverages

examples/gk-ions/2D-periodic-gk.toml

@@ -0,0 +1,87 @@
+n_ion_species = 1
+n_neutral_species = 0
+boltzmann_electron_response = true
+evolve_moments_density = false
+evolve_moments_parallel_flow = false
+evolve_moments_parallel_pressure = false
+evolve_moments_conservation = false
+gyrokinetic_ions = true
+T_e = 1.0
+T_wall = 1.0
+initial_density1 = 1.0
+initial_temperature1 = 1.0
+z_IC_option1 = "gaussian"
+z_IC_density_amplitude1 = 0.001
+z_IC_density_phase1 = 0.0
+z_IC_upar_amplitude1 = 1.0
+z_IC_upar_phase1 = 0.0
+z_IC_temperature_amplitude1 = 0.0
+z_IC_temperature_phase1 = 0.0
+vpa_IC_option1 = "gaussian"
+vpa_IC_density_amplitude1 = 1.0
+vpa_IC_density_phase1 = 0.0
+vpa_IC_upar_amplitude1 = 0.0
+vpa_IC_upar_phase1 = 0.0
+vpa_IC_temperature_amplitude1 = 0.0
+vpa_IC_temperature_phase1 = 0.0
+initial_density2 = 1.0
+initial_temperature2 = 1.0
+z_IC_option2 = "gaussian"
+z_IC_density_amplitude2 = 0.001
+z_IC_density_phase2 = 0.0
+z_IC_upar_amplitude2 = -1.0
+z_IC_upar_phase2 = 0.0
+z_IC_temperature_amplitude2 = 0.0
+z_IC_temperature_phase2 = 0.0
+vpa_IC_option2 = "gaussian"
+vpa_IC_density_amplitude2 = 1.0
+vpa_IC_density_phase2 = 0.0
+vpa_IC_upar_amplitude2 = 0.0
+vpa_IC_upar_phase2 = 0.0
+vpa_IC_temperature_amplitude2 = 0.0
+vpa_IC_temperature_phase2 = 0.0
+charge_exchange_frequency = 0.5
+ionization_frequency = 0.05
+constant_ionization_rate = true
+nstep = 50
+dt = 1.0e-3
+nwrite = 10
+nwrite_dfns = 10
+use_semi_lagrange = false
+n_rk_stages = 4
+split_operators = false
+r_ngrid = 5
+r_nelement = 2
+r_bc = "periodic"
+z_ngrid = 5
+z_nelement = 2
+z_bc = "periodic"
+z_discretization = "chebyshev_pseudospectral"
+vpa_ngrid = 5
+vpa_nelement = 2
+vpa_L = 6.0
+vpa_bc = "zero"
+vpa_discretization = "chebyshev_pseudospectral"
+vperp_ngrid = 5
+vperp_nelement = 1
+vperp_L = 3.0
+vperp_bc = "zero"
+vperp_discretization = "chebyshev_pseudospectral"
+vz_ngrid = 9
+vz_nelement = 64
+vz_L = 18.0
+vz_bc = "both_zero"
+vz_discretization = "chebyshev_pseudospectral"
+
+[numerical_dissipation]
+vpa_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+vperp_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+#r_disspipation_coefficient = 1.0e-3
+#force_minimum_pdf_value = 0.0
+
+[geometry]
+#option="1D-mirror"
+DeltaB=0.0
+option="constant-helical"
+pitch=0.1
+rhostar= 0.1

makie_post_processing/makie_post_processing/src/shared_utils.jl

@@ -88,6 +88,7 @@ function get_composition(scan_input)
     phi_wall = get(scan_input, "phi_wall", 0.0)
     # if false use true Knudsen cosine for neutral wall bc
     use_test_neutral_wall_pdf = get(scan_input, "use_test_neutral_wall_pdf", false)
+    gyrokinetic_ions = get(scan_input, "gyrokinetic_ions", false)
     # constant to be used to test nonzero Er in wall boundary condition
     Er_constant = get(scan_input, "Er_constant", 0.0)
     recycling_fraction = get(scan_input, "recycling_fraction", 1.0)
@@ -106,7 +107,7 @@ function get_composition(scan_input)
     me_over_mi = 1.0/1836.0
     composition = species_composition(n_species, n_ion_species, n_neutral_species,
         electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant,
-        mn_over_mi, me_over_mi, recycling_fraction, allocate_float(n_species))
+        mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species))
     return composition
 
 end

moment_kinetics/debug_test/gyroaverage_inputs.jl

@@ -0,0 +1,87 @@
+test_type = "gyroaverage"
+
+# default inputs for tests
+test_input = Dict(
+    "run_name" => "gyroaverage",
+    "n_ion_species" => 1,
+    "n_neutral_species" => 0,
+    "evolve_moments_density" => false,
+    "evolve_moments_parallel_flow" => false,
+    "evolve_moments_parallel_pressure" => false,
+    "evolve_moments_conservation" => false,
+    "gyrokinetic_ions" => true,
+    "T_e" => 1.0,
+    "T_wall" => 1.0,
+    "initial_density1" => 1.0,
+    "initial_temperature1" => 1.0,
+    "z_IC_option1" => "gaussian",
+    "z_IC_density_amplitude1" => 0.001,
+    "z_IC_density_phase1" => 0.0,
+    "z_IC_upar_amplitude1" => 1.0,
+    "z_IC_upar_phase1" => 0.0,
+    "z_IC_temperature_amplitude1" => 0.0,
+    "z_IC_temperature_phase1" => 0.0,
+    "vpa_IC_option1" => "gaussian",
+    "vpa_IC_density_amplitude1" => 1.0,
+    "vpa_IC_density_phase1" => 0.0,
+    "vpa_IC_upar_amplitude1" => 0.0,
+    "vpa_IC_upar_phase1" => 0.0,
+    "vpa_IC_temperature_amplitude1" => 0.0,
+    "vpa_IC_temperature_phase1" => 0.0,
+    "initial_density2" => 1.0,
+    "initial_temperature2" => 1.0,
+    "z_IC_option2" => "gaussian",
+    "z_IC_density_amplitude2" => 0.001,
+    "z_IC_density_phase2" => 0.0,
+    "z_IC_upar_amplitude2" => -1.0,
+    "z_IC_upar_phase2" => 0.0,
+    "z_IC_temperature_amplitude2" => 0.0,
+    "z_IC_temperature_phase2" => 0.0,
+    "vpa_IC_option2" => "gaussian",
+    "vpa_IC_density_amplitude2" => 1.0,
+    "vpa_IC_density_phase2" => 0.0,
+    "vpa_IC_upar_amplitude2" => 0.0,
+    "vpa_IC_upar_phase2" => 0.0,
+    "vpa_IC_temperature_amplitude2" => 0.0,
+    "vpa_IC_temperature_phase2" => 0.0,
+    "charge_exchange_frequency" => 0.5,
+    "ionization_frequency" => 0.05,
+    "constant_ionization_rate" => true,
+    "nstep" => 3,
+    "dt" => 1.0e-12,
+    "nwrite" => 2,
+    "nwrite_dfns" => 2,
+    "n_rk_stages" => 4,
+    "r_ngrid" => 5,
+    "r_nelement" => 2,
+    "r_bc" => "periodic",
+    "z_ngrid" => 5,
+    "z_nelement" => 2,
+    "z_bc" => "periodic",
+    "z_discretization" => "chebyshev_pseudospectral",
+    "vpa_ngrid" => 5,
+    "vpa_nelement" => 2,
+    "vpa_L" => 6.0,
+    "vpa_bc" => "zero",
+    "vpa_discretization" => "chebyshev_pseudospectral",
+    "vperp_ngrid" => 5,
+    "vperp_nelement" => 1,
+    "vperp_L" => 3.0,
+    "vperp_bc" => "zero",
+    "vperp_discretization" => "chebyshev_pseudospectral",
+    "vz_ngrid" => 5,
+    "vz_nelement" => 2,
+    "vz_L" => 6.0,
+    "vz_bc" => "zero",
+    "vz_discretization" => "chebyshev_pseudospectral",
+    "numerical_dissipation" => Dict("vpa_dissipation_coefficient" => 1.0e-3,
+                                    "vperp_dissipation_coefficient" => 1.0e-3),
+    "geometry" => Dict("DeltaB"=>0.0,
+                       "option"=>"constant-helical",
+                       "pitch"=>0.1,
+                       "rhostar"=> 0.1),
+)
+
+test_input_list = [
+     test_input,
+    ]

moment_kinetics/debug_test/gyroaverage_tests.jl

@@ -0,0 +1,23 @@
+module GyroaverageDebug
+
+# Debug test using wall boundary conditions.
+
+include("setup.jl")
+
+# Create a temporary directory for test output
+test_output_directory = get_MPI_tempdir()
+mkpath(test_output_directory)
+
+
+# Input parameters for the test
+include("gyroaverage_inputs.jl")
+
+# Defines the test functions, using variables defined in the *_inputs.jl file
+include("runtest_template.jl")
+
+end # GyroaverageDebug
+
+
+using .GyroaverageDebug
+
+GyroaverageDebug.runtests()

moment_kinetics/debug_test/runtests.jl

@@ -10,6 +10,7 @@ function runtests()
         include(joinpath(@__DIR__, "harrisonthompson.jl"))
         include(joinpath(@__DIR__, "mms_tests.jl"))
         include(joinpath(@__DIR__, "restart_interpolation_tests.jl"))
+        include(joinpath(@__DIR__, "gyroaverage_tests.jl"))
     end
 end
 

moment_kinetics/src/chebyshev.jl

@@ -667,7 +667,7 @@ function chebyshev_radau_weights(moments::Array{mk_float,1}, n)
     # create array for moments on extended [0,2π] domain in theta = ArcCos[z]
     fext = allocate_complex(nfft)
     # make fft plan
-    forward_transform = plan_fft!(fext, flags=FFTW.WISDOM_ONLY)
+    forward_transform = plan_fft!(fext, flags=FFTW.ESTIMATE)
     # assign values of fext from moments 
     @inbounds begin
         for j ∈ 1:n

moment_kinetics/src/coordinates.jl

@@ -98,6 +98,8 @@ struct coordinate
     element_shift::Array{mk_float,1}
     # option used to set up element spacing
     element_spacing_option::String
+    # list of element boundaries
+    element_boundaries::Array{mk_float,1}
 end
 
 """
@@ -167,7 +169,7 @@ function define_coordinate(input, parallel_io::Bool=false; run_directory=nothing
         cell_width, igrid, ielement, imin, imax, igrid_full, input.discretization, input.fd_option, input.cheb_option,
         input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch),
         scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm,
-        local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option)
+        local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option, element_boundaries)
 
     if coord.n == 1 && occursin("v", coord.name)
         spectral = null_velocity_dimension_info()

moment_kinetics/src/em_fields.jl

@@ -14,21 +14,24 @@ using ..moment_kinetics_structs: em_fields_struct
 using ..velocity_moments: update_density!
 #using ..calculus: derivative!
 using ..derivatives: derivative_r!, derivative_z!
-
+using ..gyroaverages: gyro_operators, gyroaverage_field!
 """
 """
-function setup_em_fields(nz, nr, force_phi, drive_amplitude, drive_frequency, force_Er_zero)
+function setup_em_fields(nvperp, nz, nr, n_ion_species, force_phi, drive_amplitude, drive_frequency, force_Er_zero)
     phi = allocate_shared_float(nz,nr)
     phi0 = allocate_shared_float(nz,nr)
     Er = allocate_shared_float(nz,nr)
     Ez = allocate_shared_float(nz,nr)
-    return em_fields_struct(phi, phi0, Er, Ez, force_phi, drive_amplitude, drive_frequency, force_Er_zero)
+    gphi = allocate_shared_float(nvperp,nz,nr,n_ion_species)
+    gEr = allocate_shared_float(nvperp,nz,nr,n_ion_species)
+    gEz = allocate_shared_float(nvperp,nz,nr,n_ion_species)
+    return em_fields_struct(phi, phi0, Er, Ez, gphi, gEr, gEz, force_phi, drive_amplitude, drive_frequency, force_Er_zero)
 end
 
 """
 update_phi updates the electrostatic potential, phi
 """
-function update_phi!(fields, fvec, z, r, composition, z_spectral, r_spectral, scratch_dummy)
+function update_phi!(fields, fvec, vperp, z, r, composition, z_spectral, r_spectral, scratch_dummy, gyroavs::gyro_operators)
     n_ion_species = composition.n_ion_species
     @boundscheck size(fields.phi,1) == z.n || throw(BoundsError(fields.phi))
     @boundscheck size(fields.phi,2) == r.n || throw(BoundsError(fields.phi))
@@ -137,6 +140,20 @@ function update_phi!(fields, fvec, z, r, composition, z_spectral, r_spectral, sc
             fields.Ez[:,:] .= 0.0
         end
     end
+
+    # get gyroaveraged field arrays for distribution function advance
+    gkions = composition.gyrokinetic_ions
+    if gkions
+        gyroaverage_field!(fields.gphi,fields.phi,gyroavs,vperp,z,r,composition)
+        gyroaverage_field!(fields.gEz,fields.Ez,gyroavs,vperp,z,r,composition)
+        gyroaverage_field!(fields.gEr,fields.Er,gyroavs,vperp,z,r,composition)
+    else # use the drift-kinetic form of the fields in the kinetic equation
+        @loop_s_r_z_vperp is ir iz ivperp begin
+            fields.gphi[ivperp,iz,ir,is] = fields.phi[iz,ir]
+            fields.gEz[ivperp,iz,ir,is] = fields.Ez[iz,ir]
+            fields.gEr[ivperp,iz,ir,is] = fields.Er[iz,ir]
+        end
+    end
 
 end