diff --git a/Examples/laser_acceleration/laser_acceleration_2d_PICMI.py b/Examples/laser_acceleration/laser_acceleration_2d_PICMI.py new file mode 100644 index 0000000..94e0ea8 --- /dev/null +++ b/Examples/laser_acceleration/laser_acceleration_2d_PICMI.py @@ -0,0 +1,186 @@ +#!/usr/bin/env python3 + +# 2D version of the LWFA example, provided for quick debugging +# This should be the only line that needs to be changed for different codes +# e.g. `from pywarpx import picmi` +# `from fbpic import picmi` +# `from warp import picmi` +from plasmacode import picmi + +# Create alias fror constants +cst = picmi.constants + +# Run parameters - can be in separate file +# ======================================== + +# Physics parameters +# ------------------ + +# --- laser +import math +laser_polarization = math.pi/2 # Polarization angle (in rad) +laser_a0 = 4. # Normalized potential vector +laser_wavelength = 8e-07 # Wavelength of the laser (in meters) +laser_waist = 5e-06 # Waist of the laser (in meters) +laser_duration = 15e-15 # Duration of the laser (in seconds) +laser_injection_loc = 9.e-6 # Position of injection (in meters, along z) +laser_focal_distance = 100.e-6 # Focal distance from the injection (in meters) +laser_t_peak = 30.e-15 # The time at which the laser reaches its peak + # at the antenna injection location (in seconds) +# --- plasma +plasma_density_expression = "1.e23*(1+tanh((y - 20.e-6)/10.e-6))/2." +plasma_min = [-20.e-6, 0.0e-6] +plasma_max = [ 20.e-6, 1.e-3] + +# --- electron bunch +# bunch_physical_particles = 1.e8 +# in 2D or even 1D situations where the bunch is +# infinite in the ignorable direction, then the bunch +# density is the quantity that makes more sense +# For the example where sigma is 1 micron in all 3 directions +# and N is 10^8, this corresponds to a bunch density of +# 10^8/(1e-6**3(2*pi)**3/2) = 6.35e24 +bunch_density = 6.35e24 +bunch_rms_size = [1.e-6, 1.e-6] +bunch_rms_velocity = [0.,10.*cst.c] +bunch_centroid_position = [0.,-22.e-6] +bunch_centroid_velocity = [0.,1000.*cst.c] + +# Numerics parameters +# ------------------- + +# --- Nb time steps +max_steps = 1000 + +# --- grid +nx = 64 + +ny = 480 +xmin = 1.5*plasma_min[0] +xmax = 1.5*plasma_max[0] +ymin = -38.e-6 +ymax = 10.e-6 +moving_window_velocity = [0., cst.c] +n_macroparticle_per_cell = [4, 4] + +# --- geometry and solver +em_solver_method = 'CKC' # Cole-Karkkainen-Cowan stencil +geometry = '2D' + +# Physics part - can be in separate file +# ====================================== + +# Physics components +# ------------------ + +# --- laser +laser = picmi.GaussianLaser( + wavelength = laser_wavelength, + waist = laser_waist, + duration = laser_duration, + focal_position = [0., laser_focal_distance + laser_injection_loc], + centroid_position = [0., laser_injection_loc - cst.c*laser_t_peak], + polarization_direction = [math.sin(laser_polarization), 0., math.cos(laser_polarization),], + propagation_direction = [0,1], + a0 = laser_a0) + +# --- plasma +plasma_dist = picmi.AnalyticDistribution( + density_expression = plasma_density_expression, + lower_bound = plasma_min, + upper_bound = plasma_max, + fill_in = True) +plasma = picmi.MultiSpecies( + particle_types = ['He', 'Ar', 'electron'], + names = ['He+', 'Argon', 'e-'], + charge_states = [1, 5, None], + proportions = [0.2, 0.8, 0.2 + 5*0.8], + initial_distribution=plasma_dist) + +# --- electron bunch +beam_dist = picmi.GaussianBunchDistribution( +# n_physical_particles = bunch_physical_particles, + density = bunch_density + rms_bunch_size = bunch_rms_size, + rms_velocity = bunch_rms_velocity, + centroid_position = bunch_centroid_position, + centroid_velocity = bunch_centroid_velocity ) +beam = picmi.Species( particle_type = 'electron', + name = 'beam', + initial_distribution = beam_dist) + +# Numerics components +# ------------------- + +grid = picmi.Cartesian2DGrid( + number_of_cells = [nx, ny], + lower_bound = [xmin, ymin], + upper_bound = [xmax, ymax], + lower_boundary_conditions = ['periodic', 'open'], + upper_boundary_conditions = ['periodic', 'open'], + moving_window_velocity = moving_window_velocity, + ) + + +smoother = picmi.BinomialSmoother( n_pass = 1, + compensation = True ) +solver = picmi.ElectromagneticSolver( grid = grid, + cfl = 1., + method = em_solver_method, + source_smoother = smoother) + +# Diagnostics +# ----------- +field_diag = picmi.FieldDiagnostic(name = 'diag1', + grid = grid, + period = 100,) +part_diag = picmi.ParticleDiagnostic(name = 'diag1', + period = 100, + species = [beam]) + +# Simulation setup +# ----------------- + +# Initialize the simulation object +# Note that the time step size is obtained from the solver +sim = picmi.Simulation(solver = solver, verbose = 1) + +# Inject the laser through an antenna +antenna = picmi.LaserAntenna( + position = (0, 0, 9.e-6), + normal_vector = (0, 0, 1.)) +sim.add_laser(laser, injection_method = antenna) + +# Add the plasma: continuously injected by the moving window +plasma_layout = picmi.GriddedLayout( + grid = grid, + n_macroparticle_per_cell = n_macroparticle_per_cell) +sim.add_species(species=plasma, layout=plasma_layout) + +# Add the beam +beam_layout = picmi.PseudoRandomLayout( + n_macroparticles = 10**5, + seed = 0) +initialize_self_field = True +if picmi.codename == 'warpx': + initialize_self_field = False +sim.add_species(species=beam, layout=beam_layout, + initialize_self_field=initialize_self_field) + +# Add the diagnostics +sim.add_diagnostic(field_diag) +sim.add_diagnostic(part_diag) + +# Picmi input script +# ================== + +run_python_simulation = True + +if run_python_simulation: + # `sim.step` will run the code, controlling it from Python + sim.step(max_steps) +else: + # `write_inputs` will create an input file that can be used to run + # with the compiled version. + sim.set_max_step(max_steps) + sim.write_input_file(file_name='input_script') diff --git a/PICMI_Python/particles.py b/PICMI_Python/particles.py index 6a6c291..b2a833b 100644 --- a/PICMI_Python/particles.py +++ b/PICMI_Python/particles.py @@ -19,6 +19,8 @@ class PICMI_Species(_ClassWithInit): Species - particle_type=None: A string specifying an elementary particle, atom, or other, as defined in the openPMD 2 species type extension, openPMD-standard/EXT_SpeciesType.md - name=None: Name of the species + - method=None: One of 'Boris', 'Vay', 'Higuera-Cary', 'Li' , 'free-streaming', and 'LLRK4' (Landau-Lifschitz radiation reaction formula with RK-4) (string) + code-specific method can be specified using 'other:' - charge_state=None: Charge state of the species (applies to atoms) [1] - charge=None: Particle charge, required when type is not specified, otherwise determined from type [C] - mass=None: Particle mass, required when type is not specified, otherwise determined from type [kg] @@ -27,8 +29,16 @@ class PICMI_Species(_ClassWithInit): - particle_shape: Particle shape used for deposition and gather ; if None, the default from the `Simulation` object will be used. Possible values are 'NGP', 'linear', 'quadratic', 'cubic' """ + methods_list = ['Boris' , 'Vay', 'Higuera-Cary', 'Li', 'free-streaming', 'LLRK4'] + def __init__(self, particle_type=None, name=None, charge_state=None, charge=None, mass=None, - initial_distribution=None, particle_shape=None, density_scale=None, **kw): + initial_distribution=None, particle_shape=None, density_scale=None, method=None, **kw): + + + assert method is None or method in PICMI_Species.methods_list or method.startswith('other:'), \ + Exception('method must starts with either "other:", or be one of the following '+', '.join(PICMI_Species.methods_list)) + + self.method = method self.particle_type = particle_type self.name = name self.charge = charge @@ -152,18 +162,48 @@ class PICMI_GaussianBunchDistribution(_ClassWithInit): - centroid_velocity=[0,0,0]: Velocity (gamma*V) of the bunch centroid at t=0 (vector) [m/s] - velocity_divergence=[0,0,0]: Expansion rate of the bunch at t=0 (vector) [m/s/m] """ - def __init__(self,n_physical_particles, rms_bunch_size, - rms_velocity = [0.,0.,0.], - centroid_position = [0.,0.,0.], - centroid_velocity = [0.,0.,0.], - velocity_divergence = [0.,0.,0.], + def __init__(self,n_physical_particles = None, density = None, rms_bunch_size, + rms_velocity = None, + centroid_position = None, + centroid_velocity = None, + velocity_divergence = None, **kw): + + assert n_physical_particles is not None or density is not None, 'One of n_physical_particles or density must be speficied' + + if n_physical_particles is None: + if rms_bunch_size.size==3: + n_physical_particles = np.prod(rms_bunch_size)*(np.pi*2)**1.5*density + if rms_bunch_size.size==2: + n_physical_particles = rms_bunch_size[0]**2*rms_bunch_size[1]*(np.pi*2)**1.5*density + + if density is None: + if rms_bunch_size.size==3: + density = n_physica_particles/ (np.prod(rms_bunch_size)*(np.pi*2)**1.5) + if rms_bunch_size.size==2: + density = n_physical_particles/ (rms_bunch_size[0]**2*rms_bunch_size[1]*(np.pi*2)**1.5) + + + self.n_physical_particles = n_physical_particles self.rms_bunch_size = rms_bunch_size - self.rms_velocity = rms_velocity - self.centroid_position = centroid_position - self.centroid_velocity = centroid_velocity - self.velocity_divergence = velocity_divergence + if(centroid_position is None): + self.centroid_position = np.zeros(rms_bunch_size.size) + else: + self.centroid_position = centroid_position + if(centroid_velocity is None): + self.centroid_velocity = np.zeros(rms_bunch_size.size) + else: + self.centroid_velocity = centroid_velocity + if(velocity_divergence is None): + self.velocity_divergence = np.zeros(rms_bunch_size.size) + else: + self.velocity_divergence = velocity_divergence + if(rms_velocity is None): + self.rms_velocity = np.zeros(rms_bunch_size.size) + else: + self.rms_velocity = rms_velocity + self.handle_init(kw)