From 226c74d7b73e126dbe32cfc2f91fa39bf37ef50c Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Thu, 8 Feb 2024 15:19:15 +0000 Subject: [PATCH 01/14] Attempt to use phi to modify slightly where the wall BC is imposed - likely to fail due to sqrt(-1) if potential ever reverses gradient in a fluctuation. --- moment_kinetics/src/initial_conditions.jl | 23 +++++++++++++---------- moment_kinetics/src/time_advance.jl | 4 ++-- 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index 0d8baee8f..9d3b00af3 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -23,7 +23,7 @@ using ..coordinates: coordinate using ..external_sources using ..interpolation: interpolate_to_grid_1d! using ..looping -using ..moment_kinetics_structs: scratch_pdf +using ..moment_kinetics_structs: scratch_pdf, em_fields_struct using ..velocity_moments: integrate_over_vspace, integrate_over_neutral_vspace using ..velocity_moments: integrate_over_positive_vpa, integrate_over_negative_vpa using ..velocity_moments: integrate_over_positive_vz, integrate_over_negative_vz @@ -1133,7 +1133,7 @@ end enforce boundary conditions in vpa and z on the evolved pdf; also enforce boundary conditions in z on all separately evolved velocity space moments of the pdf """ -function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, moments, vpa_bc, +function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, phi, moments, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, r_adv, composition, scratch_dummy, r_diffusion, vpa_diffusion, vperp_diffusion) if vpa.n > 1 @@ -1155,7 +1155,7 @@ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, moments, v begin_s_r_vperp_vpa_region() # enforce the z BC on the evolved velocity space moments of the pdf @views enforce_z_boundary_condition_moments!(density, moments, z_bc) - @views enforce_z_boundary_condition!(f, density, upar, ppar, moments, z_bc, z_adv, z, + @views enforce_z_boundary_condition!(f, density, upar, ppar, phi, moments, z_bc, z_adv, z, vperp, vpa, composition, scratch_dummy.buffer_vpavperprs_1, scratch_dummy.buffer_vpavperprs_2, scratch_dummy.buffer_vpavperprs_3, scratch_dummy.buffer_vpavperprs_4) @@ -1169,11 +1169,11 @@ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, moments, v r_diffusion) end end -function enforce_boundary_conditions!(fvec_out::scratch_pdf, moments, f_r_bc, vpa_bc, +function enforce_boundary_conditions!(fvec_out::scratch_pdf, moments, fields::em_fields_struct, f_r_bc, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, r_adv, composition, scratch_dummy, r_diffusion, vpa_diffusion, vperp_diffusion) enforce_boundary_conditions!(fvec_out.pdf, f_r_bc, fvec_out.density, fvec_out.upar, - fvec_out.ppar, moments, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, + fvec_out.ppar, fields.phi, moments, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, r_adv, composition, scratch_dummy, r_diffusion, vpa_diffusion, vperp_diffusion) end @@ -1231,7 +1231,7 @@ end """ enforce boundary conditions on charged particle f in z """ -function enforce_z_boundary_condition!(pdf, density, upar, ppar, moments, bc::String, adv, +function enforce_z_boundary_condition!(pdf, density, upar, ppar, phi, moments, bc::String, adv, z, vperp, vpa, composition, end1::AbstractArray{mk_float,4}, end2::AbstractArray{mk_float,4}, buffer1::AbstractArray{mk_float,4}, buffer2::AbstractArray{mk_float,4}) @@ -1293,7 +1293,7 @@ function enforce_z_boundary_condition!(pdf, density, upar, ppar, moments, bc::St else @loop_r ir begin @views enforce_zero_incoming_bc!(pdf[:,:,:,ir,is], - adv[is].speed[:,:,:,ir], z, zero) + adv[is].speed[:,:,:,ir], z, zero, phi[:,ir]) end end end @@ -1504,26 +1504,29 @@ end """ enforce a zero incoming BC in z for given species pdf at each radial location """ -function enforce_zero_incoming_bc!(pdf, speed, z, zero) +function enforce_zero_incoming_bc!(pdf, speed, z, zero, phi) nvpa = size(pdf,1) + nz = z.n # no parallel BC should be enforced for dz/dt = 0 # note that the parallel velocity coordinate vpa may be dz/dt or # some version of the peculiar velocity (dz/dt - upar), # so use advection speed below instead of vpa if z.irank == 0 + vcut = sqrt(phi[2]-phi[1]) # sqrt(-1) an option! @loop_vperp_vpa ivperp ivpa begin # for left boundary in zed (z = -Lz/2), want # f(z=-Lz/2, v_parallel > 0) = 0 - if speed[1,ivpa,ivperp] > zero + if speed[1,ivpa,ivperp] > zero - vcut pdf[ivpa,ivperp,1] = 0.0 end end end if z.irank == z.nrank - 1 + vcut = sqrt(phi[nz-1]-phi[nz]) # sqrt(-1) an option! @loop_vperp_vpa ivperp ivpa begin # for right boundary in zed (z = Lz/2), want # f(z=Lz/2, v_parallel < 0) = 0 - if speed[end,ivpa,ivperp] < -zero + if speed[end,ivpa,ivperp] < -zero + vcut pdf[ivpa,ivperp,end] = 0.0 end end diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index f3d08b899..6bb5d7f26 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -371,7 +371,7 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral, # condition enforce_boundary_conditions!( pdf.charged.norm, boundary_distributions.pdf_rboundary_charged, - moments.charged.dens, moments.charged.upar, moments.charged.ppar, moments, + moments.charged.dens, moments.charged.upar, moments.charged.ppar, fields.phi, moments, vpa.bc, z.bc, r.bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_advect, vperp_advect, z_advect, r_advect, composition, scratch_dummy, advance.r_diffusion, @@ -1268,7 +1268,7 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v # set to zero at the sheath boundary according to the final upar has a non-zero # contribution from one or more of the terms. # NB: probably need to do the same for the evolved moments - enforce_boundary_conditions!(new_scratch, moments, + enforce_boundary_conditions!(new_scratch, moments, fields, boundary_distributions.pdf_rboundary_charged, vpa.bc, z.bc, r.bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_advect, vperp_advect, z_advect, r_advect, composition, scratch_dummy, From 9027c7daf88ba78081e788cf1bf1c2db5c919fc1 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 12 Feb 2024 09:25:58 +0000 Subject: [PATCH 02/14] Use ternary operator to test the sign of delta phi and use the structure of ion phase space to justify using the vcut feature only if deltaphi has the correct sign. epsz is a parameter that would vary from 0 (vcut=0) to 1 (vcut corresponds to particles that reverse at the first grid point from the wall). Remarkably, this feature shows that the pdf at the wall is barely or not at all resolved for slow particles, because with vcut = 1, a significant part of the distribution function is set to zero. --- moment_kinetics/src/initial_conditions.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index 9d3b00af3..960b3b124 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -1511,8 +1511,13 @@ function enforce_zero_incoming_bc!(pdf, speed, z, zero, phi) # note that the parallel velocity coordinate vpa may be dz/dt or # some version of the peculiar velocity (dz/dt - upar), # so use advection speed below instead of vpa + epsz = 1.0 # the ratio |z - z_wall|/|delta z|, with delta z the grid spacing at the wall + # for epsz < 1, the cut off below would be imposed for particles travelling + # out to a distance z = epsz * delta z from the wall before returning + # epsz should really be an input parameter if z.irank == 0 - vcut = sqrt(phi[2]-phi[1]) # sqrt(-1) an option! + deltaphi = phi[2] - phi[1] + vcut = deltaphi > 0 ? sqrt(deltaphi)*(epsz^0.25) : 0.0 @loop_vperp_vpa ivperp ivpa begin # for left boundary in zed (z = -Lz/2), want # f(z=-Lz/2, v_parallel > 0) = 0 @@ -1522,7 +1527,8 @@ function enforce_zero_incoming_bc!(pdf, speed, z, zero, phi) end end if z.irank == z.nrank - 1 - vcut = sqrt(phi[nz-1]-phi[nz]) # sqrt(-1) an option! + deltaphi = phi[nz-1] - phi[nz] + vcut = deltaphi > 0 ? sqrt(deltaphi)*(epsz^0.25) : 0.0 @loop_vperp_vpa ivperp ivpa begin # for right boundary in zed (z = Lz/2), want # f(z=Lz/2, v_parallel < 0) = 0 From 92686ce25a889052cec4104642b3c56b3dffc52a Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 12 Feb 2024 09:26:28 +0000 Subject: [PATCH 03/14] A plot of f/vpa^2 to diagnose the vcut feature. --- .../src/plots_post_processing.jl | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl index f3e4aceb1..17b62354c 100644 --- a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl +++ b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl @@ -1029,7 +1029,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) composition, species, collisions, geometry, drive_input, external_source_settings, num_diss_params, manufactured_solns_input = input - if !is_1D1V + if !is_1D1V || true # make plots and animations of the phi, Ez and Er plot_charged_moments_2D(density, parallel_flow, parallel_pressure, perpendicular_pressure, thermal_speed, entropy_production, @@ -3435,6 +3435,20 @@ function plot_charged_pdf_2D_at_wall(run_name, run_name_label, r_global, z_globa @views plot(vpa.grid, pdf[:,ivperp0,ir0,is,itime0], xlabel=L"v_{\|\|}/L_{v_{\|\|}}", ylabel=L"f_i") outfile = string(run_name_label, "_pdf(vpa,vperp0,iz_"*zlabel*",ir0)"*description*"vs_vpa.pdf") trysavefig(outfile) + + # plot f(vpa,ivperp0,iz_wall,ir0,is,itime)/vpa^2 at the wall + deltavpa = minimum(vpa.grid[2:vpa.n].-vpa.grid[1:vpa.n-1]) + vpafunc = copy(vpa.grid) + for ivpa in 1:vpa.n + if abs(vpa.grid[ivpa]) > 0.5*deltavpa + vpafunc[ivpa] = 1.0/(vpa.grid[ivpa]^2) + else + vpafunc[ivpa] = 0.0 + end + end + @views plot(vpa.grid, vpafunc.*pdf[:,ivperp0,ir0,is,itime0], xlabel=L"v_{\|\|}/L_{v_{\|\|}}", ylabel=L"f_i/v_{\|\|}^2") + outfile = string(run_name_label, "_pdf(vpa,vperp0,iz_"*zlabel*",ir0)_over_vpa2"*description*"vs_vpa.pdf") + trysavefig(outfile) # plot f(vpa,vperp,iz_wall,ir0,is,itime) at the wall @views heatmap(vperp.grid, vpa.grid, pdf[:,:,ir0,is,itime0], xlabel=L"v_{\perp}", ylabel=L"v_{||}", c = :deep, interpolation = :cubic, From c84305b0156f3ef3464050018a3398c50ccf8b8c Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Tue, 13 Feb 2024 18:21:59 +0000 Subject: [PATCH 04/14] Changes to make epsz an input parameter in the newx [boundary_condition_parameters] toml list. --- moment_kinetics/src/initial_conditions.jl | 43 ++++++++++++++------ moment_kinetics/src/moment_kinetics.jl | 4 +- moment_kinetics/src/moment_kinetics_input.jl | 6 ++- moment_kinetics/src/time_advance.jl | 4 +- 4 files changed, 38 insertions(+), 19 deletions(-) diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index 960b3b124..8f890dd8e 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -6,6 +6,7 @@ export allocate_pdf_and_moments export init_pdf_and_moments! export enforce_boundary_conditions! export enforce_neutral_boundary_conditions! +export setup_boundary_parameters # functional testing export create_boundary_distributions @@ -37,6 +38,20 @@ using ..manufactured_solns: manufactured_solutions using MPI + +Base.@kwdef struct boundary_condition_parameters + # parameter controlling the cutoff of the ion distribution function + # in the vpa domain at the wall in z + epsz::mk_float = 0.0 +end + +function setup_boundary_parameters(input_section::Dict) + + input = Dict(Symbol(k)=>v for (k,v) in input_section) + + return boundary_condition_parameters(; input...) +end + """ """ struct pdf_substruct{n_distribution} @@ -76,6 +91,7 @@ struct boundary_distributions_struct pdf_rboundary_charged::MPISharedArray{mk_float,5} # neutral particle r boundary values (vz,vr,vzeta,z,r,s) pdf_rboundary_neutral::MPISharedArray{mk_float,6} + bc_parameters::boundary_condition_parameters end """ @@ -83,7 +99,7 @@ Creates the structs for the pdf and the velocity-space moments """ function allocate_pdf_and_moments(composition, r, z, vperp, vpa, vzeta, vr, vz, evolve_moments, collisions, external_source_settings, - numerical_dissipation) + numerical_dissipation, boundary_parameters::boundary_condition_parameters) pdf = create_pdf(composition, r, z, vperp, vpa, vzeta, vr, vz) # create the 'moments' struct that contains various v-space moments and other @@ -116,7 +132,7 @@ function allocate_pdf_and_moments(composition, r, z, vperp, vpa, vzeta, vr, vz, evolve_moments.parallel_pressure) boundary_distributions = create_boundary_distributions(vz, vr, vzeta, vpa, vperp, z, - composition) + composition, boundary_parameters) return pdf, moments, boundary_distributions end @@ -1105,7 +1121,7 @@ Allocate arrays for distributions to be applied as boundary conditions to the pd various boundaries. Also initialise the Knudsen cosine distribution here so it can be used when initialising the neutral pdf. """ -function create_boundary_distributions(vz, vr, vzeta, vpa, vperp, z, composition) +function create_boundary_distributions(vz, vr, vzeta, vpa, vperp, z, composition, boundary_parameters::boundary_condition_parameters) zero = 1.0e-14 #initialise knudsen distribution for neutral wall bc @@ -1119,7 +1135,7 @@ function create_boundary_distributions(vz, vr, vzeta, vpa, vperp, z, composition pdf_rboundary_neutral = allocate_shared_float(vz.n, vr.n, vzeta.n, z.n, 2, composition.n_neutral_species) - return boundary_distributions_struct(knudsen_cosine, pdf_rboundary_charged, pdf_rboundary_neutral) + return boundary_distributions_struct(knudsen_cosine, pdf_rboundary_charged, pdf_rboundary_neutral, boundary_parameters) end function init_boundary_distributions!(boundary_distributions, pdf, vz, vr, vzeta, vpa, vperp, z, r, composition) @@ -1135,7 +1151,7 @@ also enforce boundary conditions in z on all separately evolved velocity space m """ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, phi, moments, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, r_adv, composition, scratch_dummy, - r_diffusion, vpa_diffusion, vperp_diffusion) + r_diffusion, vpa_diffusion, vperp_diffusion, bc_parameters) if vpa.n > 1 begin_s_r_z_vperp_region() @loop_s_r_z_vperp is ir iz ivperp begin @@ -1156,7 +1172,7 @@ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, phi, momen # enforce the z BC on the evolved velocity space moments of the pdf @views enforce_z_boundary_condition_moments!(density, moments, z_bc) @views enforce_z_boundary_condition!(f, density, upar, ppar, phi, moments, z_bc, z_adv, z, - vperp, vpa, composition, + vperp, vpa, composition, bc_parameters, scratch_dummy.buffer_vpavperprs_1, scratch_dummy.buffer_vpavperprs_2, scratch_dummy.buffer_vpavperprs_3, scratch_dummy.buffer_vpavperprs_4) @@ -1171,11 +1187,11 @@ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, phi, momen end function enforce_boundary_conditions!(fvec_out::scratch_pdf, moments, fields::em_fields_struct, f_r_bc, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, r_adv, composition, scratch_dummy, - r_diffusion, vpa_diffusion, vperp_diffusion) + r_diffusion, vpa_diffusion, vperp_diffusion, bc_parameters) enforce_boundary_conditions!(fvec_out.pdf, f_r_bc, fvec_out.density, fvec_out.upar, fvec_out.ppar, fields.phi, moments, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, - r_adv, composition, scratch_dummy, r_diffusion, vpa_diffusion, vperp_diffusion) + r_adv, composition, scratch_dummy, r_diffusion, vpa_diffusion, vperp_diffusion, bc_parameters) end """ @@ -1232,7 +1248,7 @@ end enforce boundary conditions on charged particle f in z """ function enforce_z_boundary_condition!(pdf, density, upar, ppar, phi, moments, bc::String, adv, - z, vperp, vpa, composition, end1::AbstractArray{mk_float,4}, + z, vperp, vpa, composition, bc_parameters, end1::AbstractArray{mk_float,4}, end2::AbstractArray{mk_float,4}, buffer1::AbstractArray{mk_float,4}, buffer2::AbstractArray{mk_float,4}) # this block ensures periodic BC can be supported with distributed memory MPI @@ -1293,7 +1309,8 @@ function enforce_z_boundary_condition!(pdf, density, upar, ppar, phi, moments, b else @loop_r ir begin @views enforce_zero_incoming_bc!(pdf[:,:,:,ir,is], - adv[is].speed[:,:,:,ir], z, zero, phi[:,ir]) + adv[is].speed[:,:,:,ir], z, zero, phi[:,ir], + bc_parameters.epsz) end end end @@ -1504,17 +1521,17 @@ end """ enforce a zero incoming BC in z for given species pdf at each radial location """ -function enforce_zero_incoming_bc!(pdf, speed, z, zero, phi) +function enforce_zero_incoming_bc!(pdf, speed, z, zero, phi, epsz) nvpa = size(pdf,1) nz = z.n # no parallel BC should be enforced for dz/dt = 0 # note that the parallel velocity coordinate vpa may be dz/dt or # some version of the peculiar velocity (dz/dt - upar), # so use advection speed below instead of vpa - epsz = 1.0 # the ratio |z - z_wall|/|delta z|, with delta z the grid spacing at the wall + + # epsz is the ratio |z - z_wall|/|delta z|, with delta z the grid spacing at the wall # for epsz < 1, the cut off below would be imposed for particles travelling # out to a distance z = epsz * delta z from the wall before returning - # epsz should really be an input parameter if z.irank == 0 deltaphi = phi[2] - phi[1] vcut = deltaphi > 0 ? sqrt(deltaphi)*(epsz^0.25) : 0.0 diff --git a/moment_kinetics/src/moment_kinetics.jl b/moment_kinetics/src/moment_kinetics.jl index 970f01fae..ca193c32d 100644 --- a/moment_kinetics/src/moment_kinetics.jl +++ b/moment_kinetics/src/moment_kinetics.jl @@ -290,7 +290,7 @@ function setup_moment_kinetics(input_dict::AbstractDict; vperp, vperp_spectral, gyrophase, gyrophase_spectral, vz, vz_spectral, vr, vr_spectral, vzeta, vzeta_spectral, composition, species, collisions, geometry, drive_input, external_source_settings, num_diss_params, - manufactured_solns_input = input + manufactured_solns_input, bc_params = input # Create loop range variables for shared-memory-parallel loops if debug_loop_type === nothing @@ -316,7 +316,7 @@ function setup_moment_kinetics(input_dict::AbstractDict; pdf, moments, boundary_distributions = allocate_pdf_and_moments(composition, r, z, vperp, vpa, vzeta, vr, vz, evolve_moments, collisions, external_source_settings, - num_diss_params) + num_diss_params, bc_params) if restart === false restarting = false diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 3d387e241..0a78fe097 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -18,6 +18,7 @@ using ..krook_collisions: setup_krook_collisions using ..finite_differences: fd_check_option using ..input_structs using ..numerical_dissipation: setup_numerical_dissipation +using ..initial_conditions: setup_boundary_parameters using ..reference_parameters using ..geo: init_magnetic_geometry @@ -378,7 +379,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) vzeta.nelement_global == 1 && vr.ngrid == vr.nelement_global == 1) num_diss_params = setup_numerical_dissipation( get(scan_input, "numerical_dissipation", Dict{String,Any}()), is_1V) - + bc_params = setup_boundary_parameters( + get(scan_input, "boundary_condition_parameters", Dict{String,Any}())) # vperp.bc is set here (a bit out of place) so that we can use # num_diss_params.vperp_dissipation_coefficient to set the default. vperp.bc = get(scan_input, "vperp_bc", @@ -555,7 +557,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) vpa, vpa_spectral, vperp, vperp_spectral, gyrophase, gyrophase_spectral, vz, vz_spectral, vr, vr_spectral, vzeta, vzeta_spectral, composition, species_immutable, collisions, geometry, drive_immutable, - external_source_settings, num_diss_params, manufactured_solns_input) + external_source_settings, num_diss_params, manufactured_solns_input, bc_params) println(io, "\nAll inputs returned from mk_input():") println(io, all_inputs) close(io) diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index 6bb5d7f26..67ed9bdc6 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -375,7 +375,7 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral, vpa.bc, z.bc, r.bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_advect, vperp_advect, z_advect, r_advect, composition, scratch_dummy, advance.r_diffusion, - advance.vpa_diffusion, advance.vperp_diffusion) + advance.vpa_diffusion, advance.vperp_diffusion, boundary_distributions.bc_parameters) # Ensure normalised pdf exactly obeys integral constraints if evolving moments begin_s_r_z_region() @loop_s_r_z is ir iz begin @@ -1272,7 +1272,7 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v boundary_distributions.pdf_rboundary_charged, vpa.bc, z.bc, r.bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_advect, vperp_advect, z_advect, r_advect, composition, scratch_dummy, - advance.r_diffusion, advance.vpa_diffusion, advance.vperp_diffusion) + advance.r_diffusion, advance.vpa_diffusion, advance.vperp_diffusion, boundary_distributions.bc_parameters) if moments.evolve_density && moments.enforce_conservation begin_s_r_z_region() From 71abdbbfa93da61b5c7491119878e37a5f9f00da Mon Sep 17 00:00:00 2001 From: John Omotani Date: Thu, 5 Sep 2024 21:28:24 +0100 Subject: [PATCH 05/14] Store boundary_parameters in coordinate structs This reduces the number of function arguments that are needed overall, so keeps the code slightly simpler. --- moment_kinetics/src/coordinates.jl | 31 +++++++++++++--- moment_kinetics/src/initial_conditions.jl | 35 ++++++------------- moment_kinetics/src/load_data.jl | 12 +++---- moment_kinetics/src/moment_kinetics.jl | 4 +-- moment_kinetics/src/moment_kinetics_input.jl | 21 +++++------ moment_kinetics/src/time_advance.jl | 4 +-- moment_kinetics/test/calculus_tests.jl | 30 ++++++++-------- moment_kinetics/test/fokker_planck_tests.jl | 4 +-- moment_kinetics/test/interpolation_tests.jl | 2 +- .../test/velocity_integral_tests.jl | 8 ++--- moment_kinetics/test/wall_bc_tests.jl | 2 +- 11 files changed, 78 insertions(+), 75 deletions(-) diff --git a/moment_kinetics/src/coordinates.jl b/moment_kinetics/src/coordinates.jl index 936695a10..72850aec1 100644 --- a/moment_kinetics/src/coordinates.jl +++ b/moment_kinetics/src/coordinates.jl @@ -12,6 +12,7 @@ using ..calculus: derivative! using ..chebyshev: scaled_chebyshev_grid, scaled_chebyshev_radau_grid, setup_chebyshev_pseudospectral using ..finite_differences: finite_difference_info using ..gauss_legendre: scaled_gauss_legendre_lobatto_grid, scaled_gauss_legendre_radau_grid, setup_gausslegendre_pseudospectral +using ..input_structs using ..quadrature: composite_simpson_weights using ..input_structs: advection_input using ..moment_kinetics_structs: null_spatial_dimension_info, null_velocity_dimension_info @@ -21,7 +22,7 @@ using MPI """ structure containing basic information related to coordinates """ -struct coordinate +struct coordinate{Tbparams} # name is the name of the variable associated with this coordiante name::String # n_global is the total number of grid points associated with this coordinate @@ -62,6 +63,9 @@ struct coordinate cheb_option::String # bc is the boundary condition option for this coordinate bc::String + # struct containing some parameters that may be used for applying the boundary + # condition. + boundary_parameters::Tbparams # wgts contains the integration weights associated with each grid point wgts::Array{mk_float,1} # uniform_grid contains location of grid points mapped to a uniform grid @@ -105,7 +109,7 @@ create arrays associated with a given coordinate, setup the coordinate grid, and populate the coordinate structure containing all of this information """ -function define_coordinate(input, parallel_io::Bool=false; init_YY::Bool=true) +function define_coordinate(input, input_dict, parallel_io::Bool=false; init_YY::Bool=true) # total number of grid points is ngrid for the first element # plus ngrid-1 unique points for each additional element due # to the repetition of a point at the element boundary @@ -130,6 +134,22 @@ function define_coordinate(input, parallel_io::Bool=false; init_YY::Bool=true) imin, imax, igrid, input.discretization, input.name) # calculate the widths of the cells between neighboring grid points cell_width = grid_spacing(grid, n_local) + # Get some parameters that may be used for the boundary condition + if input_dict === nothing + boundary_parameters = nothing + else + boundary_parameters_defaults = Dict{Symbol,Any}() + if input.name == "z" + # parameter controlling the cutoff of the ion distribution function in the vpa + # domain at the wall in z + boundary_parameters_defaults[:epsz] = 0.0 + end + boundary_parameters_input = set_defaults_and_check_section!( + input_dict, "$(input.name)_boundary_condition_parameters"; + boundary_parameters_defaults... + ) + boundary_parameters = Dict_to_NamedTuple(boundary_parameters_input) + end # duniform_dgrid is the local derivative of the uniform grid with respect to # the coordinate grid duniform_dgrid = allocate_float(input.ngrid, input.nelement_local) @@ -164,9 +184,10 @@ function define_coordinate(input, parallel_io::Bool=false; init_YY::Bool=true) coord = coordinate(input.name, n_global, n_local, input.ngrid, input.nelement_global, input.nelement_local, input.nrank, input.irank, input.L, grid, cell_width, igrid, ielement, imin, imax, igrid_full, input.discretization, input.fd_option, input.cheb_option, - input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch), - scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm, - local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option) + input.bc, boundary_parameters, wgts, uniform_grid, duniform_dgrid, scratch, + copy(scratch), copy(scratch), scratch_2d, copy(scratch_2d), advection, + send_buffer, receive_buffer, input.comm, local_io_range, global_io_range, + element_scale, element_shift, input.element_spacing_option) if coord.n == 1 && occursin("v", coord.name) spectral = null_velocity_dimension_info() diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index 78f916248..0d6e4c4ed 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -6,7 +6,6 @@ export allocate_pdf_and_moments export init_pdf_and_moments! export enforce_boundary_conditions! export enforce_neutral_boundary_conditions! -export setup_boundary_parameters # functional testing export create_boundary_distributions @@ -39,19 +38,6 @@ using ..manufactured_solns: manufactured_solutions using MPI -Base.@kwdef struct boundary_condition_parameters - # parameter controlling the cutoff of the ion distribution function - # in the vpa domain at the wall in z - epsz::mk_float = 0.0 -end - -function setup_boundary_parameters(input_section::Dict) - - input = Dict(Symbol(k)=>v for (k,v) in input_section) - - return boundary_condition_parameters(; input...) -end - """ """ struct pdf_substruct{n_distribution} @@ -91,7 +77,6 @@ struct boundary_distributions_struct pdf_rboundary_charged::MPISharedArray{mk_float,5} # neutral particle r boundary values (vz,vr,vzeta,z,r,s) pdf_rboundary_neutral::MPISharedArray{mk_float,6} - bc_parameters::boundary_condition_parameters end """ @@ -99,7 +84,7 @@ Creates the structs for the pdf and the velocity-space moments """ function allocate_pdf_and_moments(composition, r, z, vperp, vpa, vzeta, vr, vz, evolve_moments, collisions, external_source_settings, - numerical_dissipation, boundary_parameters::boundary_condition_parameters) + numerical_dissipation) pdf = create_pdf(composition, r, z, vperp, vpa, vzeta, vr, vz) # create the 'moments' struct that contains various v-space moments and other @@ -132,7 +117,7 @@ function allocate_pdf_and_moments(composition, r, z, vperp, vpa, vzeta, vr, vz, evolve_moments.parallel_pressure) boundary_distributions = create_boundary_distributions(vz, vr, vzeta, vpa, vperp, z, - composition, boundary_parameters) + composition) return pdf, moments, boundary_distributions end @@ -1121,7 +1106,7 @@ Allocate arrays for distributions to be applied as boundary conditions to the pd various boundaries. Also initialise the Knudsen cosine distribution here so it can be used when initialising the neutral pdf. """ -function create_boundary_distributions(vz, vr, vzeta, vpa, vperp, z, composition, boundary_parameters::boundary_condition_parameters) +function create_boundary_distributions(vz, vr, vzeta, vpa, vperp, z, composition) zero = 1.0e-14 #initialise knudsen distribution for neutral wall bc @@ -1135,7 +1120,7 @@ function create_boundary_distributions(vz, vr, vzeta, vpa, vperp, z, composition pdf_rboundary_neutral = allocate_shared_float(vz.n, vr.n, vzeta.n, z.n, 2, composition.n_neutral_species) - return boundary_distributions_struct(knudsen_cosine, pdf_rboundary_charged, pdf_rboundary_neutral, boundary_parameters) + return boundary_distributions_struct(knudsen_cosine, pdf_rboundary_charged, pdf_rboundary_neutral) end function init_boundary_distributions!(boundary_distributions, pdf, vz, vr, vzeta, vpa, vperp, z, r, composition) @@ -1151,7 +1136,7 @@ also enforce boundary conditions in z on all separately evolved velocity space m """ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, phi, moments, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, r_adv, composition, scratch_dummy, - r_diffusion, vpa_diffusion, vperp_diffusion, bc_parameters) + r_diffusion, vpa_diffusion, vperp_diffusion) if vpa.n > 1 begin_s_r_z_vperp_region() @loop_s_r_z_vperp is ir iz ivperp begin @@ -1172,7 +1157,7 @@ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, phi, momen # enforce the z BC on the evolved velocity space moments of the pdf @views enforce_z_boundary_condition_moments!(density, moments, z_bc) @views enforce_z_boundary_condition!(f, density, upar, ppar, phi, moments, z_bc, z_adv, z, - vperp, vpa, composition, bc_parameters, + vperp, vpa, composition, scratch_dummy.buffer_vpavperprs_1, scratch_dummy.buffer_vpavperprs_2, scratch_dummy.buffer_vpavperprs_3, scratch_dummy.buffer_vpavperprs_4) @@ -1187,11 +1172,11 @@ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, phi, momen end function enforce_boundary_conditions!(fvec_out::scratch_pdf, moments, fields::em_fields_struct, f_r_bc, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, r_adv, composition, scratch_dummy, - r_diffusion, vpa_diffusion, vperp_diffusion, bc_parameters) + r_diffusion, vpa_diffusion, vperp_diffusion) enforce_boundary_conditions!(fvec_out.pdf, f_r_bc, fvec_out.density, fvec_out.upar, fvec_out.ppar, fields.phi, moments, vpa_bc, z_bc, r_bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_adv, vperp_adv, z_adv, - r_adv, composition, scratch_dummy, r_diffusion, vpa_diffusion, vperp_diffusion, bc_parameters) + r_adv, composition, scratch_dummy, r_diffusion, vpa_diffusion, vperp_diffusion) end """ @@ -1248,7 +1233,7 @@ end enforce boundary conditions on charged particle f in z """ function enforce_z_boundary_condition!(pdf, density, upar, ppar, phi, moments, bc::String, adv, - z, vperp, vpa, composition, bc_parameters, end1::AbstractArray{mk_float,4}, + z, vperp, vpa, composition, end1::AbstractArray{mk_float,4}, end2::AbstractArray{mk_float,4}, buffer1::AbstractArray{mk_float,4}, buffer2::AbstractArray{mk_float,4}) # this block ensures periodic BC can be supported with distributed memory MPI @@ -1310,7 +1295,7 @@ function enforce_z_boundary_condition!(pdf, density, upar, ppar, phi, moments, b @loop_r ir begin @views enforce_zero_incoming_bc!(pdf[:,:,:,ir,is], adv[is].speed[:,:,:,ir], z, zero, phi[:,ir], - bc_parameters.epsz) + z.boundary_parameters.epsz) end end end diff --git a/moment_kinetics/src/load_data.jl b/moment_kinetics/src/load_data.jl index 0250d667a..dcaad8b38 100644 --- a/moment_kinetics/src/load_data.jl +++ b/moment_kinetics/src/load_data.jl @@ -274,7 +274,7 @@ function load_coordinate_data(fid, name; printout=false, irank=nothing, nrank=no discretization, fd_option, cheb_option, bc, advection_input("", 0.0, 0.0, 0.0), MPI.COMM_NULL, element_spacing_option) - coord, spectral = define_coordinate(input, parallel_io) + coord, spectral = define_coordinate(input, nothing, parallel_io) return coord, spectral, chunk_size end @@ -2454,11 +2454,11 @@ function get_run_info_no_setup(run_dir::Union{AbstractString,Tuple{AbstractStrin dummy_input = grid_input("dummy", 1, 1, 1, 1, 0, 1.0, "chebyshev_pseudospectral", "", "", "periodic", dummy_adv_input, dummy_comm, "uniform") - vzeta, vzeta_spectral = define_coordinate(dummy_input) + vzeta, vzeta_spectral = define_coordinate(dummy_input, nothing) vzeta_chunk_size = 1 - vr, vr_spectral = define_coordinate(dummy_input) + vr, vr_spectral = define_coordinate(dummy_input, nothing) vr_chunk_size = 1 - vz, vz_spectral = define_coordinate(dummy_input) + vz, vz_spectral = define_coordinate(dummy_input, nothing) vz_chunk_size = 1 end end @@ -2981,8 +2981,8 @@ function construct_global_zr_coords(r_local, z_local) coord_local.advection, MPI.COMM_NULL, coord_local.element_spacing_option) end - r_global, r_global_spectral = define_coordinate(make_global_input(r_local)) - z_global, z_global_spectral = define_coordinate(make_global_input(z_local)) + r_global, r_global_spectral = define_coordinate(make_global_input(r_local), nothing) + z_global, z_global_spectral = define_coordinate(make_global_input(z_local), nothing) return r_global, r_global_spectral, z_global, z_global_spectral end diff --git a/moment_kinetics/src/moment_kinetics.jl b/moment_kinetics/src/moment_kinetics.jl index dd7427255..c305479d8 100644 --- a/moment_kinetics/src/moment_kinetics.jl +++ b/moment_kinetics/src/moment_kinetics.jl @@ -289,7 +289,7 @@ function setup_moment_kinetics(input_dict::AbstractDict; vperp, vperp_spectral, gyrophase, gyrophase_spectral, vz, vz_spectral, vr, vr_spectral, vzeta, vzeta_spectral, composition, species, collisions, geometry, drive_input, external_source_settings, num_diss_params, - manufactured_solns_input, bc_params = input + manufactured_solns_input = input # Create loop range variables for shared-memory-parallel loops if debug_loop_type === nothing @@ -315,7 +315,7 @@ function setup_moment_kinetics(input_dict::AbstractDict; pdf, moments, boundary_distributions = allocate_pdf_and_moments(composition, r, z, vperp, vpa, vzeta, vr, vz, evolve_moments, collisions, external_source_settings, - num_diss_params, bc_params) + num_diss_params) if restart === false restarting = false diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 686ea504c..2bb0d4076 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -18,7 +18,6 @@ using ..krook_collisions: setup_krook_collisions using ..finite_differences: fd_check_option using ..input_structs using ..numerical_dissipation: setup_numerical_dissipation -using ..initial_conditions: setup_boundary_parameters using ..reference_parameters using ..geo: init_magnetic_geometry, setup_geometry_input @@ -379,8 +378,6 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) vzeta.nelement_global == 1 && vr.ngrid == vr.nelement_global == 1) num_diss_params = setup_numerical_dissipation( get(scan_input, "numerical_dissipation", Dict{String,Any}()), is_1V) - bc_params = setup_boundary_parameters( - get(scan_input, "boundary_condition_parameters", Dict{String,Any}())) # vperp.bc is set here (a bit out of place) so that we can use # num_diss_params.vperp_dissipation_coefficient to set the default. vperp.bc = get(scan_input, "vperp_bc", @@ -517,21 +514,21 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) Dict(Symbol(k)=>v for (k,v) in io_settings)...) # initialize z grid and write grid point locations to file - z, z_spectral = define_coordinate(z_immutable, io_immutable.parallel_io) + z, z_spectral = define_coordinate(z_immutable, scan_input, io_immutable.parallel_io) # initialize r grid and write grid point locations to file - r, r_spectral = define_coordinate(r_immutable, io_immutable.parallel_io) + r, r_spectral = define_coordinate(r_immutable, scan_input, io_immutable.parallel_io) # initialize vpa grid and write grid point locations to file - vpa, vpa_spectral = define_coordinate(vpa_immutable, io_immutable.parallel_io) + vpa, vpa_spectral = define_coordinate(vpa_immutable, scan_input, io_immutable.parallel_io) # initialize vperp grid and write grid point locations to file - vperp, vperp_spectral = define_coordinate(vperp_immutable, io_immutable.parallel_io) + vperp, vperp_spectral = define_coordinate(vperp_immutable, scan_input, io_immutable.parallel_io) # initialize gyrophase grid and write grid point locations to file - gyrophase, gyrophase_spectral = define_coordinate(gyrophase_immutable, io_immutable.parallel_io) + gyrophase, gyrophase_spectral = define_coordinate(gyrophase_immutable, scan_input, io_immutable.parallel_io) # initialize vz grid and write grid point locations to file - vz, vz_spectral = define_coordinate(vz_immutable, io_immutable.parallel_io) + vz, vz_spectral = define_coordinate(vz_immutable, scan_input, io_immutable.parallel_io) # initialize vr grid and write grid point locations to file - vr, vr_spectral = define_coordinate(vr_immutable, io_immutable.parallel_io) + vr, vr_spectral = define_coordinate(vr_immutable, scan_input, io_immutable.parallel_io) # initialize vr grid and write grid point locations to file - vzeta, vzeta_spectral = define_coordinate(vzeta_immutable, io_immutable.parallel_io) + vzeta, vzeta_spectral = define_coordinate(vzeta_immutable, scan_input, io_immutable.parallel_io) external_source_settings = setup_external_sources!(scan_input, r, z) @@ -562,7 +559,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) vpa, vpa_spectral, vperp, vperp_spectral, gyrophase, gyrophase_spectral, vz, vz_spectral, vr, vr_spectral, vzeta, vzeta_spectral, composition, species_immutable, collisions, geometry, drive_immutable, - external_source_settings, num_diss_params, manufactured_solns_input, bc_params) + external_source_settings, num_diss_params, manufactured_solns_input) println(io, "\nAll inputs returned from mk_input():") println(io, all_inputs) close(io) diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index ee8806f9e..653c467fb 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -378,7 +378,7 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, vz_spectral, vpa.bc, z.bc, r.bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_advect, vperp_advect, z_advect, r_advect, composition, scratch_dummy, advance.r_diffusion, - advance.vpa_diffusion, advance.vperp_diffusion, boundary_distributions.bc_parameters) + advance.vpa_diffusion, advance.vperp_diffusion) # Ensure normalised pdf exactly obeys integral constraints if evolving moments begin_s_r_z_region() @loop_s_r_z is ir iz begin @@ -1275,7 +1275,7 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v boundary_distributions.pdf_rboundary_charged, vpa.bc, z.bc, r.bc, vpa, vperp, z, r, vpa_spectral, vperp_spectral, vpa_advect, vperp_advect, z_advect, r_advect, composition, scratch_dummy, - advance.r_diffusion, advance.vpa_diffusion, advance.vperp_diffusion, boundary_distributions.bc_parameters) + advance.r_diffusion, advance.vpa_diffusion, advance.vperp_diffusion) if moments.evolve_density && moments.enforce_conservation begin_s_r_z_region() diff --git a/moment_kinetics/test/calculus_tests.jl b/moment_kinetics/test/calculus_tests.jl index be8ac68ac..fddba93fd 100644 --- a/moment_kinetics/test/calculus_tests.jl +++ b/moment_kinetics/test/calculus_tests.jl @@ -43,7 +43,7 @@ function runtests() discretization, fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' - x, spectral = define_coordinate(input) + x, spectral = define_coordinate(input, nothing) # create array for the function f(x) to be differentiated/integrated f = Array{Float64,1}(undef, x.n) # create array for the derivative df/dx @@ -93,7 +93,7 @@ function runtests() "finite_difference", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' - x, spectral = define_coordinate(input) + x, spectral = define_coordinate(input, nothing) # create array for the derivative df/dx and the expected result df = Array{Float64,1}(undef, x.n) @@ -143,7 +143,7 @@ function runtests() "finite_difference", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' - x, spectral = define_coordinate(input) + x, spectral = define_coordinate(input, nothing) # create array for the derivative df/dx and the expected result df = Array{Float64,1}(undef, x.n) @@ -189,7 +189,7 @@ function runtests() "finite_difference", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' - x, spectral = define_coordinate(input) + x, spectral = define_coordinate(input, nothing) # create array for the derivative df/dx and the expected result df = Array{Float64,1}(undef, x.n) @@ -243,7 +243,7 @@ function runtests() "finite_difference", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' - x, spectral = define_coordinate(input) + x, spectral = define_coordinate(input, nothing) # create array for the derivative df/dx and the expected result df = Array{Float64,1}(undef, x.n) @@ -459,7 +459,7 @@ function runtests() "chebyshev_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. - x, spectral = define_coordinate(input) + x, spectral = define_coordinate(input, nothing) offset = randn(rng) f = @. sinpi(2.0 * x.grid / L) + offset @@ -655,7 +655,7 @@ function runtests() "chebyshev_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. - x, spectral = define_coordinate(input) + x, spectral = define_coordinate(input, nothing) offset = randn(rng) f = @. sinpi(2.0 * x.grid / L) + offset @@ -698,7 +698,7 @@ function runtests() "chebyshev_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. - x, spectral = define_coordinate(input) + x, spectral = define_coordinate(input, nothing) # test polynomials up to order ngrid-1 for n ∈ 0:ngrid-1 # create array for the function f(x) to be differentiated/integrated @@ -748,7 +748,7 @@ function runtests() "chebyshev_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. - x, spectral = define_coordinate(input) + x, spectral = define_coordinate(input, nothing) # test polynomials up to order ngrid-1 for n ∈ 0:ngrid-1 # create array for the function f(x) to be differentiated/integrated @@ -884,7 +884,7 @@ function runtests() "gausslegendre_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. - x, spectral = define_coordinate(input,init_YY=false) + x, spectral = define_coordinate(input, nothing,init_YY=false) offset = randn(rng) f = @. sinpi(2.0 * x.grid / L) + offset @@ -1001,7 +1001,7 @@ function runtests() "gausslegendre_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. - x, spectral = define_coordinate(input,init_YY=false) + x, spectral = define_coordinate(input, nothing,init_YY=false) offset = randn(rng) f = @. sinpi(2.0 * x.grid / L) + offset @@ -1045,7 +1045,7 @@ function runtests() "gausslegendre_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. - x, spectral = define_coordinate(input,init_YY=false) + x, spectral = define_coordinate(input, nothing,init_YY=false) # test polynomials up to order ngrid-1 for n ∈ 0:ngrid-1 # create array for the function f(x) to be differentiated/integrated @@ -1096,7 +1096,7 @@ function runtests() "gausslegendre_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. - x, spectral = define_coordinate(input,init_YY=false) + x, spectral = define_coordinate(input, nothing,init_YY=false) # test polynomials up to order ngrid-1 for n ∈ 0:ngrid-1 # create array for the function f(x) to be differentiated/integrated @@ -1310,7 +1310,7 @@ function runtests() "chebyshev_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. - x, spectral = define_coordinate(input) + x, spectral = define_coordinate(input, nothing) offset = randn(rng) f = @. sinpi(2.0 * x.grid / L) + offset @@ -1418,7 +1418,7 @@ function runtests() "gausslegendre_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'x' and info for derivatives, etc. - x, spectral = define_coordinate(input,init_YY=false) + x, spectral = define_coordinate(input, nothing,init_YY=false) offset = randn(rng) f = @. sinpi(2.0 * x.grid / L) + offset diff --git a/moment_kinetics/test/fokker_planck_tests.jl b/moment_kinetics/test/fokker_planck_tests.jl index 54266f518..375554adb 100644 --- a/moment_kinetics/test/fokker_planck_tests.jl +++ b/moment_kinetics/test/fokker_planck_tests.jl @@ -51,8 +51,8 @@ function create_grids(ngrid,nelement_vpa,nelement_vperp; #println("made inputs") #println("vpa: ngrid: ",ngrid," nelement: ",nelement_local_vpa, " Lvpa: ",Lvpa) #println("vperp: ngrid: ",ngrid," nelement: ",nelement_local_vperp, " Lvperp: ",Lvperp) - vpa, vpa_spectral = define_coordinate(vpa_input) - vperp, vperp_spectral = define_coordinate(vperp_input) + vpa, vpa_spectral = define_coordinate(vpa_input, nothing) + vperp, vperp_spectral = define_coordinate(vperp_input, nothing) # Set up MPI initialize_comms!() diff --git a/moment_kinetics/test/interpolation_tests.jl b/moment_kinetics/test/interpolation_tests.jl index 1951183a1..ebc349c7b 100644 --- a/moment_kinetics/test/interpolation_tests.jl +++ b/moment_kinetics/test/interpolation_tests.jl @@ -44,7 +44,7 @@ function runtests() discretization, fd_option, cheb_option, bc, adv_input, comm, element_spacing_option) # create the coordinate struct 'z' - z, spectral = define_coordinate(input) + z, spectral = define_coordinate(input, nothing) test_grid = [z for z in range(-zlim, zlim, length=ntest)] diff --git a/moment_kinetics/test/velocity_integral_tests.jl b/moment_kinetics/test/velocity_integral_tests.jl index 2ceaa7e86..d1431b6ae 100644 --- a/moment_kinetics/test/velocity_integral_tests.jl +++ b/moment_kinetics/test/velocity_integral_tests.jl @@ -45,10 +45,10 @@ function runtests() irank, Lvperp, discretization, fd_option, cheb_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) - vz, vz_spectral = define_coordinate(vz_input) - vr, vr_spectral = define_coordinate(vr_input) + vpa, vpa_spectral = define_coordinate(vpa_input, nothing) + vperp, vperp_spectral = define_coordinate(vperp_input, nothing) + vz, vz_spectral = define_coordinate(vz_input, nothing) + vr, vr_spectral = define_coordinate(vr_input, nothing) dfn = allocate_float(vpa.n,vperp.n) dfn1D = allocate_float(vz.n, vr.n) diff --git a/moment_kinetics/test/wall_bc_tests.jl b/moment_kinetics/test/wall_bc_tests.jl index b87e26fd9..c6c1988ce 100644 --- a/moment_kinetics/test/wall_bc_tests.jl +++ b/moment_kinetics/test/wall_bc_tests.jl @@ -216,7 +216,7 @@ function run_test(test_input, expected_phi, tolerance; args...) test_input["z_nelement"], nrank_per_block, irank, 1.0, test_input["z_discretization"], "", cheb_option, test_input["z_bc"], adv_input, comm, test_input["z_element_spacing_option"]) - z, z_spectral = define_coordinate(input) + z, z_spectral = define_coordinate(input, nothing) # Cross comparison of all discretizations to same benchmark if test_input["z_element_spacing_option"] == "uniform" From 8fed26a65654394e855d12e1d87b1b7be72b80e0 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Thu, 5 Sep 2024 22:04:25 +0100 Subject: [PATCH 06/14] Include modification to where wall BC is imposed in moment-kinetic bc --- moment_kinetics/src/boundary_conditions.jl | 18 ++++++++++++------ moment_kinetics/src/time_advance.jl | 3 ++- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/moment_kinetics/src/boundary_conditions.jl b/moment_kinetics/src/boundary_conditions.jl index 6597749ff..43ebfafef 100644 --- a/moment_kinetics/src/boundary_conditions.jl +++ b/moment_kinetics/src/boundary_conditions.jl @@ -177,7 +177,8 @@ function enforce_z_boundary_condition!(pdf, density, upar, ppar, phi, moments, b @loop_r ir begin @views enforce_zero_incoming_bc!( pdf[:,:,:,ir,is], z, vpa, density[:,ir,is], upar[:,ir,is], - ppar[:,ir,is], moments.evolve_upar, moments.evolve_ppar, zero) + ppar[:,ir,is], moments.evolve_upar, moments.evolve_ppar, zero, + phi[:,ir]) end else @loop_r ir begin @@ -427,27 +428,32 @@ function enforce_zero_incoming_bc!(pdf, speed, z, zero, phi, epsz) end end function get_ion_z_boundary_cutoff_indices(density, upar, ppar, evolve_upar, evolve_ppar, - z, vpa, zero) + z, vpa, zero, phi) + epsz = z.boundary_parameters.epsz if z.irank == 0 + deltaphi = phi[2] - phi[1] + vcut = deltaphi > 0 ? sqrt(deltaphi)*(epsz^0.25) : 0.0 vth = sqrt(2.0*(ppar[1]/density[1])) @. vpa.scratch = vpagrid_to_dzdt(vpa.grid, vth, upar[1], evolve_ppar, evolve_upar) - last_negative_vpa_ind = searchsortedlast(vpa.scratch, -zero) + last_negative_vpa_ind = searchsortedlast(vpa.scratch, min(-zero, -vcut)) else last_negative_vpa_ind = nothing end if z.irank == z.nrank - 1 + deltaphi = phi[end-1] - phi[end] + vcut = deltaphi > 0 ? sqrt(deltaphi)*(epsz^0.25) : 0.0 vth = sqrt(2.0*(ppar[end]/density[end])) @. vpa.scratch2 = vpagrid_to_dzdt(vpa.grid, vth, upar[end], evolve_ppar, evolve_upar) - first_positive_vpa_ind = searchsortedfirst(vpa.scratch2, zero) + first_positive_vpa_ind = searchsortedfirst(vpa.scratch2, max(zero, vcut)) else first_positive_vpa_ind = nothing end return last_negative_vpa_ind, first_positive_vpa_ind end function enforce_zero_incoming_bc!(pdf, z::coordinate, vpa::coordinate, density, upar, - ppar, evolve_upar, evolve_ppar, zero) + ppar, evolve_upar, evolve_ppar, zero, phi) if z.irank != 0 && z.irank != z.nrank - 1 # No z-boundary in this block return nothing @@ -461,7 +467,7 @@ function enforce_zero_incoming_bc!(pdf, z::coordinate, vpa::coordinate, density, # absolute velocity at left boundary last_negative_vpa_ind, first_positive_vpa_ind = get_ion_z_boundary_cutoff_indices(density, upar, ppar, evolve_upar, evolve_ppar, - z, vpa, zero) + z, vpa, zero, phi) if z.irank == 0 pdf[last_negative_vpa_ind+1:end, :, 1] .= 0.0 end diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index 3b278462e..a244b293a 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -2659,11 +2659,12 @@ function adaptive_timestep_update!(scratch, scratch_implicit, scratch_electron, density = @view scratch[t_params.n_rk_stages+1].density[:,ir,is] upar = @view scratch[t_params.n_rk_stages+1].upar[:,ir,is] ppar = @view scratch[t_params.n_rk_stages+1].ppar[:,ir,is] + phi = fields.phi[:,ir] last_negative_vpa_ind, first_positive_vpa_ind = get_ion_z_boundary_cutoff_indices(density, upar, ppar, moments.evolve_upar, moments.evolve_ppar, z, vpa, - 1.0e-14) + 1.0e-14, phi) if z.irank == 0 scratch[2].pdf[last_negative_vpa_ind,:,1,ir,is] .= scratch[t_params.n_rk_stages+1].pdf[last_negative_vpa_ind,:,1,ir,is] From ea9ebf144f6d93827b2a2222ead852190c750ca1 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 18 Sep 2024 14:16:37 +0100 Subject: [PATCH 07/14] Update boundary_parameters setup to work with new coordinate init --- moment_kinetics/src/coordinates.jl | 45 +++++++++++++++++------------- 1 file changed, 25 insertions(+), 20 deletions(-) diff --git a/moment_kinetics/src/coordinates.jl b/moment_kinetics/src/coordinates.jl index 7ab861537..b76da95f1 100644 --- a/moment_kinetics/src/coordinates.jl +++ b/moment_kinetics/src/coordinates.jl @@ -196,6 +196,27 @@ function get_coordinate_input(input_dict, name; ignore_MPI=false) # Make a copy so we do not add "name" to the global input_dict coord_input_dict = copy(coord_input_dict) coord_input_dict["name"] = name + + # Get some parameters that may be used for the boundary condition + if input_dict === nothing + boundary_parameters = nothing + else + boundary_parameters_defaults = Dict{Symbol,Any}() + if name == "z" + # parameter controlling the cutoff of the ion distribution function in the vpa + # domain at the wall in z + boundary_parameters_defaults[:epsz] = 0.0 + end + boundary_parameters_input = set_defaults_and_check_section!( + input_dict, "$(name)_boundary_condition_parameters"; + boundary_parameters_defaults... + ) + boundary_parameters = Dict_to_NamedTuple(boundary_parameters_input) + end + + coord_input_dict = deepcopy(coord_input_dict) + coord_input_dict["boundary_parameters"] = boundary_parameters + coord_input = Dict_to_NamedTuple(coord_input_dict) return coord_input @@ -277,22 +298,6 @@ function define_coordinate(coord_input::NamedTuple; parallel_io::Bool=false, coord_input.discretization, coord_input.name) # calculate the widths of the cells between neighboring grid points cell_width = grid_spacing(grid, n_local) - # Get some parameters that may be used for the boundary condition - if input_dict === nothing - boundary_parameters = nothing - else - boundary_parameters_defaults = Dict{Symbol,Any}() - if input.name == "z" - # parameter controlling the cutoff of the ion distribution function in the vpa - # domain at the wall in z - boundary_parameters_defaults[:epsz] = 0.0 - end - boundary_parameters_input = set_defaults_and_check_section!( - input_dict, "$(input.name)_boundary_condition_parameters"; - boundary_parameters_defaults... - ) - boundary_parameters = Dict_to_NamedTuple(boundary_parameters_input) - end # duniform_dgrid is the local derivative of the uniform grid with respect to # the coordinate grid duniform_dgrid = allocate_float(coord_input.ngrid, coord_input.nelement_local) @@ -372,11 +377,11 @@ function define_coordinate(coord_input::NamedTuple; parallel_io::Bool=false, coord_input.nelement, coord_input.nelement_local, nrank, irank, coord_input.L, grid, cell_width, igrid, ielement, imin, imax, igrid_full, coord_input.discretization, coord_input.finite_difference_option, - coord_input.cheb_option, coord_input.bc, boundary_parameters, wgts, uniform_grid, - duniform_dgrid, scratch, copy(scratch), copy(scratch), copy(scratch), + coord_input.cheb_option, coord_input.bc, coord_input.boundary_parameters, wgts, + uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch), copy(scratch), copy(scratch), copy(scratch), copy(scratch), copy(scratch), - scratch_shared, scratch_shared2, scratch_2d, copy(scratch_2d), advection, - send_buffer, receive_buffer, comm, local_io_range, global_io_range, + copy(scratch), scratch_shared, scratch_shared2, scratch_2d, copy(scratch_2d), + advection, send_buffer, receive_buffer, comm, local_io_range, global_io_range, element_scale, element_shift, coord_input.element_spacing_option, element_boundaries, radau_first_element, other_nodes, one_over_denominator) From ff675e507076b284080eba8eec4a671100a38a13 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 20 Sep 2024 23:24:36 +0100 Subject: [PATCH 08/14] Fix x-axis of animation of f/vpa^2 --- .../src/makie_post_processing.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl index 1b4b024e2..418b14e18 100644 --- a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl +++ b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl @@ -4161,8 +4161,8 @@ it might be useful to pass `transform=abs` (to plot the absolute value) or `axis_args` are passed as keyword arguments to `get_1d_ax()`, and from there to the `Axis` constructor. -Any extra `kwargs` are passed to [`plot_1d`](@ref) (which is used to create the plot, as -we have to handle time-varying coordinates so cannot use [`animate_1d`](@ref)). +Any extra `kwargs` are passed to `lines!()` (which is used to create the plot, as we have +to handle time-varying coordinates so cannot use [`animate_1d`](@ref)). """ function animate_f_unnorm_vs_vpa end @@ -4264,16 +4264,16 @@ function animate_f_unnorm_vs_vpa(run_info; f_over_vpa2=false, input=nothing, run_info.evolve_density, run_info.evolve_ppar) if f_over_vpa2 - dzdt = vpagrid_to_dzdt(vcoord.grid, vth[it], upar[it], - run_info.evolve_ppar, run_info.evolve_upar) - dzdt2 = dzdt.^2 - for i ∈ eachindex(dzdt2) - if dzdt2[i] == 0.0 - dzdt2[i] = 1.0 + this_dzdt = vpagrid_to_dzdt(vcoord.grid, vth[it], upar[it], + run_info.evolve_ppar, run_info.evolve_upar) + this_dzdt2 = this_dzdt.^2 + for i ∈ eachindex(this_dzdt2) + if this_dzdt2[i] == 0.0 + this_dzdt2[i] = 1.0 end end - f_unnorm = @. copy(f_unnorm) / dzdt2 + f_unnorm = @. copy(f_unnorm) / this_dzdt2 end return f_unnorm From 67b780fbadfd08f4aa1e68ff0713cf9b80571fde Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 3 Sep 2024 21:24:23 +0100 Subject: [PATCH 09/14] Fix Chodura condition diagnostic bzed is now a 2D variable, not a scalar. --- moment_kinetics/src/analysis.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/src/analysis.jl b/moment_kinetics/src/analysis.jl index fc342e050..7bc3a5037 100644 --- a/moment_kinetics/src/analysis.jl +++ b/moment_kinetics/src/analysis.jl @@ -165,7 +165,7 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, composition, for it ∈ 1:ntime, ir ∈ 1:nr v_parallel = vpagrid_to_dzdt(vpa.grid, vth[1,ir,is,it], upar[1,ir,is,it], evolve_ppar, evolve_upar) - vpabar = @. v_parallel - 0.5 * geometry.rhostar * Er[1,ir,it] / geometry.bzed + vpabar = @. v_parallel - 0.5 * geometry.rhostar * Er[1,ir,it] / geometry.bzed[1,ir] # Get rid of a zero if it is there to avoid a blow up - f should be zero at that # point anyway @@ -187,7 +187,7 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, composition, v_parallel = vpagrid_to_dzdt(vpa.grid, vth[end,ir,is,it], upar[end,ir,is,it], evolve_ppar, evolve_upar) - vpabar = @. v_parallel - 0.5 * geometry.rhostar * Er[end,ir,it] / geometry.bzed + vpabar = @. v_parallel - 0.5 * geometry.rhostar * Er[end,ir,it] / geometry.bzed[end,ir] # Get rid of a zero if it is there to avoid a blow up - f should be zero at that # point anyway From 40c67de226bafc01419c34d004bab4cc684c7a38 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 20 Sep 2024 13:23:55 +0100 Subject: [PATCH 10/14] Allow non-constant electron temperature in check_Chodura_condition() In runs that do not use a Boltzmann electron response, the electron temperature is not just a constant. --- .../src/makie_post_processing.jl | 6 ++++-- moment_kinetics/src/analysis.jl | 14 ++++++++++---- .../src/plots_post_processing.jl | 10 ++++++++-- 3 files changed, 22 insertions(+), 8 deletions(-) diff --git a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl index 418b14e18..b2926dc4e 100644 --- a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl +++ b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl @@ -5618,14 +5618,16 @@ function Chodura_condition_plots(run_info; plot_prefix=nothing, axes=nothing) density = get_variable(run_info, "density") upar = get_variable(run_info, "parallel_flow") vth = get_variable(run_info, "thermal_speed") + temp_e = get_variable(run_info, "electron_temperature") Er = get_variable(run_info, "Er") f_lower = get_variable(run_info, "f", iz=1) f_upper = get_variable(run_info, "f", iz=run_info.z.n_global) Chodura_ratio_lower, Chodura_ratio_upper = check_Chodura_condition(run_info.r_local, run_info.z_local, run_info.vperp, - run_info.vpa, density, upar, vth, run_info.composition, - Er, run_info.geometry, run_info.z.bc, nothing; + run_info.vpa, density, upar, vth, temp_e, + run_info.composition, Er, run_info.geometry, + run_info.z.bc, nothing; evolve_density=run_info.evolve_density, evolve_upar=run_info.evolve_upar, evolve_ppar=run_info.evolve_ppar, diff --git a/moment_kinetics/src/analysis.jl b/moment_kinetics/src/analysis.jl index 7bc3a5037..8179fcf6e 100644 --- a/moment_kinetics/src/analysis.jl +++ b/moment_kinetics/src/analysis.jl @@ -99,8 +99,8 @@ Note that `integrate_over_vspace()` includes the 1/sqrt(pi) factor already. If `ir0` is passed, only load the data for as single r-point (to save memory). """ -function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, composition, Er, - geometry, z_bc, nblocks, run_name=nothing, +function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, temp_e, composition, + Er, geometry, z_bc, nblocks, run_name=nothing, it0::Union{Nothing, mk_int}=nothing, ir0::Union{Nothing, mk_int}=nothing; f_lower=nothing, f_upper=nothing, @@ -127,6 +127,12 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, composition, else nr = 1 end + + if temp_e === nothing + # Assume this is from a Boltzmann electron response simulation + temp_e = fill(composition.T_e, 2, nr, ntime) + end + lower_result = zeros(nr, ntime) upper_result = zeros(nr, ntime) if f_lower !== nothing || f_upper !== nothing @@ -183,7 +189,7 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, composition, println("result lower ", lower_result[ir,it]) end - lower_result[ir,it] *= 0.5 * composition.T_e / dens[1,ir,is,it] + lower_result[ir,it] *= 0.5 * temp_e[1,ir,it] / dens[1,ir,is,it] v_parallel = vpagrid_to_dzdt(vpa.grid, vth[end,ir,is,it], upar[end,ir,is,it], evolve_ppar, evolve_upar) @@ -205,7 +211,7 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, composition, println("result upper ", upper_result[ir,it]) end - upper_result[ir,it] *= 0.5 * composition.T_e / dens[end,ir,is,it] + upper_result[ir,it] *= 0.5 * temp_e[end,ir,it] / dens[end,ir,is,it] end println("final Chodura results result ", lower_result[1,end], " ", upper_result[1,end]) diff --git a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl index caa9d1f6f..327926fb9 100644 --- a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl +++ b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl @@ -776,10 +776,13 @@ function analyze_and_plot_data(prefix...; run_index=nothing) end if pp.diagnostics_chodura_t + # At the moment electron data is not loaded, so we pass `nothing` for the `temp_e` + # argument to use the Boltzmann electron response version of + # check_Chodura_condition(). Chodura_ratio_lower, Chodura_ratio_upper = get_tuple_of_return_values(check_Chodura_condition, r, z, vperp, vpa, density_at_pdf_times, parallel_flow_at_pdf_times, - thermal_speed_at_pdf_times, composition, + thermal_speed_at_pdf_times, nothing, composition, Er_at_pdf_times, geometry, "wall", nblocks, run_names, nothing, ir0, evolve_density=evolve_density, @@ -802,10 +805,13 @@ function analyze_and_plot_data(prefix...; run_index=nothing) end if pp.diagnostics_chodura_r + # At the moment electron data is not loaded, so we pass `nothing` for the `temp_e` + # argument to use the Boltzmann electron response version of + # check_Chodura_condition(). Chodura_ratio_lower, Chodura_ratio_upper = get_tuple_of_return_values(check_Chodura_condition, r, z, vperp, vpa, density_at_pdf_times, parallel_flow_at_pdf_times, - thermal_speed_at_pdf_times, composition, + thermal_speed_at_pdf_times, nothing, composition, Er_at_pdf_times, geometry, "wall", nblocks, run_names, ntime_pdfs, nothing, evolve_density=evolve_density, From db0a1e469288a21cea7d4ff7b6f77e901c33c2a7 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 20 Sep 2024 13:37:56 +0100 Subject: [PATCH 11/14] In 'Chodura condition' check, calculate change needed to satisfy it The 'Chodura condition' is often not satisfied. By zero-ing out enough grid points near v_parallel=0, eventually we would get a distribution function at the wall that would satisfy the condition. Add an option to calculate how many grid points this would be. --- .../src/makie_post_processing.jl | 4 +- moment_kinetics/src/analysis.jl | 80 +++++++++++++++++-- 2 files changed, 76 insertions(+), 8 deletions(-) diff --git a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl index b2926dc4e..d835bc0d8 100644 --- a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl +++ b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl @@ -5623,7 +5623,7 @@ function Chodura_condition_plots(run_info; plot_prefix=nothing, axes=nothing) f_lower = get_variable(run_info, "f", iz=1) f_upper = get_variable(run_info, "f", iz=run_info.z.n_global) - Chodura_ratio_lower, Chodura_ratio_upper = + Chodura_ratio_lower, Chodura_ratio_upper, cutoff_lower, cutoff_upper = check_Chodura_condition(run_info.r_local, run_info.z_local, run_info.vperp, run_info.vpa, density, upar, vth, temp_e, run_info.composition, Er, run_info.geometry, @@ -5631,7 +5631,7 @@ function Chodura_condition_plots(run_info; plot_prefix=nothing, axes=nothing) evolve_density=run_info.evolve_density, evolve_upar=run_info.evolve_upar, evolve_ppar=run_info.evolve_ppar, - f_lower=f_lower, f_upper=f_upper) + f_lower=f_lower, f_upper=f_upper, find_extra_offset=true) if input.plot_vs_t if axes === nothing diff --git a/moment_kinetics/src/analysis.jl b/moment_kinetics/src/analysis.jl index 8179fcf6e..68e911efa 100644 --- a/moment_kinetics/src/analysis.jl +++ b/moment_kinetics/src/analysis.jl @@ -6,7 +6,7 @@ export analyze_fields_data export analyze_moments_data export analyze_pdf_data -using ..array_allocation: allocate_float +using ..array_allocation: allocate_float, allocate_int using ..calculus: integral using ..communication using ..coordinates: coordinate @@ -98,6 +98,10 @@ TeN / (2 neN sqrt(pi)) * ∫dvpaN fN / vpaN^2 ≤ 1 Note that `integrate_over_vspace()` includes the 1/sqrt(pi) factor already. If `ir0` is passed, only load the data for as single r-point (to save memory). + +If `find_extra_offset=true` is passed, calculates how many entries of `f_lower`/`f_upper` +adjacent to \$v_∥=0\$ would need to be zero-ed out in order for the condition to be +satisfied. """ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, temp_e, composition, Er, geometry, z_bc, nblocks, run_name=nothing, @@ -105,12 +109,14 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, temp_e, comp ir0::Union{Nothing, mk_int}=nothing; f_lower=nothing, f_upper=nothing, evolve_density=false, evolve_upar=false, - evolve_ppar=false) + evolve_ppar=false, find_extra_offset=false) if z_bc != "wall" return nothing, nothing end + zero = 1.0e-14 + if it0 === nothing ntime = size(Er, 3) t_range = 1:ntime @@ -127,6 +133,8 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, temp_e, comp else nr = 1 end + nvperp = size(f_lower,2) + nvpa = size(f_lower,1) if temp_e === nothing # Assume this is from a Boltzmann electron response simulation @@ -168,7 +176,18 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, temp_e, comp evolve_density, evolve_ppar) f_upper = @views get_unnormalised_f_1d(f_upper, dens[end,:,:,:], vth[end,:,:,:], evolve_density, evolve_ppar) + if find_extra_offset + # Allocate output arrays for the number of entries that would need to be zero-ed + # out. + extra_offset_lower = allocate_int(nr,ntime) + extra_offset_upper = allocate_int(nr,ntime) + cutoff_lower = allocate_float(nr,ntime) + cutoff_upper = allocate_float(nr,ntime) + end for it ∈ 1:ntime, ir ∈ 1:nr + # Lower target + ############## + v_parallel = vpagrid_to_dzdt(vpa.grid, vth[1,ir,is,it], upar[1,ir,is,it], evolve_ppar, evolve_upar) vpabar = @. v_parallel - 0.5 * geometry.rhostar * Er[1,ir,it] / geometry.bzed[1,ir] @@ -176,7 +195,7 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, temp_e, comp # Get rid of a zero if it is there to avoid a blow up - f should be zero at that # point anyway for ivpa ∈ eachindex(vpabar) - if abs(vpabar[ivpa]) < 1.e-14 + if abs(vpabar[ivpa]) < zero vpabar[ivpa] = 1.0 end end @@ -191,6 +210,28 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, temp_e, comp lower_result[ir,it] *= 0.5 * temp_e[1,ir,it] / dens[1,ir,is,it] + if find_extra_offset + if lower_result[ir,it] ≤ 1.0 + extra_offset_lower[ir,it] = 0 + else + integrand = f_lower[:,:,ir,is,it] + for ivperp ∈ 1:nvperp + # 1/sqrt(π) factor to correspond with behaviour of + # `integrate_over_vspace()`. + @. integrand[:,ivperp] *= vpabar^(-2) * vpa.wgts * vperp.wgts[ivperp] / sqrt(π) + end + vperp_integral = @view sum(integrand; dims=2)[:,1] + cumulative_vpa_integral = cumsum(vperp_integral) + cutoff_index = searchsortedfirst(cumulative_vpa_integral, 2.0 * dens[1,ir,is,it] / temp_e[1,ir,it]) - 1 + cutoff_lower[ir,it] = mean(vpabar[cutoff_index:cutoff_index+1]) + vpa_before_zero_index = searchsortedfirst(vpabar, -zero) - 1 + extra_offset_lower[ir,it] = vpa_before_zero_index - cutoff_index + end + end + + # Upper target + ############## + v_parallel = vpagrid_to_dzdt(vpa.grid, vth[end,ir,is,it], upar[end,ir,is,it], evolve_ppar, evolve_upar) vpabar = @. v_parallel - 0.5 * geometry.rhostar * Er[end,ir,it] / geometry.bzed[end,ir] @@ -198,7 +239,7 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, temp_e, comp # Get rid of a zero if it is there to avoid a blow up - f should be zero at that # point anyway for ivpa ∈ eachindex(vpabar) - if abs(vpabar[ivpa]) < 1.e-14 + if abs(vpabar[ivpa]) < zero vpabar[ivpa] = 1.0 end end @@ -212,9 +253,32 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, temp_e, comp end upper_result[ir,it] *= 0.5 * temp_e[end,ir,it] / dens[end,ir,is,it] + + if find_extra_offset + if upper_result[ir,it] ≤ 1.0 + extra_offset_upper[ir,it] = 0 + else + integrand = f_upper[:,:,ir,is,it] + for ivperp ∈ 1:nvperp + # 1/sqrt(π) factor to correspond with behaviour of + # `integrate_over_vspace()`. + @. integrand[:,ivperp] *= vpabar^(-2) * vpa.wgts * vperp.wgts[ivperp] / sqrt(π) + end + vperp_integral = @view sum(integrand; dims=2)[:,1] + cumulative_vpa_integral = reverse(cumsum(reverse(vperp_integral))) + cutoff_index = searchsortedfirst(cumulative_vpa_integral, 2.0 * dens[end,ir,is,it] / temp_e[end,ir,it]; rev=true) + cutoff_upper[ir,it] = mean(vpabar[cutoff_index-1:cutoff_index]) + vpa_after_zero_index = searchsortedlast(vpabar, zero) + 1 + extra_offset_upper[ir,it] = cutoff_index - vpa_after_zero_index + end + end end - println("final Chodura results result ", lower_result[1,end], " ", upper_result[1,end]) + if find_extra_offset + println("final Chodura results ", lower_result[1,end], " (", extra_offset_lower[1,end], ") ", upper_result[1,end], " (", extra_offset_upper[1,end], ")") + else + println("final Chodura results ", lower_result[1,end], " ", upper_result[1,end]) + end if it0 !== nothing && ir0 !== nothing lower_result = lower_result[1,1] @@ -227,7 +291,11 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, temp_e, comp upper_result = @view upper_result[1,:] end - return lower_result, upper_result + if find_extra_offset + return lower_result, upper_result, cutoff_lower, cutoff_upper + else + return lower_result, upper_result + end end """ From 62ca1d285e97af6e8af5ed6bc1c2e1e3fd54a3ae Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 20 Sep 2024 14:28:22 +0100 Subject: [PATCH 12/14] Move plots/animations of f/vpa^2 to Chodura_condition_plots() --- .../src/makie_post_processing.jl | 194 ++++++++++++++++-- 1 file changed, 177 insertions(+), 17 deletions(-) diff --git a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl index d835bc0d8..ff02a1ccf 100644 --- a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl +++ b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl @@ -758,8 +758,11 @@ function _setup_single_input!(this_input_dict::OrderedDict{String,Any}, plot_vs_t=false, plot_vs_r=false, plot_vs_r_t=false, + plot_f_over_vpa2=false, + animate_f_over_vpa2=false, it0=this_input_dict["it0"], ir0=this_input_dict["ir0"], + animation_ext=this_input_dict["animation_ext"], ) set_defaults_and_check_section!( @@ -4656,11 +4659,6 @@ function plot_charged_pdf_2D_at_wall(run_info; plot_prefix, electron=false) save(outfile, fig) end - if !electron - plot_f_unnorm_vs_vpa(run_info; f_over_vpa2=true, input=f_input, is=1, - outfile=plot_prefix * "pdf_unnorm_over_vpa2_$(label)_vs_vpa.pdf") - end - if !is_1V plot_vs_vpa_vperp(run_info, "f$electron_suffix"; is=1, input=f_input, outfile=plot_prefix * "pdf$(electron_suffix)_$(label)_vs_vpa_vperp.pdf") @@ -4759,11 +4757,6 @@ function plot_charged_pdf_2D_at_wall(run_info; plot_prefix, electron=false) save_animation(fig, frame_index, nt, outfile) end - if !electron - animate_f_unnorm_vs_vpa(run_info; f_over_vpa2=true, input=f_input, is=1, - outfile=plot_prefix * "pdf_unnorm_over_vpa2_$(label)_vs_vpa." * input.animation_ext) - end - if !is_1V animate_vs_vpa_vperp(run_info, "f$electron_suffix"; is=1, input=f_input, outfile=plot_prefix * "pdf$(electron_suffix)_$(label)_vs_vpa_vperp." * input.animation_ext) @@ -5551,6 +5544,53 @@ function Chodura_condition_plots(run_info::Tuple; plot_prefix) push!(a, nothing) end end + if input.plot_f_over_vpa2 + println("going to plot f_over_vpa2") + fig, ax = get_1d_ax(title="f/vpa^2 lower wall", xlabel="vpa", ylabel="f / vpa^2") + push!(figs, fig) + for a ∈ axes + push!(a, ax) + end + + fig, ax = get_1d_ax(title="f/vpa^2 upper wall", xlabel="vpa", ylabel="f / vpa^2") + push!(figs, fig) + for a ∈ axes + push!(a, ax) + end + else + push!(figs, nothing) + for a ∈ axes + push!(a, nothing) + end + push!(figs, nothing) + for a ∈ axes + push!(a, nothing) + end + end + if input.animate_f_over_vpa2 + fig, ax = get_1d_ax(title="f/vpa^2 lower wall", xlabel="vpa", ylabel="f / vpa^2") + frame_index = Observable(1) + push!(figs, fig) + for a ∈ axes + push!(a, (ax, frame_index)) + end + + fig, ax = get_1d_ax(title="f/vpa^2 upper wall", xlabel="vpa", ylabel="f / vpa^2") + frame_index = Observable(1) + push!(figs, fig) + for a ∈ axes + push!(a, (ax, frame_index)) + end + else + push!(figs, nothing) + for a ∈ axes + push!(a, nothing) + end + push!(figs, nothing) + for a ∈ axes + push!(a, nothing) + end + end for (ri, ax) ∈ zip(run_info, axes) Chodura_condition_plots(ri; axes=ax) @@ -5559,26 +5599,26 @@ function Chodura_condition_plots(run_info::Tuple; plot_prefix) if input.plot_vs_t fig = figs[1] ax = axes[1][1] - put_legend_right(fig, ax) + put_legend_below(fig, ax) outfile = string(plot_prefix, "Chodura_ratio_lower_vs_t.pdf") save(outfile, fig) fig = figs[2] - ax = axes[2][1] - put_legend_right(fig, ax) + ax = axes[1][2] + put_legend_below(fig, ax) outfile = string(plot_prefix, "Chodura_ratio_upper_vs_t.pdf") save(outfile, fig) end if input.plot_vs_r fig = figs[3] - ax = axes[3][1] - put_legend_right(fig, ax) + ax = axes[1][3] + put_legend_below(fig, ax) outfile = string(plot_prefix, "Chodura_ratio_lower_vs_r.pdf") save(outfile, fig) fig = figs[4] - ax = axes[4][1] - put_legend_right(fig, ax) + ax = axes[1][4] + put_legend_below(fig, ax) outfile = string(plot_prefix, "Chodura_ratio_upper_vs_r.pdf") save(outfile, fig) end @@ -5591,6 +5631,37 @@ function Chodura_condition_plots(run_info::Tuple; plot_prefix) outfile = string(plot_prefix, "Chodura_ratio_upper_vs_r_t.pdf") save(outfile, fig) end + if input.plot_f_over_vpa2 + fig = figs[7] + println("check axes ", axes) + ax = axes[1][7] + put_legend_below(fig, ax) + outfile = string(plot_prefix, "pdf_unnorm_over_vpa2_wall-_vs_vpa.pdf") + save(outfile, fig) + + fig = figs[8] + ax = axes[1][8] + put_legend_below(fig, ax) + outfile = string(plot_prefix, "pdf_unnorm_over_vpa2_wall+_vs_vpa.pdf") + save(outfile, fig) + end + if input.animate_f_over_vpa2 + nt = minimum(ri.nt for ri ∈ run_info) + + fig = figs[9] + ax = axes[1][9][1] + frame_index = axes[1][9][2] + put_legend_below(fig, ax) + outfile = string(plot_prefix, "pdf_unnorm_over_vpa2_wall-_vs_vpa." * input.animation_ext) + save_animation(fig, frame_index, nt, outfile) + + fig = figs[10] + ax = axes[1][10][1] + frame_index = axes[1][10][2] + put_legend_below(fig, ax) + outfile = string(plot_prefix, "pdf_unnorm_over_vpa2_wall+_vs_vpa." * input.animation_ext) + save_animation(fig, frame_index, nt, outfile) + end catch e return makie_post_processing_error_handler( e, @@ -5723,6 +5794,95 @@ function Chodura_condition_plots(run_info; plot_prefix=nothing, axes=nothing) end end + if input.plot_f_over_vpa2 + if axes === nothing + fig, ax, = get_1d_ax(title="f/vpa^2 lower wall", + xlabel="vpa", ylabel="f / vpa^2") + title = nothing + label = "" + else + fig = nothing + ax = axes[7] + label = run_info.run_name + end + f_input = copy(input_dict_dfns["f"]) + f_input["it0"] = input.it0 + f_input["ir0"] = input.ir0 + f_input["iz0"] = 1 + plot_f_unnorm_vs_vpa(run_info; f_over_vpa2=true, input=f_input, is=1, fig=fig, + ax=ax, label=label) + if plot_prefix !== nothing && fig !== nothing + outfile=plot_prefix * "pdf_unnorm_over_vpa2_wall-_vs_vpa.pdf" + save(outfile, fig) + end + + if axes === nothing + fig, ax, = get_1d_ax(title="f/vpa^2 upper wall", + xlabel="vpa", ylabel="f / vpa^2") + title = nothing + label = "" + else + fig = nothing + ax = axes[8] + label = run_info.run_name + end + f_input = copy(input_dict_dfns["f"]) + f_input["it0"] = input.it0 + f_input["ir0"] = input.ir0 + f_input["iz0"] = run_info.z.n + plot_f_unnorm_vs_vpa(run_info; f_over_vpa2=true, input=f_input, is=1, fig=fig, + ax=ax, label=label) + if plot_prefix !== nothing && fig !== nothing + outfile=plot_prefix * "pdf_unnorm_over_vpa2_wall+_vs_vpa.pdf" + save(outfile, fig) + end + end + + if input.animate_f_over_vpa2 + if axes === nothing + fig, ax, = get_1d_ax(title="f/vpa^2 lower wall", + xlabel="vpa", ylabel="f / vpa^2") + frame_index = Observable(1) + title = nothing + label = "" + else + fig = nothing + ax, frame_index = axes[9] + label = run_info.run_name + end + f_input = copy(input_dict_dfns["f"]) + f_input["ir0"] = input.ir0 + f_input["iz0"] = 1 + animate_f_unnorm_vs_vpa(run_info; f_over_vpa2=true, input=f_input, is=1, iz=1, + fig=fig, ax=ax, frame_index=frame_index, label=label) + if plot_prefix !== nothing && fig !== nothing + outfile=plot_prefix * "pdf_unnorm_over_vpa2_wall-_vs_vpa." * input.animation_ext + save_animation(fig, frame_index, run_info.nt, outfile) + end + + if axes === nothing + fig, ax, = get_1d_ax(title="f/vpa^2 upper wall", + xlabel="vpa", ylabel="f / vpa^2") + frame_index = Observable(1) + title = nothing + label = "" + else + fig = nothing + ax, frame_index = axes[10] + label = run_info.run_name + end + f_input = copy(input_dict_dfns["f"]) + f_input["ir0"] = input.ir0 + f_input["iz0"] = run_info.z.n + animate_f_unnorm_vs_vpa(run_info; f_over_vpa2=true, input=f_input, is=1, + iz=run_info.z.n, fig=fig, ax=ax, frame_index=frame_index, + label=label) + if plot_prefix !== nothing && fig !== nothing + outfile=plot_prefix * "pdf_unnorm_over_vpa2_wall+_vs_vpa." * input.animation_ext + save_animation(fig, frame_index, run_info.nt, outfile) + end + end + return nothing end From d0afe87198c31ece5598fb14bc5a89f1b0b395ff Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 20 Sep 2024 16:34:11 +0100 Subject: [PATCH 13/14] Plot cutoff in Chodura condition plots The line shown by the cutoff is the vpa value that, if we integrated f/vpa^2 only up to there, the result would obey the Chodura condition. --- .../makie_post_processing/src/makie_post_processing.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl index ff02a1ccf..c2aff3809 100644 --- a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl +++ b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl @@ -5811,6 +5811,7 @@ function Chodura_condition_plots(run_info; plot_prefix=nothing, axes=nothing) f_input["iz0"] = 1 plot_f_unnorm_vs_vpa(run_info; f_over_vpa2=true, input=f_input, is=1, fig=fig, ax=ax, label=label) + vlines!(ax, cutoff_lower[input.ir0,input.it0]; linestyle=:dash, color=:red) if plot_prefix !== nothing && fig !== nothing outfile=plot_prefix * "pdf_unnorm_over_vpa2_wall-_vs_vpa.pdf" save(outfile, fig) @@ -5832,6 +5833,7 @@ function Chodura_condition_plots(run_info; plot_prefix=nothing, axes=nothing) f_input["iz0"] = run_info.z.n plot_f_unnorm_vs_vpa(run_info; f_over_vpa2=true, input=f_input, is=1, fig=fig, ax=ax, label=label) + vlines!(ax, cutoff_upper[input.ir0,input.it0]; linestyle=:dash, color=:red) if plot_prefix !== nothing && fig !== nothing outfile=plot_prefix * "pdf_unnorm_over_vpa2_wall+_vs_vpa.pdf" save(outfile, fig) @@ -5855,6 +5857,7 @@ function Chodura_condition_plots(run_info; plot_prefix=nothing, axes=nothing) f_input["iz0"] = 1 animate_f_unnorm_vs_vpa(run_info; f_over_vpa2=true, input=f_input, is=1, iz=1, fig=fig, ax=ax, frame_index=frame_index, label=label) + vlines!(ax, @lift cutoff_lower[input.ir0,$frame_index]; linestyle=:dash, color=:red) if plot_prefix !== nothing && fig !== nothing outfile=plot_prefix * "pdf_unnorm_over_vpa2_wall-_vs_vpa." * input.animation_ext save_animation(fig, frame_index, run_info.nt, outfile) @@ -5877,6 +5880,7 @@ function Chodura_condition_plots(run_info; plot_prefix=nothing, axes=nothing) animate_f_unnorm_vs_vpa(run_info; f_over_vpa2=true, input=f_input, is=1, iz=run_info.z.n, fig=fig, ax=ax, frame_index=frame_index, label=label) + vlines!(ax, @lift cutoff_upper[input.ir0,$frame_index]; linestyle=:dash, color=:red) if plot_prefix !== nothing && fig !== nothing outfile=plot_prefix * "pdf_unnorm_over_vpa2_wall+_vs_vpa." * input.animation_ext save_animation(fig, frame_index, run_info.nt, outfile) From 9d40d8d98fe770dbc2c9f357b93131f6699a8ec7 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 20 Sep 2024 16:36:18 +0100 Subject: [PATCH 14/14] Implement show_element_boundaries in plot/animate_f_unnorm_vs_vpa() --- .../src/makie_post_processing.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl index c2aff3809..42991d9fb 100644 --- a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl +++ b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl @@ -3942,6 +3942,13 @@ function plot_f_unnorm_vs_vpa(run_info; f_over_vpa2=false, input=nothing, neutra l = plot_1d(dzdt, f_unnorm; ax=ax, label=run_info.run_name, kwargs...) + if input.show_element_boundaries && fig !== nothing + element_boundary_inds = + [i for i ∈ 1:run_info.vpa.ngrid-1:run_info.vpa.n_global] + element_boundary_positions = dzdt[element_boundary_inds] + vlines!(ax, element_boundary_positions, color=:black, alpha=0.3) + end + if outfile !== nothing if fig === nothing error("When ax is passed, fig must also be passed to save the plot using " @@ -4318,6 +4325,14 @@ function animate_f_unnorm_vs_vpa(run_info; f_over_vpa2=false, input=nothing, l = plot_1d(dzdt, f_unnorm; ax=ax, label=run_info.run_name, yscale=yscale, kwargs...) + if input.show_element_boundaries && fig !== nothing + element_boundary_inds = + [i for i ∈ 1:run_info.vpa.ngrid-1:run_info.vpa.n_global] + element_boundary_positions = @lift $dzdt[element_boundary_inds] + vlines!(ax, element_boundary_positions, color=:black, alpha=0.3) + end + + if outfile !== nothing if fig === nothing error("When ax is passed, fig must also be passed to save the plot using "