Add plotting for collisionality (comparing gradient scale lengths to …

…mean free path) and braginskii heat flux boundary condition at walls, following Stangeby (though I haven't checked the subtracted part for conductive heat flux yet).

makie_post_processing/makie_post_processing/src/makie_post_processing.jl

@@ -346,6 +346,8 @@ function makie_post_process(run_dir::Union{String,Tuple},
 
     sound_wave_plots(run_info; plot_prefix=plot_prefix)
 
+    collisionality_plots(run_info, plot_prefix)
+
     if all(ri === nothing for ri ∈ run_info_dfns)
         nvperp = nothing
     else
@@ -807,6 +809,24 @@ function _setup_single_input!(this_input_dict::OrderedDict{String,Any},
         plot_steady_state_residual=false,
         animate_steady_state_residual=false,
        )
+
+    set_defaults_and_check_section!(
+        this_input_dict, "collisionality_plots";
+        plot=true,
+        plot_dT_dz_vs_z=false,
+        animate_dT_dz_vs_z=false,
+        plot_mfp_vs_z=false,
+        animate_mfp_vs_z=false,
+        plot_nu_ii_vth_mfp_vs_z = false,
+        plot_LT_mfp_vs_z = false,
+        animate_LT_mfp_vs_z = false,
+        plot_LT_dT_dz_temp_vs_z = false,
+        plot_Ln_mfp_vs_z = false,
+        animate_Ln_mfp_vs_z = false,
+        plot_Lupar_mfp_vs_z = false,
+        animate_Lupar_mfp_vs_z = false,
+        animation_ext = "gif"
+       )
 
     # We allow top-level options in the post-processing input file
     check_sections!(this_input_dict; check_no_top_level_options=false)
@@ -8130,6 +8150,202 @@ function timestep_diagnostics(run_info, run_info_dfns; plot_prefix=nothing, it=n
     return steps_fig, dt_fig, CFL_fig
 end
 
+
+
+function collisionality_plots(run_info, plot_prefix=nothing)
+    # plot the mean free path of the ions at every point in z, and on the same graph
+    # plot the lengthscales of the gradients of all moments at each point in z
+    #println("run_info is ", run_info)
+    println("Making plots for collisionality")
+    if !isa(run_info, Tuple)
+        run_info = (run_info,)
+    end
+
+    input = Dict_to_NamedTuple(input_dict["collisionality_plots"])
+    temperature_input = Dict_to_NamedTuple(input_dict["temperature"])
+
+    mfp = get_variable(run_info, "mfp")
+    L_T = get_variable(run_info, "L_T")
+    L_n = get_variable(run_info, "L_n")
+    L_upar = get_variable(run_info, "L_upar")
+
+    if input.plot_mfp_vs_z
+        variable_name = "mean_free_path"
+        variable = mfp
+        variable_prefix = plot_prefix * variable_name
+        plot_vs_z(run_info, variable_name, is=1, data=variable, input=temperature_input,
+                  outfile=variable_prefix * "_vs_z.pdf")
+    end
+
+    if input.animate_mfp_vs_z
+        variable_name = "mean_free_path"
+        variable = mfp
+        variable_prefix = plot_prefix * variable_name
+        animate_vs_z(run_info, variable_name, is=1, data=variable, input=temperature_input,
+                        outfile=variable_prefix * "_vs_z." * input.animation_ext)
+    end
+    # TASKS: plot mean free path vs z, plot collision frequency vs z, make sure it agrees
+    # with original plot for it, and plot thermal speed and make sure that agrees, and then
+    # make sure that the mean free path then makes sense.
+    # Then plot the temperature and temperature gradient on the same plot
+    # Eventually plot the upar and dupar_dz and density and ddens_dz on plots too. Again,
+    # make sure they agree with the original plots.
+    if input.plot_dT_dz_vs_z
+        variable_name = "dT_dz"
+        variable = nothing
+        try
+            variable = get_variable(run_info, variable_name)
+        catch e
+            return makie_post_processing_error_handler(
+                    e,
+                    "collisionality_plots () failed for $variable_name - could not load data.")
+        end
+        #println("variable is ", variable)
+        variable_prefix = plot_prefix * variable_name
+        plot_vs_z(run_info, variable_name, is=1, data=variable, input=temperature_input,
+                  outfile=variable_prefix * "_vs_z.pdf")
+    end
+
+    if input.animate_dT_dz_vs_z
+        variable_name = "dT_dz"
+        variable = nothing
+        try
+            variable = get_variable(run_info, variable_name)
+        catch e
+            return makie_post_processing_error_handler(
+                    e,
+                    "collisionality_plots () failed for $variable_name - could not load data.")
+        end
+        #println("variable is ", variable)
+        variable_prefix = plot_prefix * variable_name
+        animate_vs_z(run_info, variable_name, is=1, data=variable, input=temperature_input,
+                        outfile=variable_prefix * "_vs_z." * input.animation_ext)
+    end
+
+    if input.plot_nu_ii_vth_mfp_vs_z
+        vth = get_variable(run_info, "thermal_speed")
+        nu_ii = get_variable(run_info, "collision_frequency_ii")
+        variable_prefix = plot_prefix * "checking_mfp_vth_and_nu_ii"
+        #println("vth[1] is ", vth[1])
+        fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+        plot_1d(run_info[1].z.grid, vth[1][:,1,1,21], xlabel="z",
+                ylabel="values", label="vth", ax=ax[1], title = "checking_mfp_vth")
+        plot_1d(run_info[1].z.grid, nu_ii[1][:,1,1,21], label="nu_ii", ax=ax[1])
+        plot_1d(run_info[1].z.grid, mfp[1][:,1,1,21], label="mfp", ax=ax[1])
+        Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+               orientation=:horizontal)
+        #plot_1d(z.grid, select_slice(error, :z; input=input), xlabel=L"z",
+        #        ylabel=norm_label, ax=ax[2])
+        outfile = variable_prefix * "_vs_z.pdf"
+        save(outfile, fig)
+    end
+
+    if input.plot_LT_dT_dz_temp_vs_z
+        dT_dz = get_variable(run_info, "dT_dz")
+        temp = get_variable(run_info, "temperature")
+        variable_prefix = plot_prefix * "LT_dT_dz_temp"
+        fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+        plot_1d(run_info[1].z.grid, L_T[1][:,1,1,21], xlabel="z",
+                ylabel="values", label="L_T", ax=ax[1], title = "LT_dT_dz_temp")
+        plot_1d(run_info[1].z.grid, dT_dz[1][:,1,1,21], label="dT_dz", ax=ax[1])
+        plot_1d(run_info[1].z.grid, temp[1][:,1,1,21], label="temp", ax=ax[1])
+        Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+               orientation=:horizontal)
+        outfile = variable_prefix * "_vs_z.pdf"
+        save(outfile, fig)
+    end
+
+    if input.plot_LT_mfp_vs_z
+        variable_prefix = plot_prefix * "LT_mfp_first_comparison"
+        fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+        plot_1d(run_info[1].z.grid, L_T[1][:,1,1,21], xlabel="z",
+                ylabel="values", label="L_T", ax=ax[1], title = "L_T and mean free path comparison")
+        plot_1d(run_info[1].z.grid, mfp[1][:,1,1,21], label="mfp", ax=ax[1])
+        Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+               orientation=:horizontal)
+        outfile = variable_prefix * "_vs_z.pdf"
+        save(outfile, fig)
+    end
+
+    if input.animate_LT_mfp_vs_z
+        nt = length(mfp[1][1,1,1,:])
+        variable_prefix = plot_prefix * "LT_mfp_first_comparison"
+
+        fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+        frame_index = Observable(1)
+        animate_1d(run_info[1].z.grid, L_T[1][:,1,1,:],
+                   frame_index=frame_index, xlabel="z", ylabel="values",
+                   label="L_T", ax=ax[1])
+        animate_1d(run_info[1].z.grid, mfp[1][:,1,1,:],
+                   frame_index=frame_index, xlabel="z", ylabel="values",
+                   label="mfp", ax=ax[1])
+        Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+               orientation=:horizontal)
+        outfile = variable_prefix * "_vs_z." * input.animation_ext
+        save_animation(fig, frame_index, nt, outfile)
+    end
+
+    if input.plot_Ln_mfp_vs_z
+        variable_prefix = plot_prefix * "Ln_mfp_first_comparison"
+        fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+        plot_1d(run_info[1].z.grid, L_n[1][:,1,1,21], xlabel="z",
+                ylabel="values", label="L_n", ax=ax[1], title = "Ln and mean free path comparison")
+        plot_1d(run_info[1].z.grid, mfp[1][:,1,1,21], label="mfp", ax=ax[1])
+        Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+               orientation=:horizontal)
+        outfile = variable_prefix * "_vs_z.pdf"
+        save(outfile, fig)
+    end
+
+    if input.animate_Ln_mfp_vs_z
+        nt = length(mfp[1][1,1,1,:])
+        variable_prefix = plot_prefix * "Ln_mfp_first_comparison"
+
+        fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+        frame_index = Observable(1)
+        animate_1d(run_info[1].z.grid, L_n[1][:,1,1,:],
+                   frame_index=frame_index, xlabel="z", ylabel="values",
+                   label="L_n", ax=ax[1])
+        animate_1d(run_info[1].z.grid, mfp[1][:,1,1,:],
+                   frame_index=frame_index, xlabel="z", ylabel="values",
+                   label="mfp", ax=ax[1])
+        Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+               orientation=:horizontal)
+        outfile = variable_prefix * "_vs_z." * input.animation_ext
+        save_animation(fig, frame_index, nt, outfile)
+    end
+
+    if input.plot_Lupar_mfp_vs_z
+        variable_prefix = plot_prefix * "Lupar_mfp_first_comparison"
+        fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+        plot_1d(run_info[1].z.grid, L_upar[1][:,1,1,21], xlabel="z",
+                ylabel="values", label="L_upar", ax=ax[1], title = "Lupar and mean free path comparison")
+        plot_1d(run_info[1].z.grid, mfp[1][:,1,1,21], label="mfp", ax=ax[1])
+        Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+               orientation=:horizontal)
+        outfile = variable_prefix * "_vs_z.pdf"
+        save(outfile, fig)
+    end
+
+    if input.animate_Lupar_mfp_vs_z
+        nt = length(mfp[1][1,1,1,:])
+        variable_prefix = plot_prefix * "Lupar_mfp_first_comparison"
+
+        fig, ax, legend_place = get_1d_ax(1; get_legend_place=:below)
+        frame_index = Observable(1)
+        animate_1d(run_info[1].z.grid, L_upar[1][:,1,1,:],
+                   frame_index=frame_index, xlabel="z", ylabel="values",
+                   label="L_upar", ax=ax[1])
+        animate_1d(run_info[1].z.grid, mfp[1][:,1,1,:],
+                   frame_index=frame_index, xlabel="z", ylabel="values",
+                   label="mfp", ax=ax[1])
+        Legend(legend_place[1], ax[1]; tellheight=true, tellwidth=false,
+               orientation=:horizontal)
+        outfile = variable_prefix * "_vs_z." * input.animation_ext
+        save_animation(fig, frame_index, nt, outfile)
+    end
+end
+
 # Utility functions
 ###################
 #

moment_kinetics/src/load_data.jl

@@ -4145,6 +4145,82 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t
     if variable_name == "temperature"
         vth = get_variable(run_info, "thermal_speed"; kwargs...)
         variable = vth.^2
+    elseif variable_name == "dT_dz"
+        T = get_variable(run_info, "temperature"; kwargs...)
+        variable = similar(T)
+        if :iz ∈ keys(kwargs) && kwargs[:iz] !== nothing
+            error("Cannot take z-derivative when iz!==nothing")
+        end
+        if :ir ∈ keys(kwargs) && isa(kwargs[:ir], mk_int)
+            for it ∈ 1:size(variable, 3)
+                @views derivative!(variable[:,:,it], T[:,:,it], run_info.z, run_info.z_spectral)
+            end
+        else
+            for it ∈ 1:size(variable, 4), ir ∈ 1:run_info.r.n
+                @views derivative!(variable[:,:,ir,it], T[:,:,ir,it], run_info.z, run_info.z_spectral)
+            end
+        end
+    elseif variable_name == "dn_dz"
+        n = get_variable(run_info, "density"; kwargs...)
+        variable = similar(n)
+        if :iz ∈ keys(kwargs) && kwargs[:iz] !== nothing
+            error("Cannot take z-derivative when iz!==nothing")
+        end
+        if :ir ∈ keys(kwargs) && isa(kwargs[:ir], mk_int)
+            for it ∈ 1:size(variable, 3)
+                @views derivative!(variable[:,:,it], n[:,:,it], run_info.z, run_info.z_spectral)
+            end
+        else
+            for it ∈ 1:size(variable, 4), ir ∈ 1:run_info.r.n
+                @views derivative!(variable[:,:,ir,it], n[:,:,ir,it], run_info.z, run_info.z_spectral)
+            end
+        end
+    elseif variable_name == "dupar_dz"
+        upar = get_variable(run_info, "parallel_flow"; kwargs...)
+        variable = similar(upar)
+        if :iz ∈ keys(kwargs) && kwargs[:iz] !== nothing
+            error("Cannot take z-derivative when iz!==nothing")
+        end
+        if :ir ∈ keys(kwargs) && isa(kwargs[:ir], mk_int)
+            for it ∈ 1:size(variable, 3)
+                @views derivative!(variable[:,:,it], upar[:,:,it], run_info.z, run_info.z_spectral)
+            end
+        else
+            for it ∈ 1:size(variable, 4), ir ∈ 1:run_info.r.n
+                @views derivative!(variable[:,:,ir,it], upar[:,:,ir,it], run_info.z, run_info.z_spectral)
+            end
+        end
+    elseif variable_name == "mfp"
+        vth = get_variable(run_info, "thermal_speed"; kwargs...)
+        nu_ii = get_variable(run_info, "collision_frequency_ii"; kwargs...)
+        variable = vth ./ nu_ii
+    elseif variable_name == "L_T"
+        dT_dz = get_variable(run_info, "dT_dz"; kwargs...)
+        temp = get_variable(run_info, "temperature"; kwargs...)
+        # We define gradient lengthscale of T as LT^-1 = dln(T)/dz (ignore negative sign
+        # tokamak convention as we're only concerned with comparing magnitudes)
+        variable = abs.(temp .* dT_dz.^(-1))
+        # flat points in temperature have diverging LT, so ignore those with NaN
+        # using a hard coded 10.0 tolerance for now
+        variable[variable .> 10.0] .= NaN
+    elseif variable_name == "L_n"
+        dn_dz = get_variable(run_info, "dn_dz"; kwargs...)
+        n = get_variable(run_info, "density"; kwargs...)
+        # We define gradient lengthscale of n as Ln^-1 = dln(n)/dz (ignore negative sign
+        # tokamak convention as we're only concerned with comparing magnitudes)
+        variable = abs.(n .* dn_dz.^(-1))
+        # flat points in temperature have diverging Ln, so ignore those with NaN
+        # using a hard coded 10.0 tolerance for now
+        variable[variable .> 10.0] .= NaN
+    elseif variable_name == "L_upar"
+        dupar_dz = get_variable(run_info, "dupar_dz"; kwargs...)
+        upar = get_variable(run_info, "parallel_flow"; kwargs...)
+        # We define gradient lengthscale of upar as Lupar^-1 = dln(upar)/dz (ignore negative sign
+        # tokamak convention as we're only concerned with comparing magnitudes)
+        variable = abs.(upar .* dupar_dz.^(-1))
+        # flat points in temperature have diverging Lupar, so ignore those with NaN
+        # using a hard coded 10.0 tolerance for now
+        variable[variable .> 10.0] .= NaN
     elseif variable_name == "collision_frequency_ii"
         n = get_variable(run_info, "density"; kwargs...)
         vth = get_variable(run_info, "thermal_speed"; kwargs...)

moment_kinetics/src/velocity_moments.jl

@@ -768,7 +768,7 @@ function update_ion_qpar_species!(qpar, density, upar, vth, dT_dz, ff, vpa, vper
         calculate_ion_qpar_from_pdf!(qpar, density, upar, vth, ff, vpa, vperp, z, r, evolve_density,
                                      evolve_upar, evolve_ppar)
     elseif ion_physics == braginskii_ions
-        calculate_ion_qpar_from_braginskii!(qpar, density, vth, dT_dz, z, r, collisions, evolve_density, 
+        calculate_ion_qpar_from_braginskii!(qpar, density, upar, vth, dT_dz, z, r, collisions, evolve_density, 
                                             evolve_upar, evolve_ppar)
     else
         throw(ArgumentError("ion model $ion_physics not implemented for qpar calculation"))
@@ -820,7 +820,7 @@ end
 """
 calculate parallel heat flux if ion composition flag is Braginskii fluid ions
 """
-function calculate_ion_qpar_from_braginskii!(qpar, density, vth, dT_dz, z, r, collisions, evolve_density, evolve_upar, evolve_ppar)
+function calculate_ion_qpar_from_braginskii!(qpar, density, upar, vth, dT_dz, z, r, collisions, evolve_density, evolve_upar, evolve_ppar)
     # Note that this is a braginskii heat flux for ions using the krook operator. The full Fokker-Planck operator
     # Braginskii heat flux is different! This also assumes one ion species, and so no friction between ions.
     @boundscheck r.n == size(qpar, 2) || throw(BoundsError(qpar))
@@ -834,7 +834,50 @@ function calculate_ion_qpar_from_braginskii!(qpar, density, vth, dT_dz, z, r, co
             qpar[iz,ir] = -(1/2) * 5/4 * density[iz,ir] * vth[iz,ir]^2 /nu_ii * dT_dz[iz,ir,1]
         end
     else
-        throw(ArgumentError("Braginskii heat flux simulation requires evolve_density, evolve_upar and evolve_ppar to be true, since it is a purely fluid simulation"))
+        throw(ArgumentError("Braginskii heat flux simulation requires evolve_density, 
+              evolve_upar and evolve_ppar to be true, since it is a purely fluid simulation"))
+    end
+
+    # add boundary condition to the heat flux, since now there is no distribution function 
+    # (in this case shape function) whose cutoff boundary condition can hold the parallel heat
+    # flux in check. See Stangeby textbook, equations (2.92) and (2.93), and the paragraph between.
+
+    if z.bc == "periodic"
+        # There's no wall boundary condition here, do nothing (qpar can be what it wants)
+        return nothing
+    end
+
+    begin_r_region()
+
+    if z.irank == 0 && (z.irank == z.nrank - 1)
+        z_indices = (1, z.n)
+    elseif z.irank == 0
+        z_indices = (1,)
+    elseif z.irank == z.nrank - 1
+        z_indices = (z.n,)
+    else
+        return nothing
+    end
+
+    @loop_r ir begin
+        for iz ∈ z_indices
+            this_ppar = vth[iz,ir]^2 * density[iz,ir]/2.0
+            this_upar = upar[iz,ir]
+            this_dens = density[iz,ir]
+            particle_flux = this_dens * this_upar
+            T_i = vth[iz,ir]^2
+
+            # Stangeby paragraph after (2.92)
+            gamma_i = 2.0
+
+            # Stangeby (2.92)
+            total_heat_flux = gamma_i * T_i * particle_flux
+
+            # E.g. Helander&Sigmar (2.14) ??????? what is it for ions?
+            conductive_heat_flux = total_heat_flux - 2.5 * this_ppar * this_upar
+
+            qpar[iz,ir] = conductive_heat_flux
+        end
     end
     return nothing
 end