aboutsummaryrefslogtreecommitdiff
path: root/Python/pywarpx
diff options
context:
space:
mode:
Diffstat (limited to 'Python/pywarpx')
-rw-r--r--Python/pywarpx/AMReX.py (renamed from Python/pywarpx/BoxLib.py)6
-rw-r--r--Python/pywarpx/PGroup.py202
-rw-r--r--Python/pywarpx/__init__.py5
-rw-r--r--Python/pywarpx/timestepper.py385
4 files changed, 594 insertions, 4 deletions
diff --git a/Python/pywarpx/BoxLib.py b/Python/pywarpx/AMReX.py
index a15119250..444df3490 100644
--- a/Python/pywarpx/BoxLib.py
+++ b/Python/pywarpx/AMReX.py
@@ -10,7 +10,7 @@ from .Particles import particles
from . import warpxC
-class BoxLib(object):
+class AMReX(object):
def init(self):
argv = []
@@ -22,7 +22,7 @@ class BoxLib(object):
argv += interpolation.attrlist()
argv += particles.attrlist()
- warpxC.boxlib_init(argv)
+ warpxC.amrex_init(argv)
def finalize(self, finalize_mpi=1):
- warpxC.boxlib_finalize(finalize_mpi)
+ warpxC.amrex_finalize(finalize_mpi)
diff --git a/Python/pywarpx/PGroup.py b/Python/pywarpx/PGroup.py
new file mode 100644
index 000000000..616141702
--- /dev/null
+++ b/Python/pywarpx/PGroup.py
@@ -0,0 +1,202 @@
+import numpy as np
+from . import WarpX
+from . import warpxC
+
+class PGroup(object):
+ """Implements a class that has the same API as a warp ParticleGroup instance.
+ """
+
+ def __init__(self):
+ self.ns = 1 # Number of species
+ self.npmax = 0 # Size of data arrays
+ self.npid = 0 # number of columns for pid.
+
+ self.gallot()
+
+ def name(self):
+ return 'WarpXParticleGroup'
+
+ def gallot(self):
+ self.lebcancel_pusher = 0 # turns on/off cancellation of E+VxB within V push (logical type)
+ self.lebcancel = 0 # turns on/off cancellation of E+VxB before V push (logical type)
+
+ self.sm = np.zeros(self.ns) # Species mass [kg]
+ self.sq = np.zeros(self.ns) # Species charge [C]
+ self.sw = np.zeros(self.ns) # Species weight, (real particles per simulation particles)
+
+ self.ins = np.ones(self.ns, dtype=int) # Index of first particle in species
+ self.nps = np.zeros(self.ns, dtype=int) # Number of particles in species
+ self.ipmax = np.zeros(self.ns+1, dtype=int) # Max extent within the arrays of each species
+
+ self.sid = np.arange(self.ns, dtype=int) # Global species index for each species
+ self.ndts = np.ones(self.ns, dtype=int) # Stride for time step advance for each species
+ self.ldts = np.ones(self.ns, dtype=int) # (logical type)
+ self.lvdts = np.ones(self.ns, dtype=int) # (logical type)
+ self.iselfb = np.zeros(self.ns, dtype=int) # Group number for particles that are affected by
+ # their own magnetic field, using the 1/gamma**2
+ # approximation. The correction is not applied to
+ # group number -1.
+ self.fselfb = np.zeros(self.ns) # The scaling factor, vz.
+ self.l_maps = np.zeros(self.ns)
+ self.dtscale = np.ones(self.ns) # Scale factor applied to time step size for each
+ # species. Only makes sense in steaday and and
+ # transverse slice modes.
+
+ self.limplicit = np.zeros(self.ns, dtype=int) # Flags implicit particle species (logical type)
+ self.iimplicit = np.full(self.ns, -1, dtype=int) # Group number for implicit particles
+ self.ldoadvance = np.ones(self.ns, dtype=int) # Flags whether particles are time advanced (logical type)
+ self.lboundaries = np.ones(self.ns, dtype=int) # Flags whether boundary conditions need to be applied (logical type)
+ self.lparaxial = np.zeros(self.ns, dtype=int) # Flags to turn on/off paraxial approximation (logical type)
+
+ self.zshift = np.zeros(self.ns)
+ self.gamma_ebcancel_max = np.ones(self.ns) # maximum value allowed for ExB cancellation
+
+ self.gaminv = np.ones(self.npmax) # inverse relativistic gamma factor
+ self._xp = np.zeros(self.npmax) # X-positions of particles [m]
+ self._yp = np.zeros(self.npmax) # Y-positions of particles [m]
+ self._zp = np.zeros(self.npmax) # Z-positions of particles [m]
+ self._uxp = np.zeros(self.npmax) # gamma * X-velocities of particles [m/s]
+ self._uyp = np.zeros(self.npmax) # gamma * Y-velocities of particles [m/s]
+ self._uzp = np.zeros(self.npmax) # gamma * Z-velocities of particles [m/s]
+ self._ex = np.zeros(self.npmax) # Ex of particles [V/m]
+ self._ey = np.zeros(self.npmax) # Ey of particles [V/m]
+ self._ez = np.zeros(self.npmax) # Ez of particles [V/m]
+ self._bx = np.zeros(self.npmax) # Bx of particles [T]
+ self._by = np.zeros(self.npmax) # By of particles [T]
+ self._bz = np.zeros(self.npmax) # Bz of particles [T]
+ self._pid = np.zeros((self.npmax, self.npid)) # Particle ID - used for various purposes
+
+ # --- Temporary fix
+ gchange = gallot
+
+ def allocated(self, name):
+ return (getattr(self, name, None) is not None)
+
+ def addspecies(self):
+ pass
+
+ def _updatelocations(self):
+ warpx = WarpX.warpx.warpx
+ mypc = warpx.GetPartContainer()
+
+ xplist = []
+ yplist = []
+ zplist = []
+ for ispecie in range(mypc.nSpecies()):
+ pc = mypc.GetParticleContainer(ispecie)
+ xx = pc.getLocations()
+ xplist.append(xx[0,:])
+ yplist.append(xx[1,:])
+ zplist.append(xx[2,:])
+ self.nps[ispecie] = len(xplist[-1])
+ if ispecie > 0:
+ self.ins[ispecie] = self.ins[ispecie-1] + self.nps[ispecie-1]
+ self.ipmax[ispecie+1] = self.ins[ispecie] + self.nps[ispecie] - 1
+
+ self._xp = np.concatenate(xplist)
+ self._yp = np.concatenate(yplist)
+ self._zp = np.concatenate(zplist)
+ self.npmax = len(self._xp)
+
+ def _updatevelocities(self):
+ warpx = WarpX.warpx.warpx
+ mypc = warpx.GetPartContainer()
+
+ uxplist = []
+ uyplist = []
+ uzplist = []
+ for ispecie in range(mypc.nSpecies()):
+ pc = mypc.GetParticleContainer(ispecie)
+ vv = pc.getData(0, 3)
+ uxplist.append(vv[0,:])
+ uyplist.append(vv[1,:])
+ uzplist.append(vv[2,:])
+ self.nps[ispecie] = len(uxplist[-1])
+ if ispecie > 0:
+ self.ins[ispecie] = self.ins[ispecie-1] + self.nps[ispecie-1]
+ self.ipmax[ispecie+1] = self.ins[ispecie] + self.nps[ispecie] - 1
+
+ self._uxp = np.concatenate(uxplist)
+ self._uyp = np.concatenate(uyplist)
+ self._uzp = np.concatenate(uzplist)
+ self.npmax = len(self._xp)
+
+ def _updatepids(self):
+ warpx = WarpX.warpx.warpx
+ mypc = warpx.GetPartContainer()
+
+ pidlist = []
+ for ispecie in range(mypc.nSpecies()):
+ pc = mypc.GetParticleContainer(ispecie)
+ self.npid = pc.nAttribs - 3
+ vv = pc.getData(3, self.npid)
+ pidlist.append(vv)
+ self.nps[ispecie] = len(uxplist[-1])
+ if ispecie > 0:
+ self.ins[ispecie] = self.ins[ispecie-1] + self.nps[ispecie-1]
+ self.ipmax[ispecie+1] = self.ins[ispecie] + self.nps[ispecie] - 1
+
+ self._pid = np.concatenate(pidlist.T, axis=0)
+ self.npmax = self._pid.shape[0]
+
+ def getxp(self):
+ self._updatelocations()
+ return self._xp
+ xp = property(getxp)
+
+ def getyp(self):
+ self._updatelocations()
+ return self._yp
+ yp = property(getyp)
+
+ def getzp(self):
+ self._updatelocations()
+ return self._zp
+ zp = property(getzp)
+
+ def getuxp(self):
+ self._updatevelocities()
+ return self._uxp
+ uxp = property(getuxp)
+
+ def getuyp(self):
+ self._updatevelocities()
+ return self._uyp
+ uyp = property(getuyp)
+
+ def getuzp(self):
+ self._updatevelocities()
+ return self._uzp
+ uzp = property(getuzp)
+
+ def getpid(self):
+ self._updatepids()
+ return self._pid
+ pid = property(getpid)
+
+ def getgaminv(self):
+ uxp = self.uxp
+ uyp = self.uyp
+ uzp = self.uzp
+ return sqrt(1. - (uxp**2 + uyp**2 + uzp**2)/warpxC.c**2)
+ gaminv = property(getgaminv)
+
+ def getex(self):
+ return np.zeros(self.npmax)
+ ex = property(getex)
+ def getey(self):
+ return np.zeros(self.npmax)
+ ey = property(getey)
+ def getez(self):
+ return np.zeros(self.npmax)
+ ez = property(getez)
+ def getbx(self):
+ return np.zeros(self.npmax)
+ bx = property(getbx)
+ def getby(self):
+ return np.zeros(self.npmax)
+ by = property(getby)
+ def getbz(self):
+ return np.zeros(self.npmax)
+ bz = property(getbz)
+
diff --git a/Python/pywarpx/__init__.py b/Python/pywarpx/__init__.py
index 98c26af36..4a509419e 100644
--- a/Python/pywarpx/__init__.py
+++ b/Python/pywarpx/__init__.py
@@ -6,5 +6,8 @@ from .Algo import algo
from .Langmuirwave import langmuirwave
from .Interpolation import interpolation
from .Particles import particles
-from .BoxLib import BoxLib
+from .AMReX import AMReX
+
+from .timestepper import TimeStepper
+from .PGroup import PGroup
diff --git a/Python/pywarpx/timestepper.py b/Python/pywarpx/timestepper.py
new file mode 100644
index 000000000..71787c7a4
--- /dev/null
+++ b/Python/pywarpx/timestepper.py
@@ -0,0 +1,385 @@
+import warp
+from warp import top
+from warp import w3d
+
+from . import WarpX
+
+class TimeStepper(object):
+
+
+ def onestep(self):
+ # --- A reference to the instance of the WarpX class
+ warpx = WarpX.warpx.warpx
+ mypc = warpx.GetPartContainer()
+
+ self.cur_time = warpx.gett_new(0)
+ self.istep = warpx.getistep(0)
+
+ #if mpi.rank == 0:
+ print "\nSTEP %d starts ..."%(self.istep + 1)
+
+ #if (ParallelDescriptor::NProcs() > 1)
+ # if (okToRegrid(step)) RegridBaseLevel();
+
+ warpx.ComputeDt()
+ dt = warpx.getdt(0)
+
+ # --- Advance level 0 by dt
+ lev = 0
+
+ # --- At the beginning, we have B^{n-1/2} and E^{n}.
+ # --- Particles have p^{n-1/2} and x^{n}.
+ warpx.EvolveB(lev, 0.5*dt) # We now B^{n}
+
+ warpx.FillBoundaryE(lev, False)
+ warpx.FillBoundaryB(lev, False)
+
+ # --- Evolve particles to p^{n+1/2} and x^{n+1}
+ # --- Depose current, j^{n+1/2}
+ warpx.PushParticlesandDepose(lev, self.cur_time)
+
+ mypc.Redistribute() # Redistribute particles
+
+ warpx.EvolveB(lev, 0.5*dt) # We now B^{n+1/2}
+
+ warpx.FillBoundaryB(lev, True)
+
+ warpx.EvolveE(lev, dt) # We now have E^{n+1}
+
+ self.istep += 1
+
+ self.cur_time += dt
+
+ warpx.MoveWindow();
+
+ #if mpi.rank == 0:
+ print "STEP %d ends. TIME = %e DT = %e"%(self.istep, self.cur_time, dt)
+
+ # --- Sync up time
+ for i in range(warpx.finestLevel()+1):
+ warpx.sett_new(i, self.cur_time)
+ warpx.setistep(i, self.istep)
+
+ max_time_reached = ((self.cur_time >= warpx.stopTime() - 1.e-6*dt) or (self.istep >= warpx.maxStep()))
+
+ if warpx.plotInt() > 0 and (self.istep+1)%warpx.plotInt() == 0 or max_time_reached:
+ warpx.WritePlotFile()
+
+ if warpx.checkInt() > 0 and (self.istep+1)%warpx.plotInt() == 0 or max_time_reached:
+ warpx.WriteCheckPointFile()
+
+
+
+
+
+
+
+
+
+
+
+
+
+class TimeStepperFromPICSAR(object):
+
+ def __init__(self, package=None, solver=None, l_debug=False):
+ self.package = package
+ self.solver = solver
+ self.l_debug = l_debug
+
+ def setpackagestepnumber(self, it):
+ if self.package is not None:
+ self.package.setstepnumber(it)
+
+ def step(self, n=1, freq_print=10, lallspecl=0):
+ """
+ This function performs a range of Particle-In-Cell iterations
+
+ Inputs:
+ - n: number of iterations
+ - freq_print: print frequency
+ """
+
+ if (self.l_debug): print("Call step")
+
+ for i in range(n):
+ if(me == 0):
+ if top.it%freq_print == 0:
+ print 'it = %g time = %g'%(top.it, top.time)
+
+ l_first = (lallspecl or (i == 0))
+ l_last = (lallspecl or (i == n-1))
+
+ self.onestep(l_first, l_last)
+
+ if (self.l_debug): print("End step")
+
+ def onestep(self, l_first, l_last):
+ """
+ Perform a single particle-in-cell step
+ """
+
+ if (self.l_debug): print("Call onestep")
+
+ # --- Iteration number
+ self.setpackagestepnumber(top.it)
+
+ # --- call beforestep functions
+ if (self.l_debug): print("Call beforestep functions")
+ warp.callbeforestepfuncs.callfuncsinlist()
+
+ # --- gather fields from grid to particles
+ if (self.l_debug): print("Call Field gathering and particle push")
+
+ # --- push
+ if l_first:
+ if self.package is None:
+ # --- Standard Warp advance
+ for specie in warp.listofallspecies:
+ for pg in specie.flatten(specie.pgroups):
+ for js in range(pg.ns):
+ self.push_velocity_second_part(js, pg)
+ self.push_positions(js, pg)
+ warp.particleboundaries3d(pg, -1, False)
+
+ else:
+ # --- Particle pusher
+ if (self.l_debug): print("Call package particle push")
+ pxr.pxrpush_particles_part2()
+
+ # --- Particle boundary consitions
+ if (self.l_debug): print("Call package particle boundary conditions")
+ pxr.particle_bcs()
+
+ if (self.l_debug): print("Call aliasparticlearrays()")
+ self.aliasparticlearrays()
+
+ else:
+ if self.package is None:
+ # --- Standard Warp advance
+
+ for specie in warp.listofallspecies:
+ for pg in specie.flatten(specie.pgroups):
+ for js in range(pg.ns):
+ self.push_velocity_full(js, pg)
+ self.push_positions(js, pg)
+
+ warp.particleboundaries3d(pg, -1, False)
+
+ else:
+ # --- Particle pusher
+ if (self.l_debug): print("Call package particle pusher")
+ pxr.field_gathering_plus_particle_pusher()
+
+ # --- Particle boundary conditions
+ if (self.l_debug): print("Call package particle boundary conditions")
+ pxr.particle_bcs()
+
+ self.aliasparticlearrays()
+
+ # --- Particle sorting
+ if (self.l_debug): print("Call Particle Sorting")
+ if self.package is not None:
+ # --- This should be a function installed before load rho
+ pxr.particle_sorting_sub()
+
+ # --- call beforeloadrho functions
+ if (self.l_debug): print("Call beforeloadrho functions")
+ warp.beforeloadrho.callfuncsinlist()
+
+ pgroups = []
+ for specie in warp.listofallspecies:
+ pgroups += specie.flatten(specie.pgroups)
+
+ self.pgroups = pgroups
+
+ # --- Call user-defined injection routines
+ if (self.l_debug): print("Call user-defined injection routines")
+ warp.userinjection.callfuncsinlist()
+
+ xgriddiff = w3d.xmmin - pxr.xmin
+ ygriddiff = w3d.ymmin - pxr.ymin
+ zgriddiff = w3d.zmmin - pxr.zmin
+
+ if (xgriddiff != 0 or ygriddiff != 0 or zgriddiff != 0):
+ pxr.pxr_move_sim_boundaries(xgriddiff, ygriddiff, zgriddiff)
+
+
+ if (self.l_debug): print("Call loadrho")
+ self.solver.loadrho(pgroups = pgroups)
+ if (self.l_debug): print("Call loadj")
+ self.solver.loadj(pgroups = pgroups)
+
+ if (self.l_debug): print("Call dosolve")
+ self.solver.dosolve()
+
+ if self.package is None:
+ for specie in warp.listofallspecies:
+ for pg in specie.flatten(specie.pgroups):
+ for js in range(pg.ns):
+ self.fetcheb(js, pg)
+ if l_last:
+ self.push_velocity_first_part(js, pg)
+
+ else:
+ if l_last:
+ if (self.l_debug): print("Call package push particles 1")
+ pxr.pxrpush_particles_part1()
+
+ # --- update time, time counter
+ top.time += top.dt
+ if top.it%top.nhist == 0:
+# zmmnt()
+ minidiag(top.it, top.time, top.lspecial)
+ top.it += 1
+
+ # --- Load balance function should be installed after step
+
+ # --- call afterstep functions
+ if (self.l_debug): print("Call callafterstepfuncs.callfuncsinlist()")
+ warp.callafterstepfuncs.callfuncsinlist()
+
+
+ # --- The following methods use the standard Warp routines
+ def fetcheb(self, js, pg=None):
+ if self.l_verbose:print me, 'enter fetcheb'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ w3d.pgroupfsapi = pg
+ w3d.ipminfsapi = pg.ins[js]
+ w3d.npfsapi = pg.nps[js]
+ pg.ex[il:iu] = 0.
+ pg.ey[il:iu] = 0.
+ pg.ez[il:iu] = 0.
+ pg.bx[il:iu] = 0.
+ pg.by[il:iu] = 0.
+ pg.bz[il:iu] = 0.
+ self.fetche()
+ self.fetchb()
+ w3d.pgroupfsapi = None
+
+ def push_velocity_full(self, js, pg=None):
+ if self.l_verbose:print me, 'enter push_velocity_full'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ if pg.lebcancel_pusher:
+ warp.ebcancelpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], top.dt, 0)
+ else:
+ # --- push velocity from electric field (half step)
+ warp.epush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt)
+ # --- update gamma
+ self.set_gamma(js, pg)
+ # --- push velocity from magnetic field
+ warp.bpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], top.dt, top.ibpush)
+ # --- push velocity from electric field (half step)
+ warp.epush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt)
+ # --- update gamma
+ self.set_gamma(js, pg)
+
+ if self.l_verbose:print me, 'exit push_velocity_first_part'
+
+ def push_velocity_first_part(self, js, pg=None):
+ if self.l_verbose:print me, 'enter push_velocity_first_part'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ if pg.lebcancel_pusher:
+ warp.ebcancelpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], top.dt, 1)
+ else:
+ # --- push velocity from electric field (half step)
+ warp.epush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt)
+ # --- update gamma
+ self.set_gamma(js,pg)
+ # --- push velocity from magnetic field
+ warp.bpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt, top.ibpush)
+
+ if self.l_verbose:print me, 'exit push_velocity_first_part'
+
+ def push_velocity_second_part(self, js, pg=None):
+ if self.l_verbose:print me, 'enter push_velocity_second_part'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ if pg.lebcancel_pusher:
+ warp.ebcancelpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], top.dt, 2)
+ else:
+ # --- push velocity from magnetic field
+ warp.bpush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], pg.gaminv[il:iu],
+ pg.bx[il:iu], pg.by[il:iu], pg.bz[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt, top.ibpush)
+ # --- push velocity from electric field (half step)
+ warp.epush3d(np, pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu],
+ pg.ex[il:iu], pg.ey[il:iu], pg.ez[il:iu],
+ pg.sq[js], pg.sm[js], 0.5*top.dt)
+ # --- update gamma
+ self.set_gamma(js, pg)
+
+ if self.l_verbose:print me, 'exit push_velocity_second_part'
+
+ def set_gamma(self, js, pg=None):
+ if self.l_verbose:print me, 'enter set_gamma'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ # --- update gamma
+ warp.gammaadv(np, pg.gaminv[il:iu], pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu], top.gamadv, top.lrelativ)
+
+ if self.l_verbose:print me, 'exit set_gamma'
+
+ def push_positions(self, js, pg=None):
+ if self.l_verbose:print me, 'enter push_positions'
+ if pg is None:
+ pg = top.pgroup
+ np = pg.nps[js]
+ if np == 0: return
+ il = pg.ins[js] - 1
+ iu = il + pg.nps[js]
+ warp.xpush3d(np, pg.xp[il:iu], pg.yp[il:iu], pg.zp[il:iu],
+ pg.uxp[il:iu], pg.uyp[il:iu], pg.uzp[il:iu],
+ pg.gaminv[il:iu], top.dt)
+
+ if self.l_verbose:print me, 'exit push_positions'
+
+
+
+
+
+
+
+