Merge pull request #225 from mabarnes/fokker-planck-vperp-bc-experiment

Fokker Planck optional features

examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml

@@ -9,7 +9,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 #force_Er_zero_at_wall = false #true
-#Er_constant = 0.0
 #epsilon_offset = 0.1
 #use_vpabar_in_mms_dfni = true
 T_e = 1.0
@@ -95,4 +94,4 @@ vperp_discretization = "gausslegendre_pseudospectral"
 [fokker_planck_collisions]
 use_fokker_planck = true
 nuii = 1.0
-frequency_option = "manual"
+frequency_option = "manual"

examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml

@@ -0,0 +1,92 @@
+# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions and collisions with fixed Maxwellian background of cold ions and electrons.
+n_ion_species = 1
+n_neutral_species = 0
+electron_physics = "boltzmann_electron_response"
+evolve_moments_density = false
+evolve_moments_parallel_flow = false
+evolve_moments_parallel_pressure = false
+evolve_moments_conservation = false
+T_e = 1.0
+T_wall = 1.0
+initial_density1 = 1.0
+initial_temperature1 = 1.0
+z_IC_option1 = "sinusoid"
+z_IC_density_amplitude1 = 0.0
+z_IC_density_phase1 = 0.0
+z_IC_upar_amplitude1 = 0.0
+z_IC_upar_phase1 = 0.0
+z_IC_temperature_amplitude1 = 0.0
+z_IC_temperature_phase1 = 0.0
+vpa_IC_option1 = "isotropic-beam"
+vpa_IC_v01 = 1.0
+vpa_IC_vth01 = 0.1
+#vpa_IC_option1 = "directed-beam"
+#vpa_IC_vpa01 = -1.5
+#vpa_IC_vperp01 = 0.0
+charge_exchange_frequency = 0.0
+ionization_frequency = 0.0
+constant_ionization_rate = false
+
+z_ngrid = 1
+z_nelement = 1
+z_nelement_local = 1
+z_bc = "wall"
+z_discretization = "chebyshev_pseudospectral"
+r_ngrid = 1
+r_nelement = 1
+r_nelement_local = 1
+r_bc = "periodic"
+r_discretization = "chebyshev_pseudospectral"
+vpa_ngrid = 5
+vpa_nelement = 32
+vpa_L = 3.0
+vpa_bc = "zero"
+vpa_discretization = "gausslegendre_pseudospectral"
+vperp_ngrid = 5
+vperp_nelement = 16
+vperp_L = 1.5
+vperp_discretization = "gausslegendre_pseudospectral"
+vperp_bc = "zero"
+# Fokker-Planck operator requires the "gausslegendre_pseudospectral
+# options for the vpa and vperp grids
+
+[fokker_planck_collisions]
+# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency
+nuii = 1.876334222e-3 #(1/nu_alphae, as computed from input diagnostic)
+Zi = 2.0
+self_collisions = true
+slowing_down_test = true
+frequency_option = "manual"
+use_fokker_planck = true
+use_conserving_corrections = true
+sd_density = 1.0
+sd_temp = 0.025 # TD/Ealpha
+sd_mi = 0.5 # mD/malpha
+sd_me = 0.000013616 # 0.25/1836.0 me/malpha
+sd_q = 1.0
+
+[ion_source]
+active=false
+source_strength=0.0
+source_T=0.005
+source_n=1.0
+r_profile="constant"
+r_width=1.0
+r_relative_minimum=0.0
+z_profile="gaussian"
+z_width=0.1
+z_relative_minimum=0.0
+source_v0 = 1.0
+#source_type="alphas"
+source_type="alphas-with-losses"
+#source_type="beam-with-losses"
+#source_vpa0 = 1.0
+#source_vperp0 = 1.0
+sink_strength = 0.0
+sink_vth=0.07071067811865475                     
+
+[timestepping]
+nstep = 50000
+dt = 2.5e-4
+nwrite = 100
+nwrite_dfns = 100

examples/fokker-planck/fokker-planck-relaxation.toml

@@ -13,7 +13,7 @@ initial_temperature1 = 1.0
 initial_density2 = 0.5
 initial_temperature2 = 1.0
 z_IC_option1 = "sinusoid"
-z_IC_density_amplitude1 = 0.001
+z_IC_density_amplitude1 = 0.0
 z_IC_density_phase1 = 0.0
 z_IC_upar_amplitude1 = 0.0
 z_IC_upar_phase1 = 0.0

examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml

@@ -0,0 +1,91 @@
+# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions and collisions with fixed Maxwellian background of cold ions and electrons.
+n_ion_species = 1
+n_neutral_species = 0
+electron_physics = "boltzmann_electron_response"
+evolve_moments_density = false
+evolve_moments_parallel_flow = false
+evolve_moments_parallel_pressure = false
+evolve_moments_conservation = false
+T_e = 1.0
+T_wall = 1.0
+initial_density1 = 1.0
+initial_temperature1 = 1.0
+z_IC_option1 = "sinusoid"
+z_IC_density_amplitude1 = 0.0
+z_IC_density_phase1 = 0.0
+z_IC_upar_amplitude1 = 0.0
+z_IC_upar_phase1 = 0.0
+z_IC_temperature_amplitude1 = 0.0
+z_IC_temperature_phase1 = 0.0
+vpa_IC_option1 = "isotropic-beam"
+vpa_IC_v01 = 1.0
+vpa_IC_vth01 = 0.1
+#vpa_IC_option1 = "directed-beam"
+#vpa_IC_vpa01 = -1.5
+#vpa_IC_vperp01 = 0.0
+charge_exchange_frequency = 0.0
+ionization_frequency = 0.0
+constant_ionization_rate = false
+
+z_ngrid = 1
+z_nelement = 1
+z_nelement_local = 1
+z_bc = "wall"
+z_discretization = "chebyshev_pseudospectral"
+r_ngrid = 1
+r_nelement = 1
+r_nelement_local = 1
+r_bc = "periodic"
+r_discretization = "chebyshev_pseudospectral"
+vpa_ngrid = 5
+vpa_nelement = 32
+vpa_L = 3.0
+vpa_bc = "zero"
+vpa_discretization = "gausslegendre_pseudospectral"
+vperp_ngrid = 5
+vperp_nelement = 16
+vperp_L = 1.5
+vperp_discretization = "gausslegendre_pseudospectral"
+vperp_bc = "zero"
+# Fokker-Planck operator requires the "gausslegendre_pseudospectral
+# options for the vpa and vperp grids
+
+[fokker_planck_collisions]
+# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency
+nuii = 1.876334222e-3 #(1/nu_alphae, as computed from input diagnostic)
+Zi = 2.0
+self_collisions = false
+slowing_down_test = true
+frequency_option = "manual"
+use_fokker_planck = true
+sd_density = 1.0
+sd_temp = 0.0025 # TD/Ealpha
+sd_mi = 0.5 # mD/malpha
+sd_me = 0.000013616 # 0.25/1836.0 me/malpha
+sd_q = 1.0
+
+[ion_source]
+active=true
+source_strength=1.0
+source_T=0.005
+source_n=1.0
+r_profile="constant"
+r_width=1.0
+r_relative_minimum=0.0
+z_profile="gaussian"
+z_width=0.1
+z_relative_minimum=0.0
+source_v0 = 1.0
+#source_type="alphas"
+source_type="alphas-with-losses"
+#source_type="beam-with-losses"
+#source_vpa0 = 1.0
+#source_vperp0 = 1.0
+sink_strength = 1.0
+sink_vth=0.07071067811865475                     
+
+[timestepping]
+nstep = 250000
+dt = 1.0e-4
+nwrite = 500
+nwrite_dfns = 500

examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml → examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml

@@ -85,7 +85,7 @@ sink_strength = 1.0
 sink_vth=0.07071067811865475                     
 
 [timestepping]
-nstep = 50000
-dt = 5.0e-4
-nwrite = 100
-nwrite_dfns = 100
+nstep = 250000
+dt = 1.0e-4
+nwrite = 500
+nwrite_dfns = 500

makie_post_processing/makie_post_processing/src/shared_utils.jl

@@ -5,15 +5,14 @@ export calculate_and_write_frequencies, get_geometry, get_composition
 using moment_kinetics.analysis: fit_delta_phi_mode
 using moment_kinetics.array_allocation: allocate_float
 using moment_kinetics.coordinates: define_coordinate
-using moment_kinetics.geo: init_magnetic_geometry
 using moment_kinetics.input_structs: boltzmann_electron_response,
                                      boltzmann_electron_response_with_simple_sheath,
                                      grid_input, geometry_input, species_composition
 using moment_kinetics.moment_kinetics_input: get_default_rhostar, setup_reference_parameters
 using moment_kinetics.type_definitions: mk_float, mk_int
 using moment_kinetics.reference_parameters: setup_reference_parameters
 using moment_kinetics.moment_kinetics_input: get_default_rhostar
-using moment_kinetics.geo: init_magnetic_geometry
+using moment_kinetics.geo: init_magnetic_geometry, setup_geometry_input
 using MPI
 
 """
@@ -90,7 +89,6 @@ function get_composition(scan_input)
     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)
     # constant to be used to control Ez divergences
     epsilon_offset = get(scan_input, "epsilon_offset", 0.001)
@@ -106,26 +104,18 @@ function get_composition(scan_input)
     # ratio of the electron particle mass to the ion particle mass
     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,
+        electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall,
         mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species))
     return composition
 
 end
 
 function get_geometry(scan_input,z,r)
     reference_params = setup_reference_parameters(scan_input)
-    # set geometry_input
-    # MRH need to get this in way that does not duplicate code
-    # MRH from moment_kinetics_input.jl
-    option = get(scan_input, "geometry_option", "constant-helical") #"1D-mirror"
-    pitch = get(scan_input, "pitch", 1.0)
-    rhostar = get(scan_input, "rhostar", get_default_rhostar(reference_params))
-    DeltaB = get(scan_input, "DeltaB", 1.0)
-    geo_in = geometry_input(rhostar,option,pitch,DeltaB)
+    reference_rhostar = get_default_rhostar(reference_params)
+    geo_in = setup_geometry_input(scan_input, reference_rhostar)
     geometry = init_magnetic_geometry(geo_in,z,r)
-
     return geometry
-
 end
 
 end # shared_utils.jl

moment_kinetics/debug_test/runtests.jl

@@ -8,7 +8,7 @@ function runtests()
         include(joinpath(@__DIR__, "fokker_planck_collisions_tests.jl"))
         include(joinpath(@__DIR__, "wall_bc_tests.jl"))
         include(joinpath(@__DIR__, "harrisonthompson.jl"))
-        include(joinpath(@__DIR__, "mms_tests.jl"))
+        #include(joinpath(@__DIR__, "mms_tests.jl"))
         include(joinpath(@__DIR__, "restart_interpolation_tests.jl"))
         include(joinpath(@__DIR__, "recycling_fraction_tests.jl"))
         include(joinpath(@__DIR__, "gyroaverage_tests.jl"))

moment_kinetics/ext/manufactured_solns_ext.jl

@@ -43,6 +43,8 @@ using IfElse
     # does not appear in a particular geometric coefficient, because 
     # that coefficient is a constant. 
     struct geometric_coefficients_sym{}
+        Er_constant::mk_float
+        Ez_constant::mk_float
         rhostar::mk_float
         Bzed::Union{mk_float,Num}
         Bzeta::Union{mk_float,Num}
@@ -62,6 +64,8 @@ using IfElse
         option = geometry_input_data.option
         rhostar = geometry_input_data.rhostar
         pitch = geometry_input_data.pitch
+        Er_constant = geometry_input_data.Er_constant
+        Ez_constant = geometry_input_data.Ez_constant
         if option == "constant-helical" || option == "default"
             bzed = pitch
             bzeta = sqrt(1 - bzed^2)
@@ -90,10 +94,20 @@ using IfElse
             end
             dBdz = expand_derivatives(Dz(Bmag))
             jacobian = 1.0
+        elseif option == "0D-Spitzer-test"
+            Bmag = 1.0
+            dBdz = geometry_input_data.dBdz_constant
+            dBdr = geometry_input_data.dBdr_constant
+            bzed = pitch
+            bzeta = sqrt(1 - bzed^2)
+            Bzed = Bmag*bzed
+            Bzeta = Bmag*bzeta
+            jacobian = 1.0
         else
             input_option_error("$option", option)
         end
-        geo_sym = geometric_coefficients_sym(rhostar,Bzed,Bzeta,Bmag,bzed,bzeta,dBdz,dBdr,jacobian)
+        geo_sym = geometric_coefficients_sym(Er_constant,Ez_constant,
+          rhostar,Bzed,Bzeta,Bmag,bzed,bzeta,dBdz,dBdr,jacobian)
         return geo_sym
     end
 
@@ -272,7 +286,7 @@ using IfElse
             upari = 0.0
         elseif z_bc == "wall"
             densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species)
-            Er, Ez, phi = electric_fields(Lr,Lz,r_bc,z_bc,composition,nr,manufactured_solns_input,species)
+            Er, Ez, phi = electric_fields(Lr,Lz,r_bc,z_bc,composition,geometry,nr,manufactured_solns_input,species)
             Bzeta = geometry.Bzeta
             Bmag = geometry.Bmag
             rhostar = geometry.rhostar
@@ -357,7 +371,7 @@ using IfElse
 
         if manufactured_solns_input.type == "default"
             # calculate the electric fields and the potential
-            Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, nr,
+            Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry, nr,
                                           manufactured_solns_input, species)
 
             # get geometric/composition data
@@ -402,7 +416,7 @@ using IfElse
         return dfni
     end
 
-    function electric_fields(Lr, Lz, r_bc, z_bc, composition, nr,
+    function electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry, nr,
                              manufactured_solns_input, species)
 
         # define derivative operators
@@ -412,7 +426,7 @@ using IfElse
         # get N_e factor for boltzmann response
         if composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr == 1 
             # so 1D MMS test with 3V neutrals where ion current can be calculated prior to knowing Er
-            jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc, composition,
+            jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc,
                                                           manufactured_solns_input)
             N_e = -2.0*sqrt(pi*composition.me_over_mi)*exp(-composition.phi_wall/composition.T_e)*jpari_into_LHS_wall
         elseif composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr > 1 
@@ -437,8 +451,8 @@ using IfElse
         # calculate the electric fields
         dense = densi # get the electron density via quasineutrality with Zi = 1
         phi = composition.T_e*log(dense/N_e) # use the adiabatic response of electrons for me/mi -> 0
-        Er = -Dr(phi)*rfac + composition.Er_constant
-        Ez = -Dz(phi)
+        Er = -Dr(phi)*rfac + geometry.Er_constant
+        Ez = -Dz(phi)      + geometry.Ez_constant
 
         Er_expanded = expand_derivatives(Er)
         Ez_expanded = expand_derivatives(Ez)
@@ -500,11 +514,13 @@ using IfElse
         return manufactured_solns_list
     end 
 
-    function manufactured_electric_fields(Lr, Lz, r_bc, z_bc, composition, nr,
+    function manufactured_electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry_input_data::geometry_input, nr,
                                           manufactured_solns_input, species)
 
+        # calculate the geometry symbolically
+        geometry = geometry_sym(geometry_input_data,Lz,Lr,nr)
         # calculate the electric fields and the potential
-        Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, nr,
+        Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry, nr,
                                       manufactured_solns_input, species)
 
         Er_func = build_function(Er, z, r, t, expression=Val{false})
@@ -603,7 +619,7 @@ using IfElse
 
         # calculate the electric fields and the potential
         Er, Ez, phi = electric_fields(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc,
-                                      composition, r_coord.n, manufactured_solns_input,
+                                      composition, geometry, r_coord.n, manufactured_solns_input,
                                       ion_species)
 
         # the adiabatic invariant (for compactness)