Wall boundary condition experiment

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!(
@@ -3939,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 "
@@ -4161,8 +4171,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 +4274,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
@@ -4315,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 "
@@ -4656,11 +4674,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 +4772,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 +5559,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 +5614,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 +5646,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,
@@ -5618,18 +5704,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
@@ -5721,6 +5809,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