diff options
Diffstat (limited to 'Example')
| -rw-r--r-- | Example/Langmuir/langmuir.py | 78 | ||||
| -rw-r--r-- | Example/Langmuir/langmuir_PICMI.py | 46 | ||||
| -rw-r--r-- | Example/gaussian_beam/gaussian_beam_PICMI.py | 59 | ||||
| -rw-r--r-- | Example/laser_acceleration/laser_acceleration_PICMI.py | 221 | ||||
| -rw-r--r-- | Example/plasma_acceleration/plasma_acceleration_PICMI.py | 62 |
5 files changed, 310 insertions, 156 deletions
diff --git a/Example/Langmuir/langmuir.py b/Example/Langmuir/langmuir.py new file mode 100644 index 000000000..da9e3d8ad --- /dev/null +++ b/Example/Langmuir/langmuir.py @@ -0,0 +1,78 @@ +import numpy as np +from pywarpx import * + +nx = 64 +ny = 64 +nz = 64 + +xmin = -20.e-6 +ymin = -20.e-6 +zmin = -20.e-6 +xmax = +20.e-6 +ymax = +20.e-6 +zmax = +20.e-6 + +dx = (xmax - xmin)/nx +dy = (ymax - ymin)/ny +dz = (zmax - zmin)/nz + +# Maximum number of time steps +max_step = 40 + +# number of grid points +amr.n_cell = [nx, ny, nz] + +# Maximum allowable size of each subdomain in the problem domain; +# this is used to decompose the domain for parallel calculations. +amr.max_grid_size = 32 + +# Maximum level in hierarchy (for now must be 0, i.e., one level in total) +amr.max_level = 0 + +amr.plot_int = 1 # How often to write plotfiles. "<= 0" means no plotfiles. + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = [1, 1, 1] # Is periodic? +geometry.prob_lo = [xmin, ymin, zmin] # physical domain +geometry.prob_hi = [xmax, ymax, zmax] + +# Verbosity +warpx.verbose = 1 +warpx.do_moving_window = 0 +warpx.moving_window_dir = 'x' +warpx.moving_window_v = 0.0 # in units of the speed of light + +# Algorithms +algo.current_deposition = 3 +algo.charge_deposition = 0 +algo.field_gathering = 0 +algo.particle_pusher = 0 + +# CFL +warpx.cfl = 1.0 + +particles.nspecies = 1 +particles.species_names = 'electrons' + +electrons.charge = '-q_e' +electrons.mass = 'm_e' +electrons.injection_style = "NUniformPerCell" +electrons.num_particles_per_cell_each_dim = [2, 2, 2] + +electrons.xmin = xmin +electrons.xmax = 0.e-6 +electrons.ymin = ymin +electrons.ymax = ymax +electrons.zmin = zmin +electrons.zmax = zmax + +electrons.profile = 'constant' +electrons.density = 1.e25 # number of electrons per m^3 + +electrons.momentum_distribution_type = "constant" +electrons.ux = 0.01 # ux = gamma*beta_x + +# --- Initialize the simulation +warpx.write_inputs('inputs_from_python', max_step=max_step) + diff --git a/Example/Langmuir/langmuir_PICMI.py b/Example/Langmuir/langmuir_PICMI.py new file mode 100644 index 000000000..0a1200ffa --- /dev/null +++ b/Example/Langmuir/langmuir_PICMI.py @@ -0,0 +1,46 @@ +import numpy as np +from pywarpx import picmi + +nx = 64 +ny = 64 +nz = 64 + +xmin = -20.e-6 +ymin = -20.e-6 +zmin = -20.e-6 +xmax = +20.e-6 +ymax = +20.e-6 +zmax = +20.e-6 + +uniform_plasma = picmi.UniformDistribution(density=1.e25, upper_bound=[0., None, None], directed_velocity=[0.1, 0., 0.]) + +electrons = picmi.Species(particle_type='electron', name='electrons', initial_distribution=uniform_plasma) + +grid = picmi.Cartesian3DGrid(number_of_cells = [nx, ny, nz], + lower_bound = [xmin, ymin, zmin], + upper_bound = [xmax, ymax, zmax], + lower_boundary_conditions = ['periodic', 'periodic', 'periodic'], + upper_boundary_conditions = ['periodic', 'periodic', 'periodic'], + moving_window_velocity = [0., 0., 0.], + warpx_max_grid_size=32, warpx_max_level=0, warpx_coord_sys=0) + +solver = picmi.ElectromagneticSolver(grid=grid, cfl=1.) + +sim = picmi.Simulation(solver = solver, + verbose = 1, + max_steps = 40, + plot_int = 1, + warpx_current_deposition_algo = 3, + warpx_charge_deposition_algo = 0, + warpx_field_gathering_algo = 0, + warpx_particle_pusher_algo = 0) + +sim.add_species(electrons, layout=picmi.GriddedLayout(n_macroparticle_per_cell=[2,2,2], grid=grid)) + +# write_inputs will create an inputs file that can be used to run +# with the compiled version. +sim.write_input_file(file_name='inputs_from_PICMI') + +# Alternatively, sim.step will run WarpX, controlling it from Python +#sim.step() + diff --git a/Example/gaussian_beam/gaussian_beam_PICMI.py b/Example/gaussian_beam/gaussian_beam_PICMI.py new file mode 100644 index 000000000..ab8d22c4b --- /dev/null +++ b/Example/gaussian_beam/gaussian_beam_PICMI.py @@ -0,0 +1,59 @@ +import numpy as np +from pywarpx import PICMI +#from warp import PICMI + +nx = 32 +ny = 32 +nz = 32 + +xmin = -2. +xmax = +2. +ymin = -2. +ymax = +2. +zmin = -2. +zmax = +2. + +number_sim_particles = 32678 +total_charge = 8.010883097437485e-07 + +beam_rms_size = 0.25 +electron_beam_divergence = -0.04 + +em_order = 3 + +grid = PICMI.Grid(nx=nx, ny=ny, nz=nz, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + bcxmin='periodic', bcxmax='periodic', bcymin='periodic', bcymax='periodic', bczmin='open', bczmax='open', + max_grid_size=16, max_level=0, coord_sys=0) + +solver = PICMI.EM_solver(current_deposition_algo = 0, + charge_deposition_algo = 0, + field_gathering_algo = 0, + particle_pusher_algo = 0, + norderx = em_order, nordery = em_order, norderz = em_order) + +electrons = PICMI.Species(type='electron', name='electrons') +protons = PICMI.Species(type='electron', name='protons') + +electron_beam = PICMI.GaussianBeam(electrons, + number_sim_particles = number_sim_particles, + number_real_particles = total_charge/PICMI.q_e, + Xrms = beam_rms_size, Yrms = beam_rms_size, Zrms = beam_rms_size, + UXdiv = electron_beam_divergence, UYdiv = electron_beam_divergence, UZdiv = electron_beam_divergence) + +proton_beam = PICMI.GaussianBeam(protons, + number_sim_particles = number_sim_particles, + number_real_particles = total_charge/PICMI.q_e, + Xrms = beam_rms_size, Yrms = beam_rms_size, Zrms = beam_rms_size) + +sim = PICMI.Simulation(plot_int = 8, + verbose = 1, + cfl = 1.0, + max_step = 1000) + +# write_inputs will create an inputs file that can be used to run +# with the compiled version. +sim.write_inputs(inputs_name = 'inputs_from_PICMI') + +# Alternatively, sim.step will run WarpX, controlling it from Python +#sim.step() + diff --git a/Example/laser_acceleration/laser_acceleration_PICMI.py b/Example/laser_acceleration/laser_acceleration_PICMI.py index 761ff7cfe..4cb906378 100644 --- a/Example/laser_acceleration/laser_acceleration_PICMI.py +++ b/Example/laser_acceleration/laser_acceleration_PICMI.py @@ -1,168 +1,77 @@ import numpy as np -from pywarpx import PICMI -#from warp import PICMI +from pywarpx import picmi +#from warp import picmi -from mpi4py import MPI -comm = MPI.COMM_WORLD +nx = 64 +ny = 64 +nz = 480 -def get_parallel_indices(Np, rank, size): - ''' +xmin = -30.e-6 +xmax = +30.e-6 +ymin = -30.e-6 +ymax = +30.e-6 +zmin = -56.e-6 +zmax = +12.e-6 - 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 load_plasma(ncells, domain_min, domain_max, injected_plasma_density, injected_plasma_ppc, plasma_min, plasma_max): - ''' - - Sets initial conditions for the plasma acceleration setup. - - ''' - - comm.Barrier() - - dx = (domain_max - domain_min) / ncells - - nplasma_cells = ((plasma_max - plasma_min)/dx + 0.5).astype('l') + 1 - iplasma_min = ((plasma_min - domain_min)/dx).astype('l')*dx + domain_min - - # Species 0 - the plasma - plasma_weight = injected_plasma_density * dx[0] * dx[1] * dx[2] / injected_plasma_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:nplasma_cells[2], 0:nplasma_cells[1], 0:nplasma_cells[0], 0:injected_plasma_ppc] - Z = Z.flatten() - Y = Y.flatten() - X = X.flatten() - P = P.flatten() - - particle_shift = (0.5 + P) / injected_plasma_ppc - - xp = (X + particle_shift)*dx[0] + iplasma_min[0] - yp = (Y + particle_shift)*dx[1] + iplasma_min[1] - zp = (Z + particle_shift)*dx[2] + iplasma_min[2] - - # 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) - - electrons.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=True) - - comm.Barrier() - -def inject_new_plasma(ncells, domain_min, domain_max, num_shift, direction, injected_plasma_density, injected_plasma_ppc, plasma_min, plasma_max): - ''' - - This injects fresh plasma into the domain after the moving window has been upated. - - ''' - - comm.Barrier() - - dx = (domain_max - domain_min) / ncells - - plasma_min_inject = plasma_min.copy() - plasma_max_inject = plasma_max.copy() - if (num_shift > 0): - plasma_min_inject[direction] = domain_max[direction] - plasma_max_inject[direction] = domain_max[direction] + num_shift*dx[direction] - else: - plasma_min_inject[direction] = domain_min[direction] - num_shift*dx[direction] - plasma_max_inject[direction] = domain_min[direction] - - load_plasma(ncells, domain_min, domain_max, injected_plasma_density, injected_plasma_ppc, plasma_min_inject, plasma_max_inject) - -ncells = np.array([64, 64, 480])/2 -domain_min = np.array([-30.e-6, -30.e-6, -56.e-6]) -domain_max = np.array([ 30.e-6, 30.e-6, 12.e-6]) -dx = (domain_max - domain_min) / ncells - -moving_window_velocity = np.array([0., 0., PICMI.clight]) +moving_window_velocity = [0., 0., picmi.c] injected_plasma_density = 1.e23 -injected_plasma_ppc = 4 -plasma_min = np.array([-20.e-6, -20.e-6, 0.0e-6]) -plasma_max = np.array([ 20.e-6, 20.e-6, domain_max[2]]) - -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) - -laser = PICMI.Gaussian_laser(antenna_z0 = 9.e-6, # This point is on the laser plane - antenna_zvec = 1., # The plane normal direction - pol_angle = PICMI.pi/2., # The main polarization vector - E0 = 16.e12, # Maximum amplitude of the laser field (in V/m) - waist = 5.e-6, # The waist of the laser (in meters) - duration = 15.e-15, # The duration of the laser (in seconds) - t_peak = 30.e-15, # The time at which the laser reaches its peak (in seconds) - focal_position = 109.e-6, # Focal position (m) - wavelength = 0.8e-6, # The wavelength of the laser (in meters) - em_solver = solver) - -electrons = PICMI.Species(type=PICMI.Electron, name='electrons') - -# Maximum number of time steps -max_step = 1000 - -sim = PICMI.Simulation(plot_int = 200, +number_per_cell_each_dim = [2, 2, 1] +plasma_min = [-20.e-6, -20.e-6, 0.0e-6] +plasma_max = [ 20.e-6, 20.e-6, zmax] + +grid = picmi.Cartesian3DGrid(number_of_cells = [nx, ny, nz], + lower_bound = [xmin, ymin, zmin], + upper_bound = [xmax, ymax, zmax], + lower_boundary_conditions = ['periodic', 'periodic', 'open'], + upper_boundary_conditions = ['periodic', 'periodic', 'open'], + moving_window_velocity = moving_window_velocity, + warpx_max_grid_size=32, warpx_max_level=0, warpx_coord_sys=0) + +solver = picmi.ElectromagneticSolver(grid=grid, cfl=1.) + +t_peak = 30.e-15 # The time at which the laser reaches its peak at the antenna (in seconds) +focal_distance = 100.e-6 # Focal distance from the antenna (in meters) +antenna_z0 = 9.e-6 # This point is on the laser plane +laser = picmi.GaussianLaser(wavelength = 0.8e-6, # The wavelength of the laser (in meters) + waist = 5.e-6, # The waist of the laser (in meters) + duration = 15.e-15, # The duration of the laser (in seconds) + polarization_angle = np.pi/2., # The main polarization vector + focal_position = [0., 0., focal_distance + antenna_z0], # Focal position (m) + E0 = 16.e12, # Maximum amplitude of the laser field (in V/m) + centroid_position = [0., 0., antenna_z0 - picmi.c*t_peak], # Position of the laser centroid in Z at time 0 + propagation_direction = [0,0,1]) + +laser_antenna = picmi.LaserAntenna(position = [0., 0., antenna_z0], # This point is on the laser plane + normal_vector = [0., 0., 1.]) # The plane normal direction + +uniform_plasma = picmi.UniformDistribution(density = injected_plasma_density, + lower_bound = plasma_min, + upper_bound = plasma_max, + fill_in = True) + +electrons = picmi.Species(particle_type = 'electron', + name = 'electrons', + initial_distribution = uniform_plasma) + +sim = picmi.Simulation(solver = solver, + plot_int = 100, verbose = 1, - cfl = 1.0) - -load_plasma(ncells, domain_min, domain_max, injected_plasma_density, injected_plasma_ppc, plasma_min, plasma_max) + cfl = 1.0, + max_steps = 1000, + warpx_current_deposition_algo = 3, + warpx_charge_deposition_algo = 0, + warpx_field_gathering_algo = 0, + warpx_particle_pusher_algo = 0) -direction = np.argmax(abs(moving_window_velocity)) -old_x = grid.getmins()[direction] -new_x = old_x -for i in range(1, max_step + 1): +sim.add_species(electrons, layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=number_per_cell_each_dim)) - # check whether the moving window has updated - num_shift = int( 0.5 + (new_x - old_x) / dx[direction]) - if (num_shift != 0): - inject_new_plasma(ncells, domain_min, domain_max, num_shift, direction, injected_plasma_density, injected_plasma_ppc, plasma_min, plasma_max) - domain_min[direction] += num_shift*dx[direction] - domain_max[direction] += num_shift*dx[direction] +sim.add_laser(laser, injection_method=laser_antenna) - sim.step(1) +# write_inputs will create an inputs file that can be used to run +# with the compiled version. +sim.write_input_file(file_name = 'inputs_from_PICMI') - old_x = new_x - new_x = grid.getmins()[direction] +# Alternatively, sim.step will run WarpX, controlling it from Python +#sim.step() -sim.finalize() diff --git a/Example/plasma_acceleration/plasma_acceleration_PICMI.py b/Example/plasma_acceleration/plasma_acceleration_PICMI.py new file mode 100644 index 000000000..0d831f3e1 --- /dev/null +++ b/Example/plasma_acceleration/plasma_acceleration_PICMI.py @@ -0,0 +1,62 @@ +import numpy as np +from pywarpx import picmi +#from warp import picmi + +nx = 64 +ny = 64 +nz = 64 + +xmin = -200.e-6 +xmax = +200.e-6 +ymin = -200.e-6 +ymax = +200.e-6 +zmin = -200.e-6 +zmax = +200.e-6 + +moving_window_velocity = [0., 0., picmi.c] + +number_per_cell_each_dim = [2, 2, 1] + +grid = picmi.Cartesian3DGrid(number_of_cells = [nx, ny, nz], + lower_bound = [xmin, ymin, zmin], + upper_bound = [xmax, ymax, zmax], + lower_boundary_conditions = ['periodic', 'periodic', 'open'], + upper_boundary_conditions = ['periodic', 'periodic', 'open'], + moving_window_velocity = moving_window_velocity, + warpx_max_grid_size=32, warpx_max_level=0, warpx_coord_sys=0) + +solver = picmi.ElectromagneticSolver(grid=grid, cfl=1) + +beam_distribution = picmi.UniformDistribution(density = 1.e23, + lower_bound = [-20.e-6, -20.e-6, -150.e-6], + upper_bound = [+20.e-6, +20.e-6, -100.e-6], + directed_velocity = [0., 0., 1.e9]) + +plasma_distribution = picmi.UniformDistribution(density = 1.e22, + lower_bound = [-200.e-6, -200.e-6, 0.], + upper_bound = [+200.e-6, +200.e-6, None], + fill_in = True) + +beam = picmi.Species(particle_type='electron', name='beam', initial_distribution=beam_distribution) +plasma = picmi.Species(particle_type='electron', name='plasma', initial_distribution=plasma_distribution) + +sim = picmi.Simulation(solver = solver, + plot_int = 2, + verbose = 1, + cfl = 1.0, + max_steps = 1000, + warpx_current_deposition_algo = 3, + warpx_charge_deposition_algo = 0, + warpx_field_gathering_algo = 0, + warpx_particle_pusher_algo = 0) + +sim.add_species(beam, layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=number_per_cell_each_dim)) +sim.add_species(plasma, layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=number_per_cell_each_dim)) + +# write_inputs will create an inputs file that can be used to run +# with the compiled version. +sim.write_input_file(file_name = 'inputs_from_PICMI') + +# Alternatively, sim.step will run WarpX, controlling it from Python +#sim.step() + |