From 6dcd85dd5a1adefc90fc186bba5768727b842a6d Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Mon, 26 Jun 2017 14:20:04 -0700 Subject: Updated to be compliant with PICMI Python standard --- Python/pywarpx/PICMI.py | 213 ++++++++++++------------------------------------ 1 file changed, 54 insertions(+), 159 deletions(-) (limited to 'Python/pywarpx/PICMI.py') diff --git a/Python/pywarpx/PICMI.py b/Python/pywarpx/PICMI.py index 3e7e7185f..95bb758dc 100644 --- a/Python/pywarpx/PICMI.py +++ b/Python/pywarpx/PICMI.py @@ -1,74 +1,14 @@ """Classes following the PICMI standard """ +from PICMI_Base import * import numpy as np from pywarpx import * -pi = 3.14159265358979323 # ratio of a circle's circumference to its diameter -euler = 0.57721566490153386 # Euler-Masceroni constant. Base of the natural logarithm. -amu = 1.660538921e-27 # Atomic Mass Unit [kg] -clight = 2.99792458e+8 # Speed of light in vacuum (exact) [m/s] -echarge = 1.602176565e-19 # Elementary charge [C] -emass = 9.10938291e-31 # Electron mass [kg] -mu0 = 4.e-7*pi # Permeability of free space [kg.m/(s.s.A.A)=H/m=T.m/A] -eps0 = 1./(mu0*clight*clight) # Permittivity of free space [F/m] -boltzmann = 1.3806488e-23 # Boltzmann's constant [J/K] -avogadro = 6.02214129e23 # Avogadro's Number [atoms/mole] -planck = 6.62606957e-34 # Planck's constant [J.s] - -class Grid(object): - """ - - `Grid` - - **type**: *object* - - `Nx=nx` - **type**: *integer* - "Number of cells along X (Nb nodes=nx+1)." - - `Ny=ny` - **type**: *integer* - "Number of cells along Y (Nb nodes=ny+1)." - - `Nr=nr` - **type**: *integer* - "Number of cells along R (Nb nodes=nr+1)." - - `Nz=nz` - **type**: *integer* - "Number of cells along Z (Nb nodes=nz+1)." - - `Nm=nm` - **type**: *integer* - "Number of azimuthal modes." - - `Xmin=xmin` - **type**: *double* - "Position of first node along X." - - `Xmax=xmax` - **type**: *double* - "Position of last node along X." - - `Ymin=ymin` - **type**: *double* - "Position of first node along Y." - - `Ymax=ymax` - **type**: *double* - "Position of last node along Y." - - `Rmax=rmax` - **type**: *double* - "Position of last node along R." - - `Zmin=zmin` - **type**: *double* - "Position of first node along Z." - - `Zmax=zmax` - **type**: *double* - "Position of last node along Z." - - `bcxmin` - **type**: *string* - "Boundary condition at min X: periodic/open/dirichlet/neumann." - - `bcxmax` - **type**: *string* - "Boundary condition at max X: periodic/open/dirichlet/neumann." - - `bcymin` - **type**: *string* - "Boundary condition at min Y: periodic/open/dirichlet/neumann." - - `bcymax` - **type**: *string* - "Boundary condition at max Y: periodic/open/dirichlet/neumann." - - `bcrmax` - **type**: *string* - "Boundary condition at max R: open/dirichlet/neumann." - - `bczmin` - **type**: *string* - "Boundary condition at min Z: periodic/open/dirichlet/neumann." - - `bczmax` - **type**: *string* - "Boundary condition at max Z: periodic/open/dirichlet/neumann." - - - max_grid_size - - max_level - - coord_sys - """ - - def __init__(self, nx=None, ny=None, nr=None, nz=None, nm=None, - xmin=None, xmax=None, ymin=None, ymax=None, rmax=None, zmin=None, zmax=None, - bcxmin=None, bcxmax=None, bcymin=None, bcymax=None, bcrmax=None, bczmin=None, bczmax=None, - **kw): - self.nx = nx - self.ny = ny - self.nr = nr - self.nz = nz - self.nm = nm - self.xmin = xmin - self.xmax = xmax - self.ymin = ymin - self.ymax = ymax - self.rmax = rmax - self.zmin = zmin - self.zmax = zmax - self.bcxmin = bcxmin - self.bcxmax = bcxmax - self.bcymin = bcymin - self.bcymax = bcymax - self.bcrmax = bcrmax - self.bczmin = bczmin - self.bczmax = bczmax - - amr.n_cell = '%d %d %d'%(nx, ny, nz) + +class Grid(PICMI_Grid): + def init(self, **kw): + + amr.n_cell = '%d %d %d'%(self.nx, self.ny, self.nz) # Maximum allowable size of each subdomain in the problem domain; # this is used to decompose the domain for parallel calculations. @@ -79,97 +19,45 @@ class Grid(object): # Geometry geometry.coord_sys = kw.get('coord_sys', 0) # 0: Cartesian - geometry.is_periodic = '%d %d %d'%(bcxmin=='periodic', bcymin=='periodic', bczmin=='periodic') # Is periodic? - geometry.prob_lo = '%7.0e %7.0e %7.0e'%(xmin, ymin, zmin) # physical domain - geometry.prob_hi = '%7.0e %7.0e %7.0e'%(xmax, ymax, zmax) - - -class EM_solver(object): - Methods_list = ['Yee', 'CK', 'CKC', 'Lehe', 'PSTD', 'PSATD', 'GPSTD'] - def __init__(self, Method=None, - norderx=None, nordery=None, norderr=None, norderz=None, - l_nodal=None, - current_deposition_algo=None, charge_deposition_algo=None, - field_gathering_algo=None, particle_pusher_algo=None, **kw): - - assert Method is None or Method in EM_solver.Methods_list, Exception('Method has incorrect value') - - if current_deposition_algo is not None: - algo.current_deposition = current_deposition_algo - if charge_deposition_algo is not None: - algo.charge_deposition = charge_deposition_algo - if field_gathering_algo is not None: - algo.field_gathering = field_gathering_algo - if particle_pusher_algo is not None: - algo.particle_pusher = particle_pusher_algo - - -class Particle(object): - def __init__(self, Charge=None, charge=None, Q=None, q=None, - Mass=None, mass=None, M=None, m=None, - Symbol=None, symbol=None, S=None, s=None, - Name=None, name=None, N=None, n=None, **kw): - # --- Accept multpiple names, but use 'charge', 'mass', 'symbol', 'name' internally. - if Charge is not None: charge = Charge - if Q is not None: charge = Q - if q is not None: charge = q - if Mass is not None: mass = Mass - if M is not None: mass = M - if m is not None: mass = m - if Symbol is not None: symbol = Symbol - if S is not None: symbol = S - if s is not None: symbol = s - if Name is not None: name = Name - if N is not None: name = N - if n is not None: name = n - self.charge = charge - self.mass = mass - self.symbol = symbol - -Electron = Particle(q=-echarge, m=emass, symbol='e-', name='Electron') -Positron = Particle(q=echarge, m=emass, symbol='e+', name='Positron') -Proton = Particle(q=echarge, m=1.6726231e-27, symbol='p', name='Proton') -AntiProton = Particle(q=-echarge, m=1.6726231e-27, symbol='p-', name='Antiproton') -Neutron = Particle(q=0. , m=1.6749286e-27, symbol='n', name='Neutron') -Muon = Particle(q=-echarge, m=1.883531475e-28, symbol='mu-', name='Muon') -Antimuon = Particle(q=echarge, m=1.883531475e-28, symbol='mu+', name='Antimuon') -Photon = Particle(q=0., m=0., symbol='gnu', name='Photon') - - -class Species(object): - def __init__(self, - Type=None, type=None, - Name=None, name=None, - Sid=None, sid=None, - Charge_state=None, charge_state=None, - Charge=None, charge=None, Q=None, q=None, - Mass=None, mass=None, M=None, m=None, - Weight=None, weight=None, W=None, w=None, **kw): - # --- Accept multpiple names, but use 'type', 'name', 'sid', 'charge_state', 'charge', 'mass', 'weight' - if Type is not None: type = Type - if Name is not None: name = Name - if Sid is not None: sid = Sid - if Charge_state is not None: charge_state = Charge_state - if Charge is not None: charge = Charge - if Q is not None: charge = Q - if q is not None: charge = q - if Mass is not None: mass = Mass - if M is not None: mass = M - if m is not None: mass = m - if Weight is not None: weight = Weight - if W is not None: weight = W - if w is not None: weight = w - self.type = type - self.name = name - self.sid = sid - self.charg_state = charg_state - self.charge = charge - self.mass = mass - self.weight = weight + 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) + + def getxmin(self): + return warpx.getProbLo(0) + def getxmax(self): + return warpx.getProbHi(0) + def getymin(self): + return warpx.getProbLo(1) + def getymax(self): + return warpx.getProbHi(1) + def getzmin(self): + return warpx.getProbLo(2) + def getzmax(self): + return warpx.getProbHi(2) + +class EM_solver(PICMI_EM_solver): + def init(self, **kw): + + if self.current_deposition_algo is not None: + algo.current_deposition = self.current_deposition_algo + if self.charge_deposition_algo is not None: + algo.charge_deposition = self.charge_deposition_algo + if self.field_gathering_algo is not None: + algo.field_gathering = self.field_gathering_algo + if self.particle_pusher_algo is not None: + algo.particle_pusher = self.particle_pusher_algo + + +class Species(PICMI_Species): + def init(self, **kw): self.species_number = particles.nspecies particles.nspecies = particles.nspecies + 1 - particles.species_names = particles.species_names + ' ' + name + particles.species_names = particles.species_names + ' ' + self.name + + self.bucket = Bucket.Bucket(self.name, mass=self.mass, charge=self.charge, injection_style = 'python') + Particles.particles_list.append(self.bucket) def add_particles(self, n=None, x=None, y=None, z=None, @@ -179,11 +67,18 @@ class Species(object): add_particles(self.species_number, x, y, z, ux, uy, uz, pid, unique_particles) -class Simulation(object): - def __init__(self, plot_int=None, verbose=None, cfl=None): - amr.plot_int = plot_int - warpx.verbose = verbose - warpx.cfl = cfl +class Simulation(PICMI_Simulation): + def set_warpx_attr(self, warpx_obj, attr, kw): + value = kw.get(attr, None) + if value is not None: + setattr(warpx_obj, attr, value) + setattr(self, attr, value) + + def init(self, **kw): + + warpx.verbose = self.verbose + warpx.cfl = self.cfl + amr.plot_int = self.plot_int self.amrex = AMReX() self.amrex.init() -- cgit v1.2.3 From f6dd7a697a15ff568bfd33abed62066bbd2639c8 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Wed, 28 Jun 2017 16:14:16 -0700 Subject: Various fixes to conform to Python PICMI standard --- Python/pywarpx/AMReX.py | 7 + Python/pywarpx/Bucket.py | 4 +- Python/pywarpx/PGroup.py | 71 ++++--- Python/pywarpx/PICMI.py | 24 ++- Python/pywarpx/_libwarpx.py | 2 + tests/plasma_acceleration/plasma_acceleration.py | 3 +- .../plasma_acceleration_PICMI.py | 229 +++++++++++++++++++++ 7 files changed, 303 insertions(+), 37 deletions(-) create mode 100644 tests/plasma_acceleration/plasma_acceleration_PICMI.py (limited to 'Python/pywarpx/PICMI.py') 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() -- cgit v1.2.3