Skip to content

Commit

Permalink
Merge pull request #133 from mabarnes/chodura-moment-kinetic
Browse files Browse the repository at this point in the history 
Make Chodura condition check compatible with moment-kinetic modes
  • Loading branch information
@johnomotani
johnomotani authored Nov 14, 2023
2 parents c35bec0 + 2572815 commit fd2e48b
Show file tree
Hide file tree
Showing 3 changed files with 106 additions and 83 deletions.
86 changes: 81 additions & 5 deletions src/analysis.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ using ..array_allocation: allocate_float
using ..calculus: integral
using ..communication
using ..coordinates: coordinate
using ..initial_conditions: vpagrid_to_dzdt
using ..interpolation: interpolate_to_grid_1d
using ..load_data: open_readonly_output_file, get_nranks, load_pdf_data, load_rank_data
using ..load_data: load_distributed_charged_pdf_slice
Expand Down Expand Up @@ -97,11 +98,13 @@ 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, composition, Er, geometry, z_bc,
nblocks, run_name=nothing,
function check_Chodura_condition(r, z, vperp, vpa, dens, upar, vth, 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)
f_lower=nothing, f_upper=nothing,
evolve_density=false, evolve_upar=false,
evolve_ppar=false)

if z_bc != "wall"
return nothing, nothing
Expand Down Expand Up @@ -153,8 +156,15 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, composition, Er, geomet
(size(f_upper, 1), size(f_upper, 2), 1, size(f_upper, 3),
size(f_upper, 4), size(f_upper, 5)))
end

f_lower = @views get_unnormalised_f_1d(f_lower, dens[1,:,:,:], vth[1,:,:,:],
evolve_density, evolve_ppar)
f_upper = @views get_unnormalised_f_1d(f_upper, dens[end,:,:,:], vth[end,:,:,:],
evolve_density, evolve_ppar)
for it ∈ 1:ntime, ir ∈ 1:nr
vpabar = @. vpa.grid - 0.5 * geometry.rhostar * Er[1,ir,it] / geometry.bzed
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

# Get rid of a zero if it is there to avoid a blow up - f should be zero at that
# point anyway
Expand All @@ -174,7 +184,9 @@ function check_Chodura_condition(r, z, vperp, vpa, dens, composition, Er, geomet

lower_result[ir,it] *= 0.5 * composition.T_e / dens[1,ir,is,it]

vpabar = @. vpa.grid - 0.5 * geometry.rhostar * Er[end,ir,it] / geometry.bzed
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

# Get rid of a zero if it is there to avoid a blow up - f should be zero at that
# point anyway
Expand Down Expand Up @@ -754,4 +766,68 @@ function _steady_state_absolute_residual(variable, variable_at_previous_time, re
return absolute_residual
end

"""
Get the unnormalised distribution function and unnormalised ('lab space') dzdt
coordinate at a point in space.
Inputs should depend only on vpa.
"""
function get_unnormalised_f_dzdt_1d(f, vpa_grid, density, upar, vth, evolve_density,
evolve_upar, evolve_ppar)

dzdt = vpagrid_to_dzdt(vpa_grid, vth, upar, evolve_ppar, evolve_upar)

f_unnorm = get_unnormalised_f_1d(f, density, vth, evolve_density, evolve_ppar)

return f_unnorm, dzdt
end
function get_unnormalised_f_1d(f, density, vth, evolve_density, evolve_ppar)
if evolve_ppar
f_unnorm = @. f * density / vth
elseif evolve_density
f_unnorm = @. f * density
else
f_unnorm = f
end
return f_unnorm
end

"""
Get the unnormalised distribution function and unnormalised ('lab space') coordinates.
Inputs should depend only on z and vpa.
"""
function get_unnormalised_f_coords_2d(f, z_grid, vpa_grid, density, upar, vth,
evolve_density, evolve_upar, evolve_ppar)

nvpa, nz = size(f)
z2d = zeros(nvpa, nz)
for iz ∈ 1:nz
z2d[:,iz] .= z_grid[iz]
end
dzdt2d = vpagrid_to_dzdt_2d(vpa_grid, vth, upar, evolve_ppar, evolve_upar)
f_unnorm = get_unnormalised_f_2d(f, density, vth, evolve_density, evolve_ppar)

return f_unnorm, z2d, dzdt2d
end
function vpagrid_to_dzdt_2d(vpa_grid, vth, upar, evolve_ppar, evolve_upar)
nvpa = length(vpa_grid)
nz = length(vth)
dzdt2d = zeros(nvpa, nz)
for iz ∈ 1:nz
@views dzdt2d[:,iz] .= vpagrid_to_dzdt(vpa_grid, vth[iz], upar[iz], evolve_ppar,
evolve_upar)
end
return dzdt2d
end
function get_unnormalised_f_2d(f, density, vth, evolve_density, evolve_ppar)
f_unnorm = similar(f)
nz = size(f, 2)
for iz ∈ 1:nz
@views f_unnorm[:,iz] .= get_unnormalised_f_1d(f[:,iz], density[iz], vth[iz],
evolve_density, evolve_ppar)
end
return f_unnorm
end

end
18 changes: 11 additions & 7 deletions src/makie_post_processing.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,9 @@ export makie_post_process, generate_example_input_file, get_variable,
irregular_heatmap!, postproc_load_variable, positive_or_nan

using ..analysis: analyze_fields_data, check_Chodura_condition, get_r_perturbation,
get_Fourier_modes_2D, get_Fourier_modes_1D, steady_state_residuals
get_Fourier_modes_2D, get_Fourier_modes_1D, steady_state_residuals,
get_unnormalised_f_dzdt_1d, get_unnormalised_f_coords_2d,
get_unnormalised_f_1d, vpagrid_to_dzdt_2d, get_unnormalised_f_2d
using ..array_allocation: allocate_float
using ..coordinates: define_coordinate
using ..input_structs: grid_input, advection_input, set_defaults_and_check_top_level!,
Expand All @@ -30,10 +32,7 @@ using ..load_data: open_readonly_output_file, get_group, load_block_data,
load_species_data, load_time_data
using ..initial_conditions: vpagrid_to_dzdt
using ..post_processing: calculate_and_write_frequencies, construct_global_zr_coords,
get_geometry_and_composition, get_unnormalised_f_dzdt_1d,
get_unnormalised_f_coords_2d, get_unnormalised_f_1d,
vpagrid_to_dzdt_2d, get_unnormalised_f_2d,
read_distributed_zr_data!
get_geometry_and_composition, read_distributed_zr_data!
using ..type_definitions: mk_float, mk_int

using Combinatorics
Expand Down Expand Up @@ -4799,14 +4798,19 @@ function Chodura_condition_plots(run_info; plot_prefix=nothing, axes=nothing)

time = run_info.time
density = postproc_load_variable(run_info, "density")
upar = postproc_load_variable(run_info, "parallel_flow")
vth = postproc_load_variable(run_info, "thermal_speed")
Er = postproc_load_variable(run_info, "Er")
f_lower = postproc_load_variable(run_info, "f", iz=1)
f_upper = postproc_load_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, run_info.composition, Er,
run_info.geometry, run_info.z.bc, nothing;
run_info.vpa, density, upar, vth, 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,
f_lower=f_lower, f_upper=f_upper)

if input.plot_vs_t
Expand Down
85 changes: 14 additions & 71 deletions src/post_processing.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,6 @@ using ..array_allocation: allocate_float
using ..coordinates: define_coordinate
using ..file_io: open_ascii_output_file
using ..type_definitions: mk_float, mk_int
using ..initial_conditions: vpagrid_to_dzdt
using ..load_data: open_readonly_output_file, get_group, load_input, load_time_data
using ..load_data: get_nranks
using ..load_data: load_fields_data, load_pdf_data
Expand All @@ -47,7 +46,8 @@ using ..load_data: load_variable
using ..load_data: load_coordinate_data, load_block_data, load_rank_data,
load_species_data, load_mk_options
using ..analysis: analyze_fields_data, analyze_moments_data, analyze_pdf_data,
check_Chodura_condition, analyze_2D_instability
check_Chodura_condition, analyze_2D_instability,
get_unnormalised_f_dzdt_1d, get_unnormalised_f_coords_2d
using ..velocity_moments: integrate_over_vspace
using ..manufactured_solns: manufactured_solutions, manufactured_electric_fields
using ..moment_kinetics_input: mk_input, get
Expand Down Expand Up @@ -853,8 +853,12 @@ function analyze_and_plot_data(prefix...; run_index=nothing)
if pp.diagnostics_chodura_t
Chodura_ratio_lower, Chodura_ratio_upper =
get_tuple_of_return_values(check_Chodura_condition, r, z, vperp, vpa,
density_at_pdf_times, composition, Er_at_pdf_times,
geometry, "wall", nblocks, run_names, nothing, ir0)
density_at_pdf_times, parallel_flow_at_pdf_times,
thermal_speed_at_pdf_times, composition,
Er_at_pdf_times, geometry, "wall", nblocks,
run_names, nothing, ir0,
evolve_density=evolve_density,
evolve_upar=evolve_upar, evolve_ppar=evolve_ppar)

plot(legend=legend)
for (t, cr, run_label) ∈ zip(time_pdfs, Chodura_ratio_lower, run_names)
Expand All @@ -875,9 +879,12 @@ function analyze_and_plot_data(prefix...; run_index=nothing)
if pp.diagnostics_chodura_r
Chodura_ratio_lower, Chodura_ratio_upper =
get_tuple_of_return_values(check_Chodura_condition, r, z, vperp, vpa,
density_at_pdf_times, composition, Er_at_pdf_times,
geometry, "wall", nblocks, run_names, ntime_pdfs,
nothing)
density_at_pdf_times, parallel_flow_at_pdf_times,
thermal_speed_at_pdf_times, composition,
Er_at_pdf_times, geometry, "wall", nblocks,
run_names, ntime_pdfs, nothing,
evolve_density=evolve_density,
evolve_upar=evolve_upar, evolve_ppar=evolve_ppar)

plot(legend=legend)
for (this_r, cr, run_label) ∈ zip(r, Chodura_ratio_lower, run_names)
Expand Down Expand Up @@ -2842,70 +2849,6 @@ function draw_v_parallel_zero!(z::AbstractVector, upar, vth, evolve_upar::Bool,
draw_v_parallel_zero!(Plots.CURRENT_PLOT, z, upar, vth, evolve_upar, evolve_ppar)
end

"""
Get the unnormalised distribution function and unnormalised ('lab space') dzdt
coordinate at a point in space.
Inputs should depend only on vpa.
"""
function get_unnormalised_f_dzdt_1d(f, vpa_grid, density, upar, vth, evolve_density,
evolve_upar, evolve_ppar)

dzdt = vpagrid_to_dzdt(vpa_grid, vth, upar, evolve_ppar, evolve_upar)

f_unnorm = get_unnormalised_f_1d(f, density, vth, evolve_density, evolve_ppar)

return f_unnorm, dzdt
end
function get_unnormalised_f_1d(f, density, vth, evolve_density, evolve_ppar)
if evolve_ppar
f_unnorm = @. f * density / vth
elseif evolve_density
f_unnorm = @. f * density
else
f_unnorm = f
end
return f_unnorm
end

"""
Get the unnormalised distribution function and unnormalised ('lab space') coordinates.
Inputs should depend only on z and vpa.
"""
function get_unnormalised_f_coords_2d(f, z_grid, vpa_grid, density, upar, vth,
evolve_density, evolve_upar, evolve_ppar)

nvpa, nz = size(f)
z2d = zeros(nvpa, nz)
for iz ∈ 1:nz
z2d[:,iz] .= z_grid[iz]
end
dzdt2d = vpagrid_to_dzdt_2d(vpa_grid, vth, upar, evolve_ppar, evolve_upar)
f_unnorm = get_unnormalised_f_2d(f, density, vth, evolve_density, evolve_ppar)

return f_unnorm, z2d, dzdt2d
end
function vpagrid_to_dzdt_2d(vpa_grid, vth, upar, evolve_ppar, evolve_upar)
nvpa = length(vpa_grid)
nz = length(vth)
dzdt2d = zeros(nvpa, nz)
for iz ∈ 1:nz
@views dzdt2d[:,iz] .= vpagrid_to_dzdt(vpa_grid, vth[iz], upar[iz], evolve_ppar,
evolve_upar)
end
return dzdt2d
end
function get_unnormalised_f_2d(f, density, vth, evolve_density, evolve_ppar)
f_unnorm = similar(f)
nz = size(f, 2)
for iz ∈ 1:nz
@views f_unnorm[:,iz] .= get_unnormalised_f_1d(f[:,iz], density[iz], vth[iz],
evolve_density, evolve_ppar)
end
return f_unnorm
end

"""
Make a 2d plot of an unnormalised f on unnormalised coordinates, as returned from
get_unnormalised_f_coords()
Expand Down

0 comments on commit fd2e48b

Please sign in to comment.