2D LWFA Example w/ GaussianBunchDistribution

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')

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:<method>'
       - 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)