search

ECP-WarpX / WarpX

WarpX is an advanced electromagnetic & electrostatic Particle-In-Cell code.
https://ecp-warpx.github.io
Other
292 stars 188 forks source link

Periodic boundaries eat particles in single precision #3347

Open KZhu-ME opened 2 years ago

KZhu-ME commented 2 years ago

We noticed that particle counts are being depleted in a single precision (fields and particles) simulation with only a reflecting and periodic boundaries. We have seen similar issues of particle depletion in our other simulation runs/plots, posted/referenced below.

We've replicated this issue using the cpp input file attached.

single_vs_double_ion_zavg_100 0_ns

Originally posted by @roelof-groenewald in https://github.com/ECP-WarpX/WarpX/issues/3231#issuecomment-1203174466

KZhu-ME commented 2 years ago 

# This is a 2D electrostatic simulation containing a circle embedded boundary at
# the center of the grid domain, pre-seeded with equal density of argon ions
# and electrons.

max_step = 25000
warpx.const_dt = 3.99e-13
warpx.do_electrostatic = labframe
warpx.self_fields_required_precision = 1e-3
warpx.self_fields_verbosity = 0
warpx.eb_potential(x,y,z,t) = -10
warpx.eb_implicit_function = -((x-0.00005)**2+(z-0.00005)**2-1e-05**2)
warpx.self_fields_absolute_tolerance = 0.02

algo.load_balance_intervals = 5
algo.load_balance_costs_update = timers
algo.load_balance_efficiency_ratio_threshold = 1.001
algo.load_balance_with_sfc = 0
algo.load_balance_knapsack_factor = 2

amr.n_cell = 128 128
amr.max_grid_size = 16
amr.max_level = 0
geometry.dims = 2
geometry.prob_lo = 0.0 0.0
geometry.prob_hi = 0.0001 0.0001
boundary.field_lo = pec periodic
boundary.field_hi = pec periodic
boundary.particle_lo = reflecting periodic
boundary.particle_hi = reflecting periodic
algo.particle_shape = 1

particles.species_names = electrons ar_ions
electrons.species_type = electron
electrons.injection_style = NRandomPerCell
electrons.initialize_self_fields = 1
electrons.num_particles_per_cell = 5
electrons.density = 10.2e16
electrons.profile = constant
electrons.momentum_distribution_type = maxwell_boltzmann
electrons.theta = (kb*30000/(m_e*clight^2))
electrons.save_particles_at_zlo = 1
electrons.save_particles_at_zhi = 1
electrons.save_particles_at_eb = 1

ar_ions.species_type = argon
ar_ions.charge = q_e
ar_ions.injection_style = NRandomPerCell
ar_ions.initialize_self_fields = 1
ar_ions.num_particles_per_cell = 5
ar_ions.profile = constant
ar_ions.density = 10.2e16
ar_ions.momentum_distribution_type = gaussian
ar_ions.ux_m = 0.0
ar_ions.uy_m = 0.0
ar_ions.uz_m = 0.0
ar_ions.ux_th = 2.6285641070031447e-06
ar_ions.uy_th = 2.6285641070031447e-06
ar_ions.uz_th = 2.6285641070031447e-06
ar_ions.save_particles_at_zlo = 1
ar_ions.save_particles_at_zhi = 1
ar_ions.save_particles_at_eb = 1

diagnostics.diags_names = diag1
diag1.diag_type = Full
diag1.format = plotfile
diag1.intervals = 1
diag1.fields_to_plot = phi rho_electrons rho_ar_ions```
PhilMiller commented 2 years ago

One thought I had in looking at the code is that the arithmetic in AMReX's enforcePeriodic feels a bit sketchy. It does iterated addition/subtraction of the domain length from the particle's position, until it's in bounds. That can potentially accumulate error, but I don't know where particles would then get lost

PhilMiller commented 2 years ago

In investigating this, there are a number of distinct issues in play. I'm digging into the one that's of most immediate concern for Modern Electron's usage, which appears with a switch from DP fields/DP particles to DP fields/SP particles. Going to SP/SP exposes much more dramatic physical differences that we're not prepared to isolate or tackle yet.

Bottom line: Work on this issue can probably be put on hold at least temporarily, if there are other things for you to focus on