Examples

Laser target

A Gaussian laser pulse interacting with a target.

Automatically plot the simulation data with a callback plot_results.

Copy/download: laser-target.py.

import numpy as np

from lambdapic import (
    Electron,
    ExtractSpeciesDensity,
    GaussianLaser2D,
    PlotFields,
    Proton,
    SaveFieldsToHDF5,
    SaveSpeciesDensityToHDF5,
    Simulation,
    Species,
    c,
    callback,
    e,
    epsilon_0,
    get_fields,
    m_e,
    pi,
)

um = 1e-6
l0 = 0.8 * um
t0 = l0 / c
omega0 = 2 * pi * c / l0
nc = epsilon_0 * m_e * omega0**2 / e**2

nx = 1024
ny = 1024
dx = l0 / 50
dy = l0 / 50

Lx = nx * dx
Ly = ny * dy


def density(n0):
    def _density(x, y):
        ne = 0.0
        if x > Lx/2 and x < Lx/2+1*um:
            ne = n0
        return ne
    return _density

laser1 = GaussianLaser2D(
    a0=10,
    w0=2e-6,
    l0=0.8e-6,
    ctau=5e-6,
    focus_position=Lx/2,
    x0=10e-6
)
laser2 = GaussianLaser2D(
    a0=10,
    w0=2e-6,
    l0=0.8e-6,
    ctau=5e-6,
    focus_position=Lx/2,
    x0=10e-6,
    pol_angle=pi/2,
    cep=pi/2
)

sim = Simulation(
    nx=nx,
    ny=ny,
    dx=dx,
    dy=dy,
    nsteps=2000,
    log_file='laser-target.log',
)

ele = Electron(density=density(10*nc), ppc=10)
proton = Proton(density=density(10*nc/8*2), ppc=10)
carbon = Species(name="C", charge=6, mass=12*1800, density=density(10*nc/8), ppc=10)

sim.add_species([ele, carbon, proton])
    
if __name__ == "__main__":
    sim.run(2001, callbacks=[
            laser1+laser2, 
            n_ele := ExtractSpeciesDensity(sim, ele, 500),
            PlotFields(
                [
                    dict(field=n_ele.density, scale=1/nc, cmap='Grays', vmin=0, vmax=20), 
                    dict(field='ey',  scale=e/(m_e*c*omega0), cmap='bwr_alpha', vmin=-laser1.a0, vmax=laser1.a0)
                ],
                prefix='laser-target/ey', interval=10e-15,
            ),
            PlotFields(
                [
                    dict(field=n_ele.density, scale=1/nc, cmap='Grays', vmin=0, vmax=20), 
                    dict(field='ez',  scale=e/(m_e*c*omega0), cmap='bwr_alpha', vmin=-laser2.a0, vmax=laser2.a0)
                ],
                prefix='laser-target/ez', interval=10e-15,
            ),
            SaveFieldsToHDF5('laser-target/fields', 500, ['ex', 'ey', 'ez', 'bx', 'by', 'bz', 'rho']),
            SaveSpeciesDensityToHDF5(carbon, 'laser-target/density', 500),
            SaveSpeciesDensityToHDF5(ele, 'laser-target/density', 500),
        ]
    )

Laser wakefield

For demostration of moving window.

Copy/download: lwfa.py.

import numpy as np
from scipy.constants import c, e, epsilon_0, m_e, mu_0, pi

from lambdapic import (
    Electron,
    ExtractSpeciesDensity,
    MovingWindow,
    PlotFields,
    Proton,
    SaveFieldsToHDF5,
    SaveSpeciesDensityToHDF5,
    SimpleLaser2D,
    Simulation,
    Species,
    c,
    callback,
    e,
    epsilon_0,
    get_fields,
    m_e,
    pi,
)

um = 1e-6
l0 = 0.8 * um
t0 = l0 / c
omega0 = 2 * pi * c / l0
nc = epsilon_0 * m_e * omega0**2 / e**2

nx = 500
ny = 800
dx = l0 / 20
dy = l0 / 20

Lx = nx * dx
Ly = ny * dy


def density(n0):
    def _density(x, y):
        ne = 0.0
        if x > 1*um:
            ne = n0
        if abs(y - Ly/2) > Ly/2 - 1*um:
            ne = 0
        return ne
    return _density

# move velocity supports constant velocity and time-dependent velocity
# here lambda t: c + (t-Lx/c)*0 is just a constant for demonstration
movingwindow = MovingWindow(velocity=lambda t: c + (t-Lx/c)*0)

laser = SimpleLaser2D(
    a0=2,
    w0=5e-6,
    l0=0.8e-6,
    ctau=5e-6,
)

ne = 0.01*nc

sim = Simulation(
    nx=nx,
    ny=ny,
    dx=dx,
    dy=dy,
    npatch_x=10,
    npatch_y=10,
    dt_cfl=0.99,
    sim_time=100e-15
)

ele = Electron(density=density(ne), ppc=10)
proton = Proton(density=density(ne/8*2), ppc=2)
carbon = Species(name="C", charge=6, mass=12*1800, density=density(ne/8), ppc=1)

sim.add_species([ele, carbon, proton])

interval = 10e-15
if __name__ == "__main__":
    sim.run(2001, callbacks=[
            movingwindow,
            laser,
            n_ele := ExtractSpeciesDensity(sim, ele, interval),
            PlotFields(
                [
                    dict(field=n_ele.density, scale=1/nc, cmap='Grays', vmin=0, vmax=ne/nc*2), 
                    dict(field='ey',  scale=e/(m_e*c*omega0), cmap='bwr_alpha', vmin=-laser.a0, vmax=laser.a0)
                ],
                prefix='lwfa', interval=interval,
            ),
            SaveFieldsToHDF5('lwfa/fields', interval, ['ex', 'ey', 'ez', 'bx', 'by', 'bz', 'jx', 'jy', 'rho']),
            SaveSpeciesDensityToHDF5(proton, 'lwfa/', interval),
        ]
    )

QED photons

For demostration of QED photons and callback turning on/off QED.

The example has a enable_radiation switch that enables radiation at desired time.

Copy/download: photons.py.

from pathlib import Path

import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap

from lambdapic import (
    Electron,
    ExtractSpeciesDensity,
    MovingWindow,
    Photon,
    PlotFields,
    Proton,
    SaveFieldsToHDF5,
    SaveSpeciesDensityToHDF5,
    SimpleLaser2D,
    Simulation,
    Species,
    c,
    callback,
    e,
    epsilon_0,
    get_fields,
    m_e,
    pi,
)
from lambdapic.core.utils.logger import logger

um = 1e-6
l0 = 0.8 * um
omega0 = 2 * pi * c / l0
nc = epsilon_0 * m_e * omega0**2 / e**2

nx = 512
ny = 512
dx = l0 / 20
dy = l0 / 20

Lx = nx * dx
Ly = ny * dy


def density(n0):
    def _density(x, y):
        ne = 0.0
        if x > 2*um:
            ne = n0
        return ne
    return _density

laser = SimpleLaser2D(
    a0=300,
    w0=2e-6,
    l0=0.8e-6,
    ctau=5e-6
)

sim = Simulation(
    nx=nx,
    ny=ny,
    dx=dx,
    dy=dy,
    sim_time=100e-15,
)

ele = Electron(density=density(5*nc), ppc=10, radiation="photons")
pho = Photon()
ele.set_photon(pho)

proton = Proton(density=density(5*nc), ppc=10)

sim.add_species([ele, proton, pho])

interval = 10e-15

@callback(interval=interval)
def npho(sim: Simulation):
    npart = 0
    for ipatch, p in enumerate(sim.patches):
        part = p.particles[pho.ispec]
        npart += part.is_alive.sum()
    
    npart = sim.mpi.comm.reduce(npart)
    if sim.mpi.rank == 0:
        logger.info(f"nphoton = {npart}")

@callback("current_deposition", interval=interval)
def prune(sim: Simulation):
    for ipatch, p in enumerate(sim.patches):
        p.particles[sim.ispec].prune()
    
if __name__ == "__main__":
    sim.run(callbacks=[
            laser, 
            n_ele := ExtractSpeciesDensity(sim, ele, interval),
            PlotFields(
                [
                    dict(field=n_ele.density, scale=1/nc, cmap='Grays', vmin=0, vmax=10), 
                    dict(field='ey',  scale=e/(m_e*c*omega0), cmap='bwr_alpha', vmin=-laser.a0, vmax=laser.a0)
                ],
                prefix='photons', interval=interval,
            ),
            npho,
            prune,
        ]
    )