Merge branch 'master' into wall-boundary-condition-experiment

.github/workflows/debug_checks.yml

@@ -31,7 +31,8 @@ jobs:
           sed -i -e "s/_debug_level = get_options.*/_debug_level = 2/" moment_kinetics/src/debugging.jl
 
           touch Project.toml
-          julia --project -O3 --check-bounds=yes -e 'using Pkg; Pkg.add(["MPI", "MPIPreferences", "PackageCompiler", "Symbolics"]); using MPIPreferences; MPIPreferences.use_system_binary()'
+          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
 
@@ -42,4 +43,4 @@ jobs:
           #   terrible performance when oversubscribing.
           ## Don't use --compiled-modules=no for now, as it currently breaks Symbolics.jl
           #mpiexec -np 4 --mca rmaps_base_oversubscribe 1 julia --project --check-bounds=yes --compiled-modules=no moment_kinetics/debug_test/sound_wave_tests.jl --debug 2
-          mpiexec -np 4 --mca rmaps_base_oversubscribe 1 julia --project -Jmoment_kinetics.so -O3 --check-bounds=yes moment_kinetics/debug_test/runtests.jl --debug 2
+          ./mpiexecjl -np 4 --mca rmaps_base_oversubscribe 1 julia --project -Jmoment_kinetics.so -O3 --check-bounds=yes moment_kinetics/debug_test/runtests.jl --debug 2

.github/workflows/parallel_test.yml

@@ -24,17 +24,18 @@ jobs:
       - uses: julia-actions/cache@v1
       - run: |
           touch Project.toml
-          julia --project -O3 --check-bounds=no -e 'import Pkg; Pkg.add(["MPI", "MPIPreferences"]); using MPIPreferences; MPIPreferences.use_system_binary()'
-          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.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.precompile()'
           # Need to use openmpi so that the following arguments work:
           # * `--mca rmaps_base_oversubscribe 1` allows oversubscription (more processes
           #   than physical cores).
           # * `--mca mpi_yield_when_idle 1` changes a setting to prevent excessively
           #   terrible performance when oversubscribing.
-          mpiexec -np 3 --mca rmaps_base_oversubscribe 1 julia --project -O3 --check-bounds=no moment_kinetics/test/runtests.jl --debug 1 --force-optional-dependencies
-          mpiexec -np 4 --mca rmaps_base_oversubscribe 1 julia --project -O3 --check-bounds=no moment_kinetics/test/runtests.jl --debug 1 --force-optional-dependencies
-          mpiexec -np 2 --mca rmaps_base_oversubscribe 1 julia --project -O3 --check-bounds=no moment_kinetics/test/runtests.jl --debug 1 --long --force-optional-dependencies
+          ./mpiexecjl -np 3 --mca rmaps_base_oversubscribe 1 julia --project -O3 --check-bounds=no moment_kinetics/test/runtests.jl --debug 1
+          ./mpiexecjl -np 4 --mca rmaps_base_oversubscribe 1 julia --project -O3 --check-bounds=no moment_kinetics/test/runtests.jl --debug 1
+          ./mpiexecjl -np 2 --mca rmaps_base_oversubscribe 1 julia --project -O3 --check-bounds=no moment_kinetics/test/runtests.jl --debug 1 --long
           # Note: MPI.jl's default implementation is mpich, which has a similar option
           # `--with-device=ch3:sock`, but that needs to be set when compiling mpich.
         shell: bash
@@ -56,15 +57,16 @@ jobs:
       - run: |
           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=\"$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.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.precompile()'
           # Need to use openmpi so that the following arguments work:
           # * `--mca rmaps_base_oversubscribe 1` allows oversubscription (more processes
           #   than physical cores).
           # * `--mca mpi_yield_when_idle 1` changes a setting to prevent excessively
           #   terrible performance when oversubscribing.
-          mpiexec -np 4 --mca rmaps_base_oversubscribe 1 julia --project -O3 --check-bounds=no moment_kinetics/test/runtests.jl --debug 1 --force-optional-dependencies
+          ./mpiexecjl -np 4 --mca rmaps_base_oversubscribe 1 julia --project -O3 --check-bounds=no moment_kinetics/test/runtests.jl --debug 1
           # Note: MPI.jl's default implementation is mpich, which has a similar option
           # `--with-device=ch3:sock`, but that needs to be set when compiling mpich.
         shell: bash

.github/workflows/test_scripts.yml

@@ -0,0 +1,25 @@
+# Based on example from https://github.com/julia-actions/julia-runtest
+name: Check test_scripts
+
+on: [push, pull_request, workflow_dispatch]
+
+jobs:
+  examples:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ubuntu-latest, macOS-latest]
+      fail-fast: false
+    timeout-minutes: 35
+
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: '1.10'
+      - uses: julia-actions/cache@v1
+      - name: Test test_scripts
+        run: |
+          touch Project.toml
+          julia -O3 --project -e 'import Pkg; Pkg.develop(path="moment_kinetics/"); Pkg.add(["FastGaussQuadrature", "LaTeXStrings", "LegendrePolynomials", "Measures", "MPI", "Plots", "SpecialFunctions"]); Pkg.precompile()'
+          julia -O3 --project -e 'include("test_scripts/2D_FEM_assembly_test.jl"); run_assembly_test(); include("test_scripts/chebyshev_radau_test.jl"); chebyshevradau_test(); include("test_scripts/fkpl_direct_integration_test.jl"); test_rosenbluth_potentials_direct_integration(); include("test_scripts/GaussLobattoLegendre_test.jl"); gausslegendre_test(); include("test_scripts/gyroaverage_test.jl"); gyroaverage_test()'

docs/make.jl

@@ -9,10 +9,19 @@ Pkg.instantiate()
 using Documenter
 using moment_kinetics, makie_post_processing, plots_post_processing
 
+if get(ENV, "CI", nothing) == "true"
+    # On the CI, run in strict mode to turn warnings into errors, so that we don't deploy
+    # documentation with formatting bugs.
+    strict = true
+else
+    strict = false
+end
+
 makedocs(
     sitename = "moment_kinetics",
     format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
-    modules = [moment_kinetics, makie_post_processing, plots_post_processing]
+    modules = [moment_kinetics, makie_post_processing, plots_post_processing],
+    strict = strict,
 )
 
 if get(ENV, "CI", nothing) == "true"

docs/src/developing.md

@@ -28,6 +28,44 @@ julia> run_moment_kinetics("input.toml")
 It might be convenient to add `using Revise` to your `startup.jl` file (`~/julia/config/startup.jl`) so it's always loaded.
 
 
+## Input options and defaults
+
+The input is read from a `.toml` file. It is also written to the output HDF5
+(or NetCDF) file, after all defaults are applied, both as a TOML-formatted
+String and as a tree of HDF5 variables.
+
+!!! warning
+    Neither TOML nor HDF5 have a 'null' type, so there is no convenient way to
+    store Julia's `nothing` when writing to TOML or HDF5.  Therefore `nothing`
+    should not be used as a default for any input option. If the code should
+    use `nothing` as a default for some setting, that is fine, but must be done
+    after the input is read, and not stored in the `input_dict`.
+
+
+## Array types
+
+Most arrays in `moment_kinetics` are declared using a custom array type
+[`moment_kinetics.communication.MPISharedArray`](@ref). Most of the time this
+type is just an alias for `Array`, and so it needs the same template parameters
+(see [Julia's Array
+documentation](https://docs.julialang.org/en/v1/manual/arrays/)) - the data
+type and the number of dimensions, e.g. `MPISharedArray{mk_float,3}`. Although
+these arrays use shared memory, Julia does not know about this. We use
+`MPI.Win_allocate_shared()` to allocate the shared memory, then wrap it in an
+`Array` in [`moment_kinetics.communication.allocate_shared`](@ref).
+
+The reason for using the alias, is that when the shared-memory debugging mode
+is activated, we instead create arrays using a type `DebugMPISharedArray`,
+which allows us to track some debugging information along with the array, see
+[Shared memory debugging](@ref), and make `MPISharedArray` an alias for
+`DebugMPISharedArray` instead. The reason for the alias is that if we declared
+our structs with just `Array` type, then when debugging is activated we would
+not be able to store `DebugMPISharedArray` instances in those structs, and if
+we declared the structs with `AbstractArray`, they would not be concretely
+typed, which could impact performance by creating code that is not 'type
+stable' (i.e. all concrete types are known at compile time).
+
+
 ## Parallelization
 
 The code is parallelized at the moment using MPI and shared-memory arrays. Arrays representing the pdf, moments, etc. are shared between all processes. Using shared memory means, for example, we can take derivatives along one dimension while parallelising the other for any dimension without having to communicate to re-distribute the arrays. Using shared memory instead of (in future as well as) distributed memory parallelism has the advantage that it is easier to split up the points within each element between processors, giving a finer-grained parallelism which should let the code use larger numbers of processors efficiently.

docs/src/zz_electron_fluid_equations.md

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

docs/src/zz_electron_kinetic_equation.md

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

docs/src/zz_electron_vpa_advection.md

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

docs/src/zz_electron_z_advection.md

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

docs/src/zz_maxwell_diffusion.md

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

docs/src/zz_nonlinear_solvers.md

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

docs/src/zz_species_input.md

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

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

@@ -1,35 +1,53 @@
-evolve_moments_density = false
-evolve_moments_parallel_flow = false
-evolve_moments_parallel_pressure = false
-evolve_moments_conservation = false
-charge_exchange_frequency = 0.0
-ionization_frequency = 1.0
-constant_ionization_rate = true
+[ion_source_1]
+z_profile = "constant"
+source_strength = 1.0
+source_T = 0.25
 
-z_ngrid = 5
-z_nelement = 16
-#z_nelement_local = 2
-z_bc = "wall"
-z_element_spacing_option = "sqrt"
-z_discretization = "chebyshev_pseudospectral"
-r_ngrid = 1
-r_nelement = 1
-r_nelement_local = 1
-r_bc = "periodic"
-r_discretization = "chebyshev_pseudospectral"
-vpa_ngrid = 6
-vpa_nelement = 15
-vpa_L = 6.0
-vpa_bc = "zero"
-#vpa_discretization = "chebyshev_pseudospectral"
-vpa_discretization = "gausslegendre_pseudospectral"
-vperp_ngrid = 5
-vperp_nelement = 4
-vperp_L = 3.0
-#vperp_bc = "periodic"
-#vperp_discretization = "finite_difference"
-#vperp_discretization = "chebyshev_pseudospectral"
-vperp_discretization = "gausslegendre_pseudospectral"
+[evolve_moments]
+density = false
+parallel_flow = false
+parallel_pressure = false
+moments_conservation = false
+
+[z]
+ngrid = 5
+nelement = 16
+#nelement_local = 2
+bc = "wall"
+element_spacing_option = "sqrt"
+discretization = "chebyshev_pseudospectral"
+
+[r]
+ngrid = 1
+nelement = 1
+nelement_local = 1
+bc = "periodic"
+discretization = "chebyshev_pseudospectral"
+
+[vpa]
+ngrid = 6
+nelement = 15
+L = 6.0
+bc = "zero"
+#discretization = "chebyshev_pseudospectral"
+discretization = "gausslegendre_pseudospectral"
+
+[vz]
+ngrid = 6
+nelement = 15
+L = 6.0
+bc = "zero"
+#discretization = "chebyshev_pseudospectral"
+discretization = "gausslegendre_pseudospectral"
+
+[vperp]
+ngrid = 5
+nelement = 4
+L = 3.0
+#bc = "periodic"
+#discretization = "finite_difference"
+#discretization = "chebyshev_pseudospectral"
+discretization = "gausslegendre_pseudospectral"
 
 [composition]
 n_ion_species = 1

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

@@ -1,32 +1,41 @@
 # cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions.
-evolve_moments_density = false
-evolve_moments_parallel_flow = false
-evolve_moments_parallel_pressure = false
-evolve_moments_conservation = false
+[reactions]
 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 = 8
-vpa_L = 6.0
-vpa_bc = "zero"
-vpa_discretization = "gausslegendre_pseudospectral"
-vperp_ngrid = 5
-vperp_nelement = 4
-vperp_L = 3.0
-vperp_discretization = "gausslegendre_pseudospectral"
-vperp_bc = "zero"
+[evolve_moments]
+density = false
+parallel_flow = false
+parallel_pressure = false
+moments_conservation = false
+
+[z]
+ngrid = 1
+nelement = 1
+nelement_local = 1
+bc = "wall"
+discretization = "chebyshev_pseudospectral"
+
+[r]
+ngrid = 1
+nelement = 1
+nelement_local = 1
+bc = "periodic"
+discretization = "chebyshev_pseudospectral"
+
+[vpa]
+ngrid = 5
+nelement = 8
+L = 6.0
+bc = "zero"
+discretization = "gausslegendre_pseudospectral"
+
+[vperp]
+ngrid = 5
+nelement = 4
+L = 3.0
+discretization = "gausslegendre_pseudospectral"
+bc = "zero"
 # Fokker-Planck operator requires the "gausslegendre_pseudospectral
 # options for the vpa and vperp grids