How to use fbpic - 10 common examples

def simulate_periodic_plasma_wave( particle_shape, show=False ):
    "Simulate a periodic plasma wave and check its fields"

    # Initialization of the simulation object
    sim = Simulation( Nz, zmax, Nr, rmax, Nm, dt,
                  p_zmin, p_zmax, p_rmin, p_rmax, p_nz, p_nr,
                  p_nt, n_e, n_order=n_order, use_cuda=use_cuda,
                  particle_shape=particle_shape )

    # Save the initial density in spectral space, and consider it
    # to be the density of the (uninitialized) ions
    sim.deposit('rho_prev', exchange=True)
    sim.fld.spect2interp('rho_prev')
    rho_ions = [ ]
    for m in range(len(sim.fld.interp)):
        rho_ions.append( -sim.fld.interp[m].rho.copy() )

    # Impart velocities to the electrons
    # (The electrons are initially homogeneous, but have an
    # intial non-zero velocity that develops into a plasma wave)
    impart_momenta( sim.ptcl[0], epsilons, k0, w0, wp )

profile with a gaussian profile
    """
    # Extract the part that has information
    if info_in_real_part:
        interp1 = interp1_complex.real
        zero_part = interp1_complex.imag
    else:
        interp1 = interp1_complex.imag
        zero_part = interp1_complex.real

    # Control that the part that has no information is 0
    assert np.allclose( 0., zero_part, atol=1.e-6*interp1.max() )

    # Get the predicted properties of the laser in the boosted frame
    if gamma_b is None:
        boost = BoostConverter(1.)
    else:
        boost = BoostConverter(gamma_b)
    ctau_b, lambda0_b, Lprop_b, z0_b = \
        boost.copropag_length([ctau, 0.8e-6, Lprop, z0])
    # Take into account whether the pulse is propagating forward or backward
    if not forward_propagating:
        Lprop_b = - Lprop_b

    # Fit the on-axis profile to extract a0
    def fit_function(z, a0, z0_phase):
        return( gaussian_laser( z, r[0], a0, z0_phase,
                                z0_b+Lprop_b, ctau_b, lambda0_b ) )
    fit_result = curve_fit( fit_function, z, interp1[:,0],
                            p0=np.array([a0, z0_b+Lprop_b]) )
    a0_fit, z0_fit = fit_result[0]

Returns
    -------
    A dictionary containing :
    - 'E' : 1d array containing the values of the amplitude
    - 'w' : 1d array containing the values of waist
    - 'fld' : the Fields object at the end of the simulation.
    - 'z_centroid' : 1d array containing the values of the centroid position
    """

    # Choose the timestep, so that the simulation advances by
    # one cell at every timestep
    dt = Lz*1./Nz * 1./c
    
    # Initialize the simulation object
    sim = Simulation( Nz, Lz, Nr, Lr, Nm, dt, p_zmin=0, p_zmax=0,
                    p_rmin=0, p_rmax=0, p_nz=2, p_nr=2, p_nt=2, n_e=0.,
                    use_cuda=True, boundaries='open' )
    # Remove the particles
    sim.ptcl = []
    # Create moving window object
    sim.set_moving_window( v=c )

    # Initialize the laser fields
    z0 = Lz/2
    init_fields( sim.fld, w0, ctau, k0, z0, E0, m )

    # Create the arrays to get the waist and amplitude
    w = np.zeros(Nt)
    E = np.zeros(Nt)
    z_center = np.zeros(Nt)

v_comoving : float
        Velocity at which the currents are assumed to move

    use_galilean: bool
        Whether to use a galilean frame that moves at the speed v_comoving

    Returns
    -------
    A dictionary containing :
    - 'E' : 1d array containing the values of the amplitude
    - 'w' : 1d array containing the values of waist
    - 'fld' : the Fields object at the end of the simulation.
    """

    # Initialize the simulation object
    sim = Simulation( Nz, zmax, Nr, Lr, Nm, dt,
                    n_order=n_order, zmin=zmin, use_cuda=use_cuda,
                    boundaries=boundaries, v_comoving=v_comoving,
                    exchange_period = 1, use_galilean=use_galilean )
    # Remove the particles
    sim.ptcl = []
    # Set the moving window object
    if v_window !=0:
        sim.set_moving_window( v=v_window )

    # Initialize the laser fields
    z0 = (zmax+zmin)/2
    init_fields( sim, w0, ctau, k0, z0, zf, E0, m )

    # Create the arrays to get the waist and amplitude
    w = np.zeros(N_diag)
    E = np.zeros(N_diag)

N_step = 300     # Number of iterations to perform
N_show = 300     # Number of timestep between every plot

# The particles
gamma0 = 25.
p_zmin = 15.e-6  # Position of the beginning of the bunch (meters)
p_zmax = 25.e-6  # Position of the end of the bunch (meters)
p_rmin = 0.      # Minimal radial position of the bunch (meters)
p_rmax = 5.e-6   # Maximal radial position of the bunch (meters)
n_e = 4.e18*1.e6 # Density (electrons.meters^-3)
p_nz = 2         # Number of particles per cell along z
p_nr = 2         # Number of particles per cell along r
p_nt = 4         # Number of particles per cell along theta

# Initialize the simulation object
sim = Simulation( Nz, zmax, Nr, rmax, Nm, dt,
    p_zmin, p_zmax, p_rmin, p_rmax, p_nz, p_nr, p_nt, n_e, n_order=n_order,
    v_comoving=-0.999999*c, use_galilean=True )

# Suppress the particles that were intialized by default and add the bunch
sim.ptcl = [ ]
add_elec_bunch( sim, gamma0, n_e, p_zmin, p_zmax,
            p_rmin, p_rmax, direction='backward' )


# Show the initial fields
plt.figure(0)
sim.fld.interp[0].show('Ez')
plt.figure(1)
sim.fld.interp[0].show('Er')
plt.show()
print( 'Done' )

v_comoving : float
        Velocity at which the currents are assumed to move

    use_galilean: bool
        Whether to use a galilean frame that moves at the speed v_comoving

    Returns
    -------
    A dictionary containing :
    - 'a' : 1d array containing the values of the amplitude
    - 'w' : 1d array containing the values of waist
    - 'fld' : the Fields object at the end of the simulation.
    """

    # Initialize the simulation object
    sim = Simulation( Nz, zmax, Nr, Lr, Nm, dt, p_zmin=0, p_zmax=0,
                    p_rmin=0, p_rmax=0, p_nz=2, p_nr=2, p_nt=2, n_e=0.,
                    n_order=n_order, zmin=zmin, use_cuda=use_cuda,
                    boundaries=boundaries, v_comoving=v_comoving,
                    exchange_period = 1, use_galilean=use_galilean,
                    use_envelope = True )
    # Remove the particles
    sim.ptcl = []
    # Set the moving window object
    if v_window !=0:
        sim.set_moving_window( v=v_window )

    # Initialize the laser fields
    z0 = (zmax+zmin)/2
    init_fields( sim, w0, ctau, k0, z0, zf, a0, m )

    # Create the arrays to get the waist and amplitude

interpolation='nearest')
    plt.xlabel('z')
    plt.ylabel('r')
    cb = plt.colorbar()
    cb.set_label('Imaginary part')

    plt.show()

Nm = 1
dt = (zmax-zmin)*1./c/Nz
dt = 0.06e-6/c
print(c*dt)
print(L_prop)
print(L_prop / c / dt)

sim = Simulation( Nz, zmax, Nr, rmax, Nm, dt,
    p_zmin=p_zmin, p_zmax=p_zmax, p_rmin=p_rmin, p_rmax=p_rmax, p_nz=p_nz,
    p_nr=p_nr, p_nt=Nm+1, n_e=n_e, n_order=n_order, zmin=zmin,
    boundaries='open', v_comoving=0.999*c, use_galilean=True, initialize_ions=True,
    use_cuda=use_cuda, use_envelope=True )
sim.set_moving_window(v=c)
tau = ctau/c
lambda0 = 2*np.pi/k0
# Create the relevant laser profile
z0 = 0
sim.diags = [ FieldDiagnostic( diag_period, sim.fld, comm=sim.comm, fieldtypes = ["a", "rho", "E"] ),
        ParticleDiagnostic(diag_period, {'electrons': sim.ptcl[0]}, sim.comm )]

profile = GaussianLaser( a0=a0, waist=w0, tau=tau,
            lambda0=lambda0, z0=z0, zf=zf )

add_laser_pulse( sim, profile, method = 'direct_envelope' )

for scheme in [ 'standard', 'galilean', 'pseudo-galilean']:

        # Choose the correct parameters for the scheme
        if scheme == 'standard':
            v_comoving = 0.
            use_galilean = False
        else:
            v_comoving = 0.9999*c
            if scheme == 'galilean':
                use_galilean = True
            else:
                use_galilean = False

        # Initialize the simulation object
        sim = Simulation( Nz, zmax, Nr, rmax, Nm, dt,
            p_zmin, p_zmax, p_rmin, p_rmax, p_nz, p_nr, p_nt, n_e,
            zmin=zmin, initialize_ions=True,
            v_comoving=v_comoving, use_galilean=use_galilean,
            boundaries='periodic', use_cuda=use_cuda )

        # Give a relativistic velocity to the particle, with some noise
        sim.ptcl[0].uz[:] = uz_m
        sim.ptcl[0].inv_gamma[:] = 1./np.sqrt( 1 + sim.ptcl[0].uz**2 )
        sim.ptcl[1].uz[:] = uz_m
        sim.ptcl[1].inv_gamma[:] = 1./np.sqrt( 1 + sim.ptcl[1].uz**2 )

        # Perform the simulation;
        # record the rms electric field every 50 timestep
        Er_rms = np.zeros(int(N_step/30)+1)
        t = np.zeros(int(N_step/30+1))
        Er_rms[0] = get_Er_rms(sim)

m: int, optional
        The mode on which to imprint the profile
        For m = 0 : gaussian profile, linearly polarized beam
        For m = 1 : annular profile
    """
    # Initialize the fields

    tau = ctau/c
    lambda0 = 2*np.pi/k0
    # Create the relevant laser profile

    if m == 0:
        profile = GaussianLaser( a0=a0, waist=w, tau=tau,
                    lambda0=lambda0, z0=z0, zf=zf )
    elif m == 1:
        profile = LaguerreGaussLaser( 0, 1, a0=a0, waist=w, tau=tau,
                    lambda0=lambda0, z0=z0, zf=zf )
    elif m == -1:
        profile = DonutLikeLaguerreGaussLaser( 0, -1, a0=a0/2**.5,
                    waist=w, tau=tau, lambda0=lambda0, z0=z0, zf=zf )

    # Add the profiles to the simulation
    add_laser_pulse( sim, profile, method = 'direct_envelope' )

The initial a0 of the pulse

    m: int, optional
        The mode on which to imprint the profile
    """
    # Initialize the fields
    tau = ctau/c
    lambda0 = 2*np.pi/k0
    # Create the relevant laser profile
    if m == 0:
        profile = GaussianLaser( a0=a0, waist=w0, tau=tau,
                    lambda0=lambda0, z0=z0, zf=zf )
    elif m == 1:
        # Put the peak of the Laguerre-Gauss at pi/4 to check that the
        # angular dependency is correctly captured
        profile = LaguerreGaussLaser( 0, 1, a0, w0, tau,
                    z0, lambda0=lambda0, zf=zf, theta0=np.pi/4 )

    # Add the profiles to the simulation
    add_laser_pulse( sim, profile, method = 'direct_envelope' )