Skip to content

Commit

Permalink
Merge pull request #181 from mabarnes/wall-boundary-condition-experiment
Browse files Browse the repository at this point in the history 
Wall boundary condition experiment
  • Loading branch information
@johnomotani
johnomotani authored Sep 21, 2024
2 parents bbc609c + 9d40d8d commit 5c1767e
Show file tree
Hide file tree
Showing 8 changed files with 377 additions and 74 deletions.
241 changes: 211 additions & 30 deletions makie_post_processing/makie_post_processing/src/makie_post_processing.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -756,8 +756,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!(
Expand Down Expand Up @@ -3937,6 +3940,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 "
Expand Down Expand Up @@ -4159,8 +4169,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

Expand Down Expand Up @@ -4262,16 +4272,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
Expand Down Expand Up @@ -4313,6 +4323,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 "
Expand Down Expand Up @@ -4654,11 +4672,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")
Expand Down Expand Up @@ -4757,11 +4770,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)
Expand Down Expand Up @@ -5549,6 +5557,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)
Expand All @@ -5557,26 +5612,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
Expand All @@ -5589,6 +5644,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,
Expand Down Expand Up @@ -5616,18 +5702,20 @@ 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 =
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, 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,
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
Expand Down Expand Up @@ -5719,6 +5807,99 @@ 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)
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)
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)
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)
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)
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)
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)
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)
end
end

return nothing
end

Expand Down
Loading

0 comments on commit 5c1767e

Please sign in to comment.