Skip to content

Commit

Permalink
Merge pull request #166 from mabarnes/geometry-upgrade
Browse files Browse the repository at this point in the history
Geometry upgrade
  • Loading branch information
mrhardman authored Mar 15, 2024
2 parents 75d4b6e + b6b7df0 commit a5752c6
Show file tree
Hide file tree
Showing 45 changed files with 2,225 additions and 376 deletions.
79 changes: 79 additions & 0 deletions docs/src/geometry.md
Original file line number Diff line number Diff line change
@@ -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}
```
6 changes: 6 additions & 0 deletions docs/src/zz_geo.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
`geo`
=====

```@autodocs
Modules = [moment_kinetics.geo]
```
Original file line number Diff line number Diff line change
@@ -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

3 changes: 0 additions & 3 deletions examples/fokker-planck/fokker-planck-relaxation.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
84 changes: 84 additions & 0 deletions examples/geometry/1D-mirror.toml
Original file line number Diff line number Diff line change
@@ -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
3 changes: 0 additions & 3 deletions examples/numerical-dissipation/num-diss-relaxation.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
6 changes: 4 additions & 2 deletions machines/shared/machine_setup.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
42 changes: 27 additions & 15 deletions makie_post_processing/makie_post_processing/src/shared_utils.jl
Original file line number Diff line number Diff line change
@@ -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

"""
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
3 changes: 0 additions & 3 deletions moment_kinetics/debug_test/fokker_planck_collisions_inputs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand All @@ -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",
Expand Down
Loading

0 comments on commit a5752c6

Please sign in to comment.