Skip to content

Commit

Permalink
Include modification to where wall BC is imposed in moment-kinetic bc
Browse files Browse the repository at this point in the history
  • Loading branch information
@johnomotani
johnomotani committed Sep 6, 2024
1 parent 39a9449 commit 8fed26a
Show file tree
Hide file tree
Showing 2 changed files with 14 additions and 7 deletions.
18 changes: 12 additions & 6 deletions moment_kinetics/src/boundary_conditions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,8 @@ function enforce_z_boundary_condition!(pdf, density, upar, ppar, phi, moments, b
@loop_r ir begin
@views enforce_zero_incoming_bc!(
pdf[:,:,:,ir,is], z, vpa, density[:,ir,is], upar[:,ir,is],
ppar[:,ir,is], moments.evolve_upar, moments.evolve_ppar, zero)
ppar[:,ir,is], moments.evolve_upar, moments.evolve_ppar, zero,
phi[:,ir])
end
else
@loop_r ir begin
Expand Down Expand Up @@ -427,27 +428,32 @@ function enforce_zero_incoming_bc!(pdf, speed, z, zero, phi, epsz)
end
end
function get_ion_z_boundary_cutoff_indices(density, upar, ppar, evolve_upar, evolve_ppar,
z, vpa, zero)
z, vpa, zero, phi)
epsz = z.boundary_parameters.epsz
if z.irank == 0
deltaphi = phi[2] - phi[1]
vcut = deltaphi > 0 ? sqrt(deltaphi)*(epsz^0.25) : 0.0
vth = sqrt(2.0*(ppar[1]/density[1]))
@. vpa.scratch = vpagrid_to_dzdt(vpa.grid, vth,
upar[1], evolve_ppar, evolve_upar)
last_negative_vpa_ind = searchsortedlast(vpa.scratch, -zero)
last_negative_vpa_ind = searchsortedlast(vpa.scratch, min(-zero, -vcut))
else
last_negative_vpa_ind = nothing
end
if z.irank == z.nrank - 1
deltaphi = phi[end-1] - phi[end]
vcut = deltaphi > 0 ? sqrt(deltaphi)*(epsz^0.25) : 0.0
vth = sqrt(2.0*(ppar[end]/density[end]))
@. vpa.scratch2 = vpagrid_to_dzdt(vpa.grid, vth,
upar[end], evolve_ppar, evolve_upar)
first_positive_vpa_ind = searchsortedfirst(vpa.scratch2, zero)
first_positive_vpa_ind = searchsortedfirst(vpa.scratch2, max(zero, vcut))
else
first_positive_vpa_ind = nothing
end
return last_negative_vpa_ind, first_positive_vpa_ind
end
function enforce_zero_incoming_bc!(pdf, z::coordinate, vpa::coordinate, density, upar,
ppar, evolve_upar, evolve_ppar, zero)
ppar, evolve_upar, evolve_ppar, zero, phi)
if z.irank != 0 && z.irank != z.nrank - 1
# No z-boundary in this block
return nothing
Expand All @@ -461,7 +467,7 @@ function enforce_zero_incoming_bc!(pdf, z::coordinate, vpa::coordinate, density,
# absolute velocity at left boundary
last_negative_vpa_ind, first_positive_vpa_ind =
get_ion_z_boundary_cutoff_indices(density, upar, ppar, evolve_upar, evolve_ppar,
z, vpa, zero)
z, vpa, zero, phi)
if z.irank == 0
pdf[last_negative_vpa_ind+1:end, :, 1] .= 0.0
end
Expand Down
3 changes: 2 additions & 1 deletion moment_kinetics/src/time_advance.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2659,11 +2659,12 @@ function adaptive_timestep_update!(scratch, scratch_implicit, scratch_electron,
density = @view scratch[t_params.n_rk_stages+1].density[:,ir,is]
upar = @view scratch[t_params.n_rk_stages+1].upar[:,ir,is]
ppar = @view scratch[t_params.n_rk_stages+1].ppar[:,ir,is]
phi = fields.phi[:,ir]
last_negative_vpa_ind, first_positive_vpa_ind =
get_ion_z_boundary_cutoff_indices(density, upar, ppar,
moments.evolve_upar,
moments.evolve_ppar, z, vpa,
1.0e-14)
1.0e-14, phi)
if z.irank == 0
scratch[2].pdf[last_negative_vpa_ind,:,1,ir,is] .=
scratch[t_params.n_rk_stages+1].pdf[last_negative_vpa_ind,:,1,ir,is]
Expand Down

0 comments on commit 8fed26a

Please sign in to comment.