Skip to content

Commit

Permalink
Fix calculation of sigma_fraction in kinetic electron bc
Browse files Browse the repository at this point in the history
  • Loading branch information
johnomotani committed Nov 19, 2024
1 parent 769fecc commit 2a2a3c7
Showing 1 changed file with 15 additions and 16 deletions.
31 changes: 15 additions & 16 deletions moment_kinetics/src/electron_kinetic_equation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2138,11 +2138,14 @@ function apply_electron_bc_and_constraints_no_r!(f_electron, phi, moments, z, vp
end

function get_cutoff_params_lower(upar, vthe, phi, me_over_mi, vpa, ir)
# Delete the upar contribution here if ignoring the 'upar shift'
vpa_unnorm = @. vpa.scratch2 = vthe[1,ir] * vpa.grid + upar[1,ir]

u_over_vt = upar[1,ir] / vthe[1,ir]

# sigma is the location we use for w_∥(v_∥=0) - set to 0 to ignore the 'upar
# shift'
sigma = -u_over_vt

vpa_unnorm = @. vpa.scratch2 = vthe[1,ir] * (vpa.grid - sigma)

# Initial guess for cut-off velocity is result from previous RK stage (which
# might be the previous timestep if this is the first stage). Recalculate this
# value from phi.
Expand All @@ -2157,10 +2160,6 @@ function get_cutoff_params_lower(upar, vthe, phi, me_over_mi, vpa, ir)
error("In lower-z electron bc, failed to find vpa=-vcut point, minus_vcut_ind=$minus_vcut_ind")
end

# sigma is the location we use for w_∥(v_∥=0) - set to 0 to ignore the 'upar
# shift'
sigma = -u_over_vt

# sigma is between sigma_ind-1 and sigma_ind
sigma_ind = searchsortedfirst(vpa_unnorm, 0.0)
if sigma_ind < 2
Expand All @@ -2172,7 +2171,7 @@ function get_cutoff_params_lower(upar, vthe, phi, me_over_mi, vpa, ir)

# sigma_fraction is the fraction of the distance between sigma_ind-1 and
# sigma_ind where sigma is.
sigma_fraction = (sigma - vpa_unnorm[sigma_ind-1]) / (vpa_unnorm[sigma_ind] - vpa_unnorm[sigma_ind-1])
sigma_fraction = -vpa_unnorm[sigma_ind-1] / (vpa_unnorm[sigma_ind] - vpa_unnorm[sigma_ind-1])

# Want the element that contains the interval on the lower side of sigma_ind. For
# points on element boundaries, the `ielement` array contains the element on the lower
Expand Down Expand Up @@ -2207,11 +2206,15 @@ function get_cutoff_params_lower(upar, vthe, phi, me_over_mi, vpa, ir)
end

function get_cutoff_params_upper(upar, vthe, phi, me_over_mi, vpa, ir)
# Delete the upar contribution here if ignoring the 'upar shift'
vpa_unnorm = @. vpa.scratch2 = vthe[end,ir] * vpa.grid + upar[end,ir]

u_over_vt = upar[end,ir] / vthe[end,ir]

# sigma is the location we use for w_∥(v_∥=0) - set to 0 to ignore the 'upar
# shift'
sigma = -u_over_vt

# Delete the upar contribution here if ignoring the 'upar shift'
vpa_unnorm = @. vpa.scratch2 = vthe[end,ir] * (vpa.grid - sigma)

# Initial guess for cut-off velocity is result from previous RK stage (which
# might be the previous timestep if this is the first stage). Recalculate this
# value from phi.
Expand All @@ -2226,10 +2229,6 @@ function get_cutoff_params_upper(upar, vthe, phi, me_over_mi, vpa, ir)
error("In upper-z electron bc, failed to find vpa=vcut point, plus_vcut_ind=$plus_vcut_ind")
end

# sigma is the location we use for w_∥(v_∥=0) - set to 0 to ignore the 'upar
# shift'
sigma = -u_over_vt

# sigma is between sigma_ind and sigma_ind+1
sigma_ind = searchsortedlast(vpa_unnorm, 0.0)
if sigma_ind < 1
Expand All @@ -2241,7 +2240,7 @@ function get_cutoff_params_upper(upar, vthe, phi, me_over_mi, vpa, ir)

# sigma_fraction is the fraction of the distance between sigma_ind+1 and
# sigma_ind where sigma is.
sigma_fraction = (sigma - vpa_unnorm[sigma_ind+1]) / (vpa_unnorm[sigma_ind] - vpa_unnorm[sigma_ind+1])
sigma_fraction = -vpa_unnorm[sigma_ind+1] / (vpa_unnorm[sigma_ind] - vpa_unnorm[sigma_ind+1])

# Want the element that contains the interval on the upper side of sigma_ind. For
# points on element boundaries, the `ielement` array contains the element on the lower
Expand Down

0 comments on commit 2a2a3c7

Please sign in to comment.