How to make this distribution symmetric?

New to Julia
1

Let’s say I have the following data file:

https://ufile.io/b2v6t45e

Using this file, I can plot a Maxwell-Boltzmann functon and see how it changes as a function of s with this code:

# using necessary packages
using LaTeXStrings, Markdown, DataFrames, Queryverse, Plots, Printf
 
df = load("data.csv", spacedelim=true, header_exists=true) |> DataFrame

m_e = 9.11E-28 # electron mass [g]
k = 1.38E-16 # boltzmann constant [erg*K^-1]

function Maxwellian(v_Max, T_Max)
    ```Maxwellian distribution function for particles moving in only one direction```
    normal = (m_e/(2*pi*k*T_Max))^1.5
    exp_term = exp(-((m_e).*v_Max.*v_Max)/(3*k*T_Max))
    return normal*exp_term
end

nu = 1 # collisional frequency [s^-1]
v_pos = range(0, 1e10, step=1e8) # positive velocity [cm/s]

f_s_pos = [Maxwellian.(v_pos, df.T[1])] # because boundary density is high at left footpoint, assume f_N = f_M

for is in 1:length(df.s) # for positive velocities, integrate from 0 to s_max
    local_F_s_pos = Maxwellian.(v_pos, df.T[is])
    push!(f_s_pos, f_s_pos[is] - (nu./v_pos).*(f_s_pos[is] - local_F_s_pos)*(df.ds[is]))
end

anim = @animate for i in 1:length(df.s)
           plot(v_pos, f_s_pos[i]; label=@sprintf("f_s at s = %.3e cm", df.s[i]))
           xlabel!("velocity (cm/s)")
           ylabel!("probability density (s/cm)")
           title!("Space-Dependent Maxwellian")
           ylims!(-1e-28, 2e-27)
       end;

gif(anim, "Simple_BGK_Spatial_Discretization_Pos.gif"; fps=20)

Simple_BGK_Spatial_Discretization_Pos

Now, I want to try to create this animation again, but instead, we have negative values for velocity (v_neg instead of v_pos) and we have to step through s values from the max value of s to the min value of s. Here is my attempt:

nu = 1 # collisional frequency [s^-1]
v_neg = range(0, -1e10, step=-1e8) # negative velocity [cm/s]

f_s_neg = [Maxwellian.(v_neg, df.T[end])] # because boundary density is high at right footpoint, assume f_N = f_M

for is in reverse(eachindex(df.s)) # for negative velocities, integrate from s_max to 0 
    local_F_s_neg = Maxwellian.(v_neg, df.T[is])
    push!(f_s_neg, f_s_neg[length(df.s)+1-is] - (nu./v_neg).*(f_s_neg[length(df.s)+1-is] - local_F_s_neg)*(df.ds[is]))
end

anim = @animate for i in 1:length(df.s)
           plot(v_neg, f_s_neg[i]; label=@sprintf("f_s_neg at s = %.3e cm", df.s[i]))
           xlabel!("velocity (cm/s)")
           ylabel!("probability density (s/cm)")
           title!("Space-Dependent Maxwellian")
           ylims!(-1e-30, 2e-27)
       end;

gif(anim, "Simple_BGK_Spatial_Discretization_Neg.gif"; fps=20)

Simple_BGK_Spatial_Discretization_Neg

I was expecting the function to look symmetric to the previous plot along the y axis but this is not the case. This might be a problem specific in my field (physics), but before investigating further, from a Julia programming perspective, are there any glaring issues with the way I set up the code for the negative velocity case that may be the problem?

2

Since you are going backwards, you probably need to change the - sign to +:

... + (nu./v_neg).*(f_s_neg[length(df.s)+1-is] - local_F_s_neg)*(df.ds[is]))
1 Like
3

@rafael.guerra This worked! Thank you so much!

Do you know how I could rewrite my code to combine the positive anfd negative velocity curves onto one graph/animation?

4

The easiest might be to concatenate the vectors to plot, say:

plot([v_neg; v_pos], [f_s_neg[i]; f_s_pos[i]]; ...)
1 Like
5

Hello again. So now, I am trying to recreate the plot but now using an implicit numerical method rather than an explicit numerical method.

# Spatial discretization 

v_pos = range(0, 1e10, step=1e8) # positive velocity [cm/s]
v_neg = range(0, -1e10, step=-1e8) # negative velocity [cm/s]

nu = 1 # constant collisional frequency [s^-1]

f_s_pos_imp = [Maxwellian.(v_pos, df.T[1])] # because boundary density is high at left footpoint, assume f_N = f_M
f_s_neg_imp = [Maxwellian.(v_neg, df.T[end])] # because boundary density is high at right footpoint, assume f_N = f_M

for is in 1:length(df.s) # for positive velocities, integrate from 0 to s_max
    local_f_s_pos = Maxwellian.(v_pos, df.T[is])
    #push!(f_s_pos_exp, f_s_pos_exp[is] - (nu ./v_pos).*(f_s_pos_exp[is] - local_f_s_pos_exp)*(df.ds[is]))   
    push!(f_s_pos_imp, (f_s_pos_imp[is] .+ (nu./v_pos).*local_f_s_pos.*df.ds[is])/(1 .+ (nu./v_pos).*(df.ds[is])))
    
end

for is in reverse(eachindex(df.s)) # for negative velocities, integrate from s_max to 0 
    local_f_s_neg = Maxwellian.(v_neg, df.T[is])
    #push!(f_s_neg_exp, f_s_neg_exp[length(df.s)+1-is] + (nu ./v_neg).*(f_s_neg_exp[length(df.s)+1-is] - local_f_s_neg_exp)
     #   *(df.ds[is]))
    push!(f_s_neg_imp, (f_s_neg_imp[is] .+ (nu./v_neg).*local_f_s_neg.*df.ds[is])/(1 .+ (nu./v_neg).*(df.ds[is])))
end

anim = @animate for i in 1:length(df.s)
           plot([v_neg; v_pos], [f_s_neg_imp[i]; f_s_pos_imp[i]]; label=@sprintf("f_s at s = %.3e cm", df.s[i]))
           xlabel!("velocity (cm/s)")
           ylabel!("probability density (s/cm)")
           title!("Space-Dependent Maxwellian")
           ylims!(-1e-28, 2e-27)
       end;

gif(anim, "TEST.gif"; fps=30)

I thought I had ensured that all my lines in the loops were performing vectorized operations by adding ., but I receive this error:

MethodError: no method matching Vector{Float64}(::Matrix{Float64})
Closest candidates are:
  Array{T, N}(::AbstractArray{S, N}) where {T, N, S} at C:\Users\Acer\AppData\Local\Programs\Julia-1.7.0\share\julia\base\array.jl:563
  Vector{T}() where T at C:\Users\Acer\AppData\Local\Programs\Julia-1.7.0\share\julia\base\boot.jl:476
  Vector{T}(::SuiteSparse.CHOLMOD.Dense{T}) where T at C:\Users\Acer\AppData\Local\Programs\Julia-1.7.0\share\julia\stdlib\v1.7\SuiteSparse\src\cholmod.jl:850

What may I be missing here?

6

This error is telling you that you can’t make a Vector{Float64} by passing the constructor a Matrix{Float64}

since you didn’t get the full error, I’m not sure what line is causing the issue.

7

You seem to have forgotten one dot for the division.
Instead of sprinkling dots all over the place and whenever possible, just use the @. macro at the beginning of the expression and leave the rest alone.

NB:

  • In nu./v_pos there is division by 0
  • As the result is symmetrical, you just need to model the positive velocities, and then mirror the resulting arrays for display
1 Like
8

Ah, thank you very much. The @. macro makes things convenient.