Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add particle weight as an explicit argument for _libwarpx.py::add_particles() #2161

Merged
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 15 additions & 12 deletions Python/pywarpx/_libwarpx.py
Original file line number Diff line number Diff line change
Expand Up @@ -370,7 +370,7 @@ def getCellSize(direction, level=0):
#
# libwarpx.warpx_ComputePMLFactors(lev, dt)

def add_particles(species_name, x=0., y=0., z=0., ux=0., uy=0., uz=0.,
def add_particles(species_name, x=0., y=0., z=0., ux=0., uy=0., uz=0., w=0.,
unique_particles=True, **kwargs):
'''

Expand All @@ -382,11 +382,12 @@ def add_particles(species_name, x=0., y=0., z=0., ux=0., uy=0., uz=0.,
species_name : the species to add the particle to
x, y, z : arrays or scalars of the particle positions (default = 0.)
ux, uy, uz : arrays or scalars of the particle momenta (default = 0.)
w : array or scalar of particle weights (default = 0.)
unique_particles : whether the particles are unique or duplicated on
several processes. (default = True)
kwargs : dictionary containing an entry for the particle weights
(with keyword 'w') and all the extra particle attribute
arrays. If an attribute is not given it will be set to 0.
kwargs : dictionary containing an entry for all the extra particle
attribute arrays. If an attribute is not given it will be
set to 0.

'''

Expand All @@ -397,18 +398,20 @@ def add_particles(species_name, x=0., y=0., z=0., ux=0., uy=0., uz=0.,
lenux = np.size(ux)
lenuy = np.size(uy)
lenuz = np.size(uz)
lenw = np.size(w)

if (lenx == 0 or leny == 0 or lenz == 0 or lenux == 0 or
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if we should make those and combined, to avoid (silent) returns where one parameter is only set to zero?
We have a check for mismatched lengths below.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cc @dpgrote , what do you think?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is tricky because of the default value of parameters. If a user calls add_particles(electrons, x=newx) each step to inject particles and on some steps the length of newx is zero, it should return silently. With the or this is what happens. Using and and relying on the length consistency check below would raise an error since the length of the default parameters is 1 and so maxlen would be 1.

Perhaps the default values should be changed to None with the argument handling taking that into account. With that it is known which arguments were supplied by the user and the consistency of those can be checked. This would require a fair amount of reworking of this code.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see. I think my confusion just lies into what happens with unprovided parameters. Do we default init those to zero or so?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, parameters not supplied default to 0., which has a length of 1. If multiple particles are being added, the zero is broadcast over all of the particles.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If an input parameter has a length of zero, using the default value instead is probably not the correct thing to do. In my example above, if newx has a length of zero, using the default value would erroneously add a particle with x=0., not the expected result.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I think I understand now. Yep, I agree that we should do a follow-up to use =None as default parameter to make things more clear/Pythonic.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can do that after this PR merges.

lenuy == 0 or lenuz == 0):
lenuy == 0 or lenuz == 0 or lenw == 0):
return

maxlen = max(lenx, leny, lenz, lenux, lenuy, lenuz)
maxlen = max(lenx, leny, lenz, lenux, lenuy, lenuz, lenw)
assert lenx==maxlen or lenx==1, "Length of x doesn't match len of others"
assert leny==maxlen or leny==1, "Length of y doesn't match len of others"
assert lenz==maxlen or lenz==1, "Length of z doesn't match len of others"
assert lenux==maxlen or lenux==1, "Length of ux doesn't match len of others"
assert lenuy==maxlen or lenuy==1, "Length of uy doesn't match len of others"
assert lenuz==maxlen or lenuz==1, "Length of uz doesn't match len of others"
assert lenw==maxlen or lenw==1, "Length of w doesn't match len of others"
for key, val in kwargs.items():
assert np.size(val)==1 or len(val)==maxlen, f"Length of {key} doesn't match len of others"

Expand All @@ -423,21 +426,21 @@ def add_particles(species_name, x=0., y=0., z=0., ux=0., uy=0., uz=0.,
if lenuy == 1:
uy = np.array(uy)*np.ones(maxlen)
if lenuz == 1:
uz = np.array(uz)*np.ones(maxlen,'d')
uz = np.array(uz)*np.ones(maxlen)
if lenw == 1:
w = np.array(w)*np.ones(maxlen)
for key, val in kwargs.items():
if np.size(val) == 1:
kwargs[key] = np.array(val)*np.ones(maxlen)

# --- The -3 is because the comps include the velocites
nattr = get_nattr_species(species_name) - 3
attr = np.zeros((maxlen, nattr))
attr[:,0] = w

for key, vals in kwargs.items():
if key == 'w':
attr[:,0] = vals
else:
# --- The -3 is because components 1 to 3 are velocities
attr[:,get_particle_comp_index(species_name, key)-3] = vals
# --- The -3 is because components 1 to 3 are velocities
attr[:,get_particle_comp_index(species_name, key)-3] = vals

libwarpx.warpx_addNParticles(
ctypes.c_char_p(species_name.encode('utf-8')), x.size,
Expand Down