Merge pull request #141 from mabarnes/Krook-collisions-MMS-port

Krook collisions manufactured solutions test port

machines/machine_setup.sh

@@ -203,10 +203,11 @@ echo
 JULIA_DIRECTORY=$JULIA_DEPOT_PATH
 echo "It can be useful or necessary to set a non-default location for the "
 echo ".julia directory. Leave this empty if the default location is OK."
-echo "Enter the current directory '.' to isolate the julia used for this "
-echo "instance of moment_kinetics - this might be useful to ensure a 'clean'"
-echo "install or to check whether some error is related to conflicting or "
-echo "corrupted dependencies or cached precompilation files, etc."
+echo "Enter a name for a subdirectory of the current directory, e.g. "
+echo "'.julia', to isolate the julia used for this instance of "
+echo "moment_kinetics - this might be useful to ensure a 'clean' install or "
+echo "to check whether some error is related to conflicting or corrupted "
+echo "dependencies or cached precompilation files, etc."
 echo "Enter location that should be used for the .julia directory [$JULIA_DIRECTORY]:"
 # Use '-e' option to get path auto-completion
 read -e -p "> "  input

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/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
 

src/initial_conditions.jl

@@ -31,7 +31,7 @@ using ..velocity_moments: create_moments_charged, create_moments_neutral, update
 using ..velocity_moments: moments_charged_substruct, moments_neutral_substruct
 using ..velocity_moments: update_neutral_density!, update_neutral_pz!, update_neutral_pr!, update_neutral_pzeta!
 using ..velocity_moments: update_neutral_uz!, update_neutral_ur!, update_neutral_uzeta!, update_neutral_qz!
-using ..velocity_moments: update_ppar!, update_upar!, update_density!, reset_moments_status!
+using ..velocity_moments: update_ppar!, update_upar!, update_density!, update_pperp!, update_vth!, reset_moments_status!
 
 using ..manufactured_solns: manufactured_solutions
 
@@ -157,8 +157,9 @@ function init_pdf_and_moments!(pdf, moments, boundary_distributions, geometry,
             init_upar!(moments.charged.upar, z, r, species.charged, n_ion_species)
             # initialise the parallel thermal speed profile
             init_vth!(moments.charged.vth, z, r, species.charged, n_ion_species)
+            # initialise pressures assuming isotropic distribution
             @. moments.charged.ppar = 0.5 * moments.charged.dens * moments.charged.vth^2
-
+            @. moments.charged.pperp = moments.charged.ppar
             if n_neutral_species > 0
                 #neutral particles
                 init_density!(moments.neutral.dens, z, r, species.neutral, n_neutral_species)
@@ -898,7 +899,7 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa,
                                               manufactured_solns_input)
     manufactured_solns_list = manufactured_solutions(manufactured_solns_input, r.L, z.L,
                                                      r.bc, z.bc, geometry, composition,
-                                                     species, r.n)
+                                                     species, r.n, vperp.n)
     dfni_func = manufactured_solns_list.dfni_func
     densi_func = manufactured_solns_list.densi_func
     dfnn_func = manufactured_solns_list.dfnn_func
@@ -924,17 +925,13 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa,
                  moments.charged.dens, moments.charged.upar, pdf.charged.norm,
                  vpa, vperp, z, r, composition, moments.evolve_density,
                  moments.evolve_upar)
+    update_pperp!(moments.charged.pperp, pdf.charged.norm, vpa, vperp, z, r, composition)
     update_qpar!(moments.charged.qpar, moments.charged.qpar_updated,
                  moments.charged.dens, moments.charged.upar,
                  moments.charged.vth, pdf.charged.norm, vpa, vperp, z, r,
                  composition, moments.evolve_density, moments.evolve_upar,
                  moments.evolve_ppar)
-    # convert from particle particle flux to parallel flow
-    begin_s_r_z_region()
-    @loop_s_r_z is ir iz begin
-        # update the thermal speed
-        moments.charged.vth[iz,ir,is] = sqrt(2.0*moments.charged.ppar[iz,ir,is]/moments.charged.dens[iz,ir,is])
-    end
+    update_vth!(moments.charged.vth, moments.charged.ppar, moments.charged.pperp, moments.charged.dens, vperp, z, r, composition)
 
     if n_neutral_species > 0
         begin_sn_r_z_region()

src/input_structs.jl

@@ -263,13 +263,7 @@ mutable struct species_composition
     phi_wall::mk_float
     # constant for testing nonzero Er
     Er_constant::mk_float
-    # constant controlling divergence at wall boundaries in MMS test
-	epsilon_offset::mk_float
-    # logical controlling whether or not dfni(vpabar,z,r) or dfni(vpa,z,r) in MMS test
-    use_vpabar_in_mms_dfni::Bool
-    # associated float controlling form of assumed potential in MMS test
-    alpha_switch::mk_float    
-	# ratio of the neutral particle mass to the ion mass
+    # ratio of the neutral particle mass to the ion mass
     mn_over_mi::mk_float
     # ratio of the electron particle mass to the ion mass
     me_over_mi::mk_float
@@ -312,6 +306,8 @@ mutable struct collisions_input
     constant_ionization_rate::Bool
     # Coulomb collision rate at the reference density and temperature
     krook_collision_frequency_prefactor::mk_float
+    # Setting to switch between different options for Krook collision operator
+    krook_collisions_option::String
 end
 
 """