Geometry upgrade

docs/src/geometry.md

@@ -0,0 +1,79 @@
+Magnetic Geometry
+===============================================
+
+We take the magnetic field $\mathbf{B}$ to have the form 
+```math
+\begin{equation}
+\mathbf{B} = B_z \hat{\mathbf{z}} + B_\zeta \hat{\mathbf{\zeta}},
+\end{equation}
+```
+with $B_\zeta = B(r,z) b_\zeta$, $B_z = B(r,z) b_z$ and $b_z$ and $b_\zeta$ the direction cosines of the magnetic field vector.
+Here the basis vectors are those of cylindrical geometry $(r,z,\zeta)$, i.e.,
+$\hat{\mathbf{r}} = \nabla r $, $\hat{\mathbf{z}} = \nabla z$, 
+and $\hat{\mathbf{\zeta}} = r \nabla \zeta$. The unit vectors $\hat{\mathbf{r}}$, $\hat{\mathbf{z}}$, and $\hat{\mathbf{\zeta}}$
+form a right-handed orthonormal basis.
+
+Supported options
+===============================================
+
+To choose the type of geometry, set the value of "option" in the geometry namelist. The namelist will have the following appearance in the TOML file.
+```
+[geometry]
+option="constant-helical" # ( or "1D-mirror" )
+pitch = 1.0
+rhostar = 1.0
+DeltaB = 0.0
+```
+If `rhostar` is not set then it is computed from reference parameters.
+
+[geometry] option = "constant-helical"
+===============================================
+Here $b_\zeta = \sqrt{1 - b_z^2}$ is a constant, $b_z$ is a constant input parameter ("pitch") and $B$ is taken to be 1 with respect to the reference value $B_{\rm ref}$.
+
+[geometry] option = "1D-mirror"
+===============================================
+Here $b_\zeta = \sqrt{1 - b_z^2}$ is a constant, $b_z$ is a constant input parameter ("pitch") and $B = B(z)$ is taken to be 
+the function 
+```math
+\begin{equation}
+\frac{B(z)}{B_{\rm ref}} = 
+    1 + \Delta B \left( 2\left(\frac{2z}{L_z}\right)^2 - \left(\frac{2z}{L_z}\right)^4\right)
+\end{equation}
+```
+where $\Delta B $ is an input parameter ("DeltaB") that must satisfy $\Delta B > -1$.
+Recalling that the coordinate $z$ runs from 
+$z = -L_z/2$ to $L_z/2$,
+if $\Delta B > 0$ than the field represents a magnetic mirror which traps particles,
+whereas if $\Delta B < 0$ then the magnetic field accelerates particles
+by the mirror force as they approach the wall.
+Note that this field does not satisfy $\nabla \cdot \mathbf{B} = 0$, and is
+only used to test the implementation of the magnetic mirror terms. 2D simulations with a radial domain and$\mathbf{E}\times\mathbf{B}$ drifts are supported in the "1D-mirror"
+geometry option.
+
+Geometric coefficients
+===============================================
+Here, we write the geometric coefficients appearing in the characteristic equations
+explicitly.
+
+The $z$ component of the $\mathbf{E}\times\mathbf{B}$ drift is given by
+```math
+\begin{equation} \frac{\mathbf{E}\times\mathbf{B}}{B^2} \cdot \nabla z = \frac{E_r B_\zeta}{B^2} \nabla r \times \hat{\mathbf{\zeta}} \cdot \nabla z 
+= - J \frac{E_r B_\zeta}{B^2},
+\end{equation}
+```
+where we have defined $J = r \nabla r \times \nabla z \cdot \nabla \zeta$.
+Note that $J$ is dimensionless.
+The $r$ component of the $\mathbf{E}\times\mathbf{B}$ drift is given by
+```math
+\begin{equation} \frac{\mathbf{E}\times\mathbf{B}}{B^2} \cdot \nabla r = \frac{E_z B_\zeta}{B^2} \nabla z \times \hat{\mathbf{\zeta}} \cdot \nabla r 
+=  J \frac{E_z B_\zeta}{B^2}.
+\end{equation}
+```
+Due to the axisymmetry of the system, the differential operator
+ $\mathbf{b} \cdot \nabla (\cdot)  = b_z \partial {(\cdot)}{\partial z}$,
+ and the convective derivative
+```math
+\begin{equation}
+\frac{d B}{d t} = \frac{d z}{d t} \frac{\partial B}{ \partial z} + \frac{dr}{dt}\frac{\partial B}{\partial r}.
+\end{equation}
+```

docs/src/zz_geo.md

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

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

@@ -0,0 +1,62 @@
+# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions.
+n_ion_species = 1
+n_neutral_species = 0
+electron_physics = "boltzmann_electron_response"
+evolve_moments_density = false
+evolve_moments_parallel_flow = false
+evolve_moments_parallel_pressure = false
+evolve_moments_conservation = false
+T_e = 1.0
+T_wall = 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
+constant_ionization_rate = false
+# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency
+nuii = 1.0
+nstep = 200000
+dt = 1.0e-3
+nwrite = 50
+nwrite_dfns = 50
+use_semi_lagrange = false
+n_rk_stages = 4
+split_operators = 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"
+# Fokker-Planck operator requires the "gausslegendre_pseudospectral
+# options for the vpa and vperp grids
+

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

@@ -8,9 +8,6 @@ evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 T_e = 1.0
 T_wall = 1.0
-rhostar = 1.0
-Bzed = 1.0
-Bmag = 1.0
 initial_density1 = 0.5
 initial_temperature1 = 1.0
 initial_density2 = 0.5

examples/geometry/1D-mirror.toml

@@ -0,0 +1,84 @@
+n_ion_species = 1
+n_neutral_species = 0
+boltzmann_electron_response = true
+evolve_moments_density = false
+evolve_moments_parallel_flow = false
+evolve_moments_parallel_pressure = false
+evolve_moments_conservation = false
+T_e = 1.0
+T_wall = 1.0
+initial_density1 = 1.0
+initial_temperature1 = 1.0
+z_IC_option1 = "gaussian"
+z_IC_density_amplitude1 = 0.001
+z_IC_density_phase1 = 0.0
+z_IC_upar_amplitude1 = 1.0
+z_IC_upar_phase1 = 0.0
+z_IC_temperature_amplitude1 = 0.0
+z_IC_temperature_phase1 = 0.0
+vpa_IC_option1 = "gaussian"
+vpa_IC_density_amplitude1 = 1.0
+vpa_IC_density_phase1 = 0.0
+vpa_IC_upar_amplitude1 = 0.0
+vpa_IC_upar_phase1 = 0.0
+vpa_IC_temperature_amplitude1 = 0.0
+vpa_IC_temperature_phase1 = 0.0
+initial_density2 = 1.0
+initial_temperature2 = 1.0
+z_IC_option2 = "gaussian"
+z_IC_density_amplitude2 = 0.001
+z_IC_density_phase2 = 0.0
+z_IC_upar_amplitude2 = -1.0
+z_IC_upar_phase2 = 0.0
+z_IC_temperature_amplitude2 = 0.0
+z_IC_temperature_phase2 = 0.0
+vpa_IC_option2 = "gaussian"
+vpa_IC_density_amplitude2 = 1.0
+vpa_IC_density_phase2 = 0.0
+vpa_IC_upar_amplitude2 = 0.0
+vpa_IC_upar_phase2 = 0.0
+vpa_IC_temperature_amplitude2 = 0.0
+vpa_IC_temperature_phase2 = 0.0
+charge_exchange_frequency = 0.5
+ionization_frequency = 0.05
+constant_ionization_rate = true
+nstep = 10000
+dt = 1.0e-3
+nwrite = 50
+nwrite_dfns = 50
+use_semi_lagrange = false
+n_rk_stages = 4
+split_operators = false
+r_ngrid = 1
+r_nelement = 1
+z_ngrid = 5
+z_nelement = 16
+z_bc = "wall"
+z_discretization = "chebyshev_pseudospectral"
+vpa_ngrid = 5
+vpa_nelement = 16
+vpa_L = 6.0
+vpa_bc = "zero"
+vpa_discretization = "chebyshev_pseudospectral"
+vperp_ngrid = 5
+vperp_nelement = 8
+vperp_L = 3.0
+vperp_bc = "zero"
+vperp_discretization = "chebyshev_pseudospectral"
+vz_ngrid = 9
+vz_nelement = 64
+vz_L = 18.0
+vz_bc = "both_zero"
+vz_discretization = "chebyshev_pseudospectral"
+
+[numerical_dissipation]
+vpa_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+vperp_dissipation_coefficient = 1.0e-3 #1.0e-2 #1.0e-1
+force_minimum_pdf_value = 0.0
+
+[geometry]
+option="1D-mirror"
+DeltaB=0.5
+#option="constant-helical"
+pitch=1.0
+rhostar= 1.0

examples/numerical-dissipation/num-diss-relaxation.toml

@@ -8,9 +8,6 @@ evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 T_e = 1.0
 T_wall = 1.0
-rhostar = 1.0
-Bzed = 1.0
-Bmag = 1.0
 initial_density1 = 0.5
 initial_temperature1 = 1.0
 initial_density2 = 0.5

machines/shared/machine_setup.jl

@@ -301,7 +301,7 @@ function machine_setup_moment_kinetics(machine::String; no_force_exit::Bool=fals
     bindir = joinpath(repo_dir, "bin")
     mkpath(bindir)
     julia_executable_name = joinpath(bindir, "julia")
-    if batch_system || julia_directory == ""
+    if batch_system || (julia_directory == "" && mk_preferences["use_plots"] == "n")
         # Make a local link to the Julia binary so scripts in the repo can find it
         println("\n** Making a symlink to the julia executable at bin/julia\n")
         islink(julia_executable_name) && rm(julia_executable_name)
@@ -311,7 +311,9 @@ function machine_setup_moment_kinetics(machine::String; no_force_exit::Bool=fals
         # needing the julia.env setup
         open(julia_executable_name, "w") do io
             println(io, "#!/usr/bin/env bash")
-            println(io, "export JULIA_DEPOT_PATH=$julia_directory")
+            if julia_directory != ""
+                println(io, "export JULIA_DEPOT_PATH=$julia_directory")
+            end
             julia_path = joinpath(Sys.BINDIR, "julia")
             println(io, "$julia_path \"\$@\"")
         end

makie_post_processing/makie_post_processing/src/makie_post_processing.jl

@@ -44,7 +44,7 @@ using moment_kinetics.load_data: close_run_info, get_run_info_no_setup, get_vari
                                  neutral_dfn_variables, all_dfn_variables, ion_variables,
                                  neutral_variables, all_variables
 using moment_kinetics.initial_conditions: vpagrid_to_dzdt
-using .shared_utils: calculate_and_write_frequencies, get_geometry_and_composition
+using .shared_utils: calculate_and_write_frequencies
 using moment_kinetics.type_definitions: mk_float, mk_int
 using moment_kinetics.velocity_moments: integrate_over_vspace,
                                         integrate_over_neutral_vspace
@@ -5397,7 +5397,7 @@ function manufactured_solutions_get_field_and_field_sym(run_info, variable_name;
                             :density_neutral, :f, :f_neutral)
         manufactured_funcs =
             manufactured_solutions(run_info.manufactured_solns_input, Lr_in, run_info.z.L,
-                                   run_info.r.bc, run_info.z.bc, run_info.geometry,
+                                   run_info.r.bc, run_info.z.bc, run_info.geometry.input,
                                    run_info.composition, run_info.species, run_info.r.n,
                                    nvperp)
     end

makie_post_processing/makie_post_processing/src/shared_utils.jl

@@ -1,15 +1,19 @@
 module shared_utils
 
-export calculate_and_write_frequencies, get_geometry_and_composition
+export calculate_and_write_frequencies, get_geometry, get_composition
 
 using moment_kinetics.analysis: fit_delta_phi_mode
 using moment_kinetics.array_allocation: allocate_float
 using moment_kinetics.coordinates: define_coordinate
+using moment_kinetics.geo: init_magnetic_geometry
 using moment_kinetics.input_structs: boltzmann_electron_response,
                                      boltzmann_electron_response_with_simple_sheath,
                                      grid_input, geometry_input, species_composition
+using moment_kinetics.moment_kinetics_input: get_default_rhostar, setup_reference_parameters
 using moment_kinetics.type_definitions: mk_float, mk_int
-
+using moment_kinetics.reference_parameters: setup_reference_parameters
+using moment_kinetics.moment_kinetics_input: get_default_rhostar
+using moment_kinetics.geo: init_magnetic_geometry
 using MPI
 
 """
@@ -61,23 +65,15 @@ end
 
 """
 """
-function get_geometry_and_composition(scan_input,n_ion_species,n_neutral_species)
-    # set geometry_input
-    # MRH need to get this in way that does not duplicate code
-    # MRH from moment_kinetics_input.jl
-    Bzed = get(scan_input, "Bzed", 1.0)
-    Bmag = get(scan_input, "Bmag", 1.0)
-    bzed = Bzed/Bmag
-    bzeta = sqrt(1.0 - bzed^2.0)
-    Bzeta = Bmag*bzeta
-    rhostar = get(scan_input, "rhostar", 0.0)
-    geometry = geometry_input(Bzed,Bmag,bzed,bzeta,Bzeta,rhostar)
-
+function get_composition(scan_input)
+    reference_params = setup_reference_parameters(scan_input)
     # set composition input
     # MRH need to get this in way that does not duplicate code
     # MRH from moment_kinetics_input.jl
     electron_physics = get(scan_input, "electron_physics", boltzmann_electron_response)
 
+    n_ion_species = get(scan_input, "n_ion_species", 1)
+    n_neutral_species = get(scan_input, "n_neutral_species", 1)
     if electron_physics ∈ (boltzmann_electron_response, boltzmann_electron_response_with_simple_sheath)
         n_species = n_ion_species + n_neutral_species
     else
@@ -111,7 +107,23 @@ function get_geometry_and_composition(scan_input,n_ion_species,n_neutral_species
     composition = species_composition(n_species, n_ion_species, n_neutral_species,
         electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant,
         mn_over_mi, me_over_mi, recycling_fraction, allocate_float(n_species))
-    return geometry, composition
+    return composition
+
+end
+
+function get_geometry(scan_input,z,r)
+    reference_params = setup_reference_parameters(scan_input)
+    # set geometry_input
+    # MRH need to get this in way that does not duplicate code
+    # MRH from moment_kinetics_input.jl
+    option = get(scan_input, "geometry_option", "constant-helical") #"1D-mirror"
+    pitch = get(scan_input, "pitch", 1.0)
+    rhostar = get(scan_input, "rhostar", get_default_rhostar(reference_params))
+    DeltaB = get(scan_input, "DeltaB", 1.0)
+    geo_in = geometry_input(rhostar,option,pitch,DeltaB)
+    geometry = init_magnetic_geometry(geo_in,z,r)
+
+    return geometry
 
 end
 

moment_kinetics/debug_test/fokker_planck_collisions_inputs.jl

@@ -7,8 +7,6 @@ test_input_full_f = Dict(
      "nstep" => 3,
      "nwrite" => 2,
      "nwrite_dfns" => 2,
-     "Bmag" => 1.0,
-     "Bzed" => 1.0,
      "T_e" => 1.0,
      "T_wall" => 1.0,
      "electron_physics" => "boltzmann_electron_response",
@@ -31,7 +29,6 @@ test_input_full_f = Dict(
      "r_discretization" => "chebyshev_pseudospectral",
      "r_nelement" => 1,
      "r_ngrid" => 3,
-     "rhostar" => 1.0,
      "split_operators" => false,
      "vpa_L" => 6.0,
      "vpa_bc" => "zero",