diff --git a/docs/src/geometry.md b/docs/src/geometry.md new file mode 100644 index 000000000..9b956e189 --- /dev/null +++ b/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} +``` diff --git a/docs/src/zz_geo.md b/docs/src/zz_geo.md new file mode 100644 index 000000000..c2608f8c5 --- /dev/null +++ b/docs/src/zz_geo.md @@ -0,0 +1,6 @@ +`geo` +===== + +```@autodocs +Modules = [moment_kinetics.geo] +``` diff --git a/examples/fokker-planck/fokker-planck-relaxation-report-resolutions.toml b/examples/fokker-planck/fokker-planck-relaxation-report-resolutions.toml new file mode 100644 index 000000000..3eb069342 --- /dev/null +++ b/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 + diff --git a/examples/fokker-planck/fokker-planck-relaxation.toml b/examples/fokker-planck/fokker-planck-relaxation.toml index b12454e6f..d3284322c 100644 --- a/examples/fokker-planck/fokker-planck-relaxation.toml +++ b/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 diff --git a/examples/geometry/1D-mirror.toml b/examples/geometry/1D-mirror.toml new file mode 100644 index 000000000..66c36fd5b --- /dev/null +++ b/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 diff --git a/examples/numerical-dissipation/num-diss-relaxation.toml b/examples/numerical-dissipation/num-diss-relaxation.toml index 29965e248..abd22e2da 100644 --- a/examples/numerical-dissipation/num-diss-relaxation.toml +++ b/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 diff --git a/machines/shared/machine_setup.jl b/machines/shared/machine_setup.jl index 75104a56a..51bb48ee9 100644 --- a/machines/shared/machine_setup.jl +++ b/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 diff --git a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl index ac5595927..4d79887a0 100644 --- a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl +++ b/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 diff --git a/makie_post_processing/makie_post_processing/src/shared_utils.jl b/makie_post_processing/makie_post_processing/src/shared_utils.jl index 001082a67..4a7b32b7b 100644 --- a/makie_post_processing/makie_post_processing/src/shared_utils.jl +++ b/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 diff --git a/moment_kinetics/debug_test/fokker_planck_collisions_inputs.jl b/moment_kinetics/debug_test/fokker_planck_collisions_inputs.jl index 5f6a8d559..199712b4c 100644 --- a/moment_kinetics/debug_test/fokker_planck_collisions_inputs.jl +++ b/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", diff --git a/moment_kinetics/debug_test/mms_inputs.jl b/moment_kinetics/debug_test/mms_inputs.jl index 328e7aa3c..5378c6602 100644 --- a/moment_kinetics/debug_test/mms_inputs.jl +++ b/moment_kinetics/debug_test/mms_inputs.jl @@ -13,9 +13,6 @@ test_input = Dict( "evolve_moments_conservation" => false, "T_e" => 1.0, "T_wall" => 1.0, - "Bzed" => 0.5, - "Bmag" => 1.0, - "rhostar" => 1.0, "initial_density1" => 0.5, "initial_temperature1" => 1.0, "initial_density2" => 0.5, @@ -74,6 +71,7 @@ test_input = Dict( "vzeta_L" => 12.0, "vzeta_bc" => "none", "vzeta_discretization" => "chebyshev_pseudospectral", + "geometry" => Dict{String,Any}("rhostar" => 1.0,), ) test_input_list = [ diff --git a/moment_kinetics/ext/manufactured_solns_ext.jl b/moment_kinetics/ext/manufactured_solns_ext.jl index 6e9318deb..2bad8dcf0 100644 --- a/moment_kinetics/ext/manufactured_solns_ext.jl +++ b/moment_kinetics/ext/manufactured_solns_ext.jl @@ -7,11 +7,12 @@ module manufactured_solns_ext using moment_kinetics.input_structs +using moment_kinetics.input_structs: geometry_input using moment_kinetics.looping -using moment_kinetics.type_definitions: mk_int +using moment_kinetics.type_definitions: mk_int, mk_float import moment_kinetics.manufactured_solns: manufactured_solutions, manufactured_sources_setup, - manufactured_electric_fields + manufactured_electric_fields, manufactured_geometry using Symbolics using IfElse @@ -34,6 +35,67 @@ using IfElse pconst = 15.0/8.0 fluxconst = 1.0 end + + # struct of symbolic functions for geometric coefficients + # Note that we restrict the types of the variables in the struct + # to be either a float or a Symbolics Num type. The Union appears + # to be required to permit geometry options where a symbolic variable + # does not appear in a particular geometric coefficient, because + # that coefficient is a constant. + struct geometric_coefficients_sym{} + rhostar::mk_float + Bzed::Union{mk_float,Num} + Bzeta::Union{mk_float,Num} + Bmag::Union{mk_float,Num} + bzed::Union{mk_float,Num} + bzeta::Union{mk_float,Num} + dBdz::Union{mk_float,Num} + dBdr::Union{mk_float,Num} + jacobian::Union{mk_float,Num} + end + + function geometry_sym(geometry_input_data::geometry_input,Lz,Lr,nr) + # define derivative operators + Dr = Differential(r) + Dz = Differential(z) + # compute symbolic geometry functions + option = geometry_input_data.option + rhostar = geometry_input_data.rhostar + pitch = geometry_input_data.pitch + if option == "constant-helical" || option == "default" + bzed = pitch + bzeta = sqrt(1 - bzed^2) + Bmag = 1.0 + Bzed = Bmag*bzed + Bzeta = Bmag*bzeta + dBdr = 0.0 + dBdz = 0.0 + jacobian = 1.0 + elseif option == "1D-mirror" + DeltaB = geometry_input_data.DeltaB + bzed = pitch + bzeta = sqrt(1 - bzed^2) + # B(z)/Bref = 1 + DeltaB*( 2(2z/L)^2 - (2z/L)^4) + # chosen so that + # B(z)/Bref = 1 + DeltaB at 2z/L = +- 1 + # d B(z)d z = 0 at 2z/L = +- 1 + zfac = 2.0*z/Lz + Bmag = 1.0 + DeltaB*( 2.0*zfac^2 - zfac^4) + Bzed = Bmag*bzed + Bzeta = Bmag*bzeta + if nr > 1 + dBdr = expand_derivatives(Dr(Bmag)) + else + dBdr = 0.0 + end + dBdz = expand_derivatives(Dz(Bmag)) + jacobian = 1.0 + else + input_option_error("$option", option) + end + geo_sym = geometric_coefficients_sym(rhostar,Bzed,Bzeta,Bmag,bzed,bzeta,dBdz,dBdr,jacobian) + return geo_sym + end #standard functions for building densities function nplus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha) @@ -211,13 +273,17 @@ using IfElse elseif z_bc == "wall" densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) Er, Ez, phi = electric_fields(Lr,Lz,r_bc,z_bc,composition,nr,manufactured_solns_input,species) + Bzeta = geometry.Bzeta + Bmag = geometry.Bmag rhostar = geometry.rhostar + jacobian = geometry.jacobian + ExBgeofac = 0.5*rhostar*Bzeta*jacobian/Bmag^2 bzed = geometry.bzed epsilon = manufactured_solns_input.epsilon_offset alpha = manufactured_solns_input.alpha_switch upari = ( (fluxconst/(sqrt(pi)*densi))*((z/Lz + 0.5)*nplus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha) - (0.5 - z/Lz)*nminus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha)) - + alpha*(rhostar/(2.0*bzed))*Er ) + + alpha*(ExBgeofac/bzed)*Er ) end return upari end @@ -296,14 +362,17 @@ using IfElse # get geometric/composition data Bzed = geometry.Bzed + Bzeta = geometry.Bzeta Bmag = geometry.Bmag rhostar = geometry.rhostar + jacobian = geometry.jacobian + ExBgeofac = 0.5*rhostar*Bzeta*jacobian/Bmag^2 epsilon = manufactured_solns_input.epsilon_offset alpha = manufactured_solns_input.alpha_switch if z_bc == "periodic" dfni = densi * exp( - vpa^2 - vperp^2) elseif z_bc == "wall" - vpabar = vpa - alpha*(rhostar/2.0)*(Bmag/Bzed)*Er # for alpha = 1.0, effective velocity in z direction * (Bmag/Bzed) + vpabar = vpa - alpha*ExBgeofac*(Bmag/Bzed)*Er # for alpha = 1.0, effective velocity in z direction * (Bmag/Bzed) Hplus = 0.5*(sign(vpabar) + 1.0) Hminus = 0.5*(sign(-vpabar) + 1.0) ffa = exp(- vperp^2) @@ -378,7 +447,9 @@ using IfElse end function manufactured_solutions(manufactured_solns_input, Lr, Lz, r_bc, z_bc, - geometry, composition, species, nr, nvperp) + geometry_input_data::geometry_input, composition, species, nr, nvperp) + # calculate the geometry symbolically + geometry = geometry_sym(geometry_input_data,Lz,Lr,nr) charged_species = species.charged[1] if composition.n_neutral_species > 0 neutral_species = species.neutral[1] @@ -445,10 +516,28 @@ using IfElse return manufactured_E_fields end + function manufactured_geometry(geometry_input_data::geometry_input,Lz,Lr,nr) + + # calculate the geometry symbolically + geosym = geometry_sym(geometry_input_data,Lz,Lr,nr) + Bmag = geosym.Bmag + bzed = geosym.bzed + dBdz = geosym.dBdz + Bmag_func = build_function(Bmag, z, r, expression=Val{false}) + bzed_func = build_function(bzed, z, r, expression=Val{false}) + dBdz_func = build_function(dBdz, z, r, expression=Val{false}) + + manufactured_geometry = (Bmag_func = Bmag_func, + bzed_func = bzed_func, + dBdz_func = dBdz_func) + return manufactured_geometry + end + function manufactured_sources_setup(manufactured_solns_input, r_coord, z_coord, vperp_coord, - vpa_coord, vzeta_coord, vr_coord, vz_coord, composition, geometry, collisions, + vpa_coord, vzeta_coord, vr_coord, vz_coord, composition, + geometry_input_data::geometry_input, collisions, num_diss_params, species) - + geometry = geometry_sym(geometry_input_data,z_coord.L,r_coord.L,r_coord.n) charged_species = species.charged[1] if composition.n_neutral_species > 0 neutral_species = species.neutral[1] @@ -491,8 +580,13 @@ using IfElse # get geometric/composition data Bzed = geometry.Bzed + Bzeta = geometry.Bzeta Bmag = geometry.Bmag + dBdz = geometry.dBdz + dBdr = geometry.dBdr rhostar = geometry.rhostar + jacobian = geometry.jacobian + ExBgeofac = 0.5*rhostar*Bzeta*jacobian/Bmag^2 #exceptions for cases with missing terms if composition.n_neutral_species > 0 cx_frequency = collisions.charge_exchange @@ -512,9 +606,21 @@ using IfElse composition, r_coord.n, manufactured_solns_input, charged_species) + # the adiabatic invariant (for compactness) + mu = 0.5*(vperp^2)/Bmag + # the ion characteristic velocities + dzdt = vpa * (Bzed/Bmag) - ExBgeofac*Er + drdt = ExBgeofac*Ez*rfac + dvpadt = 0.5*(Bzed/Bmag)*Ez - mu*(Bzed/Bmag)*dBdz + dvperpdt = (0.5*vperp/Bmag)*(dzdt*dBdz + drdt*dBdr) # the ion source to maintain the manufactured solution - Si = ( Dt(dfni) + ( vpa * (Bzed/Bmag) - 0.5*rhostar*Er ) * Dz(dfni) + ( 0.5*rhostar*Ez*rfac ) * Dr(dfni) + ( 0.5*Ez*Bzed/Bmag ) * Dvpa(dfni) - + cx_frequency*( densn*dfni - densi*gav_dfnn ) - ionization_frequency*dense*gav_dfnn) + Si = ( Dt(dfni) + + dzdt * Dz(dfni) + + drdt * Dr(dfni) + + dvpadt * Dvpa(dfni) + + dvperpdt * Dvperp(dfni) + + cx_frequency*( densn*dfni - densi*gav_dfnn ) + - ionization_frequency*dense*gav_dfnn ) nu_krook = collisions.krook_collision_frequency_prefactor if nu_krook > 0.0 Ti_over_Tref = vthi^2 diff --git a/moment_kinetics/src/advection.jl b/moment_kinetics/src/advection.jl index aad3a0df6..3a78d500c 100644 --- a/moment_kinetics/src/advection.jl +++ b/moment_kinetics/src/advection.jl @@ -64,6 +64,13 @@ function setup_advection_per_species(coords...) adv_fac = allocate_shared_float([coord.n for coord in coords]...) # create array for storing the speed along this coordinate speed = allocate_shared_float([coord.n for coord in coords]...) + # initialise speed to zero so that it can be used safely without + # introducing NaNs (if left uninitialised) when coordinate speeds + # are used but the coordinate has only a single point + # (e.g. dr/dt in dvperp/dt = (vperp/2B)dB/dt, see vperp_advection.jl) + @serial_region begin + @. speed = 0.0 + end # return advection_info struct containing necessary arrays return advection_info(rhs, df, speed, adv_fac) end diff --git a/moment_kinetics/src/chebyshev.jl b/moment_kinetics/src/chebyshev.jl index 9bac36633..4e7c25eff 100644 --- a/moment_kinetics/src/chebyshev.jl +++ b/moment_kinetics/src/chebyshev.jl @@ -794,11 +794,12 @@ function chebyshev_radau_forward_transform!(chebyf, fext, ff, transform, n) chebyshev_spectral_derivative!(cheby_df, cheby_f) # form the derivative at x = - 1 using that T_n(-1) = (-1)^n # and converting the normalisation factors to undo the normalisation in the FFT - # df = d0 + sum_n=1 (-1)^n d_n/2 with d_n the coeffs + # df = d0 + sum_n=1 (-1)^n d_n with d_n the coeffs # of the Cheb derivative in the Fourier representation + # df = sum_n=0,N-1 d_n T_n(x) df = cheby_df[1] for i in 2:coord.ngrid - df += ((-1)^(i-1))*0.5*cheby_df[i] + df += ((-1)^(i-1))*cheby_df[i] end return df end diff --git a/moment_kinetics/src/force_balance.jl b/moment_kinetics/src/force_balance.jl index 413265260..562f59972 100644 --- a/moment_kinetics/src/force_balance.jl +++ b/moment_kinetics/src/force_balance.jl @@ -26,7 +26,7 @@ function force_balance!(pflx, density_out, fvec, moments, fields, collisions, dt dt*(moments.charged.dppar_dz[iz,ir,is] + upar[iz,ir,is]*upar[iz,ir,is]*moments.charged.ddens_dz_upwind[iz,ir,is] + 2.0*density[iz,ir,is]*upar[iz,ir,is]*moments.charged.dupar_dz_upwind[iz,ir,is] - - 0.5*geometry.bzed*fields.Ez[iz,ir]*density[iz,ir,is]) + 0.5*geometry.bzed[iz,ir]*fields.Ez[iz,ir]*density[iz,ir,is]) end if ion_source_settings.active && false diff --git a/moment_kinetics/src/geo.jl b/moment_kinetics/src/geo.jl new file mode 100644 index 000000000..9ad6d4885 --- /dev/null +++ b/moment_kinetics/src/geo.jl @@ -0,0 +1,200 @@ +""" +module for including axisymmetric geometry +in coordinates (z,r), with z the vertical +coordinate and r the radial coordinate +""" +module geo + +export init_magnetic_geometry +export setup_geometry_input + +using ..input_structs: geometry_input +using ..file_io: input_option_error +using ..array_allocation: allocate_float +using ..type_definitions: mk_float, mk_int + + +""" +struct containing the geometric data necessary for +non-trivial axisymmetric geometries, to be passed +around the inside of the code, replacing the +`geometry_input` struct from input_structs.jl + +The arrays of 2 dimensions are functions of (z,r) +""" +struct geometric_coefficients +# for now include the reference parameters in `geometry_input` +input::geometry_input +# also include a copy of rhostar for ease of use +rhostar::mk_float +# the spatially varying coefficients +# Bz/Bref +Bzed::Array{mk_float,2} +# Bzeta/Bref +Bzeta::Array{mk_float,2} +# Btot/Bref +Bmag::Array{mk_float,2} +# bz -- unit vector component in z direction +bzed::Array{mk_float,2} +# bz -- unit vector component in zeta direction +bzeta::Array{mk_float,2} + + +# now the new coefficients + +# d Bmag d z +dBdz::Array{mk_float,2} +# d Bmag d r +dBdr::Array{mk_float,2} +# jacobian = r grad r x grad z . grad zeta +jacobian::Array{mk_float,2} + +# magnetic drift physics coefficients +# cvdriftr = (b/B) x (b.grad b) . grad r +cvdriftr::Array{mk_float,2} +# cvdriftz = (b/B) x (b.grad b) . grad z +cvdriftz::Array{mk_float,2} +# gbdriftr = (b/B^2) x grad B . grad r +gbdriftr::Array{mk_float,2} +# gbdriftz = (b/B^2) x grad B . grad z +gbdriftz::Array{mk_float,2} +end + +""" +function to read the geometry input data from the TOML file + +the TOML namelist should be structured like + +[geometry] +pitch = 1.0 +rhostar = 1.0 +DeltaB = 0.0 +option = "" + +""" +function setup_geometry_input(toml_input::Dict, reference_rhostar) + input_section = get(toml_input, "geometry", Dict{String,Any}()) + if !("rhostar" ∈ keys(input_section)) + # Set default rhostar with reference value + input_section["rhostar"] = get(input_section, "rhostar", reference_rhostar) + end + input = Dict(Symbol(k)=>v for (k,v) in input_section) + #println(input) + return geometry_input(; input...) +end + +""" +function to initialise the geometry coefficients +input_data -- geometry_input type +z -- coordinate type +r -- coordinate type +""" +function init_magnetic_geometry(geometry_input_data::geometry_input,z,r) + nz = z.n + nr = r.n + Bzed = allocate_float(nz,nr) + Bzeta = allocate_float(nz,nr) + Bmag = allocate_float(nz,nr) + bzed = allocate_float(nz,nr) + bzeta = allocate_float(nz,nr) + dBdr = allocate_float(nz,nr) + dBdz = allocate_float(nz,nr) + jacobian = allocate_float(nz,nr) + cvdriftr = allocate_float(nz,nr) + cvdriftz = allocate_float(nz,nr) + gbdriftr = allocate_float(nz,nr) + gbdriftz = allocate_float(nz,nr) + + option = geometry_input_data.option + rhostar = geometry_input_data.rhostar + if option == "constant-helical" || option == "default" + # \vec{B} = B ( bz \hat{z} + bzeta \hat{zeta} ) + # with B a constant and \hat{z} x \hat{r} . \hat{zeta} = 1 + pitch = geometry_input_data.pitch + for ir in 1:nr + for iz in 1:nz + bzed[iz,ir] = pitch + bzeta[iz,ir] = sqrt(1 - bzed[iz,ir]^2) + Bmag[iz,ir] = 1.0 + Bzed[iz,ir] = Bmag[iz,ir]*bzed[iz,ir] + Bzeta[iz,ir] = Bmag[iz,ir]*bzeta[iz,ir] + dBdr[iz,ir] = 0.0 + dBdz[iz,ir] = 0.0 + jacobian[iz,ir] = 1.0 + + cvdriftr[iz,ir] = 0.0 + cvdriftz[iz,ir] = 0.0 + gbdriftr[iz,ir] = 0.0 + gbdriftz[iz,ir] = 0.0 + end + end + elseif option == "1D-mirror" + # a 1D configuration for testing mirror and vperp physics + # with \vec{B} = B(z) bz \hat{z} and + # with B = B(z) a specified function + #if nr > 1 + # input_option_error("$option: You have specified nr > 1 -> set nr = 1", option) + #end + DeltaB = geometry_input_data.DeltaB + if DeltaB < -0.99999999 + input_option_error("$option: You have specified DeltaB < -1 -> set DeltaB > -1", option) + end + pitch = geometry_input_data.pitch + for ir in 1:nr + for iz in 1:nz + bzed[iz,ir] = pitch + bzeta[iz,ir] = sqrt(1 - bzed[iz,ir]^2) + # B(z)/Bref = 1 + DeltaB*( 2(2z/L)^2 - (2z/L)^4) + # chosen so that + # B(z)/Bref = 1 + DeltaB at 2z/L = +- 1 + # d B(z)d z = 0 at 2z/L = +- 1 + zfac = 2.0*z.grid[iz]/z.L + Bmag[iz,ir] = 1.0 + DeltaB*( 2.0*zfac^2 - zfac^4) + Bzed[iz,ir] = Bmag[iz,ir]*bzed[iz,ir] + Bzeta[iz,ir] = Bmag[iz,ir]*bzeta[iz,ir] + dBdr[iz,ir] = 0.0 + dBdz[iz,ir] = (2.0/z.L)*4.0*DeltaB*zfac*(1.0 - zfac^2) + jacobian[iz,ir] = 1.0 + + cvdriftr[iz,ir] = 0.0 + cvdriftz[iz,ir] = 0.0 + gbdriftr[iz,ir] = 0.0 + gbdriftz[iz,ir] = 0.0 + end + end + elseif option == "low-beta-helix" + # a 2D configuration for testing magnetic drift physics + # with \vec{B} = (B0/r) \hat{zeta} + Bz \hat{z} + # with B0 and Bz constants + pitch = geometry_input_data.pitch + B0 = 1.0 # chose reference field strength to be Bzeta at r = 1 + Bz = pitch*B0 # pitch determines ratio of Bz/B0 at r = 1 + for ir in 1:nr + rr = r.grid[ir] + for iz in 1:nz + Bmag[iz,ir] = sqrt( (B0/rr)^2 + Bz^2 ) + bzed[iz,ir] = Bz/Bmag[iz,ir] + bzeta[iz,ir] = B0/(rr*Bmag[iz,ir]) + Bzed[iz,ir] = bzed[iz,ir]*Bmag[iz,ir] + Bzeta[iz,ir] = bzeta[iz,ir]*Bmag[iz,ir] + dBdz[iz,ir] = 0.0 + dBdr[iz,ir] = -(Bmag[iz,ir]/rr)*bzeta[iz,ir]^2 + jacobian[iz,ir] = 1.0 + + cvdriftr[iz,ir] = 0.0 + cvdriftz[iz,ir] = -(bzeta[iz,ir]/Bmag[iz,ir])*(bzeta[iz,ir]^2)/rr + gbdriftr[iz,ir] = 0.0 + gbdriftz[iz,ir] = cvdriftz[iz,ir] + end + end + else + input_option_error("$option", option) + end + + geometry = geometric_coefficients(geometry_input_data, rhostar, + Bzed,Bzeta,Bmag,bzed,bzeta,dBdz,dBdr,jacobian, + cvdriftr,cvdriftz,gbdriftr,gbdriftz) + return geometry +end + +end diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index 29b57d62a..0d8baee8f 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -148,7 +148,7 @@ function init_pdf_and_moments!(pdf, moments, boundary_distributions, geometry, init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa, vperp, z, r, composition.n_ion_species, composition.n_neutral_species, - geometry, composition, species, + geometry.input, composition, species, manufactured_solns_input) else n_ion_species = composition.n_ion_species @@ -1134,8 +1134,8 @@ enforce boundary conditions in vpa and z on the evolved pdf; also enforce boundary conditions in z on all separately evolved velocity space moments of the pdf """ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, moments, vpa_bc, - z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, z_adv, r_adv, composition, scratch_dummy, - r_diffusion, vpa_diffusion) + z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, r_adv, composition, scratch_dummy, + r_diffusion, vpa_diffusion, vperp_diffusion) if vpa.n > 1 begin_s_r_z_vperp_region() @loop_s_r_z_vperp is ir iz ivperp begin @@ -1148,7 +1148,8 @@ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, moments, v end if vperp.n > 1 begin_s_r_z_vpa_region() - @views enforce_vperp_boundary_condition!(f, vperp.bc, vperp, vperp_spectral) + @views enforce_vperp_boundary_condition!(f, vperp.bc, vperp, vperp_spectral, + vperp_adv, vperp_diffusion) end if z.n > 1 begin_s_r_vperp_vpa_region() @@ -1169,12 +1170,12 @@ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, moments, v end end function enforce_boundary_conditions!(fvec_out::scratch_pdf, moments, f_r_bc, vpa_bc, - z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, z_adv, r_adv, composition, scratch_dummy, - r_diffusion, vpa_diffusion) + z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, r_adv, composition, scratch_dummy, + r_diffusion, vpa_diffusion, vperp_diffusion) enforce_boundary_conditions!(fvec_out.pdf, f_r_bc, fvec_out.density, fvec_out.upar, fvec_out.ppar, moments, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, - vpa_spectral, vperp_spectral, vpa_adv, z_adv, - r_adv, composition, scratch_dummy, r_diffusion, vpa_diffusion) + vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, + r_adv, composition, scratch_dummy, r_diffusion, vpa_diffusion, vperp_diffusion) end """ @@ -1479,12 +1480,12 @@ function enforce_neutral_z_boundary_condition!(pdf, density, uz, pz, moments, de # Note, have already calculated moments of ion distribution function(s), # so can use the moments here to get the flux if z.irank == 0 - ion_flux_0 = -density_ion[1,ir,isn] * (upar_ion[1,ir,isn]*geometry.bzed - 0.5*geometry.rhostar*Er[1,ir]) + ion_flux_0 = -density_ion[1,ir,isn] * (upar_ion[1,ir,isn]*geometry.bzed[1,ir] - 0.5*geometry.rhostar*Er[1,ir]) else ion_flux_0 = NaN end if z.irank == z.nrank - 1 - ion_flux_L = density_ion[end,ir,isn] * (upar_ion[end,ir,isn]*geometry.bzed - 0.5*geometry.rhostar*Er[end,ir]) + ion_flux_L = density_ion[end,ir,isn] * (upar_ion[end,ir,isn]*geometry.bzed[end,ir] - 0.5*geometry.rhostar*Er[end,ir]) else ion_flux_L = NaN end @@ -2173,20 +2174,24 @@ end """ enforce zero boundary condition at vperp -> infinity """ -function enforce_vperp_boundary_condition!(f, bc, vperp, vperp_spectral) +function enforce_vperp_boundary_condition!(f, bc, vperp, vperp_spectral, vperp_advect, diffusion) if bc == "zero" nvperp = vperp.n ngrid = vperp.ngrid # set zero boundary condition @loop_s_r_z_vpa is ir iz ivpa begin - f[ivpa,nvperp,iz,ir,is] = 0.0 + if diffusion || vperp_advect[is].speed[nvperp,ivpa,iz,ir] < 0.0 + f[ivpa,nvperp,iz,ir,is] = 0.0 + end end # set regularity condition d F / d vperp = 0 at vperp = 0 if vperp.discretization == "gausslegendre_pseudospectral" || vperp.discretization == "chebyshev_pseudospectral" D0 = vperp_spectral.radau.D0 @loop_s_r_z_vpa is ir iz ivpa begin - # adjust F(vperp = 0) so that d F / d vperp = 0 at vperp = 0 - f[ivpa,1,iz,ir,is] = -sum(D0[2:ngrid].*f[ivpa,2:ngrid,iz,ir,is])/D0[1] + if diffusion || vperp_advect[is].speed[1,ivpa,iz,ir] > 0.0 + # adjust F(vperp = 0) so that d F / d vperp = 0 at vperp = 0 + f[ivpa,1,iz,ir,is] = -sum(D0[2:ngrid].*f[ivpa,2:ngrid,iz,ir,is])/D0[1] + end end else println("vperp.bc=\"$bc\" not supported by discretization " diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 99c4be9eb..8e8e67bbe 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -50,6 +50,7 @@ end """ mutable struct advance_info vpa_advection::Bool + vperp_advection::Bool z_advection::Bool r_advection::Bool neutral_z_advection::Bool @@ -77,6 +78,7 @@ mutable struct advance_info manufactured_solns_test::Bool r_diffusion::Bool #flag to control how r bc is imposed when r diffusion terms are present vpa_diffusion::Bool #flag to control how vpa bc is imposed when vpa diffusion terms are present + vperp_diffusion::Bool #flag to control how vperp bc is imposed when vperp diffusion terms are present vz_diffusion::Bool #flag to control how vz bc is imposed when vz diffusion terms are present end @@ -325,19 +327,15 @@ end """ """ -mutable struct geometry_input - # Bz/Bref - Bzed::mk_float - # Btot/Bref - Bmag::mk_float - # bz -- unit vector component in z direction - bzed::mk_float - # bz -- unit vector component in zeta direction - bzeta::mk_float - # Bzeta/Bref - Bzeta::mk_float +Base.@kwdef struct geometry_input # rhostar ion (ref) - rhostar::mk_float #used to premultiply ExB drift terms + rhostar::mk_float = 0.0 #used to premultiply ExB drift terms + # magnetic geometry option + option::String = "constant-helical" # "1D-mirror" + # pitch ( = Bzed/Bmag if geometry_option == "constant-helical") + pitch::mk_float = 1.0 + # DeltaB ( = (Bzed(z=L/2) - Bzed(0))/Bref if geometry_option == "1D-mirror") + DeltaB::mk_float = 0.0 end @enum binary_format_type hdf5 netcdf diff --git a/moment_kinetics/src/load_data.jl b/moment_kinetics/src/load_data.jl index 13f8ae0c3..0250d667a 100644 --- a/moment_kinetics/src/load_data.jl +++ b/moment_kinetics/src/load_data.jl @@ -629,8 +629,8 @@ function reload_evolving_fields!(pdf, moments, boundary_distributions, restart_p neutral_1V = (vzeta.n_global == 1 && vr.n_global == 1) restart_neutral_1V = (restart_vzeta.n_global == 1 && restart_vr.n_global == 1) - if geometry.bzeta != 0.0 && ((neutral_1V && !restart_neutral_1V) || - (!neutral_1V && restart_neutral_1V)) + if any(geometry.bzeta .!= 0.0) && ((neutral_1V && !restart_neutral_1V) || + (!neutral_1V && restart_neutral_1V)) # One but not the other of the run being restarted from and this run are # 1V, but the interpolation below does not allow for vz and vpa being in # different directions. Therefore interpolation between 1V and 3V cases diff --git a/moment_kinetics/src/manufactured_solns.jl b/moment_kinetics/src/manufactured_solns.jl index fb63ac598..c5feaa503 100644 --- a/moment_kinetics/src/manufactured_solns.jl +++ b/moment_kinetics/src/manufactured_solns.jl @@ -5,6 +5,7 @@ module manufactured_solns export manufactured_solutions export manufactured_sources export manufactured_electric_fields +export manufactured_geometry using ..array_allocation: allocate_shared_float using ..looping @@ -13,6 +14,7 @@ using ..type_definitions: mk_float, mk_int function manufactured_solutions end function manufactured_sources_setup end function manufactured_electric_fields end +function manufactured_geometry end function __init__() try @@ -41,7 +43,7 @@ function manufactured_sources(manufactured_solns_input, r_coord, z_coord, vperp_ Source_i_array = allocate_shared_float(vpa_coord.n,vperp_coord.n,z_coord.n,r_coord.n) begin_s_r_z_region() - println("here loop thing ", looping.loop_ranges[].s) + #println("here loop thing ", looping.loop_ranges[].s) @loop_s is begin if is == 1 @loop_r_z_vperp_vpa ir iz ivperp ivpa begin diff --git a/moment_kinetics/src/moment_kinetics.jl b/moment_kinetics/src/moment_kinetics.jl index 80c73ace1..970f01fae 100644 --- a/moment_kinetics/src/moment_kinetics.jl +++ b/moment_kinetics/src/moment_kinetics.jl @@ -32,6 +32,7 @@ include("input_structs.jl") include("reference_parameters.jl") include("coordinates.jl") include("file_io.jl") +include("geo.jl") include("velocity_moments.jl") include("velocity_grid_transforms.jl") include("em_fields.jl") diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 47f71edfb..691571fca 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -6,6 +6,7 @@ export mk_input export performance_test #export advective_form export read_input_file +export get_default_rhostar using ..type_definitions: mk_float, mk_int using ..array_allocation: allocate_float @@ -18,6 +19,7 @@ using ..finite_differences: fd_check_option using ..input_structs using ..numerical_dissipation: setup_numerical_dissipation using ..reference_parameters +using ..geo: init_magnetic_geometry, setup_geometry_input using MPI using TOML @@ -47,6 +49,17 @@ other situations (e.g. when post-processing). """ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) + # Check for input options that used to exist, but do not any more. If these are + # present, the user probably needs to update their input file. + removed_options_list = ("Bzed", "Bmag", "rhostar", "geometry_option", "pitch", "DeltaB") + for opt in removed_options_list + if opt ∈ keys(scan_input) + error("Option '$opt' is no longer used. Please update your input file. You " + * "may need to set some new options to replicate the effect of the " + * "removed ones.") + end + end + # n_ion_species is the number of evolved ion species # currently only n_ion_species = 1 is supported n_ion_species = get(scan_input, "n_ion_species", 1) @@ -60,7 +73,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # reference value using J_||e + J_||i = 0 at z = 0 electron_physics = get(scan_input, "electron_physics", boltzmann_electron_response) - z, r, vpa, vperp, gyrophase, vz, vr, vzeta, species, composition, drive, evolve_moments, collisions, geometry = + z, r, vpa, vperp, gyrophase, vz, vr, vzeta, species, composition, drive, evolve_moments, collisions = load_defaults(n_ion_species, n_neutral_species, electron_physics) # this is the prefix for all output files associated with this run @@ -100,15 +113,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) reference_params = setup_reference_parameters(scan_input) ## set geometry_input - geometry.Bzed = get(scan_input, "Bzed", 1.0) - geometry.Bmag = get(scan_input, "Bmag", 1.0) - geometry.bzed = geometry.Bzed/geometry.Bmag - geometry.bzeta = sqrt(1.0 - geometry.bzed^2.0) - geometry.Bzeta = geometry.Bmag*geometry.bzeta - geometry.rhostar = get(scan_input, "rhostar", get_default_rhostar(reference_params)) - #println("Info: Bzed is ",geometry.Bzed) - #println("Info: Bmag is ",geometry.Bmag) - #println("Info: rhostar is ",geometry.rhostar) + geometry_in = setup_geometry_input(scan_input, get_default_rhostar(reference_params)) ispecies = 1 species.charged[1].z_IC.initialization_option = get(scan_input, "z_IC_option$ispecies", "gaussian") @@ -537,6 +542,13 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) else io = devnull end + + geometry = init_magnetic_geometry(geometry_in,z,r) + if any(geometry.dBdz .!= 0.0) && + (evolve_moments.density || evolve_moments.parallel_flow || + evolve_moments.parallel_pressure) + error("Mirror terms not yet implemented for moment-kinetic modes") + end # check input (and initialized coordinate structs) to catch errors/unsupported options check_input(io, output_dir, nstep, dt, r, z, vpa, vperp, composition, @@ -1003,15 +1015,8 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) nuii = 0.0 collisions = collisions_input(charge_exchange, ionization, constant_ionization_rate, krook_collision_frequency_prefactor,"none", nuii) - Bzed = 1.0 # magnetic field component along z - Bmag = 1.0 # magnetic field strength - bzed = 1.0 # component of b unit vector along z - bzeta = 0.0 # component of b unit vector along zeta - Bzeta = 0.0 # magnetic field component along zeta - rhostar = 0.0 #rhostar of ions for ExB drift - geometry = geometry_input(Bzed,Bmag,bzed,bzeta,Bzeta,rhostar) - return z, r, vpa, vperp, gyrophase, vz, vr, vzeta, species, composition, drive, evolve_moments, collisions, geometry + return z, r, vpa, vperp, gyrophase, vz, vr, vzeta, species, composition, drive, evolve_moments, collisions end """ diff --git a/moment_kinetics/src/r_advection.jl b/moment_kinetics/src/r_advection.jl index dcef106e6..d205b04a1 100644 --- a/moment_kinetics/src/r_advection.jl +++ b/moment_kinetics/src/r_advection.jl @@ -93,10 +93,23 @@ function update_speed_r!(advect, upar, vth, fields, evolve_upar, evolve_ppar, vp error("r_advection is not compatible with evolve_upar or evolve_ppar") end if r.advection.option == "default" && r.n > 1 - ExBfac = 0.5*geometry.rhostar + Bmag = geometry.Bmag + rhostar = geometry.rhostar + ExBfac = 0.5*rhostar + bzeta = geometry.bzeta + jacobian = geometry.jacobian + geofac = r.scratch + cvdriftr = geometry.cvdriftr + gbdriftr = geometry.gbdriftr @inbounds begin @loop_z_vperp_vpa iz ivperp ivpa begin - @views advect.speed[:,ivpa,ivperp,iz] .= ExBfac*fields.Ez[iz,:] + # ExB drift + @. geofac = bzeta[iz,:]*jacobian[iz,:]/Bmag[iz,:] + @views @. advect.speed[:,ivpa,ivperp,iz] = ExBfac*geofac*fields.Ez[iz,:] + # magnetic curvature drift + @. @views advect.speed[:,ivpa,ivperp,iz] += rhostar*(vpa.grid[ivpa]^2)*cvdriftr[iz,:] + # magnetic grad B drift + @. @views advect.speed[:,ivpa,ivperp,iz] += 0.5*rhostar*(vperp.grid[ivperp]^2)*gbdriftr[iz,:] end end elseif r.advection.option == "default" && r.n == 1 diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index 4e27925e2..796647860 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -288,11 +288,15 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral, begin_serial_region() vperp_advect = setup_advection(n_ion_species, vperp, vpa, z, r) # initialise the vperp advection speed + # Note that z_advect and r_advect are arguments of update_speed_vperp! + # This means that z_advect[is].speed and r_advect[is].speed are used to determine + # vperp_advect[is].speed, so z_advect and r_advect must always be updated before + # vperp_advect is updated and used. if vperp.n > 1 begin_serial_region() @serial_region begin for is ∈ 1:n_ion_species - @views update_speed_vperp!(vperp_advect[is], vpa, vperp, z, r) + @views update_speed_vperp!(vperp_advect[is], vpa, vperp, z, r, z_advect[is], r_advect[is], geometry) end end end @@ -355,7 +359,7 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral, if(advance.manufactured_solns_test) manufactured_source_list = manufactured_sources(manufactured_solns_input, r, z, vperp, vpa, vzeta, vr, vz, - composition, geometry, collisions, + composition, geometry.input, collisions, num_diss_params, species) else manufactured_source_list = false # dummy Bool to be passed as argument instead of list @@ -372,8 +376,9 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral, pdf.charged.norm, boundary_distributions.pdf_rboundary_charged, moments.charged.dens, moments.charged.upar, moments.charged.ppar, moments, vpa.bc, z.bc, r.bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, - vpa_advect, z_advect, r_advect, - composition, scratch_dummy, advance.r_diffusion, advance.vpa_diffusion) + vpa_advect, vperp_advect, z_advect, r_advect, + composition, scratch_dummy, advance.r_diffusion, + advance.vpa_diffusion, advance.vperp_diffusion) # Ensure normalised pdf exactly obeys integral constraints if evolving moments begin_s_r_z_region() @loop_s_r_z is ir iz begin @@ -451,6 +456,7 @@ function setup_advance_flags(moments, composition, t_input, collisions, vr, vz) # default is not to concurrently advance different operators advance_vpa_advection = false + advance_vperp_advection = false advance_z_advection = false advance_r_advection = false advance_cx_1V = false @@ -475,15 +481,17 @@ function setup_advance_flags(moments, composition, t_input, collisions, advance_neutral_energy = false r_diffusion = false vpa_diffusion = false + vperp_diffusion = false vz_diffusion = false explicit_weakform_fp_collisions = false # all advance flags remain false if using operator-splitting # otherwise, check to see if the flags need to be set to true if !t_input.split_operators # default for non-split operators is to include both vpa and z advection together - advance_vpa_advection = true && vpa.n > 1 && z.n > 1 - advance_z_advection = true && z.n > 1 - advance_r_advection = true && r.n > 1 + advance_vpa_advection = vpa.n > 1 && z.n > 1 + advance_vperp_advection = vperp.n > 1 && z.n > 1 + advance_z_advection = z.n > 1 + advance_r_advection = r.n > 1 if collisions.nuii > 0.0 && vperp.n > 1 explicit_weakform_fp_collisions = true else @@ -578,12 +586,13 @@ function setup_advance_flags(moments, composition, t_input, collisions, r_diffusion = (advance_numerical_dissipation && num_diss_params.r_dissipation_coefficient > 0.0) # flag to determine if a d^2/dvpa^2 operator is present vpa_diffusion = ((advance_numerical_dissipation && num_diss_params.vpa_dissipation_coefficient > 0.0) || explicit_weakform_fp_collisions) + vperp_diffusion = ((advance_numerical_dissipation && num_diss_params.vperp_dissipation_coefficient > 0.0) || explicit_weakform_fp_collisions) vz_diffusion = (advance_numerical_dissipation && num_diss_params.vz_dissipation_coefficient > 0.0) end manufactured_solns_test = manufactured_solns_input.use_for_advance - return advance_info(advance_vpa_advection, advance_z_advection, advance_r_advection, + return advance_info(advance_vpa_advection, advance_vperp_advection, advance_z_advection, advance_r_advection, advance_neutral_z_advection, advance_neutral_r_advection, advance_neutral_vz_advection, advance_cx, advance_cx_1V, advance_ionization, advance_ionization_1V, @@ -594,7 +603,7 @@ function setup_advance_flags(moments, composition, t_input, collisions, advance_energy, advance_neutral_external_source, advance_neutral_sources, advance_neutral_continuity, advance_neutral_force_balance, advance_neutral_energy, rk_coefs, - manufactured_solns_test, r_diffusion, vpa_diffusion, vz_diffusion) + manufactured_solns_test, r_diffusion, vpa_diffusion, vperp_diffusion, vz_diffusion) end function setup_dummy_and_buffer_arrays(nr,nz,nvpa,nvperp,nvz,nvr,nvzeta,nspecies_ion,nspecies_neutral) @@ -1236,7 +1245,7 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v z_spectral, r_spectral, vpa_spectral, vperp_spectral = spectral_objects.z_spectral, spectral_objects.r_spectral, spectral_objects.vpa_spectral, spectral_objects.vperp_spectral vzeta_spectral, vr_spectral, vz_spectral = spectral_objects.vzeta_spectral, spectral_objects.vr_spectral, spectral_objects.vz_spectral - vpa_advect, r_advect, z_advect = advect_objects.vpa_advect, advect_objects.r_advect, advect_objects.z_advect + vpa_advect, vperp_advect, r_advect, z_advect = advect_objects.vpa_advect, advect_objects.vperp_advect, advect_objects.r_advect, advect_objects.z_advect neutral_z_advect, neutral_r_advect, neutral_vz_advect = advect_objects.neutral_z_advect, advect_objects.neutral_r_advect, advect_objects.neutral_vz_advect ## @@ -1265,8 +1274,8 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v enforce_boundary_conditions!(new_scratch, moments, boundary_distributions.pdf_rboundary_charged, vpa.bc, z.bc, r.bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, - vpa_advect, z_advect, r_advect, composition, scratch_dummy, - advance.r_diffusion, advance.vpa_diffusion) + vpa_advect, vperp_advect, z_advect, r_advect, composition, scratch_dummy, + advance.r_diffusion, advance.vpa_diffusion, advance.vperp_diffusion) if moments.evolve_density && moments.enforce_conservation begin_s_r_z_region() @@ -1622,7 +1631,7 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, vpa_spectral, vperp_spectral, r_spectral, z_spectral = spectral_objects.vpa_spectral, spectral_objects.vperp_spectral, spectral_objects.r_spectral, spectral_objects.z_spectral vz_spectral, vr_spectral, vzeta_spectral = spectral_objects.vz_spectral, spectral_objects.vr_spectral, spectral_objects.vzeta_spectral - vpa_advect, r_advect, z_advect = advect_objects.vpa_advect, advect_objects.r_advect, advect_objects.z_advect + vpa_advect, vperp_advect, r_advect, z_advect = advect_objects.vpa_advect, advect_objects.vperp_advect, advect_objects.r_advect, advect_objects.z_advect neutral_z_advect, neutral_r_advect, neutral_vz_advect = advect_objects.neutral_z_advect, advect_objects.neutral_r_advect, advect_objects.neutral_vz_advect if advance.external_source @@ -1652,10 +1661,12 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, r_advection!(fvec_out.pdf, fvec_in, moments, fields, r_advect, r, z, vperp, vpa, dt, r_spectral, composition, geometry, scratch_dummy) end - - #if advance.vperp_advection - # PLACEHOLDER - #end + # vperp_advection requires information about z and r advection + # so call vperp_advection! only after z and r advection routines + if advance.vperp_advection + vperp_advection!(fvec_out.pdf, fvec_in, vperp_advect, r, z, vperp, vpa, + dt, vperp_spectral, composition, z_advect, r_advect, geometry) + end if advance.source_terms source_terms!(fvec_out.pdf, fvec_in, moments, vpa, z, r, dt, z_spectral, diff --git a/moment_kinetics/src/velocity_grid_transforms.jl b/moment_kinetics/src/velocity_grid_transforms.jl index 368ee73c2..ccb86b684 100644 --- a/moment_kinetics/src/velocity_grid_transforms.jl +++ b/moment_kinetics/src/velocity_grid_transforms.jl @@ -25,13 +25,15 @@ function vzvrvzeta_to_vpavperp!(f_out,f_in,vz,vr,vzeta,vpa,vperp,gyrophase,z,r,g begin_sn_r_z_region() @loop_sn_r_z isn ir iz begin + bzed = geometry.bzed[iz,ir] + bzeta = geometry.bzeta[iz,ir] @views vzvrvzeta_to_vpavperp_species!(f_out[:,:,iz,ir,isn],f_in[:,:,:,iz,ir,isn], - vz,vr,vzeta,vpa,vperp,gyrophase,geometry) + vz,vr,vzeta,vpa,vperp,gyrophase,bzed,bzeta) end end -function vzvrvzeta_to_vpavperp_species!(f_out,f_in,vz,vr,vzeta,vpa,vperp,gyrophase,geometry) +function vzvrvzeta_to_vpavperp_species!(f_out,f_in,vz,vr,vzeta,vpa,vperp,gyrophase,bzed,bzeta) @boundscheck vz.n == size(f_in, 1) || throw(BoundsError(f_in)) @boundscheck vr.n == size(f_in, 2) || throw(BoundsError(f_in)) @boundscheck vzeta.n == size(f_in, 3) || throw(BoundsError(f_in)) @@ -41,10 +43,7 @@ function vzvrvzeta_to_vpavperp_species!(f_out,f_in,vz,vr,vzeta,vpa,vperp,gyropha pdf_interp = linear_interpolation((vz.grid,vr.grid,vzeta.grid),f_in,extrapolation_bc = 0.0) # pdf_interp( vz_val, vr_val, vzeta_val) is interpolated value of f_in # extrapolation_bc = 0.0 makes pdf_interp = 0.0 for |vx| > vx.L/2 (x = z,r,zeta) - - bzed = geometry.bzed - bzeta = geometry.bzeta - + @loop_vperp_vpa ivperp ivpa begin # for each ivpa, ivperp, compute gyroaverage of f_in # use @@ -78,14 +77,16 @@ function vpavperp_to_vzvrvzeta!(f_out,f_in,vz,vr,vzeta,vpa,vperp,z,r,geometry,co begin_s_r_z_region() @loop_s_r_z is ir iz begin + bzed = geometry.bzed[iz,ir] + bzeta = geometry.bzeta[iz,ir] @views vpavperp_to_vzvrvzeta_species!(f_out[:,:,:,iz,ir,is],f_in[:,:,iz,ir,is], - vz,vr,vzeta,vpa,vperp,geometry) + vz,vr,vzeta,vpa,vperp,bzed,bzeta) end end -function vpavperp_to_vzvrvzeta_species!(f_out,f_in,vz,vr,vzeta,vpa,vperp,geometry) +function vpavperp_to_vzvrvzeta_species!(f_out,f_in,vz,vr,vzeta,vpa,vperp,bzed,bzeta) @boundscheck vz.n == size(f_out, 1) || throw(BoundsError(f_out)) @boundscheck vr.n == size(f_out, 2) || throw(BoundsError(f_out)) @boundscheck vzeta.n == size(f_out, 3) || throw(BoundsError(f_out)) @@ -95,10 +96,6 @@ function vpavperp_to_vzvrvzeta_species!(f_out,f_in,vz,vr,vzeta,vpa,vperp,geometr pdf_interp = linear_interpolation((vpa.grid,vperp.grid),f_in,extrapolation_bc = 0.0) # pdf_interp( vz_val, vr_val, vzeta_val) is interpolated value of f_in # extrapolation_bc = 0.0 makes pdf_interp = 0.0 for |vpa| > vpa.L/2 and vperp > vperp.L - - bzed = geometry.bzed - bzeta = geometry.bzeta - @loop_vzeta_vr_vz ivzeta ivr ivz begin # for each ivzeta, ivr, ivz interpolate f_in onto f_out # use diff --git a/moment_kinetics/src/velocity_moments.jl b/moment_kinetics/src/velocity_moments.jl index fc2f2d2ed..c3b874728 100644 --- a/moment_kinetics/src/velocity_moments.jl +++ b/moment_kinetics/src/velocity_moments.jl @@ -869,7 +869,7 @@ function update_chodura!(moments,ff,vpa,vperp,z,r,r_spectral,composition,geometr if z.irank == 0 @loop_s_r is ir begin @views moments.charged.chodura_integral_lower[ir,is] = update_chodura_integral_species!(ff[:,:,1,ir,is],dffdr[:,:,1,ir,is], - ff_dummy[:,:,1,ir,is],vpa,vperp,z,r,composition,geometry,z_advect[is].speed[1,:,:,ir],moments.charged.dens[1,ir,is],del_vpa) + ff_dummy[:,:,1,ir,is],vpa,vperp,z,r,composition,geometry,z_advect[is].speed[1,:,:,ir],moments.charged.dens[1,ir,is],del_vpa,1,ir) end else # we do not save this Chodura integral to the output file @loop_s_r is ir begin @@ -879,7 +879,7 @@ function update_chodura!(moments,ff,vpa,vperp,z,r,r_spectral,composition,geometr if z.irank == z.nrank - 1 @loop_s_r is ir begin @views moments.charged.chodura_integral_upper[ir,is] = update_chodura_integral_species!(ff[:,:,end,ir,is],dffdr[:,:,end,ir,is], - ff_dummy[:,:,end,ir,is],vpa,vperp,z,r,composition,geometry,z_advect[is].speed[end,:,:,ir],moments.charged.dens[end,ir,is],del_vpa) + ff_dummy[:,:,end,ir,is],vpa,vperp,z,r,composition,geometry,z_advect[is].speed[end,:,:,ir],moments.charged.dens[end,ir,is],del_vpa,z.n,ir) end else # we do not save this Chodura integral to the output file @loop_s_r is ir begin @@ -902,7 +902,7 @@ Chodura condition to a single species plasma with Z = 1 """ -function update_chodura_integral_species!(ff,dffdr,ff_dummy,vpa,vperp,z,r,composition,geometry,vz,dens,del_vpa) +function update_chodura_integral_species!(ff,dffdr,ff_dummy,vpa,vperp,z,r,composition,geometry,vz,dens,del_vpa,iz,ir) @boundscheck vpa.n == size(ff, 1) || throw(BoundsError(ff)) @boundscheck vperp.n == size(ff, 2) || throw(BoundsError(ff)) @boundscheck vpa.n == size(dffdr, 1) || throw(BoundsError(dffdr)) @@ -915,7 +915,7 @@ function update_chodura_integral_species!(ff,dffdr,ff_dummy,vpa,vperp,z,r,compos # we are more than a vpa mimimum grid spacing away from # the vz(vpa,r) = 0 velocity boundary if abs(vz[ivpa,ivperp]) > 0.5*del_vpa - ff_dummy[ivpa,ivperp] = (ff[ivpa,ivperp]*bzed^2/(vz[ivpa,ivperp]^2) + + ff_dummy[ivpa,ivperp] = (ff[ivpa,ivperp]*bzed[iz,ir]^2/(vz[ivpa,ivperp]^2) + geometry.rhostar*dffdr[ivpa,ivperp]/vz[ivpa,ivperp]) else ff_dummy[ivpa,ivperp] = 0.0 diff --git a/moment_kinetics/src/vpa_advection.jl b/moment_kinetics/src/vpa_advection.jl index c6dab7fb5..47da74aff 100644 --- a/moment_kinetics/src/vpa_advection.jl +++ b/moment_kinetics/src/vpa_advection.jl @@ -42,8 +42,9 @@ function update_speed_vpa!(advect, fields, fvec, moments, vpa, vperp, z, r, comp @boundscheck composition.n_ion_species == size(advect,1) || throw(BoundsError(advect)) @boundscheck vpa.n == size(advect[1].speed,1) || throw(BoundsError(speed)) if vpa.advection.option == "default" - # dvpa/dt = Ze/m ⋅ E_parallel - update_speed_default!(advect, fields, fvec, moments, vpa, z, r, composition, + # dvpa/dt = Ze/m ⋅ E_parallel - (vperp^2/2B) bz dB/dz + # magnetic mirror term only supported for standard DK implementation + update_speed_default!(advect, fields, fvec, moments, vpa, vperp, z, r, composition, collisions, ion_source_settings, t, geometry) elseif vpa.advection.option == "constant" begin_serial_region() @@ -69,7 +70,7 @@ end """ """ -function update_speed_default!(advect, fields, fvec, moments, vpa, z, r, composition, +function update_speed_default!(advect, fields, fvec, moments, vpa, vperp, z, r, composition, collisions, ion_source_settings, t, geometry) if moments.evolve_ppar && moments.evolve_upar update_speed_n_u_p_evolution!(advect, fvec, moments, vpa, z, r, composition, @@ -82,10 +83,15 @@ function update_speed_default!(advect, fields, fvec, moments, vpa, z, r, composi collisions, ion_source_settings) else bzed = geometry.bzed + dBdz = geometry.dBdz + Bmag = geometry.Bmag @inbounds @fastmath begin @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin + # mu, the adiabatic invariant + mu = 0.5*(vperp.grid[ivperp]^2)/Bmag[iz,ir] # bzed = B_z/B - advect[is].speed[ivpa,ivperp,iz,ir] = 0.5*bzed*fields.Ez[iz,ir] + advect[is].speed[ivpa,ivperp,iz,ir] = (0.5*bzed[iz,ir]*fields.Ez[iz,ir] - + mu*bzed[iz,ir]*dBdz[iz,ir]) end end end diff --git a/moment_kinetics/src/vperp_advection.jl b/moment_kinetics/src/vperp_advection.jl index 0fc738923..16c2697df 100644 --- a/moment_kinetics/src/vperp_advection.jl +++ b/moment_kinetics/src/vperp_advection.jl @@ -8,28 +8,50 @@ using ..chebyshev: chebyshev_info using ..looping # do a single stage time advance (potentially as part of a multi-stage RK scheme) -function vperp_advection!(f_out, fvec_in, advect, r, z, vperp, vpa, - dt, spectral, composition) +function vperp_advection!(f_out, fvec_in, vperp_advect, r, z, vperp, vpa, + dt, vperp_spectral, composition, z_advect, r_advect, geometry) + begin_s_r_z_vpa_region() @loop_s is begin # get the updated speed along the r direction using the current f - @views update_speed_vperp!(advect[is], vpa, vperp, z, r) + @views update_speed_vperp!(vperp_advect[is], vpa, vperp, z, r, z_advect[is], r_advect[is], geometry) @loop_r_z_vpa ir iz ivpa begin - @views advance_f_local!(f_out[ivpa,:,iz,ir,is], vperp.scratch, advect[is], ivpa, - r, dt, spectral) + @views advance_f_local!(f_out[ivpa,:,iz,ir,is], fvec_in.pdf[ivpa,:,iz,ir,is], + vperp_advect[is], ivpa, iz, ir, vperp, dt, vperp_spectral) end end end # calculate the advection speed in the vperp-direction at each grid point -function update_speed_vperp!(advect, vpa, vperp, z, r) - @boundscheck z.n == size(advect.speed,3) || throw(BoundsError(advect)) - @boundscheck vperp.n == size(advect.speed,1) || throw(BoundsError(advect)) - @boundscheck vpa.n == size(advect.speed,2) || throw(BoundsError(advect)) - @boundscheck r.n == size(advect.speed,4) || throw(BoundsError(speed)) - if vperp.advection.option == "default" || vperp.advection.option == "constant" +# note that the vperp advection speed depends on the z and r advection speeds +# It is important to ensure that z_advect and r_advect are updated before vperp_advect +function update_speed_vperp!(vperp_advect, vpa, vperp, z, r, z_advect, r_advect, geometry) + @boundscheck z.n == size(vperp_advect.speed,3) || throw(BoundsError(vperp_advect)) + @boundscheck vperp.n == size(vperp_advect.speed,1) || throw(BoundsError(vperp_advect)) + @boundscheck vpa.n == size(vperp_advect.speed,2) || throw(BoundsError(vperp_advect)) + @boundscheck r.n == size(vperp_advect.speed,4) || throw(BoundsError(vperp_speed)) + if vperp.advection.option == "default" + # advection of vperp due to conservation of + # the adiabatic invariant mu = vperp^2 / 2 B + dzdt = vperp.scratch + drdt = vperp.scratch2 + dBdr = geometry.dBdr + dBdz = geometry.dBdz + Bmag = geometry.Bmag + rfac = 0.0 + if r.n > 1 + rfac = 1.0 + end + @inbounds begin + @loop_r_z_vpa ir iz ivpa begin + @. @views dzdt = z_advect.speed[iz,ivpa,:,ir] + @. @views drdt = rfac*r_advect.speed[ir,ivpa,:,iz] + @. @views vperp_advect.speed[:,ivpa,iz,ir] = (0.5*vperp.grid[:]/Bmag[iz,ir])*(dzdt[:]*dBdz[iz,ir] + drdt[:]*dBdr[iz,ir]) + end + end + elseif vperp.advection.option == "constant" @inbounds begin @loop_r_z_vpa ir iz ivpa begin - @views advect.speed[:,ivpa,iz,ir] .= vperp.advection.constant_speed + @views vperp_advect.speed[:,ivpa,iz,ir] .= vperp.advection.constant_speed end end end diff --git a/moment_kinetics/src/z_advection.jl b/moment_kinetics/src/z_advection.jl index f67996009..8a4a3f273 100644 --- a/moment_kinetics/src/z_advection.jl +++ b/moment_kinetics/src/z_advection.jl @@ -93,10 +93,24 @@ function update_speed_z!(advect, upar, vth, evolve_upar, evolve_ppar, fields, vp if z.advection.option == "default" # bzed = B_z/B only used for z.advection.option == "default" bzed = geometry.bzed - ExBfac = -0.5*geometry.rhostar + Bmag = geometry.Bmag + bzeta = geometry.bzeta + jacobian = geometry.jacobian + rhostar = geometry.rhostar + ExBfac = -0.5*rhostar + geofac = z.scratch + cvdriftz = geometry.cvdriftz + gbdriftz = geometry.gbdriftz @inbounds begin @loop_r_vperp_vpa ir ivperp ivpa begin - @. @views advect.speed[:,ivpa,ivperp,ir] = vpa.grid[ivpa]*bzed + ExBfac*fields.Er[:,ir] + # vpa bzed + @. @views advect.speed[:,ivpa,ivperp,ir] = vpa.grid[ivpa]*bzed[:,ir] + # ExB drift + @. @views advect.speed[:,ivpa,ivperp,ir] += ExBfac*bzeta[:,ir]*jacobian[:,ir]/Bmag[:,ir]*fields.Er[:,ir] + # magnetic curvature drift + @. @views advect.speed[:,ivpa,ivperp,ir] += rhostar*(vpa.grid[ivpa]^2)*cvdriftz[:,ir] + # magnetic grad B drift + @. @views advect.speed[:,ivpa,ivperp,ir] += 0.5*rhostar*(vperp.grid[ivperp]^2)*gbdriftz[:,ir] end if evolve_ppar @loop_r_vperp_vpa ir ivperp ivpa begin diff --git a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl index a022170cb..cf324766a 100644 --- a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl +++ b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl @@ -136,9 +136,6 @@ test_input_gauss_legendre = Dict("run_name" => "gausslegendre_pseudospectral", "electron_physics" => "boltzmann_electron_response", "nuii" => 1.0, "use_semi_lagrange" => false, - "Bzed" => 1.0, - "Bmag" => 1.0, - "rhostar" => 1.0, "z_IC_upar_amplitude1" => 0.0, "z_IC_density_amplitude1" => 0.001, "z_IC_upar_amplitude2" => 0.0, diff --git a/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl b/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl index 579829503..a2c3c5d7d 100644 --- a/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl +++ b/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl @@ -2,7 +2,8 @@ """ module plot_MMS_sequence -export get_MMS_error_data +#export get_MMS_error_data +export run_mms_test # packages using Plots @@ -16,7 +17,7 @@ using ..post_processing_input: pp using ..plots_post_processing: compare_charged_pdf_symbolic_test, compare_fields_symbolic_test using ..plots_post_processing: compare_moments_symbolic_test, compare_neutral_pdf_symbolic_test using ..plots_post_processing: allocate_global_zr_neutral_moments, allocate_global_zr_charged_moments -using ..plots_post_processing: allocate_global_zr_fields, get_geometry_and_composition +using ..plots_post_processing: allocate_global_zr_fields, get_composition using ..plots_post_processing: get_coords_nelement, get_coords_ngrid using moment_kinetics.array_allocation: allocate_float using moment_kinetics.type_definitions: mk_float, mk_int @@ -28,7 +29,9 @@ using moment_kinetics.load_data: load_block_data, load_coordinate_data, load_inp using moment_kinetics.load_data: read_distributed_zr_data!, construct_global_zr_coords using moment_kinetics.velocity_moments: integrate_over_vspace using moment_kinetics.manufactured_solns: manufactured_solutions, manufactured_electric_fields -using moment_kinetics.moment_kinetics_input: mk_input, read_input_file +using moment_kinetics.moment_kinetics_input: mk_input, read_input_file, get_default_rhostar +using moment_kinetics.input_structs: geometry_input +using moment_kinetics.reference_parameters import Base: get @@ -68,12 +71,13 @@ function get_MMS_error_data(path_list,scan_type,scan_name) fid = open_readonly_output_file(run_name,"moments") scan_input = load_input(fid) - # get run-time input/composition/geometry/collisions/species info for convenience - #io_input, evolve_moments, t_input, z, z_spectral, r, r_spectral, vpa, vpa_spectral, - # vperp, vperp_spectral, gyrophase, gyrophase_spectral, vz, vz_spectral, vr, - # vr_spectral, vzeta, vzeta_spectral, composition, species, collisions, geometry, - # drive_input, external_source_settings, num_diss_params, - # manufactured_solns_input = mk_input(scan_input) + input = mk_input(scan_input) + # obtain input options from moment_kinetics_input.jl + # and check input to catch errors + io_input, evolve_moments, t_input, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, + composition, species, collisions, geometry, drive_input, external_source_settings, + num_diss_params, manufactured_solns_input = input + z_nelement, r_nelement, vpa_nelement, vperp_nelement, vz_nelement, vr_nelement, vzeta_nelement = get_coords_nelement(scan_input) z_ngrid, r_ngrid, vpa_ngrid, vperp_ngrid, @@ -132,6 +136,20 @@ function get_MMS_error_data(path_list,scan_type,scan_name) else println("ERROR: scan_type = ",scan_type," requires vpa_nelement = vperp_nelement = z_nelement") end + elseif scan_type == "vpa2vperpz_nelement" + nelement = z_nelement + if nelement == vpa_nelement && nelement == 2*vperp_nelement + nelement_sequence[isim] = nelement + else + println("ERROR: scan_type = ",scan_type," requires vpa_nelement = 2.0*vperp_nelement = z_nelement") + end + elseif scan_type == "vpa2vperpzr_nelement" + nelement = z_nelement + if nelement == vpa_nelement && nelement == 2*vperp_nelement && nelement == r_nelement + nelement_sequence[isim] = nelement + else + println("ERROR: scan_type = ",scan_type," requires vpa_nelement = 2.0*vperp_nelement = z_nelement = r_nelement") + end elseif scan_type == "vpaz_nelement" nelement = z_nelement if nelement == vpa_nelement @@ -190,24 +208,27 @@ function get_MMS_error_data(path_list,scan_type,scan_name) r_global, z_global = construct_global_zr_coords(r, z) #println("z: ",z) #println("r: ",r) - + run_names = tuple(run_name) + nbs = tuple(nblocks) + print(run_names) + iskip = 1 # stride in slice of arrays # fields - read_distributed_zr_data!(phi,"phi",run_name,"moments",nblocks,z.n,r.n) - read_distributed_zr_data!(Ez,"Ez",run_name,"moments",nblocks,z.n,r.n) - read_distributed_zr_data!(Er,"Er",run_name,"moments",nblocks,z.n,r.n) + read_distributed_zr_data!(phi,"phi",run_names,"moments",nbs,z.n,r.n,iskip) + read_distributed_zr_data!(Ez,"Ez",run_names,"moments",nbs,z.n,r.n,iskip) + read_distributed_zr_data!(Er,"Er",run_names,"moments",nbs,z.n,r.n,iskip) # charged particle moments - read_distributed_zr_data!(density,"density",run_name,"moments",nblocks,z.n,r.n) - read_distributed_zr_data!(parallel_flow,"parallel_flow",run_name,"moments",nblocks,z.n,r.n) - read_distributed_zr_data!(parallel_pressure,"parallel_pressure",run_name,"moments",nblocks,z.n,r.n) - read_distributed_zr_data!(parallel_heat_flux,"parallel_heat_flux",run_name,"moments",nblocks,z.n,r.n) - read_distributed_zr_data!(thermal_speed,"thermal_speed",run_name,"moments",nblocks,z.n,r.n) + read_distributed_zr_data!(density,"density",run_names,"moments",nbs,z.n,r.n,iskip) + read_distributed_zr_data!(parallel_flow,"parallel_flow",run_names,"moments",nbs,z.n,r.n,iskip) + read_distributed_zr_data!(parallel_pressure,"parallel_pressure",run_names,"moments",nbs,z.n,r.n,iskip) + read_distributed_zr_data!(parallel_heat_flux,"parallel_heat_flux",run_names,"moments",nbs,z.n,r.n,iskip) + read_distributed_zr_data!(thermal_speed,"thermal_speed",run_names,"moments",nbs,z.n,r.n,iskip) # neutral particle moments if n_neutral_species > 0 - read_distributed_zr_data!(neutral_density,"density_neutral",run_name,"moments",nblocks,z.n,r.n) - read_distributed_zr_data!(neutral_uz,"uz_neutral",run_name,"moments",nblocks,z.n,r.n) - read_distributed_zr_data!(neutral_pz,"pz_neutral",run_name,"moments",nblocks,z.n,r.n) - read_distributed_zr_data!(neutral_qz,"qz_neutral",run_name,"moments",nblocks,z.n,r.n) - read_distributed_zr_data!(neutral_thermal_speed,"thermal_speed_neutral",run_name,"moments",nblocks,z.n,r.n) + read_distributed_zr_data!(neutral_density,"density_neutral",run_names,"moments",nbs,z.n,r.n,iskip) + read_distributed_zr_data!(neutral_uz,"uz_neutral",run_names,"moments",nbs,z.n,r.n,iskip) + read_distributed_zr_data!(neutral_pz,"pz_neutral",run_names,"moments",nbs,z.n,r.n,iskip) + read_distributed_zr_data!(neutral_qz,"qz_neutral",run_names,"moments",nbs,z.n,r.n,iskip) + read_distributed_zr_data!(neutral_thermal_speed,"thermal_speed_neutral",run_names,"moments",nbs,z.n,r.n,iskip) end @@ -219,15 +240,22 @@ function get_MMS_error_data(path_list,scan_type,scan_name) else Lr_in = 1.0 end - geometry, composition = get_geometry_and_composition(scan_input,n_ion_species,n_neutral_species) - - manufactured_solns_list = manufactured_solutions(Lr_in,z.L,r_bc,z_bc,geometry,composition,r.n,vperp.n) + composition = get_composition(scan_input) + + reference_params = setup_reference_parameters(scan_input) + 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) + + manufactured_solns_list = manufactured_solutions(manufactured_solns_input,Lr_in,z.L,r_bc,z_bc,geo_in,composition,species,r.n,vperp.n) dfni_func = manufactured_solns_list.dfni_func densi_func = manufactured_solns_list.densi_func dfnn_func = manufactured_solns_list.dfnn_func densn_func = manufactured_solns_list.densn_func - manufactured_E_fields = manufactured_electric_fields(Lr_in,z.L,r_bc,z_bc,composition,r.n) + manufactured_E_fields = manufactured_electric_fields(Lr_in,z.L,r_bc,z_bc,composition,r.n,manufactured_solns_input,species) Er_func = manufactured_E_fields.Er_func Ez_func = manufactured_E_fields.Ez_func phi_func = manufactured_E_fields.phi_func @@ -323,6 +351,10 @@ function get_MMS_error_data(path_list,scan_type,scan_name) xlabel = L"z "*" & "*L"r "*" "*L"N_{element}" elseif scan_type == "vpavperpz_nelement" xlabel = L"N_{element}(z) = N_{element}(v_\perp) = N_{element}(v_{||})" + elseif scan_type == "vpa2vperpz_nelement" + xlabel = L"N_{element}(z) = 2N_{element}(v_\perp) = N_{element}(v_{||})" + elseif scan_type == "vpa2vperpzr_nelement" + xlabel = L"N_{element}(z) = N_{element}(r) = 2N_{element}(v_\perp) = N_{element}(v_{||})" elseif scan_type == "vpazr_nelement0.25" xlabel = L"N_{element}(z) = N_{element}(r) = N_{element}(v_{||})/4" elseif scan_type == "vpaz_nelement0.25" @@ -353,18 +385,23 @@ function get_MMS_error_data(path_list,scan_type,scan_name) ytick_sequence = Array([1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0e-0,1.0e1]) elseif scan_name == "1D-1V-wall_cheb" || scan_name == "1D-2V-wall_cheb_krook" ytick_sequence = Array([1.0e-13,1.0e-12,1.0e-11,1.0e-10,1.0e-9,1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0e-0,1.0e1]) - elseif scan_name == "1D-3V-wall_cheb-updated" || scan_name == "1D-3V-wall_cheb-new-dfni-Er" || scan_name == "1D-3V-wall_cheb-new-dfni" || scan_name == "2D-sound-wave_cheb" + elseif scan_name == "mirror_wall-1D-2V" + ytick_sequence = Array([1.0e-12,1.0e-11,1.0e-10,1.0e-9,1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0e-0,1.0e1]) + elseif scan_name == "1D-3V-wall_cheb-updated" || scan_name == "1D-3V-wall_cheb-new-dfni-Er" || scan_name == "1D-3V-wall_cheb-new-dfni" || scan_name == "2D-sound-wave_cheb" || scan_name == "mirror_wall-1D-2V" ytick_sequence = Array([1.0e-10,1.0e-9,1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0e-0,1.0e1]) elseif scan_name == "2D-1V-wall_cheb" || scan_name == "2D-1V-wall_cheb-nonzero-Er" || scan_name == "1D-1V-wall_cheb-constant-Er-ngrid-5-opt" ytick_sequence = Array([1.0e-8,1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0e-0]) #ytick_sequence = Array([1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0e-0]) + elseif scan_name == "mirror_wall-2D-2V-ngrid-5" + ytick_sequence = Array([1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0e-0,1.0e1]) else ytick_sequence = Array([1.0e-7,1.0e-6,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,1.0e-0,1.0e1]) end plot(nelement_sequence, [ion_density_error_sequence,ion_pdf_error_sequence], xlabel=xlabel, label=[ylabel_ion_density ylabel_ion_pdf], ylabel="", shape =:circle, xscale=:log10, yscale=:log10, xticks = (nelement_sequence, nelement_sequence), yticks = (ytick_sequence, ytick_sequence), markersize = 5, linewidth=2, - xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize) + xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize, + foreground_color_legend = nothing, background_color_legend = nothing, legend=:bottomleft) outfile = outprefix*".pdf" savefig(outfile) println(outfile) @@ -373,7 +410,8 @@ function get_MMS_error_data(path_list,scan_type,scan_name) plot(nelement_sequence, [ion_density_error_sequence,phi_error_sequence,Er_error_sequence,Ez_error_sequence], xlabel=xlabel, label=[ylabel_ion_density ylabel_phi ylabel_Er ylabel_Ez], ylabel="", shape =:circle, xscale=:log10, yscale=:log10, xticks = (nelement_sequence, nelement_sequence), yticks = (ytick_sequence, ytick_sequence), markersize = 5, linewidth=2, - xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize) + xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize, + foreground_color_legend = nothing, background_color_legend = nothing, legend=:bottomleft) outfile = outprefix*"_fields.pdf" savefig(outfile) println(outfile) @@ -386,7 +424,8 @@ function get_MMS_error_data(path_list,scan_type,scan_name) plot(nelement_sequence, [ion_density_error_sequence,phi_error_sequence,Ez_error_sequence], xlabel=xlabel, label=[ylabel_ion_density ylabel_phi ylabel_Ez], ylabel="", shape =:circle, xscale=:log10, yscale=:log10, xticks = (nelement_sequence, nelement_sequence), yticks = (ytick_sequence, ytick_sequence), markersize = 5, linewidth=2, - xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize) + xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize, + foreground_color_legend = nothing, background_color_legend = nothing, legend=:bottomleft) outfile = outprefix*"_fields_no_Er.pdf" savefig(outfile) println(outfile) @@ -394,16 +433,18 @@ function get_MMS_error_data(path_list,scan_type,scan_name) plot(nelement_sequence, [ion_density_error_sequence,phi_error_sequence,Ez_error_sequence,ion_pdf_error_sequence,expected_scaling], xlabel=xlabel, label=[ylabel_ion_density ylabel_phi ylabel_Ez ylabel_ion_pdf expected_label], ylabel="", shape =:circle, xscale=:log10, yscale=:log10, xticks = (nelement_sequence, nelement_sequence), yticks = (ytick_sequence, ytick_sequence), markersize = 5, linewidth=2, - xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize) + xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize, + foreground_color_legend = nothing, background_color_legend = nothing, legend=:bottomleft) outfile = outprefix*"_fields_and_ion_pdf_no_Er.pdf" savefig(outfile) println(outfile) try - plot(nelement_sequence, [ion_density_error_sequence,phi_error_sequence,Ez_error_sequence,Er_error_sequence,ion_pdf_error_sequence], xlabel=xlabel, - label=[ylabel_ion_density ylabel_phi ylabel_Ez ylabel_Er ylabel_ion_pdf], ylabel="", + plot(nelement_sequence, [ion_density_error_sequence,phi_error_sequence,Ez_error_sequence,Er_error_sequence,ion_pdf_error_sequence,expected_scaling], xlabel=xlabel, + label=[ylabel_ion_density ylabel_phi ylabel_Ez ylabel_Er ylabel_ion_pdf expected_label], ylabel="", shape =:circle, xscale=:log10, yscale=:log10, xticks = (nelement_sequence, nelement_sequence), yticks = (ytick_sequence, ytick_sequence), markersize = 5, linewidth=2, - xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize) + xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize, + foreground_color_legend = nothing, background_color_legend = nothing, legend=:bottomleft) outfile = outprefix*"_fields_and_ion_pdf.pdf" savefig(outfile) println(outfile) @@ -429,7 +470,8 @@ function get_MMS_error_data(path_list,scan_type,scan_name) plot(nelement_sequence, [ion_density_error_sequence, ion_pdf_error_sequence, neutral_density_error_sequence, neutral_pdf_error_sequence], xlabel=xlabel, label=[ylabel_ion_density ylabel_ion_pdf ylabel_neutral_density ylabel_neutral_pdf], ylabel="", shape =:circle, xscale=:log10, yscale=:log10, xticks = (nelement_sequence, nelement_sequence), yticks = (ytick_sequence, ytick_sequence), markersize = 5, linewidth=2, - xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize) + xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize, + foreground_color_legend = nothing, background_color_legend = nothing, legend=:bottomleft) outfile = outprefix*".pdf" savefig(outfile) println(outfile) @@ -437,7 +479,8 @@ function get_MMS_error_data(path_list,scan_type,scan_name) plot(nelement_sequence, [ion_density_error_sequence, neutral_density_error_sequence, phi_error_sequence, Er_error_sequence, Ez_error_sequence], xlabel=xlabel, label=[ylabel_ion_density ylabel_neutral_density ylabel_phi ylabel_Er ylabel_Ez], ylabel="", shape =:circle, xscale=:log10, yscale=:log10, xticks = (nelement_sequence, nelement_sequence), yticks = (ytick_sequence, ytick_sequence), markersize = 5, linewidth=2, - xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize) + xtickfontsize = fontsize, xguidefontsize = fontsize, ytickfontsize = fontsize, yguidefontsize = fontsize, legendfontsize = fontsize, + foreground_color_legend = nothing, background_color_legend = nothing, legend=:bottomleft) outfile = outprefix*"_fields.pdf" savefig(outfile) println(outfile) @@ -465,4 +508,213 @@ function get_MMS_error_data(path_list,scan_type,scan_name) end +function run_mms_test() + #test_option = "collisional_soundwaves" + #test_option = "collisionless_soundwaves_ions_only" + #test_option = "collisionless_soundwaves" + #test_option = "collisionless_wall-1D-3V-new-dfni-Er" + #test_option = "collisionless_wall-2D-1V-Er-nonzero-at-plate" + #test_option = "collisionless_wall-2D-1V-Er-zero-at-plate" + #test_option = "collisionless_wall-1D-1V" + #test_option = "collisionless_wall-1D-1V-constant-Er" + #test_option = "collisionless_wall-1D-1V-constant-Er-zngrid-5" + #test_option = "collisionless_wall-1D-1V-constant-Er-ngrid-5" + #test_option = "collisionless_wall-1D-1V-constant-Er-ngrid-5-opt" + #test_option = "krook_wall-1D-2V" + test_option = "mirror_wall-1D-2V" + #test_option = "mirror_wall-1D-2V-ngrid-5" + #test_option = "mirror_wall-2D-2V-ngrid-5" + #test_option = "collisionless_wall-1D-3V" + #test_option = "collisionless_wall-2D-3V" + #test_option = "collisionless_wall-2D-3V-Er-zero-at-plate" + #test_option = "collisionless_biased_wall-2D-3V" + #test_option = "collisionless_wall-1D-3V-with-sheath" + + if test_option == "collisional_soundwaves" + # collisional + path_list = ["runs/2D-sound-wave_cheb_cxiz_nel_r_2_z_2_vpa_2_vperp_2","runs/2D-sound-wave_cheb_cxiz_nel_r_2_z_2_vpa_4_vperp_4", + "runs/2D-sound-wave_cheb_cxiz_nel_r_2_z_2_vpa_8_vperp_8","runs/2D-sound-wave_cheb_cxiz_nel_r_2_z_2_vpa_16_vperp_16"] + scan_type = "velocity_nelement" + scan_name = "2D-sound-wave_cheb_cxiz" + + elseif test_option == "collisionless_soundwaves" + # collisionless + path_list = ["runs/2D-sound-wave_cheb_nel_r_2_z_2_vpa_2_vperp_2","runs/2D-sound-wave_cheb_nel_r_4_z_4_vpa_4_vperp_4",# "runs/2D-sound-wave_cheb_nel_r_6_z_6_vpa_6_vperp_6", + "runs/2D-sound-wave_cheb_nel_r_8_z_8_vpa_8_vperp_8","runs/2D-sound-wave_cheb_nel_r_16_z_16_vpa_16_vperp_16" + ] + scan_type = "nelement" + scan_name = "2D-sound-wave_cheb" + elseif test_option == "collisionless_soundwaves_ions_only" + # collisionless + path_list = ["runs/2D-sound-wave_cheb_ion_only_nel_r_2_z_2_vpa_2_vperp_2","runs/2D-sound-wave_cheb_ion_only_nel_r_4_z_4_vpa_4_vperp_4", "runs/2D-sound-wave_cheb_ion_only_nel_r_8_z_8_vpa_8_vperp_8", +# "runs/2D-sound-wave_cheb_nel_r_8_z_8_vpa_8_vperp_8" +#,"runs/2D-sound-wave_cheb_nel_r_2_z_2_vpa_16_vperp_16" + ] + scan_type = "nelement" + scan_name = "2D-sound-wave_cheb_ions_only" + elseif test_option == "collisionless_wall-2D-3V-Er-zero-at-plate" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_2_vperp_2", + "runs/2D-wall_cheb-with-neutrals_nel_r_4_z_4_vpa_4_vperp_4", + "runs/2D-wall_cheb-with-neutrals_nel_r_6_z_6_vpa_6_vperp_6", + #"runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_16_vperp_16" + ] + scan_type = "velocity_nelement" + scan_name = "2D-3V-wall_cheb" + elseif test_option == "collisionless_wall-2D-3V" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_4_vperp_4", + "runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_8_vperp_8","runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_12_vperp_12", + "runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_16_vperp_16" ] + scan_type = "velocity_nelement" + scan_name = "2D-3V-wall_cheb" + elseif test_option == "collisionless_biased_wall-2D-3V" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/2D-wall-Dirichlet_nel_r_2_z_2_vpa_2","runs/2D-wall-Dirichlet_nel_r_2_z_2_vpa_4", + "runs/2D-wall-Dirichlet_nel_r_2_z_2_vpa_8","runs/2D-wall-Dirichlet_nel_r_2_z_2_vpa_16"] + scan_type = "velocity_nelement" + scan_name = "2D-3V-biased_wall_cheb" + + elseif test_option == "collisionless_wall-1D-3V" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_2_vperp_2","runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_4_vperp_4", + "runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_8_vperp_8",#"runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_12_vperp_12", + "runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_16_vperp_16"#,"runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_24_vperp_24" + ] + scan_type = "velocity_nelement" + scan_name = "1D-3V-wall_cheb" + + elseif test_option == "collisionless_wall-1D-3V-updated" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_2_vperp_2","runs/2D-wall_cheb-with-neutrals_nel_r_1_z_4_vpa_4_vperp_4", + "runs/2D-wall_cheb-with-neutrals_nel_r_1_z_8_vpa_8_vperp_8",#"runs/2D-wall_cheb-with-neutrals_nel_r_1_z_12_vpa_12_vperp_12", + # "runs/2D-wall_cheb-with-neutrals_nel_r_1_z_16_vpa_16_vperp_16"#,"runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_24_vperp_24" + ] + scan_type = "nelement" + scan_name = "1D-3V-wall_cheb-updated" + elseif test_option == "collisionless_wall-1D-3V-new-dfni" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/2D-wall_cheb-new-MMS-dfni-ngrid-9_nel_r_1_z_2_vpa_2_vperp_2", + "runs/2D-wall_cheb-new-MMS-dfni-ngrid-9_nel_r_1_z_4_vpa_4_vperp_4", + "runs/2D-wall_cheb-new-MMS-dfni-ngrid-9_nel_r_1_z_8_vpa_8_vperp_8", + "runs/2D-wall_cheb-new-MMS-dfni-ngrid-9_nel_r_1_z_16_vpa_16_vperp_16", + "runs/2D-wall_cheb-new-MMS-dfni-ngrid-9_nel_r_1_z_32_vpa_32_vperp_32" + ] + scan_type = "nelement" + scan_name = "1D-3V-wall_cheb-new-dfni" + elseif test_option == "collisionless_wall-1D-3V-new-dfni-Er" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/2D-wall_cheb-new-MMS-dfni-Er-ngrid-9_nel_r_1_z_2_vpa_2_vperp_2", + "runs/2D-wall_cheb-new-MMS-dfni-Er-ngrid-9_nel_r_1_z_4_vpa_4_vperp_4", + "runs/2D-wall_cheb-new-MMS-dfni-Er-ngrid-9_nel_r_1_z_8_vpa_8_vperp_8", + "runs/2D-wall_cheb-new-MMS-dfni-Er-ngrid-9_nel_r_1_z_16_vpa_16_vperp_16", + "runs/2D-wall_cheb-new-MMS-dfni-Er-ngrid-9_nel_r_1_z_32_vpa_32_vperp_32" + ] + scan_type = "nelement" + scan_name = "1D-3V-wall_cheb-new-dfni-Er" + elseif test_option == "collisionless_wall-1D-3V-with-sheath" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_2_vperp_2","runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_4_vperp_4", + "runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_8_vperp_8",#"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_12_vperp_12", + "runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_16_vperp_16"#,"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_24_vperp_24" + ] + scan_type = "velocity_nelement" + scan_name = "1D-3V-wall-sheath_cheb" + elseif test_option == "collisionless_wall-2D-1V-Er-zero-at-plate" + #path_list = ["runs/2D-wall_MMS_nel_r_2_z_2_vpa_16_vperp_1_diss","runs/2D-wall_MMS_nel_r_4_z_4_vpa_16_vperp_1_diss", + # "runs/2D-wall_MMS_nel_r_8_z_8_vpa_16_vperp_1_diss",#"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_12_vperp_12", + # "runs/2D-wall_MMS_nel_r_16_z_16_vpa_16_vperp_1_diss", + # "runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss5"#,"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_24_vperp_24" + # ] + #scan_type = "zr_nelement" + path_list = ["runs/2D-wall_MMS_ngrid_5_nel_r_2_z_2_vpa_8_vperp_1_diss","runs/2D-wall_MMS_ngrid_5_nel_r_4_z_4_vpa_16_vperp_1_diss", + "runs/2D-wall_MMS_ngrid_5_nel_r_8_z_8_vpa_32_vperp_1_diss","runs/2D-wall_MMS_ngrid_5_nel_r_16_z_16_vpa_64_vperp_1_diss"] + scan_type = "vpazr_nelement0.25" + scan_name = "2D-1V-wall_cheb" + elseif test_option == "collisionless_wall-2D-1V-Er-nonzero-at-plate" + path_list = ["runs/2D-wall_MMSEr_ngrid_5_nel_r_2_z_2_vpa_8_vperp_1_diss","runs/2D-wall_MMSEr_ngrid_5_nel_r_4_z_4_vpa_16_vperp_1_diss","runs/2D-wall_MMSEr_ngrid_5_nel_r_8_z_8_vpa_32_vperp_1_diss","runs/2D-wall_MMSEr_ngrid_5_nel_r_16_z_16_vpa_64_vperp_1_diss", + ] + #path_list = ["runs/2D-wall_MMSEr_nel_r_2_z_2_vpa_2_vperp_1_diss","runs/2D-wall_MMSEr_nel_r_4_z_4_vpa_4_vperp_1_diss", + # "runs/2D-wall_MMSEr_nel_r_8_z_8_vpa_8_vperp_1_diss","runs/2D-wall_MMSEr_nel_r_16_z_16_vpa_16_vperp_1_diss" + # #"runs/2D-wall_MMSEr_nel_r_32_z_32_vpa_16_vperp_1_diss" + # ] + #path_list = ["runs/2D-wall_MMSEr_nel_r_2_z_2_vpa_16_vperp_1_diss","runs/2D-wall_MMSEr_nel_r_4_z_4_vpa_16_vperp_1_diss", + # "runs/2D-wall_MMSEr_nel_r_8_z_8_vpa_16_vperp_1_diss",#"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_12_vperp_12", + # "runs/2D-wall_MMSEr_nel_r_16_z_16_vpa_16_vperp_1_diss", + # "runs/2D-wall_MMSEr_nel_r_32_z_32_vpa_16_vperp_1_diss"#,"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_24_vperp_24" + # ] + scan_type = "vpazr_nelement0.25" + scan_name = "2D-1V-wall_cheb-nonzero-Er" + elseif test_option == "collisionless_wall-1D-1V-constant-Er-ngrid-5-opt" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_8_vpa_32_vperp_1","runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_16_vpa_64_vperp_1", + "runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_32_vpa_128_vperp_1","runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_64_vpa_256_vperp_1" + ] + scan_type = "vpaz_nelement0.25" + scan_name = "1D-1V-wall_cheb-constant-Er-ngrid-5-opt" + elseif test_option == "collisionless_wall-1D-1V-constant-Er-ngrid-5" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_8_vpa_8_vperp_1","runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_16_vpa_16_vperp_1", + "runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_32_vpa_32_vperp_1","runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_64_vpa_64_vperp_1" + ] + scan_type = "vpaz_nelement" + scan_name = "1D-1V-wall_cheb-constant-Er-ngrid-5" + elseif test_option == "collisionless_wall-1D-1V-constant-Er-zngrid-5" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/1D-wall_MMSEr_zngrid_5_nel_r_1_z_8_vpa_2_vperp_1","runs/1D-wall_MMSEr_zngrid_5_nel_r_1_z_16_vpa_4_vperp_1", + "runs/1D-wall_MMSEr_zngrid_5_nel_r_1_z_32_vpa_8_vperp_1","runs/1D-wall_MMSEr_zngrid_5_nel_r_1_z_64_vpa_16_vperp_1" + ] + scan_type = "vpaz_nelement4" + scan_name = "1D-1V-wall_cheb-constant-Er-zngrid-5" + elseif test_option == "collisionless_wall-1D-1V-constant-Er" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/1D-wall_MMSEr_nel_r_1_z_2_vpa_2_vperp_1","runs/1D-wall_MMSEr_nel_r_1_z_4_vpa_4_vperp_1", + "runs/1D-wall_MMSEr_nel_r_1_z_8_vpa_8_vperp_1","runs/1D-wall_MMSEr_nel_r_1_z_16_vpa_16_vperp_1" + ] + scan_type = "vpaz_nelement" + scan_name = "1D-1V-wall_cheb-constant-Er" + elseif test_option == "collisionless_wall-1D-1V" + # collisionless wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/1D-wall_MMS_nel_r_1_z_2_vpa_2_vperp_1","runs/1D-wall_MMS_nel_r_1_z_4_vpa_4_vperp_1", + "runs/1D-wall_MMS_nel_r_1_z_8_vpa_8_vperp_1","runs/1D-wall_MMS_nel_r_1_z_16_vpa_16_vperp_1" + ] + scan_type = "vpaz_nelement" + scan_name = "1D-1V-wall_cheb" + elseif test_option == "krook_wall-1D-2V" + # Krook wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/1D-wall_MMS_nel_r_1_z_2_vpa_2_vperp_2_krook","runs/1D-wall_MMS_nel_r_1_z_4_vpa_4_vperp_4_krook", + "runs/1D-wall_MMS_nel_r_1_z_8_vpa_8_vperp_8_krook" ] + scan_type = "vpavperpz_nelement" + scan_name = "1D-2V-wall_cheb_krook" + elseif test_option == "mirror_wall-1D-2V" + # Mirror wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss", + "runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss", + "runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss", + "runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss",] + scan_type = "vpa2vperpz_nelement" + scan_name = "mirror_wall-1D-2V" + elseif test_option == "mirror_wall-1D-2V-ngrid-5" + # Mirror wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/1D-mirror_MMS_ngrid_5_nel_r_1_z_4_vpa_4_vperp_2_diss", + "runs/1D-mirror_MMS_ngrid_5_nel_r_1_z_8_vpa_8_vperp_4_diss", + "runs/1D-mirror_MMS_ngrid_5_nel_r_1_z_16_vpa_16_vperp_8_diss", + "runs/1D-mirror_MMS_ngrid_5_nel_r_1_z_32_vpa_32_vperp_16_diss", + "runs/1D-mirror_MMS_ngrid_5_nel_r_1_z_64_vpa_64_vperp_32_diss",] + scan_type = "vpa2vperpz_nelement" + scan_name = "mirror_wall-1D-2V-ngrid-5" + elseif test_option == "mirror_wall-2D-2V-ngrid-5" + # Mirror wall test, no sheath for electrons, no radial coordinate + path_list = ["runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss", + "runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss", + "runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss", + "runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss",] + scan_type = "vpa2vperpzr_nelement" + scan_name = "mirror_wall-2D-2V-ngrid-5" + end + print(path_list) + get_MMS_error_data(path_list,scan_type,scan_name) + return nothing +end + end diff --git a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl index 149355eee..f3e4aceb1 100644 --- a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl +++ b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl @@ -57,13 +57,16 @@ using moment_kinetics.analysis: analyze_fields_data, analyze_moments_data, get_unnormalised_f_coords_2d using moment_kinetics.velocity_moments: integrate_over_vspace using moment_kinetics.manufactured_solns: manufactured_solutions, - manufactured_electric_fields -using moment_kinetics.moment_kinetics_input: mk_input, get + manufactured_electric_fields, + manufactured_geometry +using moment_kinetics.moment_kinetics_input: mk_input, get, get_default_rhostar using moment_kinetics.input_structs: geometry_input, grid_input, species_composition using moment_kinetics.input_structs: electron_physics_type, boltzmann_electron_response, boltzmann_electron_response_with_simple_sheath +using moment_kinetics.reference_parameters +using moment_kinetics.geo: init_magnetic_geometry using .post_processing_input: pp -using .shared_utils: calculate_and_write_frequencies, get_geometry_and_composition +using .shared_utils: calculate_and_write_frequencies, get_geometry, get_composition using TOML import Base: get @@ -637,9 +640,10 @@ function analyze_and_plot_data(prefix...; run_index=nothing) end end - geometry, composition = - get_tuple_of_return_values(get_geometry_and_composition, scan_input, - n_ion_species, n_neutral_species) + geometry = + get_tuple_of_return_values(get_geometry, scan_input, z, r) + composition = + get_tuple_of_return_values(get_composition, scan_input) # initialise the post-processing input options nwrite_movie, itime_min, itime_max, nwrite_movie_pdfs, itime_min_pdfs, itime_max_pdfs, @@ -1071,6 +1075,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) manufactured_solns_test = manufactured_solns_input.use_for_advance && manufactured_solns_input.use_for_init # Plots compare density and density_symbolic at last timestep #if(manufactured_solns_test && nr > 1) + println("manufactured_solns_test: ",manufactured_solns_test) if(manufactured_solns_test) # avoid passing Lr = 0 into manufactured_solns functions if r_global.n > 1 @@ -1078,10 +1083,10 @@ function analyze_and_plot_data(prefix...; run_index=nothing) else Lr_in = 1.0 end - + check_manufactured_solution_geometry(geometry,z_global,r_global) manufactured_solns_list = manufactured_solutions(manufactured_solns_input, Lr_in, z_global.L, r_global.bc, - z_global.bc, geometry, + z_global.bc, geometry.input, composition, species, r_global.n, vperp.n) dfni_func = manufactured_solns_list.dfni_func densi_func = manufactured_solns_list.densi_func @@ -3467,6 +3472,41 @@ function plot_charged_pdf_2D_at_wall(run_name, run_name_label, r_global, z_globa println("done.") end +function check_manufactured_solution_geometry(geometry,z,r) + Lz = z.L + nz = z.n + Lr = r.L + nr = r.n + + Bmag_num = geometry.Bmag + dBdz_num = geometry.dBdz + bzed_num = geometry.bzed + Bmag_sym = copy(Bmag_num) + bzed_sym = copy(bzed_num) + dBdz_sym = copy(dBdz_num) + + manufactured_geometry_list = manufactured_geometry(geometry.input,Lz,Lr,nr) + Bmag_func = manufactured_geometry_list.Bmag_func + bzed_func = manufactured_geometry_list.bzed_func + dBdz_func = manufactured_geometry_list.dBdz_func + for ir in 1:nr + for iz in 1:nz + Bmag_sym[iz,ir] = Bmag_func(z.grid[iz],r.grid[ir]) + bzed_sym[iz,ir] = bzed_func(z.grid[iz],r.grid[ir]) + dBdz_sym[iz,ir] = dBdz_func(z.grid[iz],r.grid[ir]) + end + end + println("Geometry check") + Bmag_err = sqrt(sum((Bmag_sym .- Bmag_num).^2)) + println("Bmag_err: ",Bmag_err) + bzed_err = sqrt(sum((bzed_sym .- bzed_num).^2)) + println("bzed_err: ",bzed_err) + dBdz_err = sqrt(sum((dBdz_sym .- dBdz_num).^2)) + println("dBdz_err: ",dBdz_err) + + return nothing +end + include("plot_MMS_sequence.jl") include("plot_sequence.jl") diff --git a/run_MMS_test.jl b/run_MMS_test.jl index 50f8f9a14..206b0df45 100644 --- a/run_MMS_test.jl +++ b/run_MMS_test.jl @@ -1,183 +1,7 @@ if abspath(PROGRAM_FILE) == @__FILE__ using Pkg Pkg.activate(".") - - import moment_kinetics as mk - - #test_option = "collisional_soundwaves" - #test_option = "collisionless_soundwaves_ions_only" - #test_option = "collisionless_soundwaves" - #test_option = "collisionless_wall-1D-3V-new-dfni-Er" - #test_option = "collisionless_wall-2D-1V-Er-nonzero-at-plate" - #test_option = "collisionless_wall-2D-1V-Er-zero-at-plate" - #test_option = "collisionless_wall-1D-1V" - #test_option = "collisionless_wall-1D-1V-constant-Er" - #test_option = "collisionless_wall-1D-1V-constant-Er-zngrid-5" - #test_option = "collisionless_wall-1D-1V-constant-Er-ngrid-5" - #test_option = "collisionless_wall-1D-1V-constant-Er-ngrid-5-opt" - test_option = "krook_wall-1D-2V" - #test_option = "collisionless_wall-1D-3V" - #test_option = "collisionless_wall-2D-3V" - #test_option = "collisionless_wall-2D-3V-Er-zero-at-plate" - #test_option = "collisionless_biased_wall-2D-3V" - #test_option = "collisionless_wall-1D-3V-with-sheath" - - if test_option == "collisional_soundwaves" - # collisional - path_list = ["runs/2D-sound-wave_cheb_cxiz_nel_r_2_z_2_vpa_2_vperp_2","runs/2D-sound-wave_cheb_cxiz_nel_r_2_z_2_vpa_4_vperp_4", - "runs/2D-sound-wave_cheb_cxiz_nel_r_2_z_2_vpa_8_vperp_8","runs/2D-sound-wave_cheb_cxiz_nel_r_2_z_2_vpa_16_vperp_16"] - scan_type = "velocity_nelement" - scan_name = "2D-sound-wave_cheb_cxiz" - - elseif test_option == "collisionless_soundwaves" - # collisionless - path_list = ["runs/2D-sound-wave_cheb_nel_r_2_z_2_vpa_2_vperp_2","runs/2D-sound-wave_cheb_nel_r_4_z_4_vpa_4_vperp_4",# "runs/2D-sound-wave_cheb_nel_r_6_z_6_vpa_6_vperp_6", - "runs/2D-sound-wave_cheb_nel_r_8_z_8_vpa_8_vperp_8","runs/2D-sound-wave_cheb_nel_r_16_z_16_vpa_16_vperp_16" - ] - scan_type = "nelement" - scan_name = "2D-sound-wave_cheb" - elseif test_option == "collisionless_soundwaves_ions_only" - # collisionless - path_list = ["runs/2D-sound-wave_cheb_ion_only_nel_r_2_z_2_vpa_2_vperp_2","runs/2D-sound-wave_cheb_ion_only_nel_r_4_z_4_vpa_4_vperp_4", "runs/2D-sound-wave_cheb_ion_only_nel_r_8_z_8_vpa_8_vperp_8", -# "runs/2D-sound-wave_cheb_nel_r_8_z_8_vpa_8_vperp_8" -#,"runs/2D-sound-wave_cheb_nel_r_2_z_2_vpa_16_vperp_16" - ] - scan_type = "nelement" - scan_name = "2D-sound-wave_cheb_ions_only" - elseif test_option == "collisionless_wall-2D-3V-Er-zero-at-plate" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_2_vperp_2", - "runs/2D-wall_cheb-with-neutrals_nel_r_4_z_4_vpa_4_vperp_4", - "runs/2D-wall_cheb-with-neutrals_nel_r_6_z_6_vpa_6_vperp_6", - #"runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_16_vperp_16" - ] - scan_type = "velocity_nelement" - scan_name = "2D-3V-wall_cheb" - elseif test_option == "collisionless_wall-2D-3V" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_4_vperp_4", - "runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_8_vperp_8","runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_12_vperp_12", - "runs/2D-wall_cheb-with-neutrals_nel_r_2_z_2_vpa_16_vperp_16" ] - scan_type = "velocity_nelement" - scan_name = "2D-3V-wall_cheb" - elseif test_option == "collisionless_biased_wall-2D-3V" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/2D-wall-Dirichlet_nel_r_2_z_2_vpa_2","runs/2D-wall-Dirichlet_nel_r_2_z_2_vpa_4", - "runs/2D-wall-Dirichlet_nel_r_2_z_2_vpa_8","runs/2D-wall-Dirichlet_nel_r_2_z_2_vpa_16"] - scan_type = "velocity_nelement" - scan_name = "2D-3V-biased_wall_cheb" - - elseif test_option == "collisionless_wall-1D-3V" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_2_vperp_2","runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_4_vperp_4", - "runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_8_vperp_8",#"runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_12_vperp_12", - "runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_16_vperp_16"#,"runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_24_vperp_24" - ] - scan_type = "velocity_nelement" - scan_name = "1D-3V-wall_cheb" - - elseif test_option == "collisionless_wall-1D-3V-updated" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_2_vperp_2","runs/2D-wall_cheb-with-neutrals_nel_r_1_z_4_vpa_4_vperp_4", - "runs/2D-wall_cheb-with-neutrals_nel_r_1_z_8_vpa_8_vperp_8",#"runs/2D-wall_cheb-with-neutrals_nel_r_1_z_12_vpa_12_vperp_12", - # "runs/2D-wall_cheb-with-neutrals_nel_r_1_z_16_vpa_16_vperp_16"#,"runs/2D-wall_cheb-with-neutrals_nel_r_1_z_2_vpa_24_vperp_24" - ] - scan_type = "nelement" - scan_name = "1D-3V-wall_cheb-updated" - elseif test_option == "collisionless_wall-1D-3V-new-dfni" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/2D-wall_cheb-new-MMS-dfni-ngrid-9_nel_r_1_z_2_vpa_2_vperp_2", - "runs/2D-wall_cheb-new-MMS-dfni-ngrid-9_nel_r_1_z_4_vpa_4_vperp_4", - "runs/2D-wall_cheb-new-MMS-dfni-ngrid-9_nel_r_1_z_8_vpa_8_vperp_8", - "runs/2D-wall_cheb-new-MMS-dfni-ngrid-9_nel_r_1_z_16_vpa_16_vperp_16", - "runs/2D-wall_cheb-new-MMS-dfni-ngrid-9_nel_r_1_z_32_vpa_32_vperp_32" - ] - scan_type = "nelement" - scan_name = "1D-3V-wall_cheb-new-dfni" - elseif test_option == "collisionless_wall-1D-3V-new-dfni-Er" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/2D-wall_cheb-new-MMS-dfni-Er-ngrid-9_nel_r_1_z_2_vpa_2_vperp_2", - "runs/2D-wall_cheb-new-MMS-dfni-Er-ngrid-9_nel_r_1_z_4_vpa_4_vperp_4", - "runs/2D-wall_cheb-new-MMS-dfni-Er-ngrid-9_nel_r_1_z_8_vpa_8_vperp_8", - "runs/2D-wall_cheb-new-MMS-dfni-Er-ngrid-9_nel_r_1_z_16_vpa_16_vperp_16", - "runs/2D-wall_cheb-new-MMS-dfni-Er-ngrid-9_nel_r_1_z_32_vpa_32_vperp_32" - ] - scan_type = "nelement" - scan_name = "1D-3V-wall_cheb-new-dfni-Er" - elseif test_option == "collisionless_wall-1D-3V-with-sheath" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_2_vperp_2","runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_4_vperp_4", - "runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_8_vperp_8",#"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_12_vperp_12", - "runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_16_vperp_16"#,"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_24_vperp_24" - ] - scan_type = "velocity_nelement" - scan_name = "1D-3V-wall-sheath_cheb" - elseif test_option == "collisionless_wall-2D-1V-Er-zero-at-plate" - #path_list = ["runs/2D-wall_MMS_nel_r_2_z_2_vpa_16_vperp_1_diss","runs/2D-wall_MMS_nel_r_4_z_4_vpa_16_vperp_1_diss", - # "runs/2D-wall_MMS_nel_r_8_z_8_vpa_16_vperp_1_diss",#"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_12_vperp_12", - # "runs/2D-wall_MMS_nel_r_16_z_16_vpa_16_vperp_1_diss", - # "runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss5"#,"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_24_vperp_24" - # ] - #scan_type = "zr_nelement" - path_list = ["runs/2D-wall_MMS_ngrid_5_nel_r_2_z_2_vpa_8_vperp_1_diss","runs/2D-wall_MMS_ngrid_5_nel_r_4_z_4_vpa_16_vperp_1_diss", - "runs/2D-wall_MMS_ngrid_5_nel_r_8_z_8_vpa_32_vperp_1_diss","runs/2D-wall_MMS_ngrid_5_nel_r_16_z_16_vpa_64_vperp_1_diss"] - scan_type = "vpazr_nelement0.25" - scan_name = "2D-1V-wall_cheb" - elseif test_option == "collisionless_wall-2D-1V-Er-nonzero-at-plate" - path_list = ["runs/2D-wall_MMSEr_ngrid_5_nel_r_2_z_2_vpa_8_vperp_1_diss","runs/2D-wall_MMSEr_ngrid_5_nel_r_4_z_4_vpa_16_vperp_1_diss","runs/2D-wall_MMSEr_ngrid_5_nel_r_8_z_8_vpa_32_vperp_1_diss","runs/2D-wall_MMSEr_ngrid_5_nel_r_16_z_16_vpa_64_vperp_1_diss", - ] - #path_list = ["runs/2D-wall_MMSEr_nel_r_2_z_2_vpa_2_vperp_1_diss","runs/2D-wall_MMSEr_nel_r_4_z_4_vpa_4_vperp_1_diss", - # "runs/2D-wall_MMSEr_nel_r_8_z_8_vpa_8_vperp_1_diss","runs/2D-wall_MMSEr_nel_r_16_z_16_vpa_16_vperp_1_diss" - # #"runs/2D-wall_MMSEr_nel_r_32_z_32_vpa_16_vperp_1_diss" - # ] - #path_list = ["runs/2D-wall_MMSEr_nel_r_2_z_2_vpa_16_vperp_1_diss","runs/2D-wall_MMSEr_nel_r_4_z_4_vpa_16_vperp_1_diss", - # "runs/2D-wall_MMSEr_nel_r_8_z_8_vpa_16_vperp_1_diss",#"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_12_vperp_12", - # "runs/2D-wall_MMSEr_nel_r_16_z_16_vpa_16_vperp_1_diss", - # "runs/2D-wall_MMSEr_nel_r_32_z_32_vpa_16_vperp_1_diss"#,"runs/2D-wall_cheb-with-neutrals-with-sheath_nel_r_1_z_2_vpa_24_vperp_24" - # ] - scan_type = "vpazr_nelement0.25" - scan_name = "2D-1V-wall_cheb-nonzero-Er" - elseif test_option == "collisionless_wall-1D-1V-constant-Er-ngrid-5-opt" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_8_vpa_32_vperp_1","runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_16_vpa_64_vperp_1", - "runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_32_vpa_128_vperp_1","runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_64_vpa_256_vperp_1" - ] - scan_type = "vpaz_nelement0.25" - scan_name = "1D-1V-wall_cheb-constant-Er-ngrid-5-opt" - elseif test_option == "collisionless_wall-1D-1V-constant-Er-ngrid-5" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_8_vpa_8_vperp_1","runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_16_vpa_16_vperp_1", - "runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_32_vpa_32_vperp_1","runs/1D-wall_MMSEr_ngrid_5_nel_r_1_z_64_vpa_64_vperp_1" - ] - scan_type = "vpaz_nelement" - scan_name = "1D-1V-wall_cheb-constant-Er-ngrid-5" - elseif test_option == "collisionless_wall-1D-1V-constant-Er-zngrid-5" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/1D-wall_MMSEr_zngrid_5_nel_r_1_z_8_vpa_2_vperp_1","runs/1D-wall_MMSEr_zngrid_5_nel_r_1_z_16_vpa_4_vperp_1", - "runs/1D-wall_MMSEr_zngrid_5_nel_r_1_z_32_vpa_8_vperp_1","runs/1D-wall_MMSEr_zngrid_5_nel_r_1_z_64_vpa_16_vperp_1" - ] - scan_type = "vpaz_nelement4" - scan_name = "1D-1V-wall_cheb-constant-Er-zngrid-5" - elseif test_option == "collisionless_wall-1D-1V-constant-Er" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/1D-wall_MMSEr_nel_r_1_z_2_vpa_2_vperp_1","runs/1D-wall_MMSEr_nel_r_1_z_4_vpa_4_vperp_1", - "runs/1D-wall_MMSEr_nel_r_1_z_8_vpa_8_vperp_1","runs/1D-wall_MMSEr_nel_r_1_z_16_vpa_16_vperp_1" - ] - scan_type = "vpaz_nelement" - scan_name = "1D-1V-wall_cheb-constant-Er" - elseif test_option == "collisionless_wall-1D-1V" - # collisionless wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/1D-wall_MMS_nel_r_1_z_2_vpa_2_vperp_1","runs/1D-wall_MMS_nel_r_1_z_4_vpa_4_vperp_1", - "runs/1D-wall_MMS_nel_r_1_z_8_vpa_8_vperp_1","runs/1D-wall_MMS_nel_r_1_z_16_vpa_16_vperp_1" - ] - scan_type = "vpaz_nelement" - scan_name = "1D-1V-wall_cheb" - elseif test_option == "krook_wall-1D-2V" - # Krook wall test, no sheath for electrons, no radial coordinate - path_list = ["runs/1D-wall_MMS_nel_r_1_z_2_vpa_2_vperp_2_krook","runs/1D-wall_MMS_nel_r_1_z_4_vpa_4_vperp_4_krook", - "runs/1D-wall_MMS_nel_r_1_z_8_vpa_8_vperp_8_krook" ] - scan_type = "vpavperpz_nelement" - scan_name = "1D-2V-wall_cheb_krook" - end - mk.plot_MMS_sequence.get_MMS_error_data(path_list,scan_type,scan_name) + import plots_post_processing + using plots_post_processing.plot_MMS_sequence + run_mms_test() end diff --git a/runs/1D-mirror_MMS_new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml b/runs/1D-mirror_MMS_new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml new file mode 100644 index 000000000..45d596904 --- /dev/null +++ b/runs/1D-mirror_MMS_new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml @@ -0,0 +1,98 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +geometry_option="1D-mirror" +DeltaB=10.0 +#geometry_option="constant-helical" +#pitch=1.0 +Er_constant = 0.0 +T_e = 1.0 +T_wall = 1.0 +rhostar = 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 +nstep = 2000 +dt = 0.0005 +nwrite = 200 +nwrite_dfns = 200 +use_semi_lagrange = false +n_rk_stages = 4 +split_operators = false +z_ngrid = 17 +z_nelement = 16 +z_nelement_local = 16 +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 = 17 +vpa_nelement = 16 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 17 +vperp_nelement = 8 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +vperp_discretization = "chebyshev_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.001 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.0 diff --git a/runs/1D-mirror_MMS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml b/runs/1D-mirror_MMS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml new file mode 100644 index 000000000..a465d248c --- /dev/null +++ b/runs/1D-mirror_MMS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml @@ -0,0 +1,99 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +geometry_option="1D-mirror" +DeltaB=0.5 +#geometry_option="constant-helical" +#pitch=1.0 +Er_constant = 0.0 +T_e = 1.0 +T_wall = 1.0 +rhostar = 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 +nstep = 2000 +dt = 0.0005 +nwrite = 200 +nwrite_dfns = 200 +use_semi_lagrange = false +n_rk_stages = 4 +split_operators = false +z_ngrid = 17 +z_nelement = 6 +z_nelement_local = 6 +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 = 17 +vpa_nelement = 6 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 17 +vperp_nelement = 3 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +#vperp_discretization = "chebyshev_pseudospectral" +vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.001 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.0 diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml new file mode 100644 index 000000000..6a89faab1 --- /dev/null +++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml @@ -0,0 +1,101 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +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 +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 +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 = 9 +vpa_nelement = 16 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 9 +vperp_nelement = 8 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +vperp_discretization = "chebyshev_pseudospectral" +#vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.001 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.0 + +[geometry] +option="1D-mirror" +DeltaB=0.5 +#option="constant-helical" +pitch=1.0 +rhostar = 1.0 diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml new file mode 100644 index 000000000..a90350f22 --- /dev/null +++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml @@ -0,0 +1,101 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +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 +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 +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 = 9 +vpa_nelement = 32 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 9 +vperp_nelement = 16 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +vperp_discretization = "chebyshev_pseudospectral" +#vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.001 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.0 + +[geometry] +option="1D-mirror" +DeltaB=0.5 +#option="constant-helical" +pitch=1.0 +rhostar = 1.0 diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml new file mode 100644 index 000000000..92a8aa331 --- /dev/null +++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml @@ -0,0 +1,101 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +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 +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 +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 = 9 +vpa_nelement = 4 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 9 +vperp_nelement = 2 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +#vperp_discretization = "chebyshev_pseudospectral" +vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.001 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.0 + +[geometry] +option="1D-mirror" +DeltaB=0.5 +#option="constant-helical" +pitch=1.0 +rhostar = 1.0 diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml new file mode 100644 index 000000000..5f63e0aec --- /dev/null +++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml @@ -0,0 +1,101 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +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 +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 +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 = 9 +vpa_nelement = 64 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 9 +vperp_nelement = 32 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +vperp_discretization = "chebyshev_pseudospectral" +#vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.0 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.0 + +[geometry] +option="1D-mirror" +DeltaB=0.5 +#option="constant-helical" +pitch=1.0 +rhostar = 1.0 diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml new file mode 100644 index 000000000..c8e8df9ee --- /dev/null +++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml @@ -0,0 +1,101 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +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 +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 +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 = 9 +vpa_nelement = 8 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 9 +vperp_nelement = 4 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +vperp_discretization = "chebyshev_pseudospectral" +#vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.001 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.0 + +[geometry] +option="1D-mirror" +DeltaB=0.5 +#option="constant-helical" +pitch=1.0 +rhostar = 1.0 diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml new file mode 100644 index 000000000..a300c9a56 --- /dev/null +++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml @@ -0,0 +1,101 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +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 +nstep = 2000 +dt = 0.0005 +nwrite = 200 +nwrite_dfns = 200 +use_semi_lagrange = false +n_rk_stages = 4 +split_operators = false +z_ngrid = 5 +z_nelement = 16 +z_nelement_local = 16 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 5 +r_nelement = 16 +r_nelement_local = 16 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 5 +vpa_nelement = 16 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 5 +vperp_nelement = 8 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +#vperp_discretization = "chebyshev_pseudospectral" +vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.001 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.01 + +[geometry] +option="1D-mirror" +DeltaB=0.5 +#option="constant-helical" +pitch=0.5 +rhostar = 1.0 diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml new file mode 100644 index 000000000..c98077d3a --- /dev/null +++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml @@ -0,0 +1,101 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +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 +nstep = 2000 +dt = 0.0005 +nwrite = 200 +nwrite_dfns = 200 +use_semi_lagrange = false +n_rk_stages = 4 +split_operators = false +z_ngrid = 5 +z_nelement = 32 +z_nelement_local = 32 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 5 +r_nelement = 32 +r_nelement_local = 32 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 5 +vpa_nelement = 32 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 5 +vperp_nelement = 16 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +#vperp_discretization = "chebyshev_pseudospectral" +vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.001 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.01 + +[geometry] +option="1D-mirror" +DeltaB=0.5 +#option="constant-helical" +pitch=0.5 +rhostar = 1.0 diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml new file mode 100644 index 000000000..3ec12c9fe --- /dev/null +++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml @@ -0,0 +1,101 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +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 +nstep = 2000 +dt = 0.0005 +nwrite = 200 +nwrite_dfns = 200 +use_semi_lagrange = false +n_rk_stages = 4 +split_operators = false +z_ngrid = 5 +z_nelement = 4 +z_nelement_local = 4 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 5 +r_nelement = 4 +r_nelement_local = 4 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 5 +vpa_nelement = 4 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 5 +vperp_nelement = 2 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +#vperp_discretization = "chebyshev_pseudospectral" +vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.001 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.01 + +[geometry] +option="1D-mirror" +DeltaB=0.5 +#option="constant-helical" +pitch=0.5 +rhostar = 1.0 diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml new file mode 100644 index 000000000..e936ef06e --- /dev/null +++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml @@ -0,0 +1,101 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +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 +nstep = 2000 +dt = 0.0005 +nwrite = 200 +nwrite_dfns = 200 +use_semi_lagrange = false +n_rk_stages = 4 +split_operators = false +z_ngrid = 5 +z_nelement = 8 +z_nelement_local = 8 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 5 +r_nelement = 8 +r_nelement_local = 8 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 5 +vpa_nelement = 8 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 5 +vperp_nelement = 4 +vperp_L = 6.0 +#vperp_discretization = "finite_difference" +#vperp_discretization = "chebyshev_pseudospectral" +vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = 0.001 +vperp_dissipation_coefficient = 0.001 +#z_dissipation_coefficient = 0.1 +r_dissipation_coefficient = 0.01 + +[geometry] +option="1D-mirror" +DeltaB=0.5 +#option="constant-helical" +pitch=0.5 +rhostar = 1.0 diff --git a/test_scripts/chebyshev_radau_test.jl b/test_scripts/chebyshev_radau_test.jl index 57f2a0782..a4d7d6b65 100644 --- a/test_scripts/chebyshev_radau_test.jl +++ b/test_scripts/chebyshev_radau_test.jl @@ -1,10 +1,7 @@ export chebyshevradau_test using Printf -using Plots -using LaTeXStrings using MPI -using Measures import moment_kinetics using moment_kinetics.chebyshev @@ -39,7 +36,7 @@ at the lower endpoint of the domain (-1,1] in the normalised coordinate x. Here in the tests the shifted coordinate y is used with the vperp label so that the grid runs from (0,L]. """ -function chebyshevradau_test(; ngrid=5, L_in=3.0) +function chebyshevradau_test(; ngrid=5, L_in=3.0, discretization="chebyshev_pseudospectral") # elemental grid tests #ngrid = 17 @@ -49,7 +46,6 @@ function chebyshevradau_test(; ngrid=5, L_in=3.0) y_nelement_global = y_nelement_local # total number of elements y_L = L_in bc = "zero" - discretization = "gausslegendre_pseudospectral" # fd_option and adv_input not actually used so given values unimportant fd_option = "fourth_order_centered" cheb_option = "matrix" @@ -81,6 +77,17 @@ function chebyshevradau_test(; ngrid=5, L_in=3.0) println("f(y) = exp(-4 y) test") println("exact df: ",df_exact," num df: ",df_num," abs(err): ",df_err) + for iy in 1:y.n + ff[iy] = exp(-y.grid[iy]^2) + end + df_exact = 0.0 + df_num = sum(D0.*ff)/y.element_scale[1] + df_err = abs(df_num - df_exact) + println("f(y) = exp(-y^2) test") + println("exact df: ",df_exact," num df: ",df_num," abs(err): ",df_err) + f0 = -sum(D0[2:ngrid].*ff[2:ngrid])/D0[1] + println("exact f[0]: ",ff[1]," num f[0]: ",f0," abs(err): ",abs(f0-ff[1])) + for iy in 1:y.n ff[iy] = sin(y.grid[iy]) end @@ -89,6 +96,7 @@ function chebyshevradau_test(; ngrid=5, L_in=3.0) df_err = abs(df_num - df_exact) println("f(y) = sin(y) test") println("exact df: ",df_exact," num df: ",df_num," abs(err): ",df_err) + for iy in 1:y.n ff[iy] = y.grid[iy] + (y.grid[iy])^2 + (y.grid[iy])^3 end @@ -104,4 +112,4 @@ if abspath(PROGRAM_FILE) == @__FILE__ Pkg.activate(".") chebyshevradau_test() -end \ No newline at end of file +end