[Bug] profile = constant gives wrong density in RCYLINDER / RZ when xmin > 0
Summary
With injection_style = NUniformPerCell and profile = constant, the
realized initial density of a species whose radial bounding box has
xmin > 0 is silently scaled by a geometric factor. For p = radial_numpercell_power = 0 (the default) the realized density is
-
1D
RCYLINDER: $n_\text{realized} = n_\text{input} \cdot \dfrac{r_{\max} - r_{\min}}{r_{\max}}$
-
2D
RZ: $n_\text{realized} = n_\text{input} \cdot 2\pi \cdot \dfrac{r_{\max} - r_{\min}}{r_{\max}}$
(The extra $2\pi$ in RZ comes from the "per-radian" weight convention
used there; see "Analysis" below.)
The error is per-species and local (no cross-species coupling). It
reduces to unity only when rmin = 0, which is why classic single-wire
RZ tests that start at the axis do not see it.
The parse_density_function path is not affected — it uses the
correct local Jacobian per cell.
Why this matters
Multi-shell setups (e.g. nested wire arrays) and any RZ setup that
describes a plasma starting away from the axis get the wrong initial
density. With a thin shell $\delta \ll R$ the RCYLINDER error can be
$\mathcal{O}(10^{-1})$ or worse, and the RZ error can be anywhere from
$\mathcal{O}(1)$ to $\mathcal{O}(10)$ depending on geometry.
Specifically it means that <species>.density = n0 does not set the
physical volume density to n0 in these geometries, which contradicts
the documented contract.
Environment
- WarpX branch:
development
- Build: explicit EM, Yee solver, CPU, single MPI rank
- OS: Linux
Minimal reproducers
Two tiny 0-step input files (≤ 80 lines each) and a ~60-line Python
diagnostic are included below; I am happy to attach them as a gist or a
draft PR branch on request.
1. RCYLINDER reproducer
my_constants.n0 = 1.0e20
my_constants.rmin = 1.0e-3
my_constants.rmax = 2.0e-3
my_constants.R_wall = 3.0e-3
my_constants.Nppc = 100
max_step = 0
amr.n_cell = 512
geometry.dims = RCYLINDER
geometry.prob_lo = 0.0
geometry.prob_hi = R_wall
boundary.field_lo = none
boundary.field_hi = pec
boundary.particle_lo = none
boundary.particle_hi = absorbing
algo.maxwell_solver = Yee
algo.particle_pusher = higuera
particles.species_names = ions
ions.charge = q_e
ions.mass = m_u
ions.injection_style = "NUniformPerCell"
ions.num_particles_per_cell_each_dim = Nppc 1
ions.profile = constant
ions.density = n0
ions.xmin = rmin
ions.xmax = rmax
ions.momentum_distribution_type = "at_rest"
diagnostics.diags_names = diag1
diag1.intervals = 0:0
diag1.diag_type = Full
diag1.fields_to_plot = rho
diag1.ions.variables = x w
Expected (documented contract): $n \equiv n_0 = 1.0\times 10^{20}$ m$^{-3}$ in [1, 2] mm.
Observed: $n \approx 5.0\times 10^{19}$ m$^{-3}$ (exactly
$n_0\cdot(r_\max - r_\min)/r_\max$).
2. RZ reproducer
Same idea but with a thin-$z$ periodic box and a shell at [1, 10] mm.
Here the bug predicts $n_\text{realized} \approx 5.65\times 10^{20}$,
above the target, which rules out any "visualizer-undercount"
explanation.
my_constants.n0 = 1.0e20
my_constants.rmin = 1.0e-3
my_constants.rmax = 10.0e-3
my_constants.R_wall = 12.0e-3
my_constants.Lz = 1.0e-3
my_constants.Nppc = 8
max_step = 0
amr.n_cell = 256 32
geometry.dims = RZ
geometry.prob_lo = 0.0 0.0
geometry.prob_hi = R_wall Lz
boundary.field_lo = none periodic
boundary.field_hi = pec periodic
boundary.particle_lo = none periodic
boundary.particle_hi = absorbing periodic
algo.maxwell_solver = Yee
algo.particle_pusher = higuera
particles.species_names = ions
ions.charge = q_e
ions.mass = m_u
ions.injection_style = "NUniformPerCell"
ions.num_particles_per_cell_each_dim = Nppc Nppc 1
ions.profile = constant
ions.density = n0
ions.xmin = rmin
ions.xmax = rmax
ions.momentum_distribution_type = "at_rest"
diagnostics.diags_names = diag1
diag1.intervals = 0:0
diag1.diag_type = Full
diag1.fields_to_plot = rho
diag1.ions.variables = x z w
Diagnostic script
r = ad[("ions","particle_position_x")]
w = ad[("ions","particle_weight")]
# Documented contract in RZ: weight = n0 * 2*pi*r*dr*dz / N_ppc (per radian)
# Bug formula: weight = n0 * dr*dz/N_ppc * 2*pi*(rmax-rmin) * r/rmax
# -> w/r collapses to a constant in both; only the absolute value differs.
n0_inferred = median(w/r) * Nppc / (2*pi*dr*dz) # RZNumerical check printed by check_density.py on output from
inputs_2d_rz_repro:
observed w/r : <matches bug formula to ~0.1 %>
inferred n0 : ~5.65e20 (~5.65 x target)
and on output from inputs_1d_rcylinder_repro:
observed w/r : <matches bug formula to ~0.1 %>
inferred n0 : ~5.0e19 (0.5 x target)
Analysis — source of the error
The weight assignment lives in
Source/Particles/ParticleCreation/AddParticles.cpp,
inside PhysicalParticleContainer::AddPlasma(...):
amrex::Real weight = dens;
weight *= scale_fac;
#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER)
const amrex::Real coeff = 2._rt*MathConst::pi/(1._rt + radial_numpercell_power)
*(rmax - std::pow(rmax, -radial_numpercell_power)
*std::pow(rmin, 1._rt + radial_numpercell_power));
weight *= coeff*std::pow(xb/rmax, 1._rt - radial_numpercell_power);
#elif defined(WARPX_DIM_RSPHERE)
const amrex::Real coeff = 4._rt*MathConst::pi/(1._rt + radial_numpercell_power)
*(rmax*rmax - std::pow(rmax, 1._rt - radial_numpercell_power )
*std::pow(rmin, 1._rt + radial_numpercell_power));
weight *= coeff*std::pow(xb/rmax, 2._rt - radial_numpercell_power);
#endifFor p = 0 this simplifies to
$$
w ;=; n_0 \cdot \frac{\mathrm{d}r,\mathrm{d}z}{N_\text{ppc}} \cdot 2\pi,(r_\max - r_\min) \cdot \frac{x_b}{r_\max}.
$$
The correct local per-cell Jacobian (with full-ring weight convention)
would be $2\pi,x_b$, not $2\pi,(r_\max - r_\min),x_b / r_\max$.
The ratio
$$
\frac{r_\max - r_\min}{r_\max}
$$
equals one only when rmin = 0, which is why this has not surfaced in
single-column RZ tests.
An analogous issue appears in the WARPX_DIM_RSPHERE branch just below.
The extra $2\pi$ in RZ
In RZ the stored particle weight is defined as "particles per radian" —
the full-ring count is $2\pi\cdot w$, and the rho deposition and
reductions divide by $2\pi r$ accordingly. Because the buggy coeff
above is applied identically in RZ and RCYLINDER, in RZ it happens to
combine with the per-radian convention so that the observed density
error factor is $2\pi,(r_\max-r_\min)/r_\max$ rather than
$(r_\max-r_\min)/r_\max$. We have verified this numerically in both
geometries (see diagnostic output above).
Scope: what is affected / not affected
Affected:
WARPX_DIM_RCYLINDER, WARPX_DIM_RZ, WARPX_DIM_RSPHEREinjection_style = NUniformPerCell with profile = constant and
xmin > 0.
Not affected (verified):
profile = parse_density_function — uses local per-cell Jacobian.- Any setup with
xmin = 0 (the default).
Suggested fix
For p = 0, replace the shell-integrated prefactor with the correct
local Jacobian:
// p = 0 case:
weight *= 2._rt * MathConst::pi * xb;
For general p the prefactor $2\pi/(1+p),(r_\max - r_\max^{-p},r_\min^{1+p})$
looks like it came from analytically integrating $r^{1+p}$ over
[0, rmax] assuming rmin = 0; the correct per-cell statement
independent of rmin, rmax is
$$
w ;=; \frac{n_0 \cdot 2\pi,x_b,\mathrm{d}r,\mathrm{d}z}{N_\text{ppc}(x_b)}
$$
with $N_\text{ppc}(x_b) \propto (x_b/r_\max)^p$, which evaluates to the
same final weight without any appearance of $r_\min$.
Would be happy to open a draft PR once the preferred semantics are
confirmed, including a regression test along the lines of the
reproducers above (two nested shells with profile = constant and a
target uniform density).
Workaround for current users
Pre-multiply each species .density by
$$
\frac{r_\max}{r_\max - r_\min}
$$
in RCYLINDER (or by $\dfrac{r_\max}{2\pi,(r_\max - r_\min)}$ in RZ, if
you also need to compensate the per-radian convention that surfaces in
your diagnostic). This is applied per species — no cross-species
coupling is needed.
[Bug]
profile = constantgives wrong density inRCYLINDER/RZwhenxmin > 0Summary
With
injection_style = NUniformPerCellandprofile = constant, therealized initial density of a species whose radial bounding box has
xmin > 0is silently scaled by a geometric factor. Forp = radial_numpercell_power = 0(the default) the realized density isRCYLINDER:RZ:(The extra$2\pi$ in RZ comes from the "per-radian" weight convention
used there; see "Analysis" below.)
The error is per-species and local (no cross-species coupling). It
reduces to unity only when
rmin = 0, which is why classic single-wireRZ tests that start at the axis do not see it.
The
parse_density_functionpath is not affected — it uses thecorrect local Jacobian per cell.
Why this matters
Multi-shell setups (e.g. nested wire arrays) and any RZ setup that$\delta \ll R$ the RCYLINDER error can be
$\mathcal{O}(10^{-1})$ or worse, and the RZ error can be anywhere from
$\mathcal{O}(1)$ to $\mathcal{O}(10)$ depending on geometry.
describes a plasma starting away from the axis get the wrong initial
density. With a thin shell
Specifically it means that
<species>.density = n0does not set thephysical volume density to
n0in these geometries, which contradictsthe documented contract.
Environment
developmentMinimal reproducers
Two tiny 0-step input files (≤ 80 lines each) and a ~60-line Python
diagnostic are included below; I am happy to attach them as a gist or a
draft PR branch on request.
1.
RCYLINDERreproducerExpected (documented contract):$n \equiv n_0 = 1.0\times 10^{20}$ m$^{-3}$ in
[1, 2]mm.Observed:$n \approx 5.0\times 10^{19}$ m$^{-3}$ (exactly
$n_0\cdot(r_\max - r_\min)/r_\max$ ).
2.
RZreproducerSame idea but with a thin-$z$ periodic box and a shell at$n_\text{realized} \approx 5.65\times 10^{20}$ ,
[1, 10] mm.Here the bug predicts
above the target, which rules out any "visualizer-undercount"
explanation.
Diagnostic script
Numerical check printed by
check_density.pyon output frominputs_2d_rz_repro:and on output from
inputs_1d_rcylinder_repro:Analysis — source of the error
The weight assignment lives in
Source/Particles/ParticleCreation/AddParticles.cpp,inside
PhysicalParticleContainer::AddPlasma(...):For
p = 0this simplifies toThe correct local per-cell Jacobian (with full-ring weight convention)$2\pi,x_b$ , not $2\pi,(r_\max - r_\min),x_b / r_\max$ .
would be
The ratio
equals one only when
rmin = 0, which is why this has not surfaced insingle-column RZ tests.
An analogous issue appears in the
WARPX_DIM_RSPHEREbranch just below.The extra$2\pi$ in RZ
In RZ the stored particle weight is defined as "particles per radian" —$2\pi\cdot w$ , and the rho deposition and$2\pi r$ accordingly. Because the buggy $2\pi,(r_\max-r_\min)/r_\max$ rather than
$(r_\max-r_\min)/r_\max$ . We have verified this numerically in both
the full-ring count is
reductions divide by
coeffabove is applied identically in RZ and RCYLINDER, in RZ it happens to
combine with the per-radian convention so that the observed density
error factor is
geometries (see diagnostic output above).
Scope: what is affected / not affected
Affected:
WARPX_DIM_RCYLINDER,WARPX_DIM_RZ,WARPX_DIM_RSPHEREinjection_style = NUniformPerCellwithprofile = constantandxmin > 0.Not affected (verified):
profile = parse_density_function— uses local per-cell Jacobian.xmin = 0(the default).Suggested fix
For
p = 0, replace the shell-integrated prefactor with the correctlocal Jacobian:
For general$2\pi/(1+p),(r_\max - r_\max^{-p},r_\min^{1+p})$ $r^{1+p}$ over
pthe prefactorlooks like it came from analytically integrating
[0, rmax]assumingrmin = 0; the correct per-cell statementindependent of
rmin, rmaxiswith$N_\text{ppc}(x_b) \propto (x_b/r_\max)^p$ , which evaluates to the$r_\min$ .
same final weight without any appearance of
Would be happy to open a draft PR once the preferred semantics are
confirmed, including a regression test along the lines of the
reproducers above (two nested shells with
profile = constantand atarget uniform density).
Workaround for current users
Pre-multiply each species
.densitybyin RCYLINDER (or by$\dfrac{r_\max}{2\pi,(r_\max - r_\min)}$ in RZ, if
you also need to compensate the per-radian convention that surfaces in
your diagnostic). This is applied per species — no cross-species
coupling is needed.