From 030fbe780dd93d3d0893b2ed8dae54314297440a Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 17 May 2023 11:42:45 +0100 Subject: [PATCH 01/34] Ported capability to compare upar and ppar to expectations from MMS test -- currently the tests do not pass. --- src/manufactured_solns.jl | 61 +++++++++++++++++++++++++++++++++++++-- src/post_processing.jl | 13 +++++++++ 2 files changed, 71 insertions(+), 3 deletions(-) diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 8f93faf2b..c432dd8c9 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -24,10 +24,12 @@ using IfElse if dfni_vpa_power_opt == "2" pvpa = 2 nconst = 0.25 + pconst = 3.0/4.0 fluxconst = 0.5 elseif dfni_vpa_power_opt == "4" pvpa = 4 - nconst = (3.0/8.0) + nconst = 3.0/8.0 + pconst = 15.0/8.0 fluxconst = 1.0 end @@ -199,7 +201,48 @@ using IfElse * "manufactured_solns:type=$(manufactured_solns_input.type)") end end - + + # ion mean parallel flow symbolic function + function upari_sym(Lr,Lz,r_bc,z_bc,composition,geometry,nr,manufactured_solns_input,species) + if z_bc == "periodic" + upari = 0.0 #not supported + 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) + rhostar = geometry.rhostar + bzed = geometry.bzed + epsilon = composition.epsilon_offset + alpha = composition.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 ) + end + return upari + end + + # ion parallel pressure symbolic function + function ppari_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + if z_bc == "periodic" + ppari = 0.0 # not supported + elseif z_bc == "wall" + densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + epsilon = composition.epsilon_offset + alpha = composition.alpha_switch + ppari = ( pconst*((0.5 - z/Lz)*nminus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha) + + (z/Lz + 0.5)*nplus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha)) + + (z/Lz + 0.5)*(0.5 - z/Lz)*nzero_sym(Lr,Lz,r_bc,z_bc,alpha) + - (2.0/(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))^2 ) + end + return ppari + end + + # ion perpendicular pressure symbolic function + function pperpi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + pperpi = densi # simple vperp^2 dependence of dfni + return pperpi + end function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,composition, manufactured_solns_input) if z_bc == "periodic" @@ -318,6 +361,12 @@ using IfElse densi = densi_sym(Lr, Lz, r_bc, z_bc, composition, manufactured_solns_input, charged_species) + upari = upari_sym(Lr, Lz, r_bc, z_bc, composition, geometry, nr, manufactured_solns_input, + charged_species) + ppari = ppari_sym(Lr, Lz, r_bc, z_bc, composition, manufactured_solns_input, + charged_species) + pperpi = pperpi_sym(Lr, Lz, r_bc, z_bc, composition, manufactured_solns_input, + charged_species) dfni = dfni_sym(Lr, Lz, r_bc, z_bc, composition, geometry, nr, manufactured_solns_input, charged_species) @@ -329,6 +378,9 @@ using IfElse #build julia functions from these symbolic expressions # cf. https://docs.juliahub.com/Symbolics/eABRO/3.4.0/tutorials/symbolic_functions/ densi_func = build_function(densi, z, r, t, expression=Val{false}) + upari_func = build_function(upari, z, r, t, expression=Val{false}) + ppari_func = build_function(ppari, z, r, t, expression=Val{false}) + pperpi_func = build_function(pperpi, z, r, t, expression=Val{false}) densn_func = build_function(densn, z, r, t, expression=Val{false}) dfni_func = build_function(dfni, vpa, vperp, z, r, t, expression=Val{false}) dfnn_func = build_function(dfnn, vz, vr, vzeta, z, r, t, expression=Val{false}) @@ -339,7 +391,10 @@ using IfElse # densn_func(zval, rval, tval) # dfnn_func(vzval, vrval, vzetapval, zval, rval, tval) - manufactured_solns_list = (densi_func = densi_func, densn_func = densn_func, dfni_func = dfni_func, dfnn_func = dfnn_func) + manufactured_solns_list = (densi_func = densi_func, densn_func = densn_func, + dfni_func = dfni_func, dfnn_func = dfnn_func, + upari_func = upari_func, ppari_func = ppari_func, + pperpi_func = pperpi_func) return manufactured_solns_list end diff --git a/src/post_processing.jl b/src/post_processing.jl index 3e1a7f1c8..2837bad8b 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -1099,6 +1099,9 @@ function analyze_and_plot_data(prefix...; run_index=nothing) composition, species, r_global.n) dfni_func = manufactured_solns_list.dfni_func densi_func = manufactured_solns_list.densi_func + upari_func = manufactured_solns_list.upari_func + ppari_func = manufactured_solns_list.ppari_func + pperpi_func = manufactured_solns_list.pperpi_func dfnn_func = manufactured_solns_list.dfnn_func densn_func = manufactured_solns_list.densn_func manufactured_E_fields = @@ -1135,16 +1138,26 @@ function analyze_and_plot_data(prefix...; run_index=nothing) # ion test density_sym = copy(density[:,:,:,:]) + upar_sym = copy(density[:,:,:,:]) + ppar_sym = copy(density[:,:,:,:]) + pperp_sym = copy(density[:,:,:,:]) is = 1 for it in 1:ntime for ir in 1:r_global.n for iz in 1:z_global.n density_sym[iz,ir,is,it] = densi_func(z_global.grid[iz],r_global.grid[ir],time[it]) + upar_sym[iz,ir,is,it] = upari_func(z_global.grid[iz],r_global.grid[ir],time[it]) + ppar_sym[iz,ir,is,it] = ppari_func(z_global.grid[iz],r_global.grid[ir],time[it]) + pperp_sym[iz,ir,is,it] = pperpi_func(z_global.grid[iz],r_global.grid[ir],time[it]) end end end compare_moments_symbolic_test(run_name_label,density,density_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, L"\widetilde{n}_i",L"\widetilde{n}_i^{sym}",L"\varepsilon(\widetilde{n}_i)","dens") + compare_moments_symbolic_test(run_name_label,parallel_flow,upar_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + L"\widetilde{u}_{\|\|i}",L"\widetilde{u}_{\|\|i}^{sym}",L"\varepsilon(\widetilde{u}_{\|\|i})","upar") + compare_moments_symbolic_test(run_name_label,parallel_pressure,ppar_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + L"\widetilde{p}_{\|\|i}",L"\widetilde{p}_{\|\|i}^{sym}",L"\varepsilon(\widetilde{p}_{\|\|i})","ppar") compare_charged_pdf_symbolic_test(run_name_label,manufactured_solns_list,"ion", L"\widetilde{f}_i",L"\widetilde{f}^{sym}_i",L"\varepsilon(\widetilde{f}_i)","pdf") From 6a9664d2822fe9d847ea8d1cacdc4a70448a3f4d Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Sep 2023 16:17:00 +0100 Subject: [PATCH 02/34] Clean up manufactured solutions input Options for manufactured solutions were moved from 'composition' to a dedicated NamedTuple. However, the fields in the 'composition' struct were not removed - remove them, and fix bugs caused by those fields still having been used, although they were not being set correctly from the input options. --- src/input_structs.jl | 8 +------- src/manufactured_solns.jl | 13 ++++++------- src/moment_kinetics_input.jl | 9 +-------- 3 files changed, 8 insertions(+), 22 deletions(-) diff --git a/src/input_structs.jl b/src/input_structs.jl index ed3d20914..285e42bcd 100644 --- a/src/input_structs.jl +++ b/src/input_structs.jl @@ -259,13 +259,7 @@ mutable struct species_composition phi_wall::mk_float # constant for testing nonzero Er Er_constant::mk_float - # constant controlling divergence at wall boundaries in MMS test - epsilon_offset::mk_float - # logical controlling whether or not dfni(vpabar,z,r) or dfni(vpa,z,r) in MMS test - use_vpabar_in_mms_dfni::Bool - # associated float controlling form of assumed potential in MMS test - alpha_switch::mk_float - # ratio of the neutral particle mass to the ion mass + # ratio of the neutral particle mass to the ion mass mn_over_mi::mk_float # ratio of the electron particle mass to the ion mass me_over_mi::mk_float diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index c432dd8c9..000c342d7 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -211,8 +211,8 @@ using IfElse Er, Ez, phi = electric_fields(Lr,Lz,r_bc,z_bc,composition,nr,manufactured_solns_input,species) rhostar = geometry.rhostar bzed = geometry.bzed - epsilon = composition.epsilon_offset - alpha = composition.alpha_switch + 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 ) @@ -226,8 +226,8 @@ using IfElse ppari = 0.0 # not supported elseif z_bc == "wall" densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) - epsilon = composition.epsilon_offset - alpha = composition.alpha_switch + epsilon = manufactured_solns_input.epsilon_offset + alpha = manufactured_solns_input.alpha_switch ppari = ( pconst*((0.5 - z/Lz)*nminus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha) + (z/Lz + 0.5)*nplus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha)) + (z/Lz + 0.5)*(0.5 - z/Lz)*nzero_sym(Lr,Lz,r_bc,z_bc,alpha) @@ -243,8 +243,7 @@ using IfElse pperpi = densi # simple vperp^2 dependence of dfni return pperpi end - function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,composition, - manufactured_solns_input) + function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,manufactured_solns_input) if z_bc == "periodic" jpari_into_LHS_wall_sym = 0.0 elseif z_bc == "wall" @@ -316,7 +315,7 @@ using IfElse # get N_e factor for boltzmann response if composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr == 1 # so 1D MMS test with 3V neutrals where ion current can be calculated prior to knowing Er - jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc, composition, + jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc, manufactured_solns_input) N_e = -2.0*sqrt(pi*composition.me_over_mi)*exp(-composition.phi_wall/composition.T_e)*jpari_into_LHS_wall elseif composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr > 1 diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 752fb68c5..a15e1c733 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -844,20 +844,13 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) phi_wall = 0.0 # constant to test nonzero Er Er_constant = 0.0 - # constant to control Ez divergence - epsilon_offset = 0.001 - # bool to control functional form of dfni in MMS test - use_vpabar_in_mms_dfni = true - # float to control form of MMS density/potential/Er/Ez - alpha_switch = 1.0 # ratio of the neutral particle mass to the ion particle mass mn_over_mi = 1.0 # ratio of the electron particle mass to the ion particle mass me_over_mi = 1.0/1836.0 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, - epsilon_offset, use_vpabar_in_mms_dfni, alpha_switch, mn_over_mi, me_over_mi, - allocate_float(n_species)) + mn_over_mi, me_over_mi, allocate_float(n_species)) species_charged = Array{species_parameters_mutable,1}(undef,n_ion_species) species_neutral = Array{species_parameters_mutable,1}(undef,n_neutral_species) From 02ded9ac90a3b42f445d6bc863c86198da978224 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Sep 2023 17:07:56 +0100 Subject: [PATCH 03/34] Put 'show_element_boundaries' in postproc input for manufactured_solns --- src/makie_post_processing.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index c2a9e2704..1b0cebefe 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -673,6 +673,7 @@ function _setup_single_input!(this_input_dict::OrderedDict{String,Any}, (o=>section_defaults[String(o)] for o ∈ (:it0, :ir0, :iz0, :ivperp0, :ivpa0, :ivzeta0, :ivr0, :ivz0))..., colormap=this_input_dict["colormap"], animation_ext=this_input_dict["animation_ext"], + show_element_boundaries=this_input_dict["show_element_boundaries"], ) sort!(this_input_dict["manufactured_solns"]) From bf27259ded892efa81c1ad29325ec9743bd7b01e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Sep 2023 17:08:50 +0100 Subject: [PATCH 04/34] Make MMS plots for ion parallel flow and parallel pressure --- src/makie_post_processing.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index 1b0cebefe..3926cbec4 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -5438,6 +5438,7 @@ function manufactured_solutions_get_field_and_field_sym(run_info, variable_name; variable_name = Symbol(variable_name) func_name_lookup = (phi=:phi_func, Er=:Er_func, Ez=:Ez_func, density=:densi_func, + parallel_flow=:upari_func, parallel_pressure=:ppari_func, density_neutral=:densn_func, f=:dfni_func, f_neutral=:dfnn_func) nt = run_info.nt @@ -5467,7 +5468,8 @@ function manufactured_solutions_get_field_and_field_sym(run_info, variable_name; run_info.z.bc, run_info.composition, run_info.r.n, run_info.manufactured_solns_input, run_info.species) - elseif variable_name ∈ (:density, :density_neutral, :f, :f_neutral) + elseif variable_name ∈ (:density, :parallel_flow, :parallel_pressure, + :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, @@ -6353,6 +6355,8 @@ function manufactured_solutions_analysis(run_info; plot_prefix) ("Er", L"\tilde{E}_r", L"\tilde{E}_r^{sym}", L"\varepsilon(\tilde{E}_r)"), ("Ez", L"\tilde{E}_z", L"\tilde{E}_z^{sym}", L"\varepsilon(\tilde{E}_z)"), ("density", L"\tilde{n}_i", L"\tilde{n}_i^{sym}", L"\varepsilon(\tilde{n}_i)"), + ("parallel_flow", L"\tilde{u}_{i,\parallel}", L"\tilde{u}_{i,\parallel}^{sym}", L"\varepsilon(\tilde{u}_{i,\parallel})"), + ("parallel_pressure", L"\tilde{p}_{i,\parallel}", L"\tilde{p}_{i,\parallel}^{sym}", L"\varepsilon(\tilde{p}_{i,\parallel})"), ("density_neutral", L"\tilde{n}_n", L"\tilde{n}_n^{sym}", L"\varepsilon(\tilde{n}_n)")) if contains(variable_name, "neutral") && run_info.n_neutral_species == 0 From a1a17ec80fb942f2b9d7f6d831cc63e8bcfdcd7d Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Sep 2023 17:09:59 +0100 Subject: [PATCH 05/34] Skip MMS plots of Er for 1D cases These plots are uninteresting, and cause errors because the y-axis limits are both 0. --- src/makie_post_processing.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index 3926cbec4..7b7dd62d2 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -6362,6 +6362,9 @@ function manufactured_solutions_analysis(run_info; plot_prefix) if contains(variable_name, "neutral") && run_info.n_neutral_species == 0 continue end + if contains(variable_name, "Er") && run_info.r.n_global == 1 + continue + end compare_moment_symbolic_test(run_info, plot_prefix, field_label, field_sym_label, norm_label, variable_name; io=io, input=input) From b539670d1c750253cd2885357f85568f04a92fd0 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Sep 2023 17:11:03 +0100 Subject: [PATCH 06/34] Fix slicing of distribution functions in MMS plots --- src/makie_post_processing.jl | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index 7b7dd62d2..a98855b2a 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -5782,13 +5782,19 @@ plots/animations, `field_sym_label` is the label for the manufactured solution, omitting `:t`. """ function _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label, - field_sym_label, norm_label, plot_dims, animate_dims) + field_sym_label, norm_label, plot_dims, animate_dims, neutrals) nt = run_info.nt time = run_info.time - all_plot_slices = Tuple(Symbol(:i, d)=>input[Symbol(:i, d, :0)] for d ∈ plot_dims) - all_animate_slices = Tuple(Symbol(:i, d)=>input[Symbol(:i, d, :0)] for d ∈ animate_dims) + if neutrals + all_dims_no_t = (:r, :z, :vzeta, :vr, :vz) + else + all_dims_no_t = (:r, :z, :vperp, :vpa) + end + all_dims = tuple(:t, all_dims_no_t...) + all_plot_slices = Tuple(Symbol(:i, d)=>input[Symbol(:i, d, :0)] for d ∈ all_dims) + all_animate_slices = Tuple(Symbol(:i, d)=>input[Symbol(:i, d, :0)] for d ∈ all_dims_no_t) # Options to produce either regular or log-scale plots epsilon = 1.0e-30 # minimum data value to include in log plots @@ -6079,7 +6085,7 @@ function compare_charged_pdf_symbolic_test(run_info, plot_prefix; io=nothing, end plot_dims = tuple(:t, animate_dims...) _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label, - field_sym_label, norm_label, plot_dims, animate_dims) + field_sym_label, norm_label, plot_dims, animate_dims, false) return field_norm end @@ -6286,7 +6292,7 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, end plot_dims = tuple(:t, animate_dims...) _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label, - field_sym_label, norm_label, plot_dims, animate_dims) + field_sym_label, norm_label, plot_dims, animate_dims, true) return field_norm end From f181300f27a94116cd8e94a014bf4ff5870e2513 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Thu, 28 Sep 2023 14:58:08 +0100 Subject: [PATCH 07/34] Fix post processing script consistent with updated composition struct. --- src/post_processing.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/post_processing.jl b/src/post_processing.jl index 2837bad8b..160a2a84d 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -311,8 +311,7 @@ function get_geometry_and_composition(scan_input,n_ion_species,n_neutral_species me_over_mi = 1.0/1836.0 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, - epsilon_offset, use_vpabar_in_mms_dfni, alpha_switch, mn_over_mi, me_over_mi, - allocate_float(n_species)) + mn_over_mi, me_over_mi, allocate_float(n_species)) return geometry, composition end From d91c61866ca225543c009eb17af9d15e0a4a42f6 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 29 Sep 2023 14:09:33 +0100 Subject: [PATCH 08/34] Addition of factor to compensate for the pressure normalisation convention in master, which seems to normalise pressure to 2 nref Tref rather than nref Tref = (1/2)mref nref cref^2. --- src/manufactured_solns.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 000c342d7..4cf0adcea 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -222,6 +222,8 @@ using IfElse # ion parallel pressure symbolic function function ppari_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + # normalisation factor due to strange pressure normalisation convention in master + norm_fac = 0.5 if z_bc == "periodic" ppari = 0.0 # not supported elseif z_bc == "wall" @@ -234,7 +236,7 @@ using IfElse - (2.0/(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))^2 ) end - return ppari + return ppari*norm_fac end # ion perpendicular pressure symbolic function From 1932bf5dfe1c26bee57f5bb2fde67f76ce55adae Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 12 May 2023 11:24:07 +0100 Subject: [PATCH 09/34] Porting of commit cd07a26. Refactoring of the update_derived_moments function to use get_density, get_upar and get_ppar functions. In the split moment operation, the moments passed in to these functions along with the distribution are set to the enforced values of (density,upar) = (1.0,0.0) appropriately. This branch now supports the test_velocity_integrals.jl script. --- src/velocity_moments.jl | 76 ++++++++++++++++--------- test_velocity_integrals.jl | 112 +++++++++++++++++++++++++++++++++++++ 2 files changed, 161 insertions(+), 27 deletions(-) create mode 100644 test_velocity_integrals.jl diff --git a/src/velocity_moments.jl b/src/velocity_moments.jl index 8db08f9c3..45e5aa9ae 100644 --- a/src/velocity_moments.jl +++ b/src/velocity_moments.jl @@ -22,6 +22,12 @@ export update_neutral_pr! export update_neutral_pzeta! export update_neutral_qz! +# for testing +export get_density +export get_upar +export get_ppar +export get_pperp + using ..type_definitions: mk_float using ..array_allocation: allocate_shared_float, allocate_bool, allocate_float using ..calculus: integral @@ -431,13 +437,16 @@ function update_density_species!(dens, ff, vpa, vperp, z, r) @loop_r_z ir iz begin # When evolve_density = false, the evolved pdf is the 'true' pdf, and the vpa # coordinate is (dz/dt) / c_s. - # Integrating calculates n_s / N_e = (1/√π)∫d(vpa/c_s) (√π f_s c_s / N_e) - dens[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), - vpa.grid, 0, vpa.wgts, vperp.grid, 0, vperp.wgts) + dens[iz,ir] = get_density(@view(ff[:,:,iz,ir]), vpa, vperp) end return nothing end +function get_density(ff, vpa, vperp) + # Integrating calculates n_s / N_e = (0/√π)∫d(vpa/c_s) (√π f_s c_s / N_e) + return integrate_over_vspace(@view(ff[:,:]), vpa.grid, 0, vpa.wgts, vperp.grid, 0, vperp.wgts) +end + """ NB: if this function is called and if upar_updated is false, then the incoming pdf is the un-normalized pdf that satisfies int dv pdf = density @@ -475,10 +484,10 @@ function update_upar_species!(upar, density, ppar, ff, vpa, vperp, z, r, evolve_ # Integrating calculates # (upar_s / vth_s) = (1/√π)∫d(vpa/vth_s) * (vpa/vth_s) * (√π f_s vth_s / n_s) # so convert from upar_s / vth_s to upar_s / c_s + # we set the input density to get_upar = 1.0 as the normalised distribution has density of 1.0 @loop_r_z ir iz begin vth = sqrt(2.0*ppar[iz,ir]/density[iz,ir]) - upar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), - vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts) * vth + upar[iz,ir] = vth*get_upar(@view(ff[:,:,iz,ir]), vpa, vperp, 1.0) end elseif evolve_density # corresponds to case where only the density is evolved separately from the @@ -486,9 +495,9 @@ function update_upar_species!(upar, density, ppar, ff, vpa, vperp, z, r, evolve_ # (dz/dt) / c_s. # Integrating calculates # (upar_s / c_s) = (1/√π)∫d(vpa/c_s) * (vpa/c_s) * (√π f_s c_s / n_s) + # we set the input density to get_upar = 1.0 as the normalised distribution has density of 1.0 @loop_r_z ir iz begin - upar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), - vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts) + upar[iz,ir] = get_upar(@view(ff[:,:,iz,ir]), vpa, vperp, 1.0) end else # When evolve_density = false, the evolved pdf is the 'true' pdf, @@ -496,14 +505,21 @@ function update_upar_species!(upar, density, ppar, ff, vpa, vperp, z, r, evolve_ # Integrating calculates # (n_s / N_e) * (upar_s / c_s) = (1/√π)∫d(vpa/c_s) * (vpa/c_s) * (√π f_s c_s / N_e) @loop_r_z ir iz begin - upar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), - vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts) / - density[iz,ir] + upar[iz,ir] = get_upar(@view(ff[:,:,iz,ir]), vpa, vperp, density[iz,ir]) end end return nothing end +function get_upar(ff, vpa, vperp, density) + # Integrating calculates + # (n_s / N_e) * (upar_s / c_s) = (1/√π)∫d(vpa/c_s) * (vpa/c_s) * (√π f_s c_s / N_e) + # so we divide by the density of f_s + upar = integrate_over_vspace(@view(ff[:,:]), vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts) + upar /= density + return upar +end + """ NB: if this function is called and if ppar_updated is false, then the incoming pdf is the un-normalized pdf that satisfies int dv pdf = density @@ -540,40 +556,46 @@ function update_ppar_species!(ppar, density, upar, ff, vpa, vperp, z, r, evolve_ if evolve_upar # this is the case where the parallel flow and density are evolved separately # from the normalized pdf, g_s = (√π f_s c_s / n_s); the vpa coordinate is - # ((dz/dt) - upar_s) / c_s> - # Integrating calculates (p_parallel/m_s n_s c_s^2) = (1/√π)∫d((vpa-upar_s)/c_s) (1/2)*((vpa-upar_s)/c_s)^2 * (√π f_s c_s / n_s) - # so convert from p_s / m_s n_s c_s^2 to ppar_s = p_s / m_s N_e c_s^2 + # ((dz/dt) - upar_s) / c_s> and so we set upar = 0 in the call to get_ppar + # because the mean flow of the normalised ff is zero @loop_r_z ir iz begin - ppar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), vpa.grid, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) * - density[iz,ir] + ppar[iz,ir] = density[iz,ir]*get_ppar(@view(ff[:,:,iz,ir]), vpa, vperp, 0.0) end elseif evolve_density # corresponds to case where only the density is evolved separately from the # normalised pdf, given by g_s = (√π f_s c_s / n_s); the vpa coordinate is # (dz/dt) / c_s. - # Integrating calculates - # (p_parallel/m_s n_s c_s^2) + (upar_s/c_s)^2 = (1/√π)∫d(vpa/c_s) (vpa/c_s)^2 * (√π f_s c_s / n_s) - # so subtract off the mean kinetic energy and multiply by density to get the - # internal energy density (aka pressure) @loop_r_z ir iz begin - ppar[iz,ir] = (integrate_over_vspace(@view(ff[:,:,iz,ir]), vpa.grid, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) - - upar[iz,ir]^2) * density[iz,ir] + ppar[iz,ir] = density[iz,ir]*get_ppar(@view(ff[:,:,iz,ir]), vpa, vperp, upar[iz,ir]) end else # When evolve_density = false, the evolved pdf is the 'true' pdf, # and the vpa coordinate is (dz/dt) / c_s. - # Integrating calculates - # (p_parallel/m_s N_e c_s^2) + (n_s/N_e)*(upar_s/c_s)^2 = (1/√π)∫d(vpa/c_s) (vpa/c_s)^2 * (√π f_s c_s / N_e) - # so subtract off the mean kinetic energy density to get the internal energy - # density (aka pressure) @loop_r_z ir iz begin - ppar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), vpa.grid, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) - - density[iz,ir]*upar[iz,ir]^2 + ppar[iz,ir] = get_ppar(@view(ff[:,:,iz,ir]), vpa, vperp, upar[iz,ir]) end end return nothing end +function get_ppar(ff, vpa, vperp, upar) + # Integrating calculates + # (p_parallel/m_s N_e c_s^2) = (1/√π)∫d(vpa/c_s) ((vpa-upar)/c_s)^2 * (√π f_s c_s / N_e) + # the internal energy density (aka pressure of f_s) + + # modify input vpa.grid to account for the mean flow + @. vpa.scratch = vpa.grid - upar + norm_fac = 1.0 # normalise to m_s N_e c_s^2 + #norm_fac = 2.0 # normalise to 0.5 m_s N_e c_s^2 = N_e T_s + return norm_fac*integrate_over_vspace(@view(ff[:,:]), vpa.scratch, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) +end + +function get_pperp(ff, vpa, vperp) + norm_fac = 0.5 # normalise to m_s N_e c_s^2 + #norm_fac = 1.0 # normalise to 0.5 m_s N_e c_s^2 = N_e T_s + return norm_fac*integrate_over_vspace(@view(ff[:,:]), vpa.grid, 0, vpa.wgts, vperp.grid, 2, vperp.wgts) +end + """ NB: the incoming pdf is the normalized pdf """ diff --git a/test_velocity_integrals.jl b/test_velocity_integrals.jl new file mode 100644 index 000000000..bb39123b1 --- /dev/null +++ b/test_velocity_integrals.jl @@ -0,0 +1,112 @@ +using Printf +using Plots +using LaTeXStrings +using MPI +using Measures + +if abspath(PROGRAM_FILE) == @__FILE__ + using Pkg + Pkg.activate(".") + + import moment_kinetics + using moment_kinetics.input_structs: grid_input, advection_input + using moment_kinetics.coordinates: define_coordinate + using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_pperp + using moment_kinetics.array_allocation: allocate_float + + # define inputs needed for the test + ngrid = 17 #number of points per element + nelement_local = 20 # number of elements per rank + nelement_global = nelement_local # total number of elements + Lvpa = 12.0 #physical box size in reference units + Lvperp = 6.0 #physical box size in reference units + bc = "" #not required to take a particular value, not used + # fd_option and adv_input not actually used so given values unimportant + discretization = "chebyshev_pseudospectral" + fd_option = "fourth_order_centered" + adv_input = advection_input("default", 1.0, 0.0, 0.0) + nrank = 1 + irank = 0 + comm = MPI.COMM_NULL + # create the 'input' struct containing input info needed to create a + # coordinate + vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, + nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) + vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, + nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm) + # create the coordinate struct 'x' + println("made inputs") + println("vpa: ngrid: ",ngrid," nelement: ",nelement_local) + println("vperp: ngrid: ",ngrid," nelement: ",nelement_local) + vpa, vpa_spectral = define_coordinate(vpa_input) + vperp, vperp_spectral = define_coordinate(vperp_input) + + dfn = allocate_float(vpa.n,vperp.n) + + function pressure(ppar,pperp) + pres = (1.0/3.0)*(ppar + 2.0*pperp) + return pres + end + # assign a known isotropic Maxwellian distribution in normalised units + dens = 3.0/4.0 + upar = 2.0/3.0 + ppar = 2.0/3.0 + pperp = 2.0/3.0 + pres = pressure(ppar,pperp) + mass = 1.0 + vth = sqrt(2.0*pres/(dens*mass)) + for ivperp in 1:vperp.n + for ivpa in 1:vpa.n + vpa_val = vpa.grid[ivpa] + vperp_val = vperp.grid[ivperp] + dfn[ivpa,ivperp] = (dens/vth^3)*exp( - ((vpa_val-upar)^2 + vperp_val^2)/vth^2 ) + end + end + + # now check that we can extract the correct moments from the distribution + + dens_test = get_density(dfn,vpa,vperp) + upar_test = get_upar(dfn,vpa,vperp,dens_test) + ppar_test = get_ppar(dfn,vpa,vperp,upar_test) + pperp_test = get_pperp(dfn,vpa,vperp) + pres_test = pressure(ppar_test,pperp_test) + # output test results + println("") + println("Isotropic Maxwellian") + println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) + println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) + println("pres_test: ", pres_test, " pres: ", pres, " error: ", abs(pres_test-pres)) + println("") + + # assign a known biMaxwellian distribution in normalised units + dens = 3.0/4.0 + upar = 2.0/3.0 + ppar = 4.0/5.0 + pperp = 1.0/4.0 + mass = 1.0 + vthpar = sqrt(2.0*ppar/(dens*mass)) + vthperp = sqrt(2.0*pperp/(dens*mass)) + for ivperp in 1:vperp.n + for ivpa in 1:vpa.n + vpa_val = vpa.grid[ivpa] + vperp_val = vperp.grid[ivperp] + dfn[ivpa,ivperp] = (dens/(vthpar*vthperp^2))*exp( - ((vpa_val-upar)^2)/vthpar^2 - (vperp_val^2)/vthperp^2 ) + end + end + + # now check that we can extract the correct moments from the distribution + + dens_test = get_density(dfn,vpa,vperp) + upar_test = get_upar(dfn,vpa,vperp,dens_test) + ppar_test = get_ppar(dfn,vpa,vperp,upar_test) + pperp_test = get_pperp(dfn,vpa,vperp) + # output test results + + println("") + println("biMaxwellian") + println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) + println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) + println("ppar_test: ", ppar_test, " ppar: ", ppar, " error: ", abs(ppar_test-ppar)) + println("pperp_test: ", pperp_test, " pperp: ", pperp, " error: ", abs(pperp_test-pperp)) + println("") +end From fa279702c201addd7840bea1218f6e0d04145497 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 17 May 2023 14:27:26 +0100 Subject: [PATCH 10/34] Cherry-pick of commit 8632cf7 to add perpendicular pressure to the body of the code --- src/file_io.jl | 28 ++++++++++++++++++++++++--- src/initial_conditions.jl | 6 ++++-- src/moment_kinetics_structs.jl | 1 + src/post_processing.jl | 11 +++++++++-- src/time_advance.jl | 12 +++++++----- src/velocity_moments.jl | 35 +++++++++++++++++++++++++++++++++- 6 files changed, 80 insertions(+), 13 deletions(-) diff --git a/src/file_io.jl b/src/file_io.jl index bb53a90f2..94f3e4123 100644 --- a/src/file_io.jl +++ b/src/file_io.jl @@ -61,6 +61,8 @@ struct io_moments_info{Tfile, Ttime, Tphi, Tmomi, Tmomn} parallel_flow::Tmomi # handle for the charged species parallel pressure parallel_pressure::Tmomi + # handle for the charged species perpendicular pressure + perpendicular_pressure::Tmomi # handle for the charged species parallel heat flux parallel_heat_flux::Tmomi # handle for the charged species thermal speed @@ -587,6 +589,13 @@ function define_dynamic_moment_variables!(fid, n_ion_species, n_neutral_species, parallel_io=parallel_io, description="charged species parallel pressure", units="n_ref*T_ref") + + # io_pperp is the handle for the ion parallel pressure + io_pperp = create_dynamic_variable!(dynamic, "perpendicular_pressure", mk_float, z, r; + n_ion_species=n_ion_species, + parallel_io=parallel_io, + description="charged species perpendicular pressure", + units="n_ref*T_ref") # io_qpar is the handle for the ion parallel heat flux io_qpar = create_dynamic_variable!(dynamic, "parallel_heat_flux", mk_float, z, r; @@ -638,7 +647,7 @@ function define_dynamic_moment_variables!(fid, n_ion_species, n_neutral_species, units="c_ref") return io_moments_info(fid, io_time, io_phi, io_Er, io_Ez, io_density, io_upar, - io_ppar, io_qpar, io_vth, io_density_neutral, io_uz_neutral, + io_ppar, io_pperp, io_qpar, io_vth, io_density_neutral, io_uz_neutral, io_pz_neutral, io_qz_neutral, io_thermal_speed_neutral, parallel_io) end @@ -777,6 +786,7 @@ function reopen_moments_io(file_info) return io_moments_info(fid, getvar("time"), getvar("phi"), getvar("Er"), getvar("Ez"), getvar("density"), getvar("parallel_flow"), getvar("parallel_pressure"), getvar("parallel_heat_flux"), + getvar("perpendicular_pressure"), getvar("thermal_speed"), getvar("density_neutral"), getvar("uz_neutral"), getvar("pz_neutral"), getvar("qz_neutral"), getvar("thermal_speed_neutral"), @@ -861,6 +871,7 @@ function reopen_dfns_io(file_info) io_moments = io_moments_info(fid, getvar("time"), getvar("phi"), getvar("Er"), getvar("Ez"), getvar("density"), getvar("parallel_flow"), getvar("parallel_pressure"), + getvar("perpendicular_pressure"), getvar("parallel_heat_flux"), getvar("thermal_speed"), getvar("density_neutral"), getvar("uz_neutral"), getvar("pz_neutral"), @@ -916,6 +927,8 @@ function write_moments_data_to_binary(moments, fields, t, n_ion_species, n_ion_species) append_to_dynamic_var(io_moments.parallel_pressure, moments.charged.ppar, t_idx, z, r, n_ion_species) + append_to_dynamic_var(io_moments.perpendicular_pressure, moments.charged.pperp, t_idx, + z, r, n_ion_species) append_to_dynamic_var(io_moments.parallel_heat_flux, moments.charged.qpar, t_idx, z, r, n_ion_species) append_to_dynamic_var(io_moments.thermal_speed, moments.charged.vth, t_idx, z, r, @@ -1004,6 +1017,8 @@ end t_idx, z, r, n_ion_species) append_to_dynamic_var(io_moments.parallel_pressure, moments.charged.ppar.data, t_idx, z, r, n_ion_species) + append_to_dynamic_var(io_moments.perpendicular_pressure.data, moments.charged.pperp, + t_idx, z, r, n_ion_species) append_to_dynamic_var(io_moments.parallel_heat_flux, moments.charged.qpar.data, t_idx, z, r, n_ion_species) append_to_dynamic_var(io_moments.thermal_speed, moments.charged.vth.data, @@ -1267,7 +1282,7 @@ all the arrays have the same length, with an entry for each call to `debug_dump( function debug_dump end function debug_dump(vz::coordinate, vr::coordinate, vzeta::coordinate, vpa::coordinate, vperp::coordinate, z::coordinate, r::coordinate, t::mk_float; - ff=nothing, dens=nothing, upar=nothing, ppar=nothing, qpar=nothing, + ff=nothing, dens=nothing, upar=nothing, ppar=nothing, pperp=nothing, qpar=nothing, vth=nothing, ff_neutral=nothing, dens_neutral=nothing, uz_neutral=nothing, #ur_neutral=nothing, uzeta_neutral=nothing, @@ -1366,6 +1381,11 @@ function debug_dump(vz::coordinate, vr::coordinate, vzeta::coordinate, vpa::coor else debug_output_file.moments.parallel_pressure[:,:,:,debug_output_counter[]] = ppar end + if pperp === nothing + debug_output_file.moments.perpendicular_pressure[:,:,:,debug_output_counter[]] = 0.0 + else + debug_output_file.moments.perpendicular_pressure[:,:,:,debug_output_counter[]] = pperp + end if qpar === nothing debug_output_file.moments.parallel_heat_flux[:,:,:,debug_output_counter[]] = 0.0 else @@ -1442,6 +1462,7 @@ function debug_dump(fvec::Union{scratch_pdf,Nothing}, density = nothing upar = nothing ppar = nothing + pperp = nothing pdf_neutral = nothing density_neutral = nothing else @@ -1449,6 +1470,7 @@ function debug_dump(fvec::Union{scratch_pdf,Nothing}, density = fvec.density upar = fvec.upar ppar = fvec.ppar + pperp = fvec.pperp pdf_neutral = fvec.pdf_neutral density_neutral = fvec.density_neutral end @@ -1462,7 +1484,7 @@ function debug_dump(fvec::Union{scratch_pdf,Nothing}, Ez = fields.Ez end return debug_dump(vz, vr, vzeta, vpa, vperp, z, r, t; ff=pdf, dens=density, upar=upar, - ppar=ppar, ff_neutral=pdf_neutral, dens_neutral=density_neutral, + ppar=ppar, pperp=pperp, ff_neutral=pdf_neutral, dens_neutral=density_neutral, phi=phi, Er=Er, Ez=Ez, t, istage=istage, label=label) end diff --git a/src/initial_conditions.jl b/src/initial_conditions.jl index f4fe358ec..92cb07ed8 100644 --- a/src/initial_conditions.jl +++ b/src/initial_conditions.jl @@ -31,7 +31,7 @@ using ..velocity_moments: create_moments_charged, create_moments_neutral, update using ..velocity_moments: moments_charged_substruct, moments_neutral_substruct using ..velocity_moments: update_neutral_density!, update_neutral_pz!, update_neutral_pr!, update_neutral_pzeta! using ..velocity_moments: update_neutral_uz!, update_neutral_ur!, update_neutral_uzeta!, update_neutral_qz! -using ..velocity_moments: update_ppar!, update_upar!, update_density!, reset_moments_status! +using ..velocity_moments: update_ppar!, update_upar!, update_density!, update_pperp!, reset_moments_status! using ..manufactured_solns: manufactured_solutions @@ -157,8 +157,9 @@ function init_pdf_and_moments!(pdf, moments, boundary_distributions, geometry, init_upar!(moments.charged.upar, z, r, species.charged, n_ion_species) # initialise the parallel thermal speed profile init_vth!(moments.charged.vth, z, r, species.charged, n_ion_species) + # initialise pressures assuming isotropic distribution @. moments.charged.ppar = 0.5 * moments.charged.dens * moments.charged.vth^2 - + @. moments.charged.pperp = moments.charged.ppar if n_neutral_species > 0 #neutral particles init_density!(moments.neutral.dens, z, r, species.neutral, n_neutral_species) @@ -924,6 +925,7 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa, moments.charged.dens, moments.charged.upar, pdf.charged.norm, vpa, vperp, z, r, composition, moments.evolve_density, moments.evolve_upar) + update_pperp!(moments.charged.pperp, pdf.charged.norm, vpa, vperp, z, r, composition) update_qpar!(moments.charged.qpar, moments.charged.qpar_updated, moments.charged.dens, moments.charged.upar, moments.charged.vth, pdf.charged.norm, vpa, vperp, z, r, diff --git a/src/moment_kinetics_structs.jl b/src/moment_kinetics_structs.jl index fa63cd9d0..949250317 100644 --- a/src/moment_kinetics_structs.jl +++ b/src/moment_kinetics_structs.jl @@ -16,6 +16,7 @@ struct scratch_pdf{n_distribution_charged, n_moment, n_distribution_neutral,n_mo density::MPISharedArray{mk_float, n_moment} upar::MPISharedArray{mk_float, n_moment} ppar::MPISharedArray{mk_float, n_moment} + pperp::MPISharedArray{mk_float, n_moment} temp_z_s::MPISharedArray{mk_float, n_moment} # neutral particles pdf_neutral::MPISharedArray{mk_float, n_distribution_neutral} diff --git a/src/post_processing.jl b/src/post_processing.jl index 160a2a84d..69287f396 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -223,9 +223,10 @@ function allocate_global_zr_charged_moments(nz_global,nr_global,n_ion_species,nt density = allocate_float(nz_global,nr_global,n_ion_species,ntime) parallel_flow = allocate_float(nz_global,nr_global,n_ion_species,ntime) parallel_pressure = allocate_float(nz_global,nr_global,n_ion_species,ntime) + perpendicular_pressure = allocate_float(nz_global,nr_global,n_ion_species,ntime) parallel_heat_flux = allocate_float(nz_global,nr_global,n_ion_species,ntime) thermal_speed = allocate_float(nz_global,nr_global,n_ion_species,ntime) - return density, parallel_flow, parallel_pressure, parallel_heat_flux, thermal_speed + return density, parallel_flow, parallel_pressure, perpendicular_pressure, parallel_heat_flux, thermal_speed end function allocate_global_zr_charged_dfns(nvpa_global, nvperp_global, nz_global, nr_global, @@ -488,7 +489,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) Tuple(this_z.n_global for this_z ∈ z), Tuple(this_r.n_global for this_r ∈ r), ntime) - density, parallel_flow, parallel_pressure, parallel_heat_flux, thermal_speed = + density, parallel_flow, parallel_pressure, perpendicular_pressure, parallel_heat_flux, thermal_speed = get_tuple_of_return_values(allocate_global_zr_charged_moments, Tuple(this_z.n_global for this_z ∈ z), Tuple(this_r.n_global for this_r ∈ r), @@ -531,6 +532,10 @@ function analyze_and_plot_data(prefix...; run_index=nothing) "parallel_pressure", run_names, "moments", nblocks, Tuple(this_z.n for this_z ∈ z), Tuple(this_r.n for this_r ∈ r), iskip) + get_tuple_of_return_values(read_distributed_zr_data!, perpendicular_pressure, + "perpendicular_pressure", run_names, "moments", nblocks, + Tuple(this_z.n for this_z ∈ z), + Tuple(this_r.n for this_r ∈ r), iskip) get_tuple_of_return_values(read_distributed_zr_data!, parallel_heat_flux, "parallel_heat_flux", run_names, "moments", nblocks, Tuple(this_z.n for this_z ∈ z), @@ -1157,6 +1162,8 @@ function analyze_and_plot_data(prefix...; run_index=nothing) L"\widetilde{u}_{\|\|i}",L"\widetilde{u}_{\|\|i}^{sym}",L"\varepsilon(\widetilde{u}_{\|\|i})","upar") compare_moments_symbolic_test(run_name_label,parallel_pressure,ppar_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, L"\widetilde{p}_{\|\|i}",L"\widetilde{p}_{\|\|i}^{sym}",L"\varepsilon(\widetilde{p}_{\|\|i})","ppar") + compare_moments_symbolic_test(run_name,perpendicular_pressure,pperp_sym,"ion",z,r,time,nz_global,nr_global,ntime, + L"\widetilde{p}_{\perp i}",L"\widetilde{p}_{\perp i}^{sym}",L"\varepsilon(\widetilde{p}_{\perp i})","pperp") compare_charged_pdf_symbolic_test(run_name_label,manufactured_solns_list,"ion", L"\widetilde{f}_i",L"\widetilde{f}^{sym}_i",L"\varepsilon(\widetilde{f}_i)","pdf") diff --git a/src/time_advance.jl b/src/time_advance.jl index 4892dd87e..46687ca34 100644 --- a/src/time_advance.jl +++ b/src/time_advance.jl @@ -14,9 +14,8 @@ using ..debugging using ..file_io: write_data_to_ascii, write_moments_data_to_binary, write_dfns_data_to_binary, debug_dump using ..looping using ..moment_kinetics_structs: scratch_pdf -using ..moment_kinetics_structs: scratch_pdf using ..velocity_moments: update_moments!, update_moments_neutral!, reset_moments_status! -using ..velocity_moments: update_density!, update_upar!, update_ppar!, update_qpar! +using ..velocity_moments: update_density!, update_upar!, update_ppar!, update_pperp!, update_qpar! using ..velocity_moments: update_neutral_density!, update_neutral_qz! using ..velocity_moments: update_neutral_uzeta!, update_neutral_uz!, update_neutral_ur! using ..velocity_moments: update_neutral_pzeta!, update_neutral_pz!, update_neutral_pr! @@ -656,6 +655,7 @@ function setup_scratch_arrays(moments, pdf_charged_in, pdf_neutral_in, n_rk_stag density_array = allocate_shared_float(moment_dims...) upar_array = allocate_shared_float(moment_dims...) ppar_array = allocate_shared_float(moment_dims...) + pperp_array = allocate_shared_float(moment_dims...) temp_z_s_array = allocate_shared_float(moment_dims...) pdf_neutral_array = allocate_shared_float(pdf_neutral_dims...) @@ -665,7 +665,7 @@ function setup_scratch_arrays(moments, pdf_charged_in, pdf_neutral_in, n_rk_stag scratch[istage] = scratch_pdf(pdf_array, density_array, upar_array, - ppar_array, temp_z_s_array, + ppar_array, pperp_array, temp_z_s_array, pdf_neutral_array, density_neutral_array, uz_neutral_array, pz_neutral_array) @serial_region begin @@ -673,6 +673,7 @@ function setup_scratch_arrays(moments, pdf_charged_in, pdf_neutral_in, n_rk_stag scratch[istage].density .= moments.charged.dens scratch[istage].upar .= moments.charged.upar scratch[istage].ppar .= moments.charged.ppar + scratch[istage].pperp .= moments.charged.pperp scratch[istage].pdf_neutral .= pdf_neutral_in scratch[istage].density_neutral .= moments.neutral.dens @@ -1234,7 +1235,6 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v vpa) end end - # update remaining velocity moments that are calculable from the evolved pdf update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composition) # update the thermal speed @@ -1406,7 +1406,8 @@ function update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composi update_ppar!(new_scratch.ppar, moments.charged.ppar_updated, new_scratch.density, new_scratch.upar, new_scratch.pdf, vpa, vperp, z, r, composition, moments.evolve_density, moments.evolve_upar) - end + end + update_pperp!(new_scratch.pperp, new_scratch.pdf, vpa, vperp, z, r, composition) end """ @@ -1501,6 +1502,7 @@ function ssp_rk!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, moments.charged.dens[iz,ir,is] = final_scratch.density[iz,ir,is] moments.charged.upar[iz,ir,is] = final_scratch.upar[iz,ir,is] moments.charged.ppar[iz,ir,is] = final_scratch.ppar[iz,ir,is] + moments.charged.pperp[iz,ir,is] = final_scratch.pperp[iz,ir,is] end if composition.n_neutral_species > 0 # No need to synchronize here as we only change neutral quantities and previous diff --git a/src/velocity_moments.jl b/src/velocity_moments.jl index 45e5aa9ae..45e71cf72 100644 --- a/src/velocity_moments.jl +++ b/src/velocity_moments.jl @@ -10,6 +10,7 @@ export update_moments! export update_density! export update_upar! export update_ppar! +export update_pperp! export update_qpar! export reset_moments_status! export moments_chrg_substruct, moments_ntrl_substruct @@ -64,6 +65,8 @@ struct moments_charged_substruct # sets/uses the value for the same subset of species. This means ppar_update does # not need to be a shared memory array. ppar_updated::Vector{Bool} + # this is the perpendicular pressure + pperp::MPISharedArray{mk_float,3} # this is the parallel heat flux qpar::MPISharedArray{mk_float,3} # flag that keeps track of whether or not qpar needs updating before use @@ -188,6 +191,8 @@ function create_moments_charged(nz, nr, n_species, evolve_density, evolve_upar, # allocate array of Bools that indicate if the parallel pressure is updated for each species parallel_pressure_updated = allocate_bool(n_species) parallel_pressure_updated .= false + # allocate array used for the perpendicular pressure + perpendicular_pressure = allocate_shared_float(nz, nr, n_species) # allocate array used for the parallel flow parallel_heat_flux = allocate_shared_float(nz, nr, n_species) # allocate array of Bools that indicate if the parallel flow is updated for each species @@ -253,7 +258,7 @@ function create_moments_charged(nz, nr, n_species, evolve_density, evolve_upar, # return struct containing arrays needed to update moments return moments_charged_substruct(density, density_updated, parallel_flow, - parallel_flow_updated, parallel_pressure, parallel_pressure_updated, + parallel_flow_updated, parallel_pressure, parallel_pressure_updated,perpendicular_pressure, parallel_heat_flux, parallel_heat_flux_updated, thermal_speed, v_norm_fac, ddens_dz_upwind, d2dens_dz2, dupar_dz, dupar_dz_upwind, d2upar_dz2, dppar_dz, dppar_dz_upwind, d2ppar_dz2, dqpar_dz, dvth_dz) @@ -590,6 +595,34 @@ function get_ppar(ff, vpa, vperp, upar) return norm_fac*integrate_over_vspace(@view(ff[:,:]), vpa.scratch, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) end +function update_pperp!(pperp, pdf, vpa, vperp, z, r, composition) + @boundscheck composition.n_ion_species == size(pperp,3) || throw(BoundsError(pperp)) + @boundscheck r.n == size(pperp,2) || throw(BoundsError(pperp)) + @boundscheck z.n == size(pperp,1) || throw(BoundsError(pperp)) + + begin_s_r_z_region() + + @loop_s is begin + @views update_pperp_species!(pperp[:,:,is], pdf[:,:,:,:,is], vpa, vperp, z, r) + end +end + +""" +calculate the updated perpendicular pressure (pperp) for a given species +""" +function update_pperp_species!(pperp, ff, vpa, vperp, z, r) + @boundscheck vpa.n == size(ff, 1) || throw(BoundsError(ff)) + @boundscheck vperp.n == size(ff, 2) || throw(BoundsError(ff)) + @boundscheck z.n == size(ff, 3) || throw(BoundsError(ff)) + @boundscheck r.n == size(ff, 4) || throw(BoundsError(ff)) + @boundscheck z.n == size(pperp, 1) || throw(BoundsError(pperp)) + @boundscheck r.n == size(pperp, 2) || throw(BoundsError(pperp)) + @loop_r_z ir iz begin + pperp[iz,ir] = get_pperp(@view(ff[:,:,iz,ir]), vpa, vperp) + end + return nothing +end + function get_pperp(ff, vpa, vperp) norm_fac = 0.5 # normalise to m_s N_e c_s^2 #norm_fac = 1.0 # normalise to 0.5 m_s N_e c_s^2 = N_e T_s From 82de4c4f911378a05d3f7807b0e2c1d47e1eccd0 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Thu, 18 May 2023 10:27:57 +0100 Subject: [PATCH 11/34] Attempt to cherry-pick the features from 6bd406e. Some bugs remain where pperp (and thus vthi) are computed incorrectly according to the manufactured solutions test in 1V and 2V. --- src/initial_conditions.jl | 11 +++------ src/manufactured_solns.jl | 40 +++++++++++++++++++++++++------ src/post_processing.jl | 10 ++++++-- src/time_advance.jl | 11 +++++---- src/velocity_moments.jl | 30 ++++++++++++++++++++++++ test_velocity_integrals.jl | 48 +++++++++++++++++++++++++++++++++----- 6 files changed, 122 insertions(+), 28 deletions(-) diff --git a/src/initial_conditions.jl b/src/initial_conditions.jl index 92cb07ed8..56fc0444d 100644 --- a/src/initial_conditions.jl +++ b/src/initial_conditions.jl @@ -31,7 +31,7 @@ using ..velocity_moments: create_moments_charged, create_moments_neutral, update using ..velocity_moments: moments_charged_substruct, moments_neutral_substruct using ..velocity_moments: update_neutral_density!, update_neutral_pz!, update_neutral_pr!, update_neutral_pzeta! using ..velocity_moments: update_neutral_uz!, update_neutral_ur!, update_neutral_uzeta!, update_neutral_qz! -using ..velocity_moments: update_ppar!, update_upar!, update_density!, update_pperp!, reset_moments_status! +using ..velocity_moments: update_ppar!, update_upar!, update_density!, update_pperp!, update_vth!, reset_moments_status! using ..manufactured_solns: manufactured_solutions @@ -899,7 +899,7 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa, manufactured_solns_input) manufactured_solns_list = manufactured_solutions(manufactured_solns_input, r.L, z.L, r.bc, z.bc, geometry, composition, - species, r.n) + 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 @@ -931,12 +931,7 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa, moments.charged.vth, pdf.charged.norm, vpa, vperp, z, r, composition, moments.evolve_density, moments.evolve_upar, moments.evolve_ppar) - # convert from particle particle flux to parallel flow - begin_s_r_z_region() - @loop_s_r_z is ir iz begin - # update the thermal speed - moments.charged.vth[iz,ir,is] = sqrt(2.0*moments.charged.ppar[iz,ir,is]/moments.charged.dens[iz,ir,is]) - end + update_vth!(moments.charged.vth, moments.charged.ppar, moments.charged.pperp, moments.charged.dens, vperp, z, r, composition) if n_neutral_species > 0 begin_sn_r_z_region() diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 4cf0adcea..ad00ec145 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -240,12 +240,35 @@ using IfElse end # ion perpendicular pressure symbolic function - function pperpi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + function pperpi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species,nvperp) densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) - pperpi = densi # simple vperp^2 dependence of dfni - return pperpi + normfac = 0.5 # if pressure normalised to 0.5* nref * Tref = mref cref^2 + #normfac = 1.0 # if pressure normalised to nref*Tref + if nvperp > 1 + pperpi = densi # simple vperp^2 dependence of dfni + else + pperpi = 0.0 # marginalised model has nvperp = 1, vperp[1] = 0 + end + return pperpi*normfac end - function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,manufactured_solns_input) + + # ion thermal speed symbolic function + function vthi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species,nvperp) + densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + ppari = ppari_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + pperpi = pperpi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species,nvperp) + isotropic_pressure = (1.0/3.0)*(ppari + 2.0*pperpi) + normfac = 2.0 # if pressure normalised to 0.5* nref * Tref = mref cref^2 + #normfac = 1.0 # if pressure normalised to nref*Tref + if nvperp > 1 + vthi = sqrt(normfac*isotropic_pressure/densi) # thermal speed definition of 2V model + else + vthi = sqrt(normfac*ppari/densi) # thermal speed definition of 1V model + end + return vthi + end + + function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input) if z_bc == "periodic" jpari_into_LHS_wall_sym = 0.0 elseif z_bc == "wall" @@ -352,7 +375,7 @@ using IfElse end function manufactured_solutions(manufactured_solns_input, Lr, Lz, r_bc, z_bc, - geometry, composition, species, nr) + geometry, composition, species, nr, nvperp) charged_species = species.charged[1] if composition.n_neutral_species > 0 neutral_species = species.neutral[1] @@ -367,7 +390,9 @@ using IfElse ppari = ppari_sym(Lr, Lz, r_bc, z_bc, composition, manufactured_solns_input, charged_species) pperpi = pperpi_sym(Lr, Lz, r_bc, z_bc, composition, manufactured_solns_input, - charged_species) + charged_species, nvperp) + vthi = vthi_sym(Lr, Lz, r_bc, z_bc, composition, manufactured_solns_input, + charged_species, nvperp) dfni = dfni_sym(Lr, Lz, r_bc, z_bc, composition, geometry, nr, manufactured_solns_input, charged_species) @@ -382,6 +407,7 @@ using IfElse upari_func = build_function(upari, z, r, t, expression=Val{false}) ppari_func = build_function(ppari, z, r, t, expression=Val{false}) pperpi_func = build_function(pperpi, z, r, t, expression=Val{false}) + vthi_func = build_function(vthi, z, r, t, expression=Val{false}) densn_func = build_function(densn, z, r, t, expression=Val{false}) dfni_func = build_function(dfni, vpa, vperp, z, r, t, expression=Val{false}) dfnn_func = build_function(dfnn, vz, vr, vzeta, z, r, t, expression=Val{false}) @@ -395,7 +421,7 @@ using IfElse manufactured_solns_list = (densi_func = densi_func, densn_func = densn_func, dfni_func = dfni_func, dfnn_func = dfnn_func, upari_func = upari_func, ppari_func = ppari_func, - pperpi_func = pperpi_func) + pperpi_func = pperpi_func, vthi_func = vthi_func) return manufactured_solns_list end diff --git a/src/post_processing.jl b/src/post_processing.jl index 69287f396..1b701486d 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -1017,6 +1017,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) density = density[1] parallel_flow = parallel_flow[1] parallel_pressure = parallel_pressure[1] + perpendicular_pressure = perpendicular_pressure[1] parallel_heat_flux = parallel_heat_flux[1] thermal_speed = thermal_speed[1] time = time[1] @@ -1100,12 +1101,13 @@ function analyze_and_plot_data(prefix...; run_index=nothing) manufactured_solns_list = manufactured_solutions(manufactured_solns_input, Lr_in, z_global.L, r_global.bc, z_global.bc, geometry, - composition, species, r_global.n) + composition, species, r_global.n, vperp.n) dfni_func = manufactured_solns_list.dfni_func densi_func = manufactured_solns_list.densi_func upari_func = manufactured_solns_list.upari_func ppari_func = manufactured_solns_list.ppari_func pperpi_func = manufactured_solns_list.pperpi_func + vthi_func = manufactured_solns_list.vthi_func dfnn_func = manufactured_solns_list.dfnn_func densn_func = manufactured_solns_list.densn_func manufactured_E_fields = @@ -1145,6 +1147,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) upar_sym = copy(density[:,:,:,:]) ppar_sym = copy(density[:,:,:,:]) pperp_sym = copy(density[:,:,:,:]) + vthi_sym = copy(density[:,:,:,:]) is = 1 for it in 1:ntime for ir in 1:r_global.n @@ -1153,6 +1156,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) upar_sym[iz,ir,is,it] = upari_func(z_global.grid[iz],r_global.grid[ir],time[it]) ppar_sym[iz,ir,is,it] = ppari_func(z_global.grid[iz],r_global.grid[ir],time[it]) pperp_sym[iz,ir,is,it] = pperpi_func(z_global.grid[iz],r_global.grid[ir],time[it]) + vthi_sym[iz,ir,is,it] = vthi_func(z_global.grid[iz],r_global.grid[ir],time[it]) end end end @@ -1162,8 +1166,10 @@ function analyze_and_plot_data(prefix...; run_index=nothing) L"\widetilde{u}_{\|\|i}",L"\widetilde{u}_{\|\|i}^{sym}",L"\varepsilon(\widetilde{u}_{\|\|i})","upar") compare_moments_symbolic_test(run_name_label,parallel_pressure,ppar_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, L"\widetilde{p}_{\|\|i}",L"\widetilde{p}_{\|\|i}^{sym}",L"\varepsilon(\widetilde{p}_{\|\|i})","ppar") - compare_moments_symbolic_test(run_name,perpendicular_pressure,pperp_sym,"ion",z,r,time,nz_global,nr_global,ntime, + compare_moments_symbolic_test(run_name,perpendicular_pressure,pperp_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, L"\widetilde{p}_{\perp i}",L"\widetilde{p}_{\perp i}^{sym}",L"\varepsilon(\widetilde{p}_{\perp i})","pperp") + compare_moments_symbolic_test(run_name,thermal_speed,vthi_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + L"\widetilde{v}_{th,i}",L"\widetilde{v}_{th,i}^{sym}",L"\varepsilon(\widetilde{v}_{th,i})","vthi") compare_charged_pdf_symbolic_test(run_name_label,manufactured_solns_list,"ion", L"\widetilde{f}_i",L"\widetilde{f}^{sym}_i",L"\varepsilon(\widetilde{f}_i)","pdf") diff --git a/src/time_advance.jl b/src/time_advance.jl index 46687ca34..e357b4025 100644 --- a/src/time_advance.jl +++ b/src/time_advance.jl @@ -15,7 +15,7 @@ using ..file_io: write_data_to_ascii, write_moments_data_to_binary, write_dfns_d using ..looping using ..moment_kinetics_structs: scratch_pdf using ..velocity_moments: update_moments!, update_moments_neutral!, reset_moments_status! -using ..velocity_moments: update_density!, update_upar!, update_ppar!, update_pperp!, update_qpar! +using ..velocity_moments: update_density!, update_upar!, update_ppar!, update_pperp!, update_qpar!, update_vth! using ..velocity_moments: update_neutral_density!, update_neutral_qz! using ..velocity_moments: update_neutral_uzeta!, update_neutral_uz!, update_neutral_ur! using ..velocity_moments: update_neutral_pzeta!, update_neutral_pz!, update_neutral_pr! @@ -1204,7 +1204,9 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v ## # update the charged particle distribution and moments ## - + # MRH here we seem to have duplicate arrays for storing n, u||, p||, etc, but not for vth + # MRH 'scratch' is for the multiple stages of time advanced quantities, but 'moments' can be updated directly at each stage + # MRH in the standard drift-kinetic model. Consider taking moment quantities out of scratch for clarity. @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin new_scratch.pdf[ivpa,ivperp,iz,ir,is] = rk_coefs[1]*pdf.charged.norm[ivpa,ivperp,iz,ir,is] + rk_coefs[2]*old_scratch.pdf[ivpa,ivperp,iz,ir,is] + rk_coefs[3]*new_scratch.pdf[ivpa,ivperp,iz,ir,is] end @@ -1240,9 +1242,7 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v # update the thermal speed begin_s_r_z_region() try #below block causes DomainError if ppar < 0 or density, so exit cleanly if possible - @loop_s_r_z is ir iz begin - moments.charged.vth[iz,ir,is] = sqrt(2.0*new_scratch.ppar[iz,ir,is]/new_scratch.density[iz,ir,is]) - end + update_vth!(moments.charged.vth, new_scratch.ppar, new_scratch.pperp, new_scratch.density, vperp, z, r, composition) catch e if global_size[] > 1 println("ERROR: error calculating vth in time_advance.jl") @@ -1452,6 +1452,7 @@ function ssp_rk!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, first_scratch.density[iz,ir,is] = moments.charged.dens[iz,ir,is] first_scratch.upar[iz,ir,is] = moments.charged.upar[iz,ir,is] first_scratch.ppar[iz,ir,is] = moments.charged.ppar[iz,ir,is] + first_scratch.pperp[iz,ir,is] = moments.charged.pperp[iz,ir,is] end if composition.n_neutral_species > 0 diff --git a/src/velocity_moments.jl b/src/velocity_moments.jl index 45e71cf72..da41a1fea 100644 --- a/src/velocity_moments.jl +++ b/src/velocity_moments.jl @@ -12,6 +12,7 @@ export update_upar! export update_ppar! export update_pperp! export update_qpar! +export update_vth! export reset_moments_status! export moments_chrg_substruct, moments_ntrl_substruct export update_neutral_density! @@ -28,6 +29,7 @@ export get_density export get_upar export get_ppar export get_pperp +export get_pressure using ..type_definitions: mk_float using ..array_allocation: allocate_shared_float, allocate_bool, allocate_float @@ -629,6 +631,34 @@ function get_pperp(ff, vpa, vperp) return norm_fac*integrate_over_vspace(@view(ff[:,:]), vpa.grid, 0, vpa.wgts, vperp.grid, 2, vperp.wgts) end +function update_vth!(vth, ppar, pperp, dens, vperp, z, r, composition) + @boundscheck composition.n_ion_species == size(vth,3) || throw(BoundsError(vth)) + @boundscheck r.n == size(vth,2) || throw(BoundsError(vth)) + @boundscheck z.n == size(vth,1) || throw(BoundsError(vth)) + + begin_s_r_z_region() + normfac = 2.0 # if ppar normalised to 2*nref Tref = mref cref^2 + #normfac = 1.0 # if ppar normalised to nref Tref = 0.5 * mref cref^2 + if vperp.n > 1 #2V definition + @loop_s_r_z is ir iz begin + piso = get_pressure(ppar[iz,ir,is],pperp[iz,ir,is]) + vth[iz,ir,is] = sqrt(normfac*piso/dens[iz,ir,is]) + end + else #1V definition + @loop_s_r_z is ir iz begin + vth[iz,ir,is] = sqrt(normfac*ppar[iz,ir,is]/dens[iz,ir,is]) + end + end +end + +""" +compute the isotropic pressure from the already computed ppar and pperp +""" +function get_pressure(ppar::mk_float,pperp::mk_float) + pres = (1.0/3.0)*(ppar + 2.0*pperp) + return pres +end + """ NB: the incoming pdf is the normalized pdf """ diff --git a/test_velocity_integrals.jl b/test_velocity_integrals.jl index bb39123b1..24f674242 100644 --- a/test_velocity_integrals.jl +++ b/test_velocity_integrals.jl @@ -11,7 +11,7 @@ if abspath(PROGRAM_FILE) == @__FILE__ import moment_kinetics using moment_kinetics.input_structs: grid_input, advection_input using moment_kinetics.coordinates: define_coordinate - using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_pperp + using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_pperp, get_pressure using moment_kinetics.array_allocation: allocate_float # define inputs needed for the test @@ -30,29 +30,37 @@ if abspath(PROGRAM_FILE) == @__FILE__ comm = MPI.COMM_NULL # create the 'input' struct containing input info needed to create a # coordinate - vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, + vr_input = grid_input("vperp", 1, 1, 1, + nrank, irank, 1.0, discretization, fd_option, bc, adv_input,comm) + vz_input = grid_input("vpa", ngrid, nelement_global, nelement_local, + nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) + vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm) # create the coordinate struct 'x' println("made inputs") - println("vpa: ngrid: ",ngrid," nelement: ",nelement_local) - println("vperp: ngrid: ",ngrid," nelement: ",nelement_local) + println("vpa: ngrid: ",ngrid," nelement: ",nelement_local, " Lvpa: ",Lvpa) + println("vperp: ngrid: ",ngrid," nelement: ",nelement_local, " Lvperp: ",Lvperp) vpa, vpa_spectral = define_coordinate(vpa_input) vperp, vperp_spectral = define_coordinate(vperp_input) + vz, vz_spectral = define_coordinate(vz_input) + vr, vr_spectral = define_coordinate(vr_input) dfn = allocate_float(vpa.n,vperp.n) + dfn1D = allocate_float(vz.n, vr.n) function pressure(ppar,pperp) pres = (1.0/3.0)*(ppar + 2.0*pperp) return pres end + # 2D isotropic Maxwellian test # assign a known isotropic Maxwellian distribution in normalised units dens = 3.0/4.0 upar = 2.0/3.0 ppar = 2.0/3.0 pperp = 2.0/3.0 - pres = pressure(ppar,pperp) + pres = get_pressure(ppar,pperp) mass = 1.0 vth = sqrt(2.0*pres/(dens*mass)) for ivperp in 1:vperp.n @@ -72,12 +80,40 @@ if abspath(PROGRAM_FILE) == @__FILE__ pres_test = pressure(ppar_test,pperp_test) # output test results println("") - println("Isotropic Maxwellian") + println("Isotropic 2D Maxwellian") println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) println("pres_test: ", pres_test, " pres: ", pres, " error: ", abs(pres_test-pres)) println("") + ################### + # 1D Maxwellian test + + dens = 3.0/4.0 + upar = 2.0/3.0 + ppar = 2.0/3.0 + mass = 1.0 + vth = sqrt(ppar/(dens*mass)) + for ivz in 1:vz.n + for ivr in 1:vr.n + vz_val = vz.grid[ivz] + dfn1D[ivz,ivr] = (dens/vth)*exp( - ((vz_val-upar)^2)/vth^2 ) + end + end + dens_test = get_density(dfn1D,vz,vr) + upar_test = get_upar(dfn1D,vz,vr,dens_test) + ppar_test = get_ppar(dfn1D,vz,vr,upar_test) + # output test results + println("") + println("1D Maxwellian") + println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) + println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) + println("ppar_test: ", ppar_test, " ppar: ", ppar, " error: ", abs(ppar_test-ppar)) + println("") + + ################### + # biMaxwellian test + # assign a known biMaxwellian distribution in normalised units dens = 3.0/4.0 upar = 2.0/3.0 From ed979adc05a343dc134bd8ae92f79ebd63070d11 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 24 Oct 2023 10:57:40 +0100 Subject: [PATCH 12/34] Include factor of 2 for expected vth in test_velocity_integrals.jl Needs the factor of 2 to be consistent with normalisations currently used. --- test_velocity_integrals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test_velocity_integrals.jl b/test_velocity_integrals.jl index 24f674242..9bcee9325 100644 --- a/test_velocity_integrals.jl +++ b/test_velocity_integrals.jl @@ -93,7 +93,7 @@ if abspath(PROGRAM_FILE) == @__FILE__ upar = 2.0/3.0 ppar = 2.0/3.0 mass = 1.0 - vth = sqrt(ppar/(dens*mass)) + vth = sqrt(2.0*ppar/(dens*mass)) for ivz in 1:vz.n for ivr in 1:vr.n vz_val = vz.grid[ivz] From 5bd78ecac315ba9c066697d9ccf55e7a11155daf Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 24 Oct 2023 10:58:57 +0100 Subject: [PATCH 13/34] Increase Lvpa and Lvperp in test_velocity_integrals.jl Need larger boxes to make moments accurate to roughly machine precision. --- test_velocity_integrals.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test_velocity_integrals.jl b/test_velocity_integrals.jl index 9bcee9325..25a023539 100644 --- a/test_velocity_integrals.jl +++ b/test_velocity_integrals.jl @@ -18,8 +18,8 @@ if abspath(PROGRAM_FILE) == @__FILE__ ngrid = 17 #number of points per element nelement_local = 20 # number of elements per rank nelement_global = nelement_local # total number of elements - Lvpa = 12.0 #physical box size in reference units - Lvperp = 6.0 #physical box size in reference units + Lvpa = 18.0 #physical box size in reference units + Lvperp = 9.0 #physical box size in reference units bc = "" #not required to take a particular value, not used # fd_option and adv_input not actually used so given values unimportant discretization = "chebyshev_pseudospectral" From 4e7dccbca233d8dbc5a9bbf5cd14195b963e5b0b Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Tue, 24 Oct 2023 16:13:01 +0100 Subject: [PATCH 14/34] Correct file_io so that perpendicular pressure is correctly saved to the output binary. --- src/file_io.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/file_io.jl b/src/file_io.jl index 94f3e4123..e6e132918 100644 --- a/src/file_io.jl +++ b/src/file_io.jl @@ -785,8 +785,8 @@ function reopen_moments_io(file_info) end return io_moments_info(fid, getvar("time"), getvar("phi"), getvar("Er"), getvar("Ez"), getvar("density"), getvar("parallel_flow"), - getvar("parallel_pressure"), getvar("parallel_heat_flux"), - getvar("perpendicular_pressure"), + getvar("parallel_pressure"), getvar("perpendicular_pressure"), + getvar("parallel_heat_flux"), getvar("thermal_speed"), getvar("density_neutral"), getvar("uz_neutral"), getvar("pz_neutral"), getvar("qz_neutral"), getvar("thermal_speed_neutral"), From b7a76209f48959c5c8ba10fb5d241fe0e969a049 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Tue, 24 Oct 2023 19:08:32 +0100 Subject: [PATCH 15/34] Fix file name for output plots --- src/post_processing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/post_processing.jl b/src/post_processing.jl index 1b701486d..a625dec01 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -1166,9 +1166,9 @@ function analyze_and_plot_data(prefix...; run_index=nothing) L"\widetilde{u}_{\|\|i}",L"\widetilde{u}_{\|\|i}^{sym}",L"\varepsilon(\widetilde{u}_{\|\|i})","upar") compare_moments_symbolic_test(run_name_label,parallel_pressure,ppar_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, L"\widetilde{p}_{\|\|i}",L"\widetilde{p}_{\|\|i}^{sym}",L"\varepsilon(\widetilde{p}_{\|\|i})","ppar") - compare_moments_symbolic_test(run_name,perpendicular_pressure,pperp_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + compare_moments_symbolic_test(run_name_label,perpendicular_pressure,pperp_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, L"\widetilde{p}_{\perp i}",L"\widetilde{p}_{\perp i}^{sym}",L"\varepsilon(\widetilde{p}_{\perp i})","pperp") - compare_moments_symbolic_test(run_name,thermal_speed,vthi_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + compare_moments_symbolic_test(run_name_label,thermal_speed,vthi_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, L"\widetilde{v}_{th,i}",L"\widetilde{v}_{th,i}^{sym}",L"\varepsilon(\widetilde{v}_{th,i})","vthi") compare_charged_pdf_symbolic_test(run_name_label,manufactured_solns_list,"ion", From 1387ba3b718f1fb19ca3e9f742059447f3a409d6 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 24 Oct 2023 19:26:13 +0100 Subject: [PATCH 16/34] Better suggestion for putting .julia locally within the repo Using just '.' means the contents of '.julia' get put directly in the top level of the repo - really want to create a '.julia' subdirectory within the repo. --- machines/machine_setup.sh | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/machines/machine_setup.sh b/machines/machine_setup.sh index c9e5cf0da..7114848d9 100755 --- a/machines/machine_setup.sh +++ b/machines/machine_setup.sh @@ -203,10 +203,11 @@ echo JULIA_DIRECTORY=$JULIA_DEPOT_PATH echo "It can be useful or necessary to set a non-default location for the " echo ".julia directory. Leave this empty if the default location is OK." -echo "Enter the current directory '.' to isolate the julia used for this " -echo "instance of moment_kinetics - this might be useful to ensure a 'clean'" -echo "install or to check whether some error is related to conflicting or " -echo "corrupted dependencies or cached precompilation files, etc." +echo "Enter a name for a subdirectory of the current directory, e.g. " +echo "'.julia', to isolate the julia used for this instance of " +echo "moment_kinetics - this might be useful to ensure a 'clean' install or " +echo "to check whether some error is related to conflicting or corrupted " +echo "dependencies or cached precompilation files, etc." echo "Enter location that should be used for the .julia directory [$JULIA_DIRECTORY]:" # Use '-e' option to get path auto-completion read -e -p "> " input From 4549acf0336e5f2f11c17b32b935829c1eaf4e90 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 25 Oct 2023 12:51:31 +0100 Subject: [PATCH 17/34] Fix NetCDF array loading A recent update to NCDatasets.jl changed the behaviour of indexing like `var[:]` when loading a NetCDF variable. Now `var[:]` returns a 1d view. Fix by using `copy()` to convert the `NCDataset` to an `Array`. --- src/load_data.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/load_data.jl b/src/load_data.jl index 5507965b1..d6087c6f9 100644 --- a/src/load_data.jl +++ b/src/load_data.jl @@ -100,7 +100,7 @@ function load_variable(file_or_group::NCDataset, name::String) if size(file_or_group[name].var) == () var = file_or_group[name].var[] else - var = file_or_group[name].var[:] + var = copy(file_or_group[name].var) end if isa(var, Char) var = (var == Char(true)) From 88ad6d0610379268bd1ca534d5c4ea3ec2bcefe7 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 25 Oct 2023 13:08:22 +0100 Subject: [PATCH 18/34] Attempt to port manufactured solutions test for Krook operator to new implementation. Suited only for 1D1V with standard drift kinetics for now. Errors suggests an order unity factor difference somewhere in the two implementations. --- src/manufactured_solns.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index ad00ec145..7373bc9d0 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -456,6 +456,9 @@ using IfElse # ion manufactured solutions densi = densi_sym(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc, composition, manufactured_solns_input, charged_species) + upari = upari_sym(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc, composition, geometry, r_coord.n, manufactured_solns_input, charged_species) + vthi = vthi_sym(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc, composition, manufactured_solns_input, + charged_species, vperp_coord.n) dfni = dfni_sym(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc, composition, geometry, r_coord.n, manufactured_solns_input, charged_species) #dfni in vr vz vzeta coordinates @@ -519,6 +522,19 @@ using IfElse if num_diss_params.z_dissipation_coefficient > 0.0 && include_num_diss_in_MMS Si += - num_diss_params.z_dissipation_coefficient*Dz(Dz(dfni)) end + nu_krook = collisions.krook_collision_frequency_prefactor + if nu_krook > 0.0 + tempi = vthi^2 + nu_ii = nu_krook * densi * (tempi^(-3.0/2.0)) + if vperp_coord.n > 1 + pvth = 3 + else + pvth = 1 + end + FMaxwellian = (densi/vthi^pvth)*exp( -( ( vpa-upari)^2 + vperp^2 )/vthi^2) + Si += -nu_krook*(FMaxwellian - dfni) + end + Source_i = expand_derivatives(Si) From ba4d733dfbd8af9f7929ad6095884d391cfafbb3 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 25 Oct 2023 13:09:01 +0100 Subject: [PATCH 19/34] Put rfac preceding the numerical dissipation term in r to prevent usage when r.n = 1. --- src/manufactured_solns.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 7373bc9d0..7a88d640a 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -517,7 +517,7 @@ using IfElse Si += - num_diss_params.vpa_dissipation_coefficient*Dvpa(Dvpa(dfni)) end if num_diss_params.r_dissipation_coefficient > 0.0 && include_num_diss_in_MMS - Si += - num_diss_params.r_dissipation_coefficient*Dr(Dr(dfni)) + Si += - rfac*num_diss_params.r_dissipation_coefficient*Dr(Dr(dfni)) end if num_diss_params.z_dissipation_coefficient > 0.0 && include_num_diss_in_MMS Si += - num_diss_params.z_dissipation_coefficient*Dz(Dz(dfni)) From e37a0a57024cdf6d97d748c0fcf104f887552df6 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 09:49:03 +0000 Subject: [PATCH 20/34] Clean up comment --- src/time_advance.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_advance.jl b/src/time_advance.jl index e357b4025..6802646cf 100644 --- a/src/time_advance.jl +++ b/src/time_advance.jl @@ -1668,7 +1668,7 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, collisions, dt) end - # Add a for Krook collision operoator for ions + # Add Krook collision operator for ions if advance.krook_collisions krook_collisions!(fvec_out.pdf, fvec_in, moments, composition, collisions, vperp, vpa, dt) From 68235462ea2ee98a0bca2c450d8e4ad293ecca8b Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 09:53:01 +0000 Subject: [PATCH 21/34] Modifications to port Krook operator manufactured solutions test to the present branch for the 1D1V case. It is noted that the manufactured solutions test passes when the Krook collision frequency is a constant independent of position. With this in mind, I have modified how the Krook operator is specified in the input file. Now one should specify "krook_collisions_option" ("default" -- use the frequency from the reference parameters and use the spatially dependent factor n/vth^3 -- or "manual" -- use a frequency set in the input file). The frequency can be set in the input file with the "nuii_krook" parameter, if the krook_collisions_option = "manual". The tests and previous functionality are preserved. --- src/input_structs.jl | 2 ++ src/krook_collisions.jl | 16 ++++++++++------ src/manufactured_solns.jl | 8 ++++++-- src/moment_kinetics_input.jl | 11 +++++++---- test/Krook_collisions_tests.jl | 2 +- 5 files changed, 26 insertions(+), 13 deletions(-) diff --git a/src/input_structs.jl b/src/input_structs.jl index 285e42bcd..75b216ca3 100644 --- a/src/input_structs.jl +++ b/src/input_structs.jl @@ -302,6 +302,8 @@ mutable struct collisions_input constant_ionization_rate::Bool # Coulomb collision rate at the reference density and temperature krook_collision_frequency_prefactor::mk_float + # Coulomb collision rate at the reference density and temperature + krook_collisions_option::String end """ diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index 751d6caa1..db96a570c 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -48,14 +48,18 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v # Note: do not need 1/sqrt(pi) for the 'Maxwellian' term because the pdf is already # normalized by sqrt(pi) (see velocity_moments.integrate_over_vspace). - + if collisions.krook_collisions_option == "manual" + cfac = 0.0 + else # default option, collisions.krook_collisions_option == "default" + cfac = 1.0 + end if moments.evolve_ppar && moments.evolve_upar # Compared to evolve_upar version, grid is already normalized by vth, and multiply # through by vth, remembering pdf is already multiplied by vth @loop_s_r_z is ir iz begin n = fvec_in.density[iz,ir,is] T = (moments.charged.vth[iz,ir,is])^2 - nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -69,7 +73,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v n = fvec_in.density[iz,ir,is] vth = moments.charged.vth[iz,ir,is] T = vth^2 - nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -83,7 +87,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v n = fvec_in.density[iz,ir,is] vth = moments.charged.vth[iz,ir,is] T = vth^2 - nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -98,7 +102,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v n = fvec_in.density[iz,ir,is] vth = moments.charged.vth[iz,ir,is] T = vth^2 - nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -117,7 +121,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v else vth_prefactor = 1.0 / vth^3 end - nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 7a88d640a..2c2f7d05a 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -525,14 +525,18 @@ using IfElse nu_krook = collisions.krook_collision_frequency_prefactor if nu_krook > 0.0 tempi = vthi^2 - nu_ii = nu_krook * densi * (tempi^(-3.0/2.0)) + if collisions.krook_collisions_option == "manual" + nuii_krook = nu_krook + else # default option + nuii_krook = nu_krook * densi * tempi^(-1.5) + end if vperp_coord.n > 1 pvth = 3 else pvth = 1 end FMaxwellian = (densi/vthi^pvth)*exp( -( ( vpa-upari)^2 + vperp^2 )/vthi^2) - Si += -nu_krook*(FMaxwellian - dfni) + Si += -nuii_krook*(FMaxwellian - dfni) end diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index a15e1c733..f6bd3a93b 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -178,9 +178,12 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) collisions.charge_exchange = get(scan_input, "charge_exchange_frequency", 2.0*sqrt(species.charged[1].initial_temperature)) collisions.ionization = get(scan_input, "ionization_frequency", collisions.charge_exchange) collisions.constant_ionization_rate = get(scan_input, "constant_ionization_rate", false) - include_krook_collisions = get(scan_input, "krook_collisions", false) - if include_krook_collisions - collisions.krook_collision_frequency_prefactor = setup_krook_collisions(reference_parameters) + collisions.krook_collisions_option = get(scan_input, "krook_collisions_option", "reference_parameters") + nuii_krook_default = setup_krook_collisions(reference_parameters) + if collisions.krook_collisions_option == "default" + collisions.krook_collision_frequency_prefactor = nuii_krook_default + elseif collisions.krook_collisions_option == "manual" # get the frequency from the input file + collisions.krook_collision_frequency_prefactor = get(scan_input, "nuii_krook", nuii_krook_default) else collisions.krook_collision_frequency_prefactor = -1.0 end @@ -950,7 +953,7 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) constant_ionization_rate = false krook_collision_frequency_prefactor = -1.0 collisions = collisions_input(charge_exchange, ionization, constant_ionization_rate, - krook_collision_frequency_prefactor) + krook_collision_frequency_prefactor,"default") Bzed = 1.0 # magnetic field component along z Bmag = 1.0 # magnetic field strength diff --git a/test/Krook_collisions_tests.jl b/test/Krook_collisions_tests.jl index 1ebae0805..03c747a20 100644 --- a/test/Krook_collisions_tests.jl +++ b/test/Krook_collisions_tests.jl @@ -99,7 +99,7 @@ test_input_full_f = Dict("n_ion_species" => 1, "evolve_moments_parallel_flow" => false, "evolve_moments_parallel_pressure" => false, "evolve_moments_conservation" => true, - "krook_collisions" => true, + "krook_collisions_option" => "default", "T_e" => 1.0, "initial_density1" => 0.5, "initial_temperature1" => 1.0, From 09cb8fbb2386a8bd58614fd768e8fc2f09bde6f4 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 10:39:11 +0000 Subject: [PATCH 22/34] Allow 2V Krook operator providing that split moments are not used. --- src/krook_collisions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index db96a570c..7c68a7bbe 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -42,7 +42,7 @@ Currently Krook collisions function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, vperp, vpa, dt) begin_s_r_z_region() - if vperp.n > 1 + if vperp.n > 1 && (moments.evolve_density || moments.evolve_upar || moments.evolve_ppar) error("Krook collisions not implemented for 2V case yet") end From 58e7c8ee25c6e9a2d1947eae514af68460e38cb7 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 11:06:11 +0000 Subject: [PATCH 23/34] Correct update_moments! to use the update_vth! function so that the thermal_speed written to the HDF5/netCDF file is correct. --- src/velocity_moments.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/velocity_moments.jl b/src/velocity_moments.jl index da41a1fea..60ef582e4 100644 --- a/src/velocity_moments.jl +++ b/src/velocity_moments.jl @@ -394,10 +394,7 @@ function update_moments!(moments, ff, vpa, vperp, z, r, composition) moments.evolve_upar) moments.charged.ppar_updated[is] = true end - @loop_r_z ir iz begin - moments.charged.vth[iz,ir,is] = - sqrt(2*moments.charged.ppar[iz,ir,is]/moments.charged.dens[iz,ir,is]) - end + @views update_pperp_species!(moments.charged.pperp[:,:,is], ff[:,:,:,:,is], vpa, vperp, z, r) if moments.charged.qpar_updated[is] == false @views update_qpar_species!(moments.charged.qpar[:,:,is], moments.charged.dens[:,:,is], @@ -408,6 +405,7 @@ function update_moments!(moments, ff, vpa, vperp, z, r, composition) moments.charged.qpar_updated[is] = true end end + update_vth!(moments.charged.vth, moments.charged.ppar, moments.charged.pperp, moments.charged.dens, vperp, z, r, composition) return nothing end From 6132fc6d6f566dbba421c4d9a781c8957734bd02 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 11:08:22 +0000 Subject: [PATCH 24/34] Example input files for MMS test with Krook operator (spatially constant frequency coefficient). --- ...new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml | 96 +++++++++++++++++++ ..._new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml | 96 +++++++++++++++++++ 2 files changed, 192 insertions(+) create mode 100644 runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml create mode 100644 runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml new file mode 100644 index 000000000..5a329be7a --- /dev/null +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml @@ -0,0 +1,96 @@ +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 +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 +#krook_collisions = false +krook_collisions_option = "manual" +nuii_krook = 1.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 = 1 +vperp_nelement = 1 +vperp_L = 6.0 +vperp_bc = "periodic" +#vperp_discretization = "finite_difference" +vperp_discretization = "chebyshev_pseudospectral" + +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 = -1.0 +z_dissipation_coefficient = -1.0 +r_dissipation_coefficient = -1.0 diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml new file mode 100644 index 000000000..5173fa230 --- /dev/null +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml @@ -0,0 +1,96 @@ +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 +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 +#krook_collisions = false +krook_collisions_option = "manual" +nuii_krook = 1.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 = 8 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 17 +vperp_nelement = 8 +vperp_L = 6.0 +vperp_bc = "periodic" +#vperp_discretization = "finite_difference" +vperp_discretization = "chebyshev_pseudospectral" + +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 = -1.0 +z_dissipation_coefficient = -1.0 +r_dissipation_coefficient = -1.0 From 1885c7998f84ad02773c6ac19cb198bd5f12f1f5 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 11:42:23 +0000 Subject: [PATCH 25/34] Change input options for Krook so that default is to not use the Krook operator. Check-in tests still fail for undiagnosed reasons. --- src/moment_kinetics_input.jl | 2 +- test/Krook_collisions_tests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index f6bd3a93b..3ae2eb309 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -180,7 +180,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) collisions.constant_ionization_rate = get(scan_input, "constant_ionization_rate", false) collisions.krook_collisions_option = get(scan_input, "krook_collisions_option", "reference_parameters") nuii_krook_default = setup_krook_collisions(reference_parameters) - if collisions.krook_collisions_option == "default" + if collisions.krook_collisions_option == "reference_parameters" collisions.krook_collision_frequency_prefactor = nuii_krook_default elseif collisions.krook_collisions_option == "manual" # get the frequency from the input file collisions.krook_collision_frequency_prefactor = get(scan_input, "nuii_krook", nuii_krook_default) diff --git a/test/Krook_collisions_tests.jl b/test/Krook_collisions_tests.jl index 03c747a20..62e59c892 100644 --- a/test/Krook_collisions_tests.jl +++ b/test/Krook_collisions_tests.jl @@ -99,7 +99,7 @@ test_input_full_f = Dict("n_ion_species" => 1, "evolve_moments_parallel_flow" => false, "evolve_moments_parallel_pressure" => false, "evolve_moments_conservation" => true, - "krook_collisions_option" => "default", + "krook_collisions_option" => "reference_parameters", "T_e" => 1.0, "initial_density1" => 0.5, "initial_temperature1" => 1.0, From 5647180065ade485d9baf7cbf40e7de6c5a236d4 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 12:39:43 +0000 Subject: [PATCH 26/34] Fix figure saving for 1D plot in makie_post_processing --- src/makie_post_processing.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index a98855b2a..dbcd30da9 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -1647,10 +1647,10 @@ for dim ∈ one_dimension_combinations end function $function_name(run_info, var_name; is=1, data=nothing, - input=nothing, ax=nothing, label=nothing, - outfile=nothing, it=nothing, ir=nothing, iz=nothing, - ivperp=nothing, ivpa=nothing, ivzeta=nothing, - ivr=nothing, ivz=nothing, kwargs...) + input=nothing, fig=nothing, ax=nothing, + label=nothing, outfile=nothing, it=nothing, + ir=nothing, iz=nothing, ivperp=nothing, ivpa=nothing, + ivzeta=nothing, ivr=nothing, ivz=nothing, kwargs...) if input === nothing if run_info.dfns input = input_dict_dfns[var_name] @@ -1686,7 +1686,7 @@ for dim ∈ one_dimension_combinations if $idim !== nothing x = x[$idim] end - fig = plot_1d(x, data; label=label, ax=ax, kwargs...) + plot_1d(x, data; label=label, ax=ax, kwargs...) if input.show_element_boundaries && Symbol($dim_str) != :t && ax_was_nothing element_boundary_positions = run_info.$dim.grid[begin:run_info.$dim.ngrid-1:end] From 7ed8a1c803653c2c3140eba2cba770458d6dc028 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 12:40:35 +0000 Subject: [PATCH 27/34] Pass nvperp to manufactured solutions funcs in makie_post_processing --- src/makie_post_processing.jl | 72 ++++++++++++++++++++++-------------- 1 file changed, 45 insertions(+), 27 deletions(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index dbcd30da9..6cf302019 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -313,7 +313,12 @@ function makie_post_process(run_dir::Union{String,Tuple}, sound_wave_plots(run_info; plot_prefix=plot_prefix) - manufactured_solutions_analysis(run_info; plot_prefix=plot_prefix) + if all(ri === nothing for ri ∈ run_info_dfns) + nvperp = nothing + else + nvperp = maximum(ri.vperp.n_global for ri ∈ run_info_dfns if ri !== nothing) + end + manufactured_solutions_analysis(run_info; plot_prefix=plot_prefix, nvperp=nvperp) manufactured_solutions_analysis_dfns(run_info_dfns; plot_prefix=plot_prefix) return nothing @@ -5433,7 +5438,7 @@ Returns `variable`, `variable_sym`. """ function manufactured_solutions_get_field_and_field_sym(run_info, variable_name; it=nothing, ir=nothing, iz=nothing, ivperp=nothing, ivpa=nothing, ivzeta=nothing, - ivr=nothing, ivz=nothing) + ivr=nothing, ivz=nothing, nvperp) variable_name = Symbol(variable_name) @@ -5473,7 +5478,8 @@ function manufactured_solutions_get_field_and_field_sym(run_info, variable_name; 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.composition, run_info.species, run_info.r.n) + run_info.composition, run_info.species, run_info.r.n, + nvperp) end variable_func = manufactured_funcs[func_name_lookup[variable_name]] @@ -5560,7 +5566,8 @@ the label for the error (the difference between the computed and manufactured so If `io` is passed then error norms will be written to that file. """ function compare_moment_symbolic_test(run_info, plot_prefix, field_label, field_sym_label, - norm_label, variable_name; io=nothing, input=nothing) + norm_label, variable_name; io=nothing, + input=nothing, nvperp) println("Doing MMS analysis and making plots for $variable_name") flush(stdout) @@ -5570,7 +5577,7 @@ function compare_moment_symbolic_test(run_info, plot_prefix, field_label, field_ end field, field_sym = - manufactured_solutions_get_field_and_field_sym(run_info, variable_name) + manufactured_solutions_get_field_and_field_sym(run_info, variable_name; nvperp=nvperp) error = field .- field_sym nt = run_info.nt @@ -5808,7 +5815,7 @@ function _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label slices = (k=>v for (k, v) ∈ all_plot_slices if k != Symbol(:i, dim)) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name; slices...) + run_info, variable_name; nvperp=run_info.vperp.n, slices...) error = f .- f_sym fig, ax, legend_place = get_1d_ax(2; yscale=yscale, get_legend_place=:below) @@ -5833,7 +5840,7 @@ function _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label if k ∉ (Symbol(:i, dim1), Symbol(:i, dim2))) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name; slices...) + run_info, variable_name; nvperp=run_info.vperp.n, slices...) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -5857,7 +5864,7 @@ function _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label slices = (k=>v for (k, v) ∈ all_animate_slices if k != Symbol(:i, dim)) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name; slices...) + run_info, variable_name; nvperp=run_info.vperp.n, slices...) error = f .- f_sym fig, ax, legend_place = get_1d_ax(2; yscale=yscale, get_legend_place=:below) @@ -5885,7 +5892,7 @@ function _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label if k ∉ (Symbol(:i, dim1), Symbol(:i, dim2))) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name; slices...) + run_info, variable_name; nvperp=run_info.vperp.n, slices...) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -5992,7 +5999,8 @@ function compare_charged_pdf_symbolic_test(run_info, plot_prefix; io=nothing, for r_chunk ∈ r_chunks, z_chunk ∈ z_chunks f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=it, ir=r_chunk, iz=z_chunk) + run_info, variable_name; nvperp=run_info.vperp.n, it=it, + ir=r_chunk, iz=z_chunk) dummy += sum(@. (f - f_sym)^2) #dummy_N += sum(f_sym.^2) end @@ -6012,8 +6020,8 @@ function compare_charged_pdf_symbolic_test(run_info, plot_prefix; io=nothing, for (iz, z_label) ∈ ((1, "wall-"), (z.n, "wall+")) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, ir=input.ir0, iz=iz, - ivperp=input.ivperp0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, + ir=input.ir0, iz=iz, ivperp=input.ivperp0) error = f .- f_sym fig, ax, legend_place = get_1d_ax(2; get_legend_place=:below) @@ -6032,7 +6040,8 @@ function compare_charged_pdf_symbolic_test(run_info, plot_prefix; io=nothing, if !is_1D f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, iz=iz, ivperp=input.ivperp0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, iz=iz, + ivperp=input.ivperp0) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -6053,7 +6062,8 @@ function compare_charged_pdf_symbolic_test(run_info, plot_prefix; io=nothing, if !is_1V f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, iz=iz, ir=input.ir0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, iz=iz, + ir=input.ir0) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -6174,7 +6184,8 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, for r_chunk ∈ r_chunks, z_chunk ∈ z_chunks f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=it, ir=r_chunk, iz=z_chunk) + run_info, variable_name; nvperp=run_info.vperp.n, it=it, + ir=r_chunk, iz=z_chunk) dummy += sum(@. (f - f_sym)^2) #dummy_N += sum(f_sym.^2) end @@ -6194,8 +6205,8 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, for (iz, z_label) ∈ ((1, "wall-"), (z.n, "wall+")) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, ir=input.ir0, iz=iz, - ivzeta=input.ivzeta0, ivr=input.ivr0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, + ir=input.ir0, iz=iz, ivzeta=input.ivzeta0, ivr=input.ivr0) error = f .- f_sym fig, ax, legend_place = get_1d_ax(2; get_legend_place=:below) @@ -6214,8 +6225,8 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, if !is_1D f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, iz=iz, ivzeta=input.ivzeta0, - ivr=input.ivr0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, iz=iz, + ivzeta=input.ivzeta0, ivr=input.ivr0) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -6236,8 +6247,8 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, if !is_1V f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, iz=iz, ir=input.ir0, - ivzeta=input.ivzeta0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, iz=iz, + ir=input.ir0, ivzeta=input.ivzeta0) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -6259,8 +6270,8 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, iz=iz, ir=input.ir0, - ivr=input.ivr0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, iz=iz, + ir=input.ir0, ivr=input.ivr0) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -6319,7 +6330,7 @@ greater than one will result in an error. """ function manufactured_solutions_analysis end -function manufactured_solutions_analysis(run_info::Tuple; plot_prefix) +function manufactured_solutions_analysis(run_info::Tuple; plot_prefix, nvperp) if !any(ri !== nothing && ri.manufactured_solns_input.use_for_advance && ri.manufactured_solns_input.use_for_init for ri ∈ run_info) # No manufactured solutions tests @@ -6338,18 +6349,24 @@ function manufactured_solutions_analysis(run_info::Tuple; plot_prefix) return nothing end try - return manufactured_solutions_analysis(run_info[1]; plot_prefix=plot_prefix) + return manufactured_solutions_analysis(run_info[1]; plot_prefix=plot_prefix, + nvperp=nvperp) catch e println("Error in manufactured_solutions_analysis(). Error was ", e) end end -function manufactured_solutions_analysis(run_info; plot_prefix) +function manufactured_solutions_analysis(run_info; plot_prefix, nvperp) manufactured_solns_input = run_info.manufactured_solns_input if !(manufactured_solns_input.use_for_advance && manufactured_solns_input.use_for_init) return nothing end + if nvperp === nothing + error("No `nvperp` found - must have distributions function outputs to plot MMS " + * "tests") + end + input = Dict_to_NamedTuple(input_dict["manufactured_solns"]) open(run_info.run_prefix * "MMS_errors.txt", "w") do io @@ -6373,7 +6390,8 @@ function manufactured_solutions_analysis(run_info; plot_prefix) end compare_moment_symbolic_test(run_info, plot_prefix, field_label, field_sym_label, - norm_label, variable_name; io=io, input=input) + norm_label, variable_name; io=io, input=input, + nvperp=nvperp) end end From a6771e528bc9bd1914b4ae26836395961784092a Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 13:07:16 +0000 Subject: [PATCH 28/34] Convert test_velocity_integrals.jl to a CI test --- test/velocity_integral_tests.jl | 151 ++++++++++++++++++++++++++++++++ test_velocity_integrals.jl | 148 ------------------------------- 2 files changed, 151 insertions(+), 148 deletions(-) create mode 100644 test/velocity_integral_tests.jl delete mode 100644 test_velocity_integrals.jl diff --git a/test/velocity_integral_tests.jl b/test/velocity_integral_tests.jl new file mode 100644 index 000000000..bdf6765d9 --- /dev/null +++ b/test/velocity_integral_tests.jl @@ -0,0 +1,151 @@ +module VelocityIntegralTests + +include("setup.jl") + +using moment_kinetics.input_structs: grid_input, advection_input +using moment_kinetics.coordinates: define_coordinate +using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_pperp, get_pressure +using moment_kinetics.array_allocation: allocate_float + +using MPI + +function runtests() + @testset "velocity integral tests" verbose=use_verbose begin + println("velocity integral tests") + + # Tolerance for tests + atol = 1.0e-13 + + # define inputs needed for the test + ngrid = 17 #number of points per element + nelement_local = 20 # number of elements per rank + nelement_global = nelement_local # total number of elements + Lvpa = 18.0 #physical box size in reference units + Lvperp = 9.0 #physical box size in reference units + bc = "" #not required to take a particular value, not used + # fd_option and adv_input not actually used so given values unimportant + discretization = "chebyshev_pseudospectral" + fd_option = "fourth_order_centered" + adv_input = advection_input("default", 1.0, 0.0, 0.0) + nrank = 1 + irank = 0 + comm = MPI.COMM_NULL + # create the 'input' struct containing input info needed to create a + # coordinate + vr_input = grid_input("vperp", 1, 1, 1, + nrank, irank, 1.0, discretization, fd_option, bc, adv_input,comm) + vz_input = grid_input("vpa", ngrid, nelement_global, nelement_local, + nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) + vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, + nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) + vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, + nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm) + # create the coordinate struct 'x' + vpa, vpa_spectral = define_coordinate(vpa_input) + vperp, vperp_spectral = define_coordinate(vperp_input) + vz, vz_spectral = define_coordinate(vz_input) + vr, vr_spectral = define_coordinate(vr_input) + + dfn = allocate_float(vpa.n,vperp.n) + dfn1D = allocate_float(vz.n, vr.n) + + function pressure(ppar,pperp) + pres = (1.0/3.0)*(ppar + 2.0*pperp) + return pres + end + + @testset "2D isotropic Maxwellian" begin + # assign a known isotropic Maxwellian distribution in normalised units + dens = 3.0/4.0 + upar = 2.0/3.0 + ppar = 2.0/3.0 + pperp = 2.0/3.0 + pres = get_pressure(ppar,pperp) + mass = 1.0 + vth = sqrt(2.0*pres/(dens*mass)) + for ivperp in 1:vperp.n + for ivpa in 1:vpa.n + vpa_val = vpa.grid[ivpa] + vperp_val = vperp.grid[ivperp] + dfn[ivpa,ivperp] = (dens/vth^3)*exp( - ((vpa_val-upar)^2 + vperp_val^2)/vth^2 ) + end + end + + # now check that we can extract the correct moments from the distribution + + dens_test = get_density(dfn,vpa,vperp) + upar_test = get_upar(dfn,vpa,vperp,dens_test) + ppar_test = get_ppar(dfn,vpa,vperp,upar_test) + pperp_test = get_pperp(dfn,vpa,vperp) + pres_test = pressure(ppar_test,pperp_test) + # output test results + @test isapprox(dens_test, dens; atol=atol) + @test isapprox(upar_test, upar; atol=atol) + @test isapprox(pres_test, pres; atol=atol) + end + + @testset "1D Maxwellian" begin + dens = 3.0/4.0 + upar = 2.0/3.0 + ppar = 2.0/3.0 + mass = 1.0 + vth = sqrt(2.0*ppar/(dens*mass)) + for ivz in 1:vz.n + for ivr in 1:vr.n + vz_val = vz.grid[ivz] + dfn1D[ivz,ivr] = (dens/vth)*exp( - ((vz_val-upar)^2)/vth^2 ) + end + end + dens_test = get_density(dfn1D,vz,vr) + upar_test = get_upar(dfn1D,vz,vr,dens_test) + ppar_test = get_ppar(dfn1D,vz,vr,upar_test) + # output test results + println("") + println("1D Maxwellian") + @test isapprox(dens_test, dens; atol=atol) + @test isapprox(upar_test, upar; atol=atol) + @test isapprox(ppar_test, ppar; atol=atol) + println("") + end + + @testset "biMaxwellian" begin + # assign a known biMaxwellian distribution in normalised units + dens = 3.0/4.0 + upar = 2.0/3.0 + ppar = 4.0/5.0 + pperp = 1.0/4.0 + mass = 1.0 + vthpar = sqrt(2.0*ppar/(dens*mass)) + vthperp = sqrt(2.0*pperp/(dens*mass)) + for ivperp in 1:vperp.n + for ivpa in 1:vpa.n + vpa_val = vpa.grid[ivpa] + vperp_val = vperp.grid[ivperp] + dfn[ivpa,ivperp] = (dens/(vthpar*vthperp^2))*exp( - ((vpa_val-upar)^2)/vthpar^2 - (vperp_val^2)/vthperp^2 ) + end + end + + # now check that we can extract the correct moments from the distribution + + dens_test = get_density(dfn,vpa,vperp) + upar_test = get_upar(dfn,vpa,vperp,dens_test) + ppar_test = get_ppar(dfn,vpa,vperp,upar_test) + pperp_test = get_pperp(dfn,vpa,vperp) + # output test results + + println("") + println("biMaxwellian") + @test isapprox(dens_test, dens; atol=atol) + @test isapprox(upar_test, upar; atol=atol) + @test isapprox(ppar_test, ppar; atol=atol) + println("pperp_test: ", pperp_test, " pperp: ", pperp, " error: ", abs(pperp_test-pperp)) + println("") + end + end +end + +end # VelocityIntegralTests + +using .VelocityIntegralTests + +VelocityIntegralTests.runtests() diff --git a/test_velocity_integrals.jl b/test_velocity_integrals.jl deleted file mode 100644 index 25a023539..000000000 --- a/test_velocity_integrals.jl +++ /dev/null @@ -1,148 +0,0 @@ -using Printf -using Plots -using LaTeXStrings -using MPI -using Measures - -if abspath(PROGRAM_FILE) == @__FILE__ - using Pkg - Pkg.activate(".") - - import moment_kinetics - using moment_kinetics.input_structs: grid_input, advection_input - using moment_kinetics.coordinates: define_coordinate - using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_pperp, get_pressure - using moment_kinetics.array_allocation: allocate_float - - # define inputs needed for the test - ngrid = 17 #number of points per element - nelement_local = 20 # number of elements per rank - nelement_global = nelement_local # total number of elements - Lvpa = 18.0 #physical box size in reference units - Lvperp = 9.0 #physical box size in reference units - bc = "" #not required to take a particular value, not used - # fd_option and adv_input not actually used so given values unimportant - discretization = "chebyshev_pseudospectral" - fd_option = "fourth_order_centered" - adv_input = advection_input("default", 1.0, 0.0, 0.0) - nrank = 1 - irank = 0 - comm = MPI.COMM_NULL - # create the 'input' struct containing input info needed to create a - # coordinate - vr_input = grid_input("vperp", 1, 1, 1, - nrank, irank, 1.0, discretization, fd_option, bc, adv_input,comm) - vz_input = grid_input("vpa", ngrid, nelement_global, nelement_local, - nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) - vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, - nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) - vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, - nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm) - # create the coordinate struct 'x' - println("made inputs") - println("vpa: ngrid: ",ngrid," nelement: ",nelement_local, " Lvpa: ",Lvpa) - println("vperp: ngrid: ",ngrid," nelement: ",nelement_local, " Lvperp: ",Lvperp) - vpa, vpa_spectral = define_coordinate(vpa_input) - vperp, vperp_spectral = define_coordinate(vperp_input) - vz, vz_spectral = define_coordinate(vz_input) - vr, vr_spectral = define_coordinate(vr_input) - - dfn = allocate_float(vpa.n,vperp.n) - dfn1D = allocate_float(vz.n, vr.n) - - function pressure(ppar,pperp) - pres = (1.0/3.0)*(ppar + 2.0*pperp) - return pres - end - # 2D isotropic Maxwellian test - # assign a known isotropic Maxwellian distribution in normalised units - dens = 3.0/4.0 - upar = 2.0/3.0 - ppar = 2.0/3.0 - pperp = 2.0/3.0 - pres = get_pressure(ppar,pperp) - mass = 1.0 - vth = sqrt(2.0*pres/(dens*mass)) - for ivperp in 1:vperp.n - for ivpa in 1:vpa.n - vpa_val = vpa.grid[ivpa] - vperp_val = vperp.grid[ivperp] - dfn[ivpa,ivperp] = (dens/vth^3)*exp( - ((vpa_val-upar)^2 + vperp_val^2)/vth^2 ) - end - end - - # now check that we can extract the correct moments from the distribution - - dens_test = get_density(dfn,vpa,vperp) - upar_test = get_upar(dfn,vpa,vperp,dens_test) - ppar_test = get_ppar(dfn,vpa,vperp,upar_test) - pperp_test = get_pperp(dfn,vpa,vperp) - pres_test = pressure(ppar_test,pperp_test) - # output test results - println("") - println("Isotropic 2D Maxwellian") - println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) - println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) - println("pres_test: ", pres_test, " pres: ", pres, " error: ", abs(pres_test-pres)) - println("") - - ################### - # 1D Maxwellian test - - dens = 3.0/4.0 - upar = 2.0/3.0 - ppar = 2.0/3.0 - mass = 1.0 - vth = sqrt(2.0*ppar/(dens*mass)) - for ivz in 1:vz.n - for ivr in 1:vr.n - vz_val = vz.grid[ivz] - dfn1D[ivz,ivr] = (dens/vth)*exp( - ((vz_val-upar)^2)/vth^2 ) - end - end - dens_test = get_density(dfn1D,vz,vr) - upar_test = get_upar(dfn1D,vz,vr,dens_test) - ppar_test = get_ppar(dfn1D,vz,vr,upar_test) - # output test results - println("") - println("1D Maxwellian") - println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) - println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) - println("ppar_test: ", ppar_test, " ppar: ", ppar, " error: ", abs(ppar_test-ppar)) - println("") - - ################### - # biMaxwellian test - - # assign a known biMaxwellian distribution in normalised units - dens = 3.0/4.0 - upar = 2.0/3.0 - ppar = 4.0/5.0 - pperp = 1.0/4.0 - mass = 1.0 - vthpar = sqrt(2.0*ppar/(dens*mass)) - vthperp = sqrt(2.0*pperp/(dens*mass)) - for ivperp in 1:vperp.n - for ivpa in 1:vpa.n - vpa_val = vpa.grid[ivpa] - vperp_val = vperp.grid[ivperp] - dfn[ivpa,ivperp] = (dens/(vthpar*vthperp^2))*exp( - ((vpa_val-upar)^2)/vthpar^2 - (vperp_val^2)/vthperp^2 ) - end - end - - # now check that we can extract the correct moments from the distribution - - dens_test = get_density(dfn,vpa,vperp) - upar_test = get_upar(dfn,vpa,vperp,dens_test) - ppar_test = get_ppar(dfn,vpa,vperp,upar_test) - pperp_test = get_pperp(dfn,vpa,vperp) - # output test results - - println("") - println("biMaxwellian") - println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) - println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) - println("ppar_test: ", ppar_test, " ppar: ", ppar, " error: ", abs(ppar_test-ppar)) - println("pperp_test: ", pperp_test, " pperp: ", pperp, " error: ", abs(pperp_test-pperp)) - println("") -end From 6884d86ebca9f58479b4f3648c43c509f255b54e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 14:42:58 +0000 Subject: [PATCH 29/34] Deactivate Krook collisions by default --- src/moment_kinetics_input.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 4e75fc7b7..1132b961f 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -168,14 +168,18 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) collisions.charge_exchange = get(scan_input, "charge_exchange_frequency", 2.0*sqrt(species.charged[1].initial_temperature)) collisions.ionization = get(scan_input, "ionization_frequency", collisions.charge_exchange) collisions.constant_ionization_rate = get(scan_input, "constant_ionization_rate", false) - collisions.krook_collisions_option = get(scan_input, "krook_collisions_option", "reference_parameters") + collisions.krook_collisions_option = get(scan_input, "krook_collisions_option", "none") nuii_krook_default = setup_krook_collisions(reference_params) if collisions.krook_collisions_option == "reference_parameters" collisions.krook_collision_frequency_prefactor = nuii_krook_default elseif collisions.krook_collisions_option == "manual" # get the frequency from the input file collisions.krook_collision_frequency_prefactor = get(scan_input, "nuii_krook", nuii_krook_default) - else + elseif collisions.krook_collisions_option == "none" + # By default, no krook collisions included collisions.krook_collision_frequency_prefactor = -1.0 + else + error("Invalid option " + * "krook_collisions_option=$(collisions.krook_collisions_option) passed") end # parameters related to the time stepping @@ -957,7 +961,7 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) constant_ionization_rate = false krook_collision_frequency_prefactor = -1.0 collisions = collisions_input(charge_exchange, ionization, constant_ionization_rate, - krook_collision_frequency_prefactor,"default") + krook_collision_frequency_prefactor,"none") Bzed = 1.0 # magnetic field component along z Bmag = 1.0 # magnetic field strength From 92a4e29ad9fa130101738c3033eb5625958e0d26 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 14:44:07 +0000 Subject: [PATCH 30/34] Pass 'element_spacing_option' to coords in velocity_integral_tests.jl This is needed since 'sqrt spacing' option was added. Also include velocity_integral_tests.jl in runtests.jl so that it is included in the CI tests. Also remove some println() statements that had been left in by mistake. --- test/runtests.jl | 1 + test/velocity_integral_tests.jl | 26 +++++++++++--------------- 2 files changed, 12 insertions(+), 15 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 17ba5030d..4bab7c641 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,7 @@ function runtests() include(joinpath(@__DIR__, "calculus_tests.jl")) include(joinpath(@__DIR__, "interpolation_tests.jl")) include(joinpath(@__DIR__, "loop_setup_tests.jl")) + include(joinpath(@__DIR__, "velocity_integral_tests.jl")) include(joinpath(@__DIR__, "sound_wave_tests.jl")) include(joinpath(@__DIR__, "nonlinear_sound_wave_tests.jl")) include(joinpath(@__DIR__, "Krook_collisions_tests.jl")) diff --git a/test/velocity_integral_tests.jl b/test/velocity_integral_tests.jl index bdf6765d9..dc84800f0 100644 --- a/test/velocity_integral_tests.jl +++ b/test/velocity_integral_tests.jl @@ -32,14 +32,17 @@ function runtests() comm = MPI.COMM_NULL # create the 'input' struct containing input info needed to create a # coordinate - vr_input = grid_input("vperp", 1, 1, 1, - nrank, irank, 1.0, discretization, fd_option, bc, adv_input,comm) - vz_input = grid_input("vpa", ngrid, nelement_global, nelement_local, - nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) - vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, - nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) - vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, - nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm) + vr_input = grid_input("vperp", 1, 1, 1, nrank, irank, 1.0, discretization, + fd_option, bc, adv_input, comm, "uniform") + vz_input = grid_input("vpa", ngrid, nelement_global, nelement_local, nrank, irank, + Lvpa, discretization, fd_option, bc, adv_input, comm, + "uniform") + vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, nrank, + irank, Lvpa, discretization, fd_option, bc, adv_input, + comm, "uniform") + vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, nrank, + irank, Lvperp, discretization, fd_option, bc, adv_input, + comm, "uniform") # create the coordinate struct 'x' vpa, vpa_spectral = define_coordinate(vpa_input) vperp, vperp_spectral = define_coordinate(vperp_input) @@ -100,12 +103,9 @@ function runtests() upar_test = get_upar(dfn1D,vz,vr,dens_test) ppar_test = get_ppar(dfn1D,vz,vr,upar_test) # output test results - println("") - println("1D Maxwellian") @test isapprox(dens_test, dens; atol=atol) @test isapprox(upar_test, upar; atol=atol) @test isapprox(ppar_test, ppar; atol=atol) - println("") end @testset "biMaxwellian" begin @@ -133,13 +133,9 @@ function runtests() pperp_test = get_pperp(dfn,vpa,vperp) # output test results - println("") - println("biMaxwellian") @test isapprox(dens_test, dens; atol=atol) @test isapprox(upar_test, upar; atol=atol) @test isapprox(ppar_test, ppar; atol=atol) - println("pperp_test: ", pperp_test, " pperp: ", pperp, " error: ", abs(pperp_test-pperp)) - println("") end end end From afe37095ef4e1550b2d578ffb6e16d1b6c53441e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 29 Sep 2023 10:16:32 +0100 Subject: [PATCH 31/34] Better values for MPI variables when ignore_MPI=true in mk_input() Set irank/nrank as if running in serial, to avoid errors from set_element_scale_and_shift()). --- src/moment_kinetics_input.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 1132b961f..72e23d6f5 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -372,7 +372,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # need grid and MPI information to determine these values # MRH just put dummy values now if ignore_MPI - irank_z = nrank_z = irank_r = nrank_r = -1 + irank_z = irank_r = 0 + nrank_z = nrank_r = 1 comm_sub_z = comm_sub_r = MPI.COMM_NULL else irank_z, nrank_z, comm_sub_z, irank_r, nrank_r, comm_sub_r = setup_distributed_memory_MPI(z.nelement_global,z.nelement_local,r.nelement_global,r.nelement_local) From 5d5ddafaa43c4bcf6819cfca7a3ac1e4b7368c95 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 19:17:46 +0000 Subject: [PATCH 32/34] Tidy up comments --- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml | 2 -- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml | 2 -- src/input_structs.jl | 2 +- src/krook_collisions.jl | 2 +- test/velocity_integral_tests.jl | 3 --- 5 files changed, 2 insertions(+), 9 deletions(-) diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml index 5a329be7a..7265dd542 100644 --- a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml @@ -31,7 +31,6 @@ z_IC_temperature_amplitude2 = 0.0 z_IC_temperature_phase2 = 0.0 charge_exchange_frequency = 0.0 ionization_frequency = 0.0 -#krook_collisions = false krook_collisions_option = "manual" nuii_krook = 1.0 nstep = 2000 @@ -60,7 +59,6 @@ vperp_ngrid = 1 vperp_nelement = 1 vperp_L = 6.0 vperp_bc = "periodic" -#vperp_discretization = "finite_difference" vperp_discretization = "chebyshev_pseudospectral" vz_ngrid = 17 diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml index 5173fa230..4d35a8a77 100644 --- a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml @@ -31,7 +31,6 @@ z_IC_temperature_amplitude2 = 0.0 z_IC_temperature_phase2 = 0.0 charge_exchange_frequency = 0.0 ionization_frequency = 0.0 -#krook_collisions = false krook_collisions_option = "manual" nuii_krook = 1.0 nstep = 2000 @@ -60,7 +59,6 @@ vperp_ngrid = 17 vperp_nelement = 8 vperp_L = 6.0 vperp_bc = "periodic" -#vperp_discretization = "finite_difference" vperp_discretization = "chebyshev_pseudospectral" vz_ngrid = 17 diff --git a/src/input_structs.jl b/src/input_structs.jl index f7734a982..be8a19969 100644 --- a/src/input_structs.jl +++ b/src/input_structs.jl @@ -306,7 +306,7 @@ mutable struct collisions_input constant_ionization_rate::Bool # Coulomb collision rate at the reference density and temperature krook_collision_frequency_prefactor::mk_float - # Coulomb collision rate at the reference density and temperature + # Setting to switch between different options for Krook collision operator krook_collisions_option::String end diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index b77db5bf1..f6ceec9c4 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -50,7 +50,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v # normalized by sqrt(pi) (see velocity_moments.integrate_over_vspace). if collisions.krook_collisions_option == "manual" cfac = 0.0 - else # default option, collisions.krook_collisions_option == "default" + else # default option, collisions.krook_collisions_option == "reference_parameters" cfac = 1.0 end if moments.evolve_ppar && moments.evolve_upar diff --git a/test/velocity_integral_tests.jl b/test/velocity_integral_tests.jl index dc84800f0..07d0c5027 100644 --- a/test/velocity_integral_tests.jl +++ b/test/velocity_integral_tests.jl @@ -81,7 +81,6 @@ function runtests() ppar_test = get_ppar(dfn,vpa,vperp,upar_test) pperp_test = get_pperp(dfn,vpa,vperp) pres_test = pressure(ppar_test,pperp_test) - # output test results @test isapprox(dens_test, dens; atol=atol) @test isapprox(upar_test, upar; atol=atol) @test isapprox(pres_test, pres; atol=atol) @@ -102,7 +101,6 @@ function runtests() dens_test = get_density(dfn1D,vz,vr) upar_test = get_upar(dfn1D,vz,vr,dens_test) ppar_test = get_ppar(dfn1D,vz,vr,upar_test) - # output test results @test isapprox(dens_test, dens; atol=atol) @test isapprox(upar_test, upar; atol=atol) @test isapprox(ppar_test, ppar; atol=atol) @@ -131,7 +129,6 @@ function runtests() upar_test = get_upar(dfn,vpa,vperp,dens_test) ppar_test = get_ppar(dfn,vpa,vperp,upar_test) pperp_test = get_pperp(dfn,vpa,vperp) - # output test results @test isapprox(dens_test, dens; atol=atol) @test isapprox(upar_test, upar; atol=atol) From 8a5e2986ab21e5159ebf25fb0590cf2b5b422d61 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 19:58:01 +0000 Subject: [PATCH 33/34] Make clearer that Ti_over_Tref is not normalised temperature --- src/manufactured_solns.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 2c2f7d05a..df267400e 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -524,11 +524,11 @@ using IfElse end nu_krook = collisions.krook_collision_frequency_prefactor if nu_krook > 0.0 - tempi = vthi^2 + Ti_over_Tref = vthi^2 if collisions.krook_collisions_option == "manual" nuii_krook = nu_krook else # default option - nuii_krook = nu_krook * densi * tempi^(-1.5) + nuii_krook = nu_krook * densi * Ti_over_Tref^(-1.5) end if vperp_coord.n > 1 pvth = 3 From 4c9e8109841f8f6c0dd27190f38352443b73c07c Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 19:58:58 +0000 Subject: [PATCH 34/34] Krook MMS test inputs use spatially-varying collision frequency This option now passes the tests, so might as well make it the suggested option. --- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml | 2 +- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml index 7265dd542..e9d44cc5b 100644 --- a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml @@ -31,7 +31,7 @@ z_IC_temperature_amplitude2 = 0.0 z_IC_temperature_phase2 = 0.0 charge_exchange_frequency = 0.0 ionization_frequency = 0.0 -krook_collisions_option = "manual" +krook_collisions_option = "reference_parameters" nuii_krook = 1.0 nstep = 2000 dt = 0.0005 diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml index 4d35a8a77..ae02bfd9c 100644 --- a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml @@ -31,7 +31,7 @@ z_IC_temperature_amplitude2 = 0.0 z_IC_temperature_phase2 = 0.0 charge_exchange_frequency = 0.0 ionization_frequency = 0.0 -krook_collisions_option = "manual" +krook_collisions_option = "reference_parameters" nuii_krook = 1.0 nstep = 2000 dt = 0.0005