Skip to content

Commit

Permalink
Switch from reconcile_element_boundaries_MPI!() using centred averagi…
Browse files Browse the repository at this point in the history 
…ng of element boundaries to upwind choice for element boundaries. This requires extra dummy arrays to be passed to the boundary condition functions.
  • Loading branch information
@mrhardman
mrhardman committed Dec 12, 2023
1 parent ee12e1a commit 816e058
Showing 1 changed file with 36 additions and 21 deletions.
57 changes: 36 additions & 21 deletions src/initial_conditions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1157,14 +1157,15 @@ function enforce_boundary_conditions!(f, f_r_bc, density, upar, ppar, moments, v
@views enforce_z_boundary_condition!(f, density, upar, ppar, moments, z_bc, z_adv, z,
vperp, vpa, composition,
scratch_dummy.buffer_vpavperprs_1, scratch_dummy.buffer_vpavperprs_2,
scratch_dummy.buffer_vpavperprs_3, scratch_dummy.buffer_vpavperprs_4)
scratch_dummy.buffer_vpavperprs_3, scratch_dummy.buffer_vpavperprs_4,
scratch_dummy.buffer_vpavperprs_5, scratch_dummy.buffer_vpavperprs_6)

end
if r.n > 1
begin_s_z_vperp_vpa_region()
@views enforce_r_boundary_condition!(f, f_r_bc, r_bc, r_adv, vpa, vperp, z, r, composition,
scratch_dummy.buffer_vpavperpzs_1, scratch_dummy.buffer_vpavperpzs_2,
scratch_dummy.buffer_vpavperpzs_3, scratch_dummy.buffer_vpavperpzs_4,
scratch_dummy.buffer_vpavperpzs_3, scratch_dummy.buffer_vpavperpzs_4, scratch_dummy.buffer_vpavperpzs_5, scratch_dummy.buffer_vpavperpzs_6,
r_diffusion)
end
end
Expand All @@ -1181,7 +1182,8 @@ end
enforce boundary conditions on f in r
"""
function enforce_r_boundary_condition!(f::AbstractArray{mk_float,5}, f_r_bc, bc::String,
adv, vpa, vperp, z, r, composition, end1::AbstractArray{mk_float,4},
r_advect, vpa, vperp, z, r, composition, adv1::AbstractArray{mk_float,4},
adv2::AbstractArray{mk_float,4}, end1::AbstractArray{mk_float,4},
end2::AbstractArray{mk_float,4}, buffer1::AbstractArray{mk_float,4},
buffer2::AbstractArray{mk_float,4}, r_diffusion::Bool)

Expand All @@ -1191,10 +1193,12 @@ function enforce_r_boundary_condition!(f::AbstractArray{mk_float,5}, f_r_bc, bc:
# reconcile internal element boundaries across processes
# & enforce periodicity and external boundaries if needed
@loop_s_z_vperp_vpa is iz ivperp ivpa begin
adv1[ivpa,ivperp,iz,is] = r_advect[is].adv_fac[1,ivpa,ivperp,iz]
adv2[ivpa,ivperp,iz,is] = r_advect[is].adv_fac[end,ivpa,ivperp,iz]
end1[ivpa,ivperp,iz,is] = f[ivpa,ivperp,iz,1,is]
end2[ivpa,ivperp,iz,is] = f[ivpa,ivperp,iz,nr,is]
end
@views reconcile_element_boundaries_MPI!(f,
@views reconcile_element_boundaries_MPI!(f, adv1, adv2,
end1, end2, buffer1, buffer2, r)
end

Expand All @@ -1216,11 +1220,11 @@ function enforce_r_boundary_condition!(f::AbstractArray{mk_float,5}, f_r_bc, bc:
# impose bc on both sides of the domain to accomodate a diffusion operator d^2 / d r^2
@loop_s_z_vperp_vpa is iz ivperp ivpa begin
ir = 1 # r = -L/2 -- check that the point is on lowest rank
if r.irank == 0 && (r_diffusion || adv[is].speed[ir,ivpa,ivperp,iz] > zero)
if r.irank == 0 && (r_diffusion || r_advect[is].speed[ir,ivpa,ivperp,iz] > zero)
f[ivpa,ivperp,iz,ir,is] = f_r_bc[ivpa,ivperp,iz,1,is]
end
ir = r.n # r = L/2 -- check that the point is on highest rank
if r.irank == r.nrank - 1 && (r_diffusion || adv[is].speed[ir,ivpa,ivperp,iz] < -zero)
if r.irank == r.nrank - 1 && (r_diffusion || r_advect[is].speed[ir,ivpa,ivperp,iz] < -zero)
f[ivpa,ivperp,iz,ir,is] = f_r_bc[ivpa,ivperp,iz,end,is]
end
end
Expand All @@ -1230,8 +1234,9 @@ end
"""
enforce boundary conditions on charged particle f in z
"""
function enforce_z_boundary_condition!(pdf, density, upar, ppar, moments, bc::String, adv,
z, vperp, vpa, composition, end1::AbstractArray{mk_float,4},
function enforce_z_boundary_condition!(pdf, density, upar, ppar, moments, bc::String, z_advect,
z, vperp, vpa, composition, adv1::AbstractArray{mk_float,4},
adv2::AbstractArray{mk_float,4}, end1::AbstractArray{mk_float,4},
end2::AbstractArray{mk_float,4}, buffer1::AbstractArray{mk_float,4},
buffer2::AbstractArray{mk_float,4})
# this block ensures periodic BC can be supported with distributed memory MPI
Expand All @@ -1240,11 +1245,13 @@ function enforce_z_boundary_condition!(pdf, density, upar, ppar, moments, bc::St
# & enforce periodicity and external boundaries if needed
nz = z.n
@loop_s_r_vperp_vpa is ir ivperp ivpa begin
adv1[ivpa,ivperp,ir,is] = z_advect[is].adv_fac[1,ivpa,ivperp,ir]
adv2[ivpa,ivperp,ir,is] = z_advect[is].adv_fac[end,ivpa,ivperp,ir]
end1[ivpa,ivperp,ir,is] = pdf[ivpa,ivperp,1,ir,is]
end2[ivpa,ivperp,ir,is] = pdf[ivpa,ivperp,nz,ir,is]
end
# check on periodic bc happens inside this call below
@views reconcile_element_boundaries_MPI!(pdf,
@views reconcile_element_boundaries_MPI!(pdf, adv1, adv2,
end1, end2, buffer1, buffer2, z)
end
# define a zero that accounts for finite precision
Expand All @@ -1257,14 +1264,14 @@ function enforce_z_boundary_condition!(pdf, density, upar, ppar, moments, bc::St
vwidth = 1.0
if z.irank == 0
@loop_s_r_vperp_vpa is ir ivperp ivpa begin
if adv[is].speed[ivpa,1,ir] > 0.0
if z_advect[is].speed[ivpa,1,ir] > 0.0
pdf[ivpa,ivperp,1,ir,is] = density_offset * exp(-(vpa.grid[ivpa]^2 + vperp.grid[ivperp]^2)/vwidth^2) / sqrt(pi)
end
end
end
if z.irank == z.nrank - 1
@loop_s_r_vperp_vpa is ir ivperp ivpa begin
if adv[is].speed[ivpa,end,ir] > 0.0
if z_advect[is].speed[ivpa,end,ir] > 0.0
pdf[ivpa,ivperp,end,ir,is] = density_offset * exp(-(vpa.grid[ivpa]^2 + vperp.grid[ivperp]^2)/vwidth^2) / sqrt(pi)
end
end
Expand Down Expand Up @@ -1292,7 +1299,7 @@ function enforce_z_boundary_condition!(pdf, density, upar, ppar, moments, bc::St
else
@loop_r ir begin
@views enforce_zero_incoming_bc!(pdf[:,:,:,ir,is],
adv[is].speed[:,:,:,ir], z, zero)
z_advect[is].speed[:,:,:,ir], z, zero)
end
end
end
Expand Down Expand Up @@ -1348,20 +1355,23 @@ function enforce_neutral_boundary_conditions!(f_neutral, f_charged,
pz_neutral, moments, density_ion, upar_ion, Er, boundary_distributions,
z_adv, z, vzeta, vr, vz, composition, geometry,
scratch_dummy.buffer_vzvrvzetarsn_1, scratch_dummy.buffer_vzvrvzetarsn_2,
scratch_dummy.buffer_vzvrvzetarsn_3, scratch_dummy.buffer_vzvrvzetarsn_4)
scratch_dummy.buffer_vzvrvzetarsn_3, scratch_dummy.buffer_vzvrvzetarsn_4,
scratch_dummy.buffer_vzvrvzetarsn_5, scratch_dummy.buffer_vzvrvzetarsn_6)
end
if r.n > 1
begin_sn_z_vzeta_vr_vz_region()
@views enforce_neutral_r_boundary_condition!(f_neutral, boundary_distributions.pdf_rboundary_neutral,
r_adv, vz, vr, vzeta, z, r, composition,
scratch_dummy.buffer_vzvrvzetazsn_1, scratch_dummy.buffer_vzvrvzetazsn_2,
scratch_dummy.buffer_vzvrvzetazsn_3, scratch_dummy.buffer_vzvrvzetazsn_4,
scratch_dummy.buffer_vzvrvzetazsn_5, scratch_dummy.buffer_vzvrvzetazsn_6,
r_diffusion)
end
end

function enforce_neutral_r_boundary_condition!(f::AbstractArray{mk_float,6},
f_r_bc::AbstractArray{mk_float,6}, adv, vz, vr, vzeta, z, r, composition,
f_r_bc::AbstractArray{mk_float,6}, r_advect, vz, vr, vzeta, z, r, composition,
adv1::AbstractArray{mk_float,5}, adv2::AbstractArray{mk_float,5},
end1::AbstractArray{mk_float,5}, end2::AbstractArray{mk_float,5},
buffer1::AbstractArray{mk_float,5}, buffer2::AbstractArray{mk_float,5},
r_diffusion) #f_initial,
Expand All @@ -1373,10 +1383,12 @@ function enforce_neutral_r_boundary_condition!(f::AbstractArray{mk_float,6},
# reconcile internal element boundaries across processes
# & enforce periodicity and external boundaries if needed
@loop_sn_z_vzeta_vr_vz isn iz ivzeta ivr ivz begin
adv1[ivz,ivr,ivzeta,iz,isn] = r_advect[isn].speed[1,ivz,ivr,ivzeta,iz]
adv2[ivz,ivr,ivzeta,iz,isn] = r_advect[isn].speed[nr,ivz,ivr,ivzeta,iz]
end1[ivz,ivr,ivzeta,iz,isn] = f[ivz,ivr,ivzeta,iz,1,isn]
end2[ivz,ivr,ivzeta,iz,isn] = f[ivz,ivr,ivzeta,iz,nr,isn]
end
@views reconcile_element_boundaries_MPI!(f,
@views reconcile_element_boundaries_MPI!(f, adv1, adv2,
end1, end2, buffer1, buffer2, r)
end
# 'periodic' BC enforces periodicity by taking the average of the boundary points
Expand All @@ -1396,12 +1408,12 @@ function enforce_neutral_r_boundary_condition!(f::AbstractArray{mk_float,6},
@loop_sn_z_vzeta_vr_vz isn iz ivzeta ivr ivz begin
ir = 1 # r = -L/2
# incoming particles and on lowest rank
if r.irank == 0 && (r_diffusion || adv[isn].speed[ir,ivz,ivr,ivzeta,iz] > zero)
if r.irank == 0 && (r_diffusion || r_advect[isn].speed[ir,ivz,ivr,ivzeta,iz] > zero)
f[ivz,ivr,ivzeta,iz,ir,isn] = f_r_bc[ivz,ivr,ivzeta,iz,1,isn]
end
ir = nr # r = L/2
# incoming particles and on highest rank
if r.irank == r.nrank - 1 && (r_diffusion || adv[isn].speed[ir,ivz,ivr,ivzeta,iz] < -zero)
if r.irank == r.nrank - 1 && (r_diffusion || r_advect[isn].speed[ir,ivz,ivr,ivzeta,iz] < -zero)
f[ivz,ivr,ivzeta,iz,ir,isn] = f_r_bc[ivz,ivr,ivzeta,iz,end,isn]
end
end
Expand All @@ -1412,8 +1424,9 @@ end
enforce boundary conditions on neutral particle f in z
"""
function enforce_neutral_z_boundary_condition!(pdf, density, uz, pz, moments, density_ion,
upar_ion, Er, boundary_distributions, adv,
upar_ion, Er, boundary_distributions, z_advect,
z, vzeta, vr, vz, composition, geometry,
adv1::AbstractArray{mk_float,5}, adv2::AbstractArray{mk_float,5},
end1::AbstractArray{mk_float,5}, end2::AbstractArray{mk_float,5},
buffer1::AbstractArray{mk_float,5}, buffer2::AbstractArray{mk_float,5})

Expand All @@ -1423,11 +1436,13 @@ function enforce_neutral_z_boundary_condition!(pdf, density, uz, pz, moments, de
# & enforce periodicity and external boundaries if needed
nz = z.n
@loop_sn_r_vzeta_vr_vz isn ir ivzeta ivr ivz begin
adv1[ivz,ivr,ivzeta,ir,isn] = z_advect[isn].speed[1,ivz,ivr,ivzeta,ir]
adv2[ivz,ivr,ivzeta,ir,isn] = z_advect[isn].speed[nz,ivz,ivr,ivzeta,ir]
end1[ivz,ivr,ivzeta,ir,isn] = pdf[ivz,ivr,ivzeta,1,ir,isn]
end2[ivz,ivr,ivzeta,ir,isn] = pdf[ivz,ivr,ivzeta,nz,ir,isn]
end
# check on periodic bc occurs within this call below
@views reconcile_element_boundaries_MPI!(pdf,
@views reconcile_element_boundaries_MPI!(pdf, adv1, adv2,
end1, end2, buffer1, buffer2, z)
end

Expand All @@ -1440,7 +1455,7 @@ function enforce_neutral_z_boundary_condition!(pdf, density, uz, pz, moments, de
vwidth = 1.0
if z.irank == 0
@loop_sn_r_vzeta_vr_vz isn ir ivzeta ivr ivz begin
if adv[isn].speed[ivz,ivr,ivzeta,1,ir] > 0.0
if z_advect[isn].speed[1,ivz,ivr,ivzeta,ir] > 0.0
pdf[ivz,ivr,ivzeta,1,ir,is] = density_offset *
exp(-(vzeta.grid[ivzeta]^2 + vr.grid[ivr] + vz.grid[ivz])/vwidth^2) /
sqrt(pi)
Expand All @@ -1449,7 +1464,7 @@ function enforce_neutral_z_boundary_condition!(pdf, density, uz, pz, moments, de
end
if z.irank == z.nrank - 1
@loop_sn_r_vzeta_vr_vz isn ir ivzeta ivr ivz begin
if adv[isn].speed[ivz,ivr,ivzeta,end,ir] > 0.0
if z_advect[isn].speed[end,ivz,ivr,ivzeta,ir] > 0.0
pdf[ivz,ivr,ivzeta,end,ir,is] = density_offset *
exp(-(vzeta.grid[ivzeta]^2 + vr.grid[ivr] + vz.grid[ivz])/vwidth^2) /
sqrt(pi)
Expand Down

0 comments on commit 816e058

Please sign in to comment.