Implicit pseudo-timestep for electrons

.github/workflows/debug_checks.yml

@@ -31,7 +31,7 @@ jobs:
           julia --project -O3 --check-bounds=yes -e 'using Pkg; Pkg.add(["MPI", "MPIPreferences", "PackageCompiler", "Symbolics"]); using MPIPreferences; MPIPreferences.use_jll_binary("OpenMPI_jll")'
           julia --project -O3 --check-bounds=no -e 'using MPI; MPI.install_mpiexecjl(; destdir=".")'
           julia --project -O3 --check-bounds=yes -e 'using Pkg; Pkg.develop(path="moment_kinetics/"); Pkg.precompile()'
-          julia --project -O3 --check-bounds=yes precompile.jl --debug 2
+          julia --project -O3 --check-bounds=yes precompile-with-check-bounds.jl --debug 2
 
           # Need to use openmpi so that we can use `--oversubscribe` to allow using more MPI ranks than physical cores
           ## Don't use --compiled-modules=no for now, as it currently breaks Symbolics.jl

.github/workflows/examples.yml

@@ -23,8 +23,4 @@ jobs:
           touch Project.toml
           julia -O3 --project -e 'import Pkg; Pkg.develop(path="moment_kinetics/"); Pkg.add("NCDatasets"); Pkg.precompile()'
           # Reduce nstep for each example to 10 to avoid the CI job taking too long
-          # Note we skip the example `if (occursin("ARK", get(t_input, "type", "") && Sys.isapple())`
-          # because the way we use MINPACK.jl (needed for nonlinear solvers
-          # used for implicit parts of timestep) doesn't currently work on
-          # macOS.
-          julia -O3 --project -e 'using moment_kinetics; for (root, dirs, files) in walkdir("examples") for file in files if endswith(file, ".toml") filename = joinpath(root, file); println(filename); input = moment_kinetics.moment_kinetics_input.read_input_file(filename); t_input = get(input, "timestepping", Dict{String,Any}()); if (occursin("ARK", get(t_input, "type", "")) && Sys.isapple()) continue end; t_input["nstep"] = 10; t_input["dt"] = 1.0e-12; input["timestepping"] = t_input; pop!(input, "z_nelement_local", ""); pop!(input, "r_nelement_local", ""); electron_t_input = get(input, "electron_timestepping", Dict{String,Any}()); electron_t_input["initialization_residual_value"] = 1.0e8; electron_t_input["converged_residual_value"] = 1.0e8; input["electron_timestepping"] = electron_t_input; nl_solver_input = get(input, "nonlinear_solver", Dict{String,Any}()); nl_solver_input["rtol"] = 1.0e6; nl_solver_input["atol"] = 1.0e6; input["nonlinear_solver"] = nl_solver_input; run_moment_kinetics(input) end end end'
+          julia -O3 --project -e 'using moment_kinetics; for (root, dirs, files) in walkdir("examples") for file in files if endswith(file, ".toml") filename = joinpath(root, file); println(filename); input = moment_kinetics.moment_kinetics_input.read_input_file(filename); t_input = get(input, "timestepping", Dict{String,Any}()); t_input["nstep"] = 10; t_input["dt"] = 1.0e-12; input["timestepping"] = t_input; pop!(get(input, "z", Dict{String,Any}()), "nelement_local", ""); pop!(get(input, "r", Dict{String,Any}()), "nelement_local", ""); electron_t_input = get(input, "electron_timestepping", Dict{String,Any}()); electron_t_input["initialization_residual_value"] = 1.0e8; electron_t_input["converged_residual_value"] = 1.0e8; input["electron_timestepping"] = electron_t_input; nl_solver_input = get(input, "nonlinear_solver", Dict{String,Any}()); nl_solver_input["rtol"] = 1.0e6; nl_solver_input["atol"] = 1.0e6; input["nonlinear_solver"] = nl_solver_input; run_moment_kinetics(input) end end end'

.github/workflows/parallel_test.yml

@@ -23,7 +23,7 @@ jobs:
           touch Project.toml
           julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.add(["MPI", "MPIPreferences"]); using MPIPreferences; MPIPreferences.use_jll_binary("OpenMPI_jll")'
           julia --project -O3 --check-bounds=no -e 'using MPI; MPI.install_mpiexecjl(; destdir=".")'
-          julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.add(["Random", "SpecialFunctions", "Test"]); Pkg.develop(path="moment_kinetics/")'
+          julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.add(["Random", "SpecialFunctions", "StatsBase", "Test"]); Pkg.develop(path="moment_kinetics/")'
           julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.precompile()'
           # Need to use openmpi so that we can use `--oversubscribe` to allow using more MPI ranks than physical cores
           ./mpiexecjl -np 3 --oversubscribe julia --project -O3 --check-bounds=no moment_kinetics/test/runtests.jl --debug 1
@@ -48,7 +48,7 @@ jobs:
           touch Project.toml
           julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.add(["MPI", "MPIPreferences"]); using MPIPreferences; MPIPreferences.use_jll_binary("OpenMPI_jll")'
           julia --project -O3 --check-bounds=no -e 'using MPI; MPI.install_mpiexecjl(; destdir=".")'
-          julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.add(["Random", "SpecialFunctions", "Test"]); Pkg.develop(path="moment_kinetics/")'
+          julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.add(["Random", "SpecialFunctions", "StatsBase", "Test"]); Pkg.develop(path="moment_kinetics/")'
           julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.precompile()'
           # Need to use openmpi so that we can use `--oversubscribe` to allow using more MPI ranks than physical cores
           ./mpiexecjl -np 4 --oversubscribe julia --project -O3 --check-bounds=no moment_kinetics/test/runtests.jl --debug 1

examples/kinetic-electrons/README.md

@@ -0,0 +1,10 @@
+This directory contains input files for some kinetic electron simulations that
+are known to run (and probably some other experimental input files too). Inputs
+that are expected to work:
+* Wall bc with uniform grid. First converge a Boltzmann-electron simulation to
+  steady state, then restart kinetic electron simulation from that, e.g.
+  ```julia
+  run_moment_kinetics("wall-bc_recyclefraction0.5_split3_boltzmann-coarse_tails-uniform-z-init.toml")
+  run_moment_kinetics("wall-bc_recyclefraction0.5_split3_boltzmann-coarse_tails-uniform-z.toml; restart="runs/wall-bc_recyclefraction0.5_split3_boltzmann-coarse_tails-uniform-z-init/wall-bc_recyclefraction0.5_split3_boltzmann-coarse_tails-uniform-z-init.dfns.h5")
+  run_moment_kinetics("wall-bc_recyclefraction0.5_split3_kinetic-coarse_tails-uniform-z-PareschiRusso2222.toml"; restart="runs/wall-bc_recyclefraction0.5_split3_boltzmann-coarse_tails-uniform-z/wall-bc_recyclefraction0.5_split3_boltzmann-coarse_tails-uniform-z.dfns.h5")
+  ```

examples/kinetic-electrons/periodic_split3_boltzmann.toml

@@ -2,9 +2,6 @@
 charge_exchange_frequency = 0.75
 ionization_frequency = 0.0
 
-[electron_fluid_collisions]
-nu_ei = 1000.0
-
 [evolve_moments]
 density = true
 parallel_flow = true

examples/kinetic-electrons/periodic_split3_kinetic-IMEX.toml

@@ -2,9 +2,6 @@
 charge_exchange_frequency = 0.75
 ionization_frequency = 0.0
 
-[electron_fluid_collisions]
-nu_ei = 1000.0
-
 [evolve_moments]
 density = true
 parallel_flow = true
@@ -16,25 +13,25 @@ ngrid = 1
 nelement = 1
 
 [z]
-ngrid = 17
-nelement = 16
+ngrid = 5
+nelement = 32
 #nelement_local = 16
 bc = "periodic"
-discretization = "chebyshev_pseudospectral"
+discretization = "gausslegendre_pseudospectral"
 
 [vpa]
 ngrid = 6
 nelement = 31
 L = 12.0
 bc = "zero"
-discretization = "chebyshev_pseudospectral"
+discretization = "gausslegendre_pseudospectral"
 
 [vz]
 ngrid = 6
 nelement = 31
 L = 12.0
 bc = "zero"
-discretization = "chebyshev_pseudospectral"
+discretization = "gausslegendre_pseudospectral"
 
 [composition]
 n_ion_species = 1
@@ -89,6 +86,7 @@ upar_phase = 0.0
 [timestepping]
 type = "KennedyCarpenterARK324"
 implicit_electron_advance = true
+implicit_electron_ppar = false
 implicit_ion_advance = false
 implicit_vpa_advection = false
 nstep = 1000000
@@ -116,6 +114,9 @@ type = "Fekete4(3)"
 rtol = 1.0e-6
 atol = 1.0e-14
 minimum_dt = 1.0e-10
+decrease_dt_iteration_threshold = 100
+increase_dt_iteration_threshold = 20
+cap_factor_ion_dt = 10.0
 initialization_residual_value = 2.5
 converged_residual_value = 0.1 #1.0e-3
 #debug_io = 10000

examples/kinetic-electrons/periodic_split3_kinetic-implicit-electron_ppar-loworder-PareschiRusso2222.toml

@@ -0,0 +1,152 @@
+[r]
+ngrid = 1
+nelement = 1
+
+[evolve_moments]
+parallel_pressure = true
+density = true
+moments_conservation = true
+parallel_flow = true
+
+[reactions]
+ionization_frequency = 0.0
+charge_exchange_frequency = 0.75
+
+[vz]
+ngrid = 6
+discretization = "gausslegendre_pseudospectral"
+nelement = 31
+L = 12.0
+element_spacing_option = "coarse_tails"
+bc = "zero"
+
+[ion_species_1]
+initial_temperature = 1.0
+initial_density = 1.0
+
+[krook_collisions]
+nuee0 = 1000.0
+use_krook = true
+frequency_option = "reference_parameters"
+nuei0 = 1000.0
+
+[vpa]
+ngrid = 6
+discretization = "gausslegendre_pseudospectral"
+nelement = 31
+L = 12.0
+element_spacing_option = "coarse_tails"
+bc = "zero"
+
+[z]
+ngrid = 5
+discretization = "gausslegendre_pseudospectral"
+nelement = 16
+nelement_local = 2
+bc = "periodic"
+
+[vpa_IC_ion_species_1]
+initialization_option = "gaussian"
+density_amplitude = 1.0
+temperature_amplitude = 0.0
+density_phase = 0.0
+upar_amplitude = 0.0
+temperature_phase = 0.0
+upar_phase = 0.0
+
+[z_IC_neutral_species_1]
+initialization_option = "sinusoid"
+temperature_amplitude = 0.0
+density_amplitude = 0.001
+density_phase = 0.0
+upar_amplitude = 0.0
+temperature_phase = 0.0
+upar_phase = 0.0
+
+[composition]
+T_wall = 0.1
+T_e = 1.0
+electron_physics = "kinetic_electrons"
+recycling_fraction = 0.5
+n_ion_species = 1
+n_neutral_species = 1
+
+[vz_IC_neutral_species_1]
+initialization_option = "gaussian"
+density_amplitude = 1.0
+temperature_amplitude = 0.0
+density_phase = 0.0
+upar_amplitude = 0.0
+temperature_phase = 0.0
+upar_phase = 0.0
+
+[z_IC_ion_species_1]
+initialization_option = "sinusoid"
+density_amplitude = 0.1
+temperature_amplitude = 0.1
+density_phase = 0.0
+upar_amplitude = 0.1
+temperature_phase = 1.0
+upar_phase = 0.0
+
+[neutral_species_1]
+initial_temperature = 1.0
+initial_density = 1.0
+
+[timestepping]
+type = "PareschiRusso2(2,2,2)"
+implicit_electron_advance = false
+implicit_electron_ppar = true
+implicit_ion_advance = false
+implicit_vpa_advection = false
+nstep = 500000
+dt = 2.0e-4
+#nwrite = 50
+#nwrite_dfns = 50
+nwrite = 5
+nwrite_dfns = 5
+steady_state_residual = true
+converged_residual_value = 1.0e-3
+
+[electron_timestepping]
+nstep = 5000000
+#nstep = 1
+#dt = 2.0e-8
+dt = 5.0e-5
+maximum_dt = 1.0
+nwrite = 10000
+nwrite_dfns = 100000
+#type = "SSPRK4"
+type = "Fekete4(3)"
+rtol = 1.0e-6
+atol = 1.0e-14
+minimum_dt = 1.0e-10
+decrease_dt_iteration_threshold = 100
+increase_dt_iteration_threshold = 20
+cap_factor_ion_dt = 10.0
+initialization_residual_value = 2.5
+#converged_residual_value = 0.1 #1.0e-3
+converged_residual_value = 1.0e-2
+#debug_io = 10000
+constraint_forcing_rate = 1.0e-4
+
+[nonlinear_solver]
+#nonlinear_max_iterations = 20 #100
+nonlinear_max_iterations = 1000
+rtol = 1.0e-8 #1.0e-5
+atol = 1.0e-16
+linear_restart = 5
+preconditioner_update_interval = 100 #1000
+
+[ion_numerical_dissipation]
+vpa_dissipation_coefficient = 1.0e0
+force_minimum_pdf_value = 0.0
+
+[electron_numerical_dissipation]
+vpa_dissipation_coefficient = 2.0
+#vpa_dissipation_coefficient = 2.0e2
+force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
+vz_dissipation_coefficient = 1.0e-1
+force_minimum_pdf_value = 0.0