Krook 'collisions'

Project.toml

@@ -43,6 +43,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 julia = "1.7.0"

docs/src/utils.md

@@ -0,0 +1,6 @@
+`utils`
+===============
+
+```@autodocs
+Modules = [moment_kinetics.utils]
+```

runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml

@@ -0,0 +1,94 @@
+n_ion_species = 1
+n_neutral_species = 0
+electron_physics = "boltzmann_electron_response"
+#electron_physics = "boltzmann_electron_response_with_simple_sheath"
+evolve_moments_density = false
+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
+T_e = 1.0
+T_wall = 1.0
+rhostar = 1.0
+initial_density1 = 0.5
+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_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
+z_IC_option2 = "sinusoid"
+z_IC_density_amplitude2 = 0.001
+z_IC_density_phase2 = 0.0
+z_IC_upar_amplitude2 = 0.0
+z_IC_upar_phase2 = 0.0
+z_IC_temperature_amplitude2 = 0.0
+z_IC_temperature_phase2 = 0.0
+charge_exchange_frequency = 0.0
+ionization_frequency = 0.0
+krook_collisions_option = "reference_parameters"
+nuii_krook = 1.0
+nstep = 2000
+dt = 0.0005
+nwrite = 200
+nwrite_dfns = 200
+use_semi_lagrange = false
+n_rk_stages = 4
+split_operators = false
+z_ngrid = 17
+z_nelement = 16
+z_nelement_local = 16
+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 = 17
+vpa_nelement = 16
+vpa_L = 12.0
+vpa_bc = "zero"
+vpa_discretization = "chebyshev_pseudospectral"
+vperp_ngrid = 1
+vperp_nelement = 1
+vperp_L = 6.0
+vperp_bc = "periodic"
+vperp_discretization = "chebyshev_pseudospectral"
+
+vz_ngrid = 17
+vz_nelement = 4
+vz_L = 12.0
+vz_bc = "periodic"
+vz_discretization = "chebyshev_pseudospectral"
+
+vr_ngrid = 17
+vr_nelement = 4
+vr_L = 12.0
+vr_bc = "periodic"
+vr_discretization = "chebyshev_pseudospectral"
+
+vzeta_ngrid = 17
+vzeta_nelement = 4
+vzeta_L = 12.0
+vzeta_bc = "periodic"
+vzeta_discretization = "chebyshev_pseudospectral"
+
+[manufactured_solns]
+ use_for_advance=true
+ use_for_init=true
+ # constant to be used to control Ez divergence in MMS tests
+ epsilon_offset=0.1
+ # bool to control if dfni is a function of vpa or vpabar in MMS test
+ use_vpabar_in_mms_dfni=true
+ alpha_switch=1.0
+ type="default"
+[numerical_dissipation]
+vpa_dissipation_coefficient = -1.0
+z_dissipation_coefficient = -1.0
+r_dissipation_coefficient = -1.0

runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml

@@ -0,0 +1,94 @@
+n_ion_species = 1
+n_neutral_species = 0
+electron_physics = "boltzmann_electron_response"
+#electron_physics = "boltzmann_electron_response_with_simple_sheath"
+evolve_moments_density = false
+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
+T_e = 1.0
+T_wall = 1.0
+rhostar = 1.0
+initial_density1 = 0.5
+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_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
+z_IC_option2 = "sinusoid"
+z_IC_density_amplitude2 = 0.001
+z_IC_density_phase2 = 0.0
+z_IC_upar_amplitude2 = 0.0
+z_IC_upar_phase2 = 0.0
+z_IC_temperature_amplitude2 = 0.0
+z_IC_temperature_phase2 = 0.0
+charge_exchange_frequency = 0.0
+ionization_frequency = 0.0
+krook_collisions_option = "reference_parameters"
+nuii_krook = 1.0
+nstep = 2000
+dt = 0.0005
+nwrite = 200
+nwrite_dfns = 200
+use_semi_lagrange = false
+n_rk_stages = 4
+split_operators = false
+z_ngrid = 17
+z_nelement = 16
+z_nelement_local = 16
+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 = 17
+vpa_nelement = 8
+vpa_L = 12.0
+vpa_bc = "zero"
+vpa_discretization = "chebyshev_pseudospectral"
+vperp_ngrid = 17
+vperp_nelement = 8
+vperp_L = 6.0
+vperp_bc = "periodic"
+vperp_discretization = "chebyshev_pseudospectral"
+
+vz_ngrid = 17
+vz_nelement = 4
+vz_L = 12.0
+vz_bc = "periodic"
+vz_discretization = "chebyshev_pseudospectral"
+
+vr_ngrid = 17
+vr_nelement = 4
+vr_L = 12.0
+vr_bc = "periodic"
+vr_discretization = "chebyshev_pseudospectral"
+
+vzeta_ngrid = 17
+vzeta_nelement = 4
+vzeta_L = 12.0
+vzeta_bc = "periodic"
+vzeta_discretization = "chebyshev_pseudospectral"
+
+[manufactured_solns]
+ use_for_advance=true
+ use_for_init=true
+ # constant to be used to control Ez divergence in MMS tests
+ epsilon_offset=0.1
+ # bool to control if dfni is a function of vpa or vpabar in MMS test
+ use_vpabar_in_mms_dfni=true
+ alpha_switch=1.0
+ type="default"
+[numerical_dissipation]
+vpa_dissipation_coefficient = -1.0
+z_dissipation_coefficient = -1.0
+r_dissipation_coefficient = -1.0

src/constants.jl

@@ -0,0 +1,33 @@
+"""
+Some physical constants
+"""
+module constants
+
+export epsilon0, mu0
+export electron_mass
+export proton_charge, proton_mass
+export deuteron_mass
+export amu
+
+# https://physics.nist.gov/cgi-bin/cuu/Value?ep0
+const epsilon0 = 8.8541878128e-12 # F m^-1
+
+# https://physics.nist.gov/cgi-bin/cuu/Value?mu0
+const mu0 = 1.25663706212e-6 # N A^-2
+
+# https://physics.nist.gov/cgi-bin/cuu/Value?me
+const electron_mass = 9.109383701e-31 # kg
+
+# https://physics.nist.gov/cgi-bin/cuu/Value?e
+const proton_charge = 1.602176634e-19 # C
+
+# https://physics.nist.gov/cgi-bin/cuu/Value?mp
+const proton_mass = 1.67262192369e-27 # kg
+
+# https://physics.nist.gov/cgi-bin/cuu/Value?md
+const deuteron_mass = 3.3435837724e-27
+
+# https://physics.nist.gov/cgi-bin/cuu/Value?ukg
+const amu = 1.66053906660e-27
+
+end

src/file_io.jl

@@ -61,6 +61,8 @@ struct io_moments_info{Tfile, Ttime, Tphi, Tmomi, Tmomn}
     parallel_flow::Tmomi
     # handle for the charged species parallel pressure
     parallel_pressure::Tmomi
+    # handle for the charged species perpendicular pressure
+    perpendicular_pressure::Tmomi
     # handle for the charged species parallel heat flux
     parallel_heat_flux::Tmomi
     # handle for the charged species thermal speed
@@ -591,6 +593,13 @@ function define_dynamic_moment_variables!(fid, n_ion_species, n_neutral_species,
                                            parallel_io=parallel_io,
                                            description="charged species parallel pressure",
                                            units="n_ref*T_ref")
+
+        # io_pperp is the handle for the ion parallel pressure
+        io_pperp = create_dynamic_variable!(dynamic, "perpendicular_pressure", mk_float, z, r;
+                                           n_ion_species=n_ion_species,
+                                           parallel_io=parallel_io,
+                                           description="charged species perpendicular pressure",
+                                           units="n_ref*T_ref")
 
         # io_qpar is the handle for the ion parallel heat flux
         io_qpar = create_dynamic_variable!(dynamic, "parallel_heat_flux", mk_float, z, r;
@@ -642,7 +651,7 @@ function define_dynamic_moment_variables!(fid, n_ion_species, n_neutral_species,
             units="c_ref")
 
         return io_moments_info(fid, io_time, io_phi, io_Er, io_Ez, io_density, io_upar,
-                               io_ppar, io_qpar, io_vth, io_density_neutral, io_uz_neutral,
+                               io_ppar, io_pperp, io_qpar, io_vth, io_density_neutral, io_uz_neutral,
                                io_pz_neutral, io_qz_neutral, io_thermal_speed_neutral,
                                parallel_io)
     end
@@ -780,7 +789,8 @@ function reopen_moments_io(file_info)
         end
         return io_moments_info(fid, getvar("time"), getvar("phi"), getvar("Er"),
                                getvar("Ez"), getvar("density"), getvar("parallel_flow"),
-                               getvar("parallel_pressure"), getvar("parallel_heat_flux"),
+                               getvar("parallel_pressure"), getvar("perpendicular_pressure"),
+                               getvar("parallel_heat_flux"),
                                getvar("thermal_speed"), getvar("density_neutral"),
                                getvar("uz_neutral"), getvar("pz_neutral"),
                                getvar("qz_neutral"), getvar("thermal_speed_neutral"),
@@ -865,6 +875,7 @@ function reopen_dfns_io(file_info)
         io_moments = io_moments_info(fid, getvar("time"), getvar("phi"), getvar("Er"),
                                      getvar("Ez"), getvar("density"),
                                      getvar("parallel_flow"), getvar("parallel_pressure"),
+                                     getvar("perpendicular_pressure"),
                                      getvar("parallel_heat_flux"),
                                      getvar("thermal_speed"), getvar("density_neutral"),
                                      getvar("uz_neutral"), getvar("pz_neutral"),
@@ -920,6 +931,8 @@ function write_moments_data_to_binary(moments, fields, t, n_ion_species,
                               n_ion_species)
         append_to_dynamic_var(io_moments.parallel_pressure, moments.charged.ppar, t_idx,
                               z, r, n_ion_species)
+        append_to_dynamic_var(io_moments.perpendicular_pressure, moments.charged.pperp, t_idx,
+                              z, r, n_ion_species)
         append_to_dynamic_var(io_moments.parallel_heat_flux, moments.charged.qpar, t_idx,
                               z, r, n_ion_species)
         append_to_dynamic_var(io_moments.thermal_speed, moments.charged.vth, t_idx, z, r,
@@ -1008,6 +1021,8 @@ end
                                   t_idx, z, r, n_ion_species)
             append_to_dynamic_var(io_moments.parallel_pressure, moments.charged.ppar.data,
                                   t_idx, z, r, n_ion_species)
+            append_to_dynamic_var(io_moments.perpendicular_pressure.data, moments.charged.pperp,
+                                  t_idx, z, r, n_ion_species)
             append_to_dynamic_var(io_moments.parallel_heat_flux,
                                   moments.charged.qpar.data, t_idx, z, r, n_ion_species)
             append_to_dynamic_var(io_moments.thermal_speed, moments.charged.vth.data,
@@ -1271,7 +1286,7 @@ all the arrays have the same length, with an entry for each call to `debug_dump(
 function debug_dump end
 function debug_dump(vz::coordinate, vr::coordinate, vzeta::coordinate, vpa::coordinate,
                     vperp::coordinate, z::coordinate, r::coordinate, t::mk_float;
-                    ff=nothing, dens=nothing, upar=nothing, ppar=nothing, qpar=nothing,
+                    ff=nothing, dens=nothing, upar=nothing, ppar=nothing, pperp=nothing, qpar=nothing,
                     vth=nothing,
                     ff_neutral=nothing, dens_neutral=nothing, uz_neutral=nothing,
                     #ur_neutral=nothing, uzeta_neutral=nothing,
@@ -1370,6 +1385,11 @@ function debug_dump(vz::coordinate, vr::coordinate, vzeta::coordinate, vpa::coor
         else
             debug_output_file.moments.parallel_pressure[:,:,:,debug_output_counter[]] = ppar
         end
+        if pperp === nothing
+            debug_output_file.moments.perpendicular_pressure[:,:,:,debug_output_counter[]] = 0.0
+        else
+            debug_output_file.moments.perpendicular_pressure[:,:,:,debug_output_counter[]] = pperp
+        end
         if qpar === nothing
             debug_output_file.moments.parallel_heat_flux[:,:,:,debug_output_counter[]] = 0.0
         else
@@ -1446,13 +1466,15 @@ function debug_dump(fvec::Union{scratch_pdf,Nothing},
         density = nothing
         upar = nothing
         ppar = nothing
+        pperp = nothing
         pdf_neutral = nothing
         density_neutral = nothing
     else
         pdf = fvec.pdf
         density = fvec.density
         upar = fvec.upar
         ppar = fvec.ppar
+        pperp = fvec.pperp
         pdf_neutral = fvec.pdf_neutral
         density_neutral = fvec.density_neutral
     end
@@ -1466,7 +1488,7 @@ function debug_dump(fvec::Union{scratch_pdf,Nothing},
         Ez = fields.Ez
     end
     return debug_dump(vz, vr, vzeta, vpa, vperp, z, r, t; ff=pdf, dens=density, upar=upar,
-                      ppar=ppar, ff_neutral=pdf_neutral, dens_neutral=density_neutral,
+                      ppar=ppar, pperp=pperp, ff_neutral=pdf_neutral, dens_neutral=density_neutral,
                       phi=phi, Er=Er, Ez=Ez, t, istage=istage, label=label)
 end