Fokker Planck optional features

.github/workflows/parallel_test.yml

@@ -54,9 +54,9 @@ jobs:
           version: '1.10'
       - uses: julia-actions/cache@v1
       - run: |
-          MPILIBPATH=$(find /opt/homebrew/Cellar/open-mpi/ -name libmpi.dylib)
+          export MPILIBPATH=$(find /opt/homebrew/Cellar/open-mpi/ -name libmpi.dylib)
           touch Project.toml
-          julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.add(["MPI", "MPIPreferences"]); using MPIPreferences; MPIPreferences.use_system_binary(library_names="/opt/homebrew/Cellar/open-mpi/5.0.3/lib/libmpi.dylib")'
+          julia --project -O3 --check-bounds=no -e "import Pkg; Pkg.add([\"MPI\", \"MPIPreferences\"]); using MPIPreferences; MPIPreferences.use_system_binary(library_names=\"$MPILIBPATH\")"
           julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.add(["NCDatasets", "Random", "SpecialFunctions", "Test"]); Pkg.develop(path="moment_kinetics/")'
           julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.precompile()'
           # Need to use openmpi so that the following arguments work:

docs/src/moment_kinetic_equations.md

@@ -89,7 +89,7 @@ energy equation over $\tilde{v}_{\|}$ instead of $w_{\|}$,
   \tilde{p}_{\|,s}
   & = \frac{1}{\sqrt{\pi}}\int d\tilde{v}_{\|}\left(\tilde{v}_{\|}
       - \tilde{u}_{s}\right)^{2}\tilde{f}_{s}
-    = \int d\tilde{v}_{\|}\tilde{v}_{\|}^{2}\tilde{f}_{s}
+    = \frac{1}{\sqrt{\pi}}\int d\tilde{v}_{\|}\tilde{v}_{\|}^{2}\tilde{f}_{s}
       - \tilde{n}_{s}\tilde{u}_{s}^{2}\\
 %
   \tilde{q}_{\|,s}
@@ -184,7 +184,7 @@ tildes from here on)
   \frac{\partial u_{i}}{\partial t} + \frac{1}{n_{i}}\frac{\partial p_{\|,i}}{\partial z}
   + u_{i}\frac{\partial u_{i}}{\partial z} + \frac{1}{2}\frac{\partial\phi}{\partial z}
   &= -R_{in}n_{n}\left(u_{i}-u_{n}\right)
-    + R_{\mathrm{ion}}\frac{n_{i}n_{n}}{n_{s}}\left(u_{n}-u_{i}\right)
+    + R_{\mathrm{ion}}n_{i}n_{n}\left(u_{n}-u_{i}\right)
     - \frac{u_{i}}{n_{i}} \int dv_\parallel S_{i}
 \end{align}
 ```
@@ -440,16 +440,16 @@ The derivatives transform as
 \begin{align}
   \left.\frac{\partial f_{s}}{\partial t}\right|_{z,v\|}
   & \rightarrow\left.\frac{\partial f_{s}}{\partial t}\right|_{z,w\|}
-               - \frac{1}{v_{\mathrm{th},s}}\frac{\partial u_{s}}{\partial t}\left.\frac{\partial f_{s}}{\partial w_{\|,s}}\right|_{z,w\|}
-               - \frac{w_{\|,s}}{v_{\mathrm{th},s}}\frac{\partial v_{\mathrm{th},s}}{\partial t}\left.\frac{\partial f_{s}}{\partial w_{\|,s}}\right|_{z,w\|}\\
+               - \frac{1}{v_{\mathrm{th},s}}\frac{\partial u_{s}}{\partial t}\left.\frac{\partial f_{s}}{\partial w_{\|,s}}\right|_{z,t}
+               - \frac{w_{\|,s}}{v_{\mathrm{th},s}}\frac{\partial v_{\mathrm{th},s}}{\partial t}\left.\frac{\partial f_{s}}{\partial w_{\|,s}}\right|_{z,t}\\
 %
   \left.\frac{\partial f_{s}}{\partial z}\right|_{z,v\|}
-  & \rightarrow\left.\frac{\partial f_{s}}{\partial z}\right|_{z,w\|}
-               - \frac{1}{v_{\mathrm{th},s}}\frac{\partial u_{s}}{\partial z}\left.\frac{\partial f_{s}}{\partial w_{\|,s}}\right|_{z,w\|}
-               - \frac{w_{\|,s}}{v_{\mathrm{th},s}}\frac{\partial v_{\mathrm{th},s}}{\partial z}\left.\frac{\partial f_{s}}{\partial w_{\|,s}}\right|_{z,w\|}\\
+  & \rightarrow\left.\frac{\partial f_{s}}{\partial z}\right|_{t,w\|}
+               - \frac{1}{v_{\mathrm{th},s}}\frac{\partial u_{s}}{\partial z}\left.\frac{\partial f_{s}}{\partial w_{\|,s}}\right|_{z,t}
+               - \frac{w_{\|,s}}{v_{\mathrm{th},s}}\frac{\partial v_{\mathrm{th},s}}{\partial z}\left.\frac{\partial f_{s}}{\partial w_{\|,s}}\right|_{z,wt}\\
 %
   \left.\frac{\partial f_{s}}{\partial v_{\|}}\right|_{z,v\|}
-  & \rightarrow\frac{1}{v_{\mathrm{th},s}}\left.\frac{\partial f_{s}}{\partial w_{\|,s}}\right|_{z,w\|}
+  & \rightarrow\frac{1}{v_{\mathrm{th},s}}\left.\frac{\partial f_{s}}{\partial w_{\|,s}}\right|_{z,t}
 \end{align}
 ```
 

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-report-resolutions.toml

@@ -49,6 +49,7 @@ vperp_ngrid = 5
 vperp_nelement = 4
 vperp_L = 3.0
 vperp_discretization = "gausslegendre_pseudospectral"
+vperp_bc = "zero"
 # Fokker-Planck operator requires the "gausslegendre_pseudospectral
 # options for the vpa and vperp grids
 

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-no-regularity"
+# 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-no-regularity"
+# 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

@@ -46,7 +46,7 @@ vperp_ngrid = 5
 vperp_nelement = 16
 vperp_L = 1.5
 vperp_discretization = "gausslegendre_pseudospectral"
-vperp_bc = "zero"
+vperp_bc = "zero-no-regularity"
 # Fokker-Planck operator requires the "gausslegendre_pseudospectral
 # options for the vpa and vperp grids
 
@@ -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,7 @@ 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)
+    #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 +105,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, #Er_constant,
         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