diff options
| -rw-r--r-- | Python/pywarpx/AMReX.py | 7 | ||||
| -rw-r--r-- | Python/pywarpx/Bucket.py | 4 | ||||
| -rw-r--r-- | Python/pywarpx/PGroup.py | 71 | ||||
| -rw-r--r-- | Python/pywarpx/PICMI.py | 24 | ||||
| -rwxr-xr-x | Python/pywarpx/_libwarpx.py | 2 | ||||
| -rw-r--r-- | tests/plasma_acceleration/plasma_acceleration.py | 3 | ||||
| -rw-r--r-- | tests/plasma_acceleration/plasma_acceleration_PICMI.py | 229 |
7 files changed, 303 insertions, 37 deletions
diff --git a/Python/pywarpx/AMReX.py b/Python/pywarpx/AMReX.py index 581133695..bdadc49bd 100644 --- a/Python/pywarpx/AMReX.py +++ b/Python/pywarpx/AMReX.py @@ -6,6 +6,7 @@ from .Geometry import geometry from .Algo import algo from .Langmuirwave import langmuirwave from .Interpolation import interpolation +from . import Particles from .Particles import particles, particles_list import ctypes @@ -24,6 +25,12 @@ class AMReX(object): argv += interpolation.attrlist() argv += particles.attrlist() + if not particles_list: + # --- This is needed in case only species_names has been set, + # --- assuming that only the built in particle types are being used. + for pstring in particles.species_names.split(' '): + particles_list.append(getattr(Particles, pstring)) + for particle in particles_list: argv += particle.attrlist() diff --git a/Python/pywarpx/Bucket.py b/Python/pywarpx/Bucket.py index 2f4a1704c..916b030d1 100644 --- a/Python/pywarpx/Bucket.py +++ b/Python/pywarpx/Bucket.py @@ -25,7 +25,9 @@ class Bucket(object): "Concatenate the attributes into a string" result = [] for attr, value in self.argvattrs.iteritems(): - attrstring = '{0}.{1}={2} '.format(self.instancename, attr, value) + # --- repr is applied to value so that for floats, all of the digits are included. + # --- The strip is then needed when value is a string. + attrstring = '{0}.{1}={2} '.format(self.instancename, attr, repr(value).strip("'\"")) result += [attrstring] return result diff --git a/Python/pywarpx/PGroup.py b/Python/pywarpx/PGroup.py index cbc30943b..fa8572c0a 100644 --- a/Python/pywarpx/PGroup.py +++ b/Python/pywarpx/PGroup.py @@ -6,8 +6,9 @@ class PGroup(object): """Implements a class that has the same API as a warp ParticleGroup instance. """ - def __init__(self, igroup): + def __init__(self, igroup, ispecie): self.igroup = igroup + self.ispecie = ispecie self.ns = 1 # Number of species self.gallot() @@ -75,81 +76,85 @@ class PGroup(object): return self.nps.sum() npmax = property(getnpmax) - def getxp(self, js=0): - return _libwarpx.get_particle_x(js)[self.igroup] + def getxp(self): + return _libwarpx.get_particle_x(self.ispecie)[self.igroup] xp = property(getxp) - def getyp(self, js=0): - return _libwarpx.get_particle_y(js)[self.igroup] + def getyp(self): + return _libwarpx.get_particle_y(self.ispecie)[self.igroup] yp = property(getyp) - def getzp(self, js=0): - return _libwarpx.get_particle_z(js)[self.igroup] + def getzp(self): + return _libwarpx.get_particle_z(self.ispecie)[self.igroup] zp = property(getzp) - def getuxp(self, js=0): - return _libwarpx.get_particle_ux(js)[self.igroup] + def getuxp(self): + return _libwarpx.get_particle_ux(self.ispecie)[self.igroup] uxp = property(getuxp) - def getuyp(self, js=0): - return _libwarpx.get_particle_uy(js)[self.igroup] + def getuyp(self): + return _libwarpx.get_particle_uy(self.ispecie)[self.igroup] uyp = property(getuyp) - def getuzp(self, js=0): - return _libwarpx.get_particle_uz(js)[self.igroup] + def getuzp(self): + return _libwarpx.get_particle_uz(self.ispecie)[self.igroup] uzp = property(getuzp) - def getpid(self, js=0): - id = _libwarpx.get_particle_id(js)[self.igroup] + def getpid(self): + id = _libwarpx.get_particle_id(self.ispecie)[self.igroup] pid = np.array([id]).T return pid pid = property(getpid) - def getgaminv(self, js=0): - uxp = self.getuxp(js) - uyp = self.getuyp(js) - uzp = self.getuzp(js) + def getgaminv(self): + uxp = self.getuxp() + uyp = self.getuyp() + uzp = self.getuzp() return np.sqrt(1. - (uxp**2 + uyp**2 + uzp**2)/_libwarpx.clight**2) gaminv = property(getgaminv) - def getex(self, js=0): - return _libwarpx.get_particle_Ex(js)[self.igroup] + def getex(self): + return _libwarpx.get_particle_Ex(self.ispecie)[self.igroup] ex = property(getex) - def getey(self, js=0): - return _libwarpx.get_particle_Ey(js)[self.igroup] + def getey(self): + return _libwarpx.get_particle_Ey(self.ispecie)[self.igroup] ey = property(getey) - def getez(self, js=0): - return _libwarpx.get_particle_Ez(js)[self.igroup] + def getez(self): + return _libwarpx.get_particle_Ez(self.ispecie)[self.igroup] ez = property(getez) - def getbx(self, js=0): - return _libwarpx.get_particle_Bx(js)[self.igroup] + def getbx(self): + return _libwarpx.get_particle_Bx(self.ispecie)[self.igroup] bx = property(getbx) - def getby(self, js=0): - return _libwarpx.get_particle_By(js)[self.igroup] + def getby(self): + return _libwarpx.get_particle_By(self.ispecie)[self.igroup] by = property(getby) - def getbz(self, js=0): - return _libwarpx.get_particle_Bz(js)[self.igroup] + def getbz(self): + return _libwarpx.get_particle_Bz(self.ispecie)[self.igroup] bz = property(getbz) class PGroups(object): def __init__(self, ispecie=0): self.ispecie = ispecie - xall = _libwarpx.get_particle_x(ispecie) + + def setuppgroups(self): + xall = _libwarpx.get_particle_x(self.ispecie) self.ngroups = len(xall) self._pgroups = [] for igroup in range(self.ngroups): - self._pgroups.append(PGroup(igroup)) + self._pgroups.append(PGroup(igroup, self.ispecie)) def __iter__(self): + self.setuppgroups() for igroup in range(self.ngroups): yield self._pgroups[igroup] def __getitem__(self, key): + self.setuppgroups() return self._pgroups[key] diff --git a/Python/pywarpx/PICMI.py b/Python/pywarpx/PICMI.py index 95bb758dc..30957970e 100644 --- a/Python/pywarpx/PICMI.py +++ b/Python/pywarpx/PICMI.py @@ -4,6 +4,8 @@ from PICMI_Base import * import numpy as np from pywarpx import * +codename = 'WarpX' + class Grid(PICMI_Grid): def init(self, **kw): @@ -20,8 +22,26 @@ class Grid(PICMI_Grid): # Geometry geometry.coord_sys = kw.get('coord_sys', 0) # 0: Cartesian geometry.is_periodic = '%d %d %d'%(self.bcxmin=='periodic', self.bcymin=='periodic', self.bczmin=='periodic') # Is periodic? - geometry.prob_lo = '%7.0e %7.0e %7.0e'%(self.xmin, self.ymin, self.zmin) # physical domain - geometry.prob_hi = '%7.0e %7.0e %7.0e'%(self.xmax, self.ymax, self.zmax) + geometry.prob_lo = '%s %s %s'%(repr(self.xmin), repr(self.ymin), repr(self.zmin)) # physical domain + geometry.prob_hi = '%s %s %s'%(repr(self.xmax), repr(self.ymax), repr(self.zmax)) + + if self.moving_window_velocity is not None and np.any(self.moving_window_velocity != 0): + warpx.do_moving_window = 1 + if self.moving_window_velocity[0] != 0.: + warpx.moving_window_dir = 'x' + warpx.moving_window_v = self.moving_window_velocity[0]/clight # in units of the speed of light + if self.moving_window_velocity[1] != 0.: + warpx.moving_window_dir = 'y' + warpx.moving_window_v = self.moving_window_velocity[1]/clight # in units of the speed of light + if self.moving_window_velocity[2] != 0.: + warpx.moving_window_dir = 'z' + warpx.moving_window_v = self.moving_window_velocity[2]/clight # in units of the speed of light + + def getmins(self, **kw): + return np.array([warpx.getProbLo(0), warpx.getProbLo(1), warpx.getProbLo(2)]) + + def getmaxs(self, **kw): + return np.array([warpx.getProbHi(0), warpx.getProbHi(1), warpx.getProbHi(2)]) def getxmin(self): return warpx.getProbLo(0) diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py index e552a2f58..9dcc879c5 100755 --- a/Python/pywarpx/_libwarpx.py +++ b/Python/pywarpx/_libwarpx.py @@ -107,6 +107,8 @@ f.argtypes = (ctypes.c_int, ctypes.c_int, ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"), ctypes.c_int) +libwarpx.warpx_getProbLo.restype = ctypes.c_double +libwarpx.warpx_getProbHi.restype = ctypes.c_double libwarpx.warpx_getistep.restype = ctypes.c_int libwarpx.warpx_gett_new.restype = ctypes.c_double libwarpx.warpx_getdt.restype = ctypes.c_double diff --git a/tests/plasma_acceleration/plasma_acceleration.py b/tests/plasma_acceleration/plasma_acceleration.py index 51d540017..5bd6bff16 100644 --- a/tests/plasma_acceleration/plasma_acceleration.py +++ b/tests/plasma_acceleration/plasma_acceleration.py @@ -209,6 +209,7 @@ algo.particle_pusher = 0 # Particles particles.nspecies = 2 +particles.species_names = "electrons protons" warpx.verbose = 1 warpx.cfl = 1.0 @@ -241,7 +242,7 @@ for i in range(1, max_step + 1): domain_min[direction] += num_shift*dx[direction] domain_max[direction] += num_shift*dx[direction] - warpx.evolve(i) + warpx.evolve(1) old_x = new_x new_x = warpx.getProbLo(direction) diff --git a/tests/plasma_acceleration/plasma_acceleration_PICMI.py b/tests/plasma_acceleration/plasma_acceleration_PICMI.py new file mode 100644 index 000000000..e309773fc --- /dev/null +++ b/tests/plasma_acceleration/plasma_acceleration_PICMI.py @@ -0,0 +1,229 @@ +import numpy as np +from pywarpx import PICMI +#from warp import PICMI +if PICMI.codename == 'WarpX': + from pywarpx import PGroups + +from mpi4py import MPI +comm = MPI.COMM_WORLD + +def get_parallel_indices(Np, rank, size): + ''' + + This decomposes the arrays of particle attributes into subarrays + for parallel processing. + + ''' + Navg = Np / size + Nleft = Np - Navg * size + if (rank < Nleft): + lo = rank*(Navg + 1) + hi = lo + Navg + 1 + else: + lo = rank * Navg + Nleft + hi = lo + Navg + return lo, hi + + +def set_initial_conditions(ncells, domain_min, domain_max): + ''' + + Sets initial conditions for the plasma acceleration setup. + + + ''' + + comm.Barrier() + + num_ppc = 4 + c = 299792458.0 + dx = (domain_max - domain_min) / ncells + + # Species 0 - the beam + beam_min = np.array([ -20e-6, -20e-6, -150e-6]) + beam_max = np.array([ 20e-6, 20e-6, -100e-6]) + beam_density = 1.e23 + beam_gamma = 1.e9 + beam_uz = np.sqrt(beam_gamma**2 - 1)*c + beam_weight = beam_density * dx[0] * dx[1] * dx[2] / num_ppc + + # Species 1 - the plasma + plasma_min = np.array([-200e-6, -200e-6, 0.0e-6]) + plasma_max = np.array([ 200e-6, 200e-6, 200e-6]) + plasma_density = 1.e22 + plasma_weight = plasma_density * dx[0] * dx[1] * dx[2] / num_ppc + + # Fill the entire domain with particles. Only a subset of these points + # will be selected for each species + Z, Y, X, P = np.mgrid[0:ncells[0], 0:ncells[1], 0:ncells[2], 0:num_ppc] + Z = Z.flatten() + Y = Y.flatten() + X = X.flatten() + P = P.flatten() + + particle_shift = (0.5 + P) / num_ppc + + xp = (X + particle_shift)*dx[0] + domain_min[0] + yp = (Y + particle_shift)*dx[1] + domain_min[1] + zp = (Z + particle_shift)*dx[2] + domain_min[2] + + # Do the beam species first + beam_locs = np.logical_and(xp >= beam_min[0], xp < beam_max[0]) + beam_locs = np.logical_and(beam_locs, yp >= beam_min[1]) + beam_locs = np.logical_and(beam_locs, zp >= beam_min[2]) + beam_locs = np.logical_and(beam_locs, yp < beam_max[1]) + beam_locs = np.logical_and(beam_locs, zp < beam_max[2]) + + beam_xp = xp[beam_locs] + beam_yp = yp[beam_locs] + beam_zp = zp[beam_locs] + + beam_uxp = np.zeros_like(beam_xp) + beam_uyp = np.zeros_like(beam_xp) + beam_uzp = beam_uz * np.ones_like(beam_xp) + + # --- and weights + Np = beam_xp.shape[0] + beam_wp = np.full(Np, beam_weight) + + lo, hi = get_parallel_indices(Np, comm.rank, comm.size) + + beam.add_particles(x=beam_xp[lo:hi], y=beam_yp[lo:hi], z=beam_zp[lo:hi], + ux=beam_uxp[lo:hi], uy=beam_uyp[lo:hi], uz=beam_uzp[lo:hi], + w=beam_wp[lo:hi], unique_particles=1) + + # now do the plasma species + plasma_locs = np.logical_and(xp >= plasma_min[0], xp < plasma_max[0]) + plasma_locs = np.logical_and(plasma_locs, yp >= plasma_min[1]) + plasma_locs = np.logical_and(plasma_locs, zp >= plasma_min[2]) + plasma_locs = np.logical_and(plasma_locs, yp < plasma_max[1]) + plasma_locs = np.logical_and(plasma_locs, zp < plasma_max[2]) + + plasma_xp = xp[plasma_locs] + plasma_yp = yp[plasma_locs] + plasma_zp = zp[plasma_locs] + + plasma_uxp = np.zeros_like(plasma_xp) + plasma_uyp = np.zeros_like(plasma_xp) + plasma_uzp = np.zeros_like(plasma_xp) + + # --- and weights + Np = plasma_xp.shape[0] + plasma_wp = np.full(Np, plasma_weight) + + lo, hi = get_parallel_indices(Np, comm.rank, comm.size) + + plasma.add_particles(x=plasma_xp[lo:hi], y=plasma_yp[lo:hi], z=plasma_zp[lo:hi], + ux=plasma_uxp[lo:hi], uy=plasma_uyp[lo:hi], uz=plasma_uzp[lo:hi], + w=plasma_wp[lo:hi], unique_particles=1) + + comm.Barrier() + + +def inject_plasma(num_shift, direction): + ''' + + This injects fresh plasma into the domain after the moving window has been upated. + + ''' + + comm.Barrier() + + num_ppc = 4 + density = 1.e22 + weight = density * dx[0]*dx[1]*dx[2] / num_ppc + + shift_box = np.array(ncells) + shift_box[direction] = abs(num_shift) + + Z, Y, X, P = np.mgrid[0:shift_box[2], 0:shift_box[1], 0:shift_box[0], 0:num_ppc] + X = X.flatten() + Y = Y.flatten() + Z = Z.flatten() + P = P.flatten() + particle_shift = (0.5 + P) / num_ppc + + pos = [[],[],[]] + + if (num_shift > 0): + pos[0] = (X + particle_shift)*dx[0] + domain_min[0] + pos[1] = (Y + particle_shift)*dx[1] + domain_min[1] + pos[2] = (Z + particle_shift)*dx[2] + domain_min[2] + pos[direction] -= domain_min[direction] + pos[direction] += domain_max[direction] + + if (num_shift < 0): + pos[0] = (X + particle_shift)*dx[0] + domain_min[0] + pos[1] = (Y + particle_shift)*dx[1] + domain_min[1] + pos[2] = (Z + particle_shift)*dx[2] + domain_min[2] + pos[direction] += num_shift*dx[direction] + + xp = pos[0] + yp = pos[1] + zp = pos[2] + + uxp = np.zeros_like(xp) + uyp = np.zeros_like(xp) + uzp = np.zeros_like(xp) + + Np = xp.shape[0] + wp = np.full(Np, weight) + + lo, hi = get_parallel_indices(Np, comm.rank, comm.size) + + plasma.add_particles(x=xp[lo:hi], y=yp[lo:hi], z=zp[lo:hi], + ux=uxp[lo:hi], uy=uyp[lo:hi], uz=uzp[lo:hi], + w=wp[lo:hi], unique_particles=1) + + comm.Barrier() + +ncells = np.array([64, 64, 64]) +domain_min = np.array([-200e-6, -200e-6, -200e-6]) +domain_max = np.array([ 200e-6, 200e-6, 200e-6]) +dx = (domain_max - domain_min) / ncells + +moving_window_velocity = np.array([0., 0., PICMI.clight]) + +grid = PICMI.Grid(nx=ncells[0], ny=ncells[1], nz=ncells[2], + xmin=domain_min[0], xmax=domain_max[0], ymin=domain_min[1], ymax=domain_max[1], zmin=domain_min[2], zmax=domain_max[2], + bcxmin='periodic', bcxmax='periodic', bcymin='periodic', bcymax='periodic', bczmin='open', bczmax='open', + moving_window_velocity = moving_window_velocity, + max_grid_size=32, max_level=0, coord_sys=0) + + +solver = PICMI.EM_solver(current_deposition = 3, + charge_deposition = 0, + field_gathering = 0, + particle_pusher = 0) + +beam = PICMI.Species(type=PICMI.Electron, name='beam') +plasma = PICMI.Species(type=PICMI.Electron, name='plasma') + + +# Maximum number of time steps +max_step = 160 + +sim = PICMI.Simulation(plot_int = 2, + verbose = 1, + cfl = 1.0) + +set_initial_conditions(ncells, domain_min, domain_max) + +direction = np.argmax(abs(moving_window_velocity)) +old_x = grid.getmins()[direction] +new_x = old_x +for i in range(1, max_step + 1): + + # check whether the moving window has updated + num_shift = int( 0.5 + (new_x - old_x) / dx[direction]) + if (num_shift != 0): + inject_plasma(num_shift, direction) + domain_min[direction] += num_shift*dx[direction] + domain_max[direction] += num_shift*dx[direction] + + sim.step(1) + + old_x = new_x + new_x = grid.getmins()[direction] + +sim.finalize() |