MMS test bugfixes and docs

README.md

@@ -393,3 +393,55 @@ so, it skips many cases. To run more comprehensive tests, you can activate the
 To get more output on what tests were successful, an option `--verbose` (or
 `-v`) can be passed in a similar way to `--long` (if any tests fail, the output
 is printed by default).
+
+## Manufactured Solutions Tests
+In addition to the test suite in the `test/` subdirectory, the `moment_kinetics` project
+utilises the method of manufactured solutions to test more complicated models in 1D1V, 
+and 2D2V or 2D3V (for neutral particles). To run these tests we run a normal `moment_kinetics`
+simulation, making use of the manufacted solutions test TOML options. We describe how to use 
+the existing tests below. To set up `moment_kinetics` to use the manufactured solutions features,
+ take the following steps:
+* Install `moment_kinetics` using the setup instructions above ([Setup](https://github.com/mabarnes/moment_kinetics/tree/mms_bugfixes_and_docs#setup)), 
+  using the `plots_post_processing` project and make sure that the `Symbolics` package is installed, e.g., if following
+  the manual setup instructions ([Manual setup](https://mabarnes.github.io/moment_kinetics/dev/manual_setup/)), these commands would be
+    ```
+    $ julia -O3 --project
+    julia> ]
+    develop ./moment_kinetics
+    develop ./plots_post_processing/plots_post_processing
+    add Symbolics
+    ```
+    if you will run the tests with MPI, make sure that MPI is also installed at this step.
+* Select an input file representing the desired test. For example, we can pick from the list 
+  [MMS input TOML list](https://mabarnes.github.io/moment_kinetics/dev/manufactured_solution_test_examples/).
+* Run the input file using the usual command.
+    ```
+    julia> using moment_kinetics
+    julia> run_moment_kinetics("runs/your_MMS_test_input.toml")
+    ```
+* Use the post processing module to test the error norms for the simulation of interest.
+    ```
+    julia> using plots_post_processing
+    julia> analyze_and_plot_data("runs/your_MMS_test_input")
+    ```
+  This will print out a series of numbers to the terminal which represent the error norms
+  for each field and distribution function compared to the exact analytical solution, at 
+  each time step in the simulation. This error data can be computed for different resolutions.
+
+* Finally, to partially automate this last step when a resolution scan is performed, we provide
+  functions for generating plots of the error data versus resolutions in the file `plot_MMS_sequence.jl`
+  in the `plots_post_processing` project. This can be accessed by using the `run_MMS_test.jl` 
+  script from the command line
+  ```
+  $ julia -O3 --project run_MMS_test.kl
+  ```
+  or by using the underlying functions in the REPL
+  ```
+  import plots_post_processing
+  using plots_post_processing.plot_MMS_sequence
+  run_mms_test()
+  ```
+  Note that currently the lists of files used as input for the plotting functions
+  are hardcoded for the purposes of self-documenting the tests -- these lists could be made
+  input parameters to improve these scripts.
+

docs/src/index.md

@@ -18,6 +18,7 @@ Pages = ["getting_started.md",
          "manual_setup.md",
          "machine_setup_notes.md",
          "parameter_scans.md",
+         "manufactured_solution_test_examples.md",
         ]
 ```
 

docs/src/manufactured_solution_test_examples.md

@@ -0,0 +1,117 @@
+# List of Manufactured Solutions Test TOML inputs
+
+Here we list the existing manufactured solution test inputs.
+These inputs are examples only, and in most cases we only
+keep the lowest resolution examples. The user should copy
+these inputs and make a series of TOML with increasing resolutions
+to generate a series of simulations on which the numerical errors
+can be tested and compared to the expected scaling of the 
+numerical method employed.
+
+# 1D1V tests
+
+There are 1D1V tests which complement the check-in testing suite.
+The example input files in this category are as follows:
+
+* 1D1V simulation of kinetic ions (no neutrals) and numerical
+  velocity dissipation.
+```
+runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_diss.toml
+```
+* 1D1V simulation of kinetic ions (no neutrals) and a krook
+  collision operator.
+```
+runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml
+```
+
+# 1D2V tests
+
+* 1D2V simulation of kinetic ions (no neutrals) and a krook
+  collision operator.
+```
+runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml
+```
+* 1D2V simulation of a open field lines in 1D magnetic mirror
+  (no neutrals)
+```
+runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml
+runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml
+runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml
+runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml
+```
+
+# 1D2V/1D3V test (with neutrals)
+
+* A test with periodic boundary conditions in 1D, Boltzmann
+  electrons, neutrals, and ions.
+```
+runs/1D-sound-wave_cheb_nel_r_1_z_2_vpa_4_vperp_4.toml
+```
+
+* A test with ions and neutral species and wall boundary conditions, 
+  using the Boltzmann electron response to model the electron species.
+```
+runs/1D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_4_vperp_4.toml
+```
+
+* A test with neutral species and wall boundary conditions, using a simple
+  sheath model for electrons based on the Boltzmann electron response.
+```
+runs/1D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_4_vperp_4.toml
+```
+
+# 2D1V tests
+
+These tests are used to test the spatial advection in simple cases
+with wall boundary conditions.
+
+* 2D1V test of wall boundary conditions and the E x B drift.
+  Numerical dissipation in the radial domain is imposed to stabilise
+  an instability that otherwise appears at the grid scale.
+```
+runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss.toml
+```
+
+# 2D2V tests
+
+These tests are used to test the spatial advection in cases
+with wall boundary conditions or geometrical features where
+two velocity dimensions are necessary.
+
+* 2D2V simulation with periodic boundary conditions, 
+  Boltzmann electrons and ions
+```
+runs/2D-sound-wave_cheb_ion_only_nel_r_2_z_2_vpa_4_vperp_4.toml
+```
+
+* 2D2V simulation of a open field lines in 1D magnetic mirror
+  (no neutrals)
+```
+runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml
+runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml
+runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_16_vperp_16_diss.toml
+```
+
+# 2D2V/2D3V tests
+
+These tests include a two-dimensional domain, ions, and neutrals 
+(which have three velocity dimensions).
+
+* 2D2V/2D3V simulation on a domain with wall boundaries, with
+  helical geometry, Boltzmann electrons, neutrals, and ions.
+```
+runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_4_vperp_4.toml
+```
+
+* 2D2V/2D3V simulation on a periodic domain, with
+  helical geometry, Boltzmann electrons, neutrals, and ions.
+```
+runs/2D-sound-wave_cheb_nel_r_2_z_2_vpa_4_vperp_4.toml
+```
+
+* 2D2V/2D3V simulation on a periodic domain, with
+  model charge exchange and ionisation collisions,
+  helical geometry, Boltzmann electrons, neutrals, and ions.
+```
+runs/2D-sound-wave_cheb_cxiz_nel_r_2_z_2_vpa_4_vperp_4.toml
+```

moment_kinetics/ext/manufactured_solns_ext.jl

@@ -412,7 +412,7 @@ using IfElse
         # get N_e factor for boltzmann response
         if composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr == 1 
             # so 1D MMS test with 3V neutrals where ion current can be calculated prior to knowing Er
-            jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc,
+            jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc, composition,
                                                           manufactured_solns_input)
             N_e = -2.0*sqrt(pi*composition.me_over_mi)*exp(-composition.phi_wall/composition.T_e)*jpari_into_LHS_wall
         elseif composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr > 1 

runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml

@@ -30,13 +30,6 @@ z_IC_temperature_amplitude2 = 0.0
 z_IC_temperature_phase2 = 0.0
 charge_exchange_frequency = 0.0
 ionization_frequency = 0.0
-nstep = 2000
-dt = 0.0005
-nwrite = 200
-nwrite_dfns = 200
-use_semi_lagrange = false
-n_rk_stages = 4
-split_operators = false
 z_ngrid = 9
 z_nelement = 16
 z_nelement_local = 16
@@ -78,6 +71,13 @@ vzeta_L = 12.0
 vzeta_bc = "periodic"
 vzeta_discretization = "chebyshev_pseudospectral"
 
+[timestepping]
+nstep = 2000
+dt = 0.0005
+nwrite = 200
+nwrite_dfns = 200
+split_operators = false
+
 [manufactured_solns]
  use_for_advance=true
  use_for_init=true

runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml

@@ -30,13 +30,6 @@ z_IC_temperature_amplitude2 = 0.0
 z_IC_temperature_phase2 = 0.0
 charge_exchange_frequency = 0.0
 ionization_frequency = 0.0
-nstep = 2000
-dt = 0.0005
-nwrite = 200
-nwrite_dfns = 200
-use_semi_lagrange = false
-n_rk_stages = 4
-split_operators = false
 z_ngrid = 9
 z_nelement = 32
 z_nelement_local = 32
@@ -78,6 +71,13 @@ vzeta_L = 12.0
 vzeta_bc = "periodic"
 vzeta_discretization = "chebyshev_pseudospectral"
 
+[timestepping]
+nstep = 2000
+dt = 0.0005
+nwrite = 200
+nwrite_dfns = 200
+split_operators = false
+
 [manufactured_solns]
  use_for_advance=true
  use_for_init=true

runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml

@@ -30,13 +30,6 @@ z_IC_temperature_amplitude2 = 0.0
 z_IC_temperature_phase2 = 0.0
 charge_exchange_frequency = 0.0
 ionization_frequency = 0.0
-nstep = 2000
-dt = 0.0005
-nwrite = 200
-nwrite_dfns = 200
-use_semi_lagrange = false
-n_rk_stages = 4
-split_operators = false
 z_ngrid = 9
 z_nelement = 4
 z_nelement_local = 4
@@ -78,6 +71,13 @@ vzeta_L = 12.0
 vzeta_bc = "periodic"
 vzeta_discretization = "chebyshev_pseudospectral"
 
+[timestepping]
+nstep = 2000
+dt = 0.0005
+nwrite = 200
+nwrite_dfns = 200
+split_operators = false
+
 [manufactured_solns]
  use_for_advance=true
  use_for_init=true

runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml

@@ -30,13 +30,6 @@ z_IC_temperature_amplitude2 = 0.0
 z_IC_temperature_phase2 = 0.0
 charge_exchange_frequency = 0.0
 ionization_frequency = 0.0
-nstep = 2000
-dt = 0.0005
-nwrite = 200
-nwrite_dfns = 200
-use_semi_lagrange = false
-n_rk_stages = 4
-split_operators = false
 z_ngrid = 9
 z_nelement = 64
 z_nelement_local = 64
@@ -78,6 +71,13 @@ vzeta_L = 12.0
 vzeta_bc = "periodic"
 vzeta_discretization = "chebyshev_pseudospectral"
 
+[timestepping]
+nstep = 2000
+dt = 0.0005
+nwrite = 200
+nwrite_dfns = 200
+split_operators = false
+
 [manufactured_solns]
  use_for_advance=true
  use_for_init=true

runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml

@@ -30,13 +30,6 @@ z_IC_temperature_amplitude2 = 0.0
 z_IC_temperature_phase2 = 0.0
 charge_exchange_frequency = 0.0
 ionization_frequency = 0.0
-nstep = 2000
-dt = 0.0005
-nwrite = 200
-nwrite_dfns = 200
-use_semi_lagrange = false
-n_rk_stages = 4
-split_operators = false
 z_ngrid = 9
 z_nelement = 8
 z_nelement_local = 8
@@ -78,6 +71,13 @@ vzeta_L = 12.0
 vzeta_bc = "periodic"
 vzeta_discretization = "chebyshev_pseudospectral"
 
+[timestepping]
+nstep = 2000
+dt = 0.0005
+nwrite = 200
+nwrite_dfns = 200
+split_operators = false
+
 [manufactured_solns]
  use_for_advance=true
  use_for_init=true

runs/2D-sound-wave_cheb_nel_r_2_z_2_vpa_8_vperp_8.toml → runs/1D-sound-wave_cheb_nel_r_1_z_2_vpa_4_vperp_4.toml

@@ -6,9 +6,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 T_e = 1.0
-Bzed = 0.5
-Bmag = 1.0
-rhostar = 1.0
 initial_density1 = 0.5
 initial_temperature1 = 1.0
 initial_density2 = 0.5
@@ -33,36 +30,36 @@ z_ngrid = 5
 z_nelement = 2
 z_bc = "periodic"
 z_discretization = "chebyshev_pseudospectral"
-r_ngrid = 5
-r_nelement = 2
+r_ngrid = 1
+r_nelement = 1
 r_bc = "periodic"
 r_discretization = "chebyshev_pseudospectral"
 vpa_ngrid = 5
-vpa_nelement = 8
+vpa_nelement = 4
 vpa_L = 12.0
 vpa_bc = "periodic"
 vpa_discretization = "chebyshev_pseudospectral"
 
 vperp_ngrid = 5
-vperp_nelement = 8
+vperp_nelement = 4
 vperp_L = 6.0
 #vperp_discretization = "finite_difference"
 vperp_discretization = "chebyshev_pseudospectral"
 
 vz_ngrid = 5
-vz_nelement = 8
+vz_nelement = 4
 vz_L = 12.0
 vz_bc = "periodic"
 vz_discretization = "chebyshev_pseudospectral"
 
 vr_ngrid = 5
-vr_nelement = 8
+vr_nelement = 4
 vr_L = 12.0
 vr_bc = "periodic"
 vr_discretization = "chebyshev_pseudospectral"
 
 vzeta_ngrid = 5
-vzeta_nelement = 8
+vzeta_nelement = 4
 vzeta_L = 12.0
 vzeta_bc = "periodic"
 vzeta_discretization = "chebyshev_pseudospectral"
@@ -71,8 +68,12 @@ vzeta_discretization = "chebyshev_pseudospectral"
 nstep = 160
 dt = 0.002
 nwrite = 40
-n_rk_stages = 4
 split_operators = false
 
 [manufactured_solns]
 use_for_advance = true
+
+[geometry]
+option="constant-helical"
+pitch=1.0
+rhostar = 0.0