Bad termination for 1d smilei laser wakefield simulation

[YCAOsz]

Hi. I experienced a problem using smilei for 1d laser wakefield simulation. When I used the script below for a long 1d simulation, the process only produces about 25 outputs and shows this error:

Species creation summary
 --------------------------------------------------------------------------------
         Species 0 (electron) created with 0 particles
         Species 1 (proton) created with 0 particles
(this part is weird as I clearly reads some good rho fields in the outputs produced)

Time-Loop started: number of time-steps n_time = 6000000
 --------------------------------------------------------------------------------
          timestep       sim time   cpu time [s]   (    diff [s] )
 >>> Window starts moving
      5000/6000000     7.6896e+02     2.0917e+01   (  2.0917e+01 )
     10000/6000000     1.5378e+03     2.0693e+02   (  1.8602e+02 )
     15000/6000000     2.3067e+03     4.0574e+02   (  1.9881e+02 )
     20000/6000000     3.0756e+03     5.9944e+02   (  1.9370e+02 )

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 26310 RUNNING AT cx1-136-10-3.cx1.hpc.ic.ac.uk
=   EXIT CODE: 11
=   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES

I'm not sure why this happens. Maybe the density function? (The script is shown below) Thank you

---

################### 1D Laser Wakefield with envelope
import math
import numpy as np
from scipy.constants import c, e, m_e, m_p
#Import basic units
e0=8.854187817e-12
me=m_e
qe=e

lambda_laser=800*10**(-9)
w_laser=2*math.pi*c/lambda_laser
Lr=c/w_laser
Tr=1/w_laser
Nr=e0*me*w_laser**2/(qe**2)

#Fbpic code
ne_plasma=7e23/Nr
Nz = 5760
zmax = 7.287716e-5
zmin = -3.989999e-5
dt_fbpic = (zmax-zmin)/Nz/c   
N_step =6000000

a0 = 2         # Laser amplitude
w0 = 3.8e-5   # Laser waist, 3 times radius, 4 cells raidally, 4 cells , 5 asumuthal
ctau = 1.65e-5  # Laser duration,
z0 = 3.99e-5   # Laser centroid

# The density profile
ramp_start = 2.e-4/Lr
ramp_length = 5.e-4/Lr

def dens_func( x ) :
    """Returns relative density at position z and r"""
    # Allocate relative density
    n = np.ones_like(x)
    # Make linear ramp
    n = np.where( x<ramp_start+ramp_length, (x-ramp_start)/ramp_length, n )
    # Supress density before the ramp
    n = np.where( x<ramp_start, 0., n )
    return (n*ne_plasma)

#Covert to Smilei Variables
dx = (zmax-zmin)/Nz/Lr
#dtrans = 2*rmax/Nr/Lr
dt = dt_fbpic/Tr
nx = Nz
Lx = nx * dx
npatch_x = 64
laser_fwhm = ctau/Lr
center_laser = Lx-2*laser_fwhm # the temporal center here is the same as waist position, but in principle they can differ
time_start_moving_window = 0.

Main(
    geometry = "1Dcartesian",
    interpolation_order = 2,
    timestep = dt,
    simulation_time = N_step*dt,
    cell_length  = [dx],
    grid_length = [ Lx],
    number_of_patches =[npatch_x], 
    clrw = 6,
    EM_boundary_conditions = [ ["silver-muller"] ],
    solve_poisson = False,
    print_every = 5000,
    random_seed = smilei_mpi_rank
)

MovingWindow(
    time_start = time_start_moving_window,
    velocity_x = 0.9995 
)

LoadBalancing(
    initial_balance = False,
        every = 150,
    cell_load = 1.,
    frozen_particle_load = 0.1
)

Species(
    name = "electron",
    position_initialization = "regular",
    momentum_initialization = "cold",
    particles_per_cell = 80,
    c_part_max = 1.0,
    ponderomotive_dynamics = True, # = this species interacts with laser envelope
    mass = 1.0,
    charge = -1.0,
    number_density = dens_func,
    mean_velocity = [0.0, 0.0, 0.0],
    temperature = [0.0],
    pusher = "ponderomotive_boris", # pusher to interact with envelope
    #pusher = "boris",
    time_frozen = 0.0,
    boundary_conditions =  [
        ["remove", "remove"],
    ],
)

Species(
    name = "proton",
    position_initialization = "regular",
    momentum_initialization = "cold",
    particles_per_cell = 80,
    c_part_max = 1.0,
    ponderomotive_dynamics = True, # = this species interacts with laser envelope
    mass = 1836.152,
    charge = 1.0,
    number_density = dens_func,
    mean_velocity = [0.0, 0.0, 0.0],
    temperature = [0.0],
    pusher = "ponderomotive_boris", # pusher to interact with envelope
    #pusher = "boris",
    time_frozen = 0.0,
    boundary_conditions =  [
        ["remove", "remove"],
    ],
)

LaserEnvelopePlanar1D( # linear regime of LWFA
    a0              = 3.12,     
    time_envelope   = tgaussian(center=center_laser, fwhm=laser_fwhm),
    envelope_solver = 'explicit',
     Envelope_boundary_conditions = [ ["reflective", "reflective"],
     ],
)

Checkpoints(
    dump_step = 0,
    dump_minutes = 0.0,
    exit_after_dump = False,
)

list_fields = ['Ex','Ey','Rho','Rho_electron','Env_A_abs','Env_Chi','Env_E_abs']

DiagFields(
   every = 1000,
        fields = list_fields
)

DiagScalar(every = 1000, vars=['Utot','Env_A_absMax','Env_E_absMax'])

[jderouillat]

Hi, Could you provide the whole namelist within a file that we can run ? Here, some characters are missing.

Julien

[YCAOsz]

This is the python file (just change the .txt to .py would work) If I change the density function to constant(), it would run longer, but still stops at some point. Also, if I extract the envelope (Envabs), it gets blended into the field and the wakefield wasn't generated some time into the simulation. What should I modify to get a good 1d simulation of laser wakefield? Thank you (by the way is there a forum to discuss these issues of smiles?) Smilei_envelope_111519.txt

[Z10Frank]

Hello, I see that your timestep may be too high for an envelope simulation. The Courant-Friedrichs-Levy condition of the envelope solver is even more restrictive than the one of the "standard" Yee solver. Your envelope values maybe are increasing exponentially, blocking the simulation at some point. Normally I use a timestep dt=0.8*dx, could you try it?

With this modification I seem to arrive to a longer distance, with a well-formed wakefield:

PS since you are not interested in the ion motion, you can comment the ion Species block, the Esirkepov method will take care of adding an implicit layer of neutralising ions (see the section "Why ions are not present" in this tutorial https://smileipic.github.io/tutorials/advanced_wakefield_electron_beam.html)

[YCAOsz]

Thank you. The simulation can last to about 3cm in real unit before the envelope blows up. However, is it possible to make the simulation last even longer (to about 6cm)? Would it be better to use a gaussian laser wave? (though I'm not sure if i can put the laser at some offset in x coordinate so that I don't need to start at x_min and time the window movement.)

[Z10Frank]

Following the envelope evolution, what do you mean that it blows up? If it increases in an uncontrolled way, maybe you should decrease even more the dt, like 0.7*dx or 0.6*dx.

If you are seeing strange effects at the borders, since the envelope has only reflective boundary conditions, probably you will need a larger window.

Besides, you are using a dx that is very small (dx = 0.153 in code units), which should be used for a standard laser simulation where you need to sample the laser high frequency oscillations. With the envelope model you can increase it, to sample only the envelope of the laser (see Figure in http://www.maisondelasimulation.fr/smilei/laser_envelope.html) and reduce your simulation time.

In your case, your laser fwhm in field (not intensity fwhm) is laser_fwhm = 129 in code units (~21 optical cycles), so you can increase your dx even up to dx = 1 in code units. This way, your simulation will be much faster. If however you think that you will have to resolve finer envelope modulation effects (e.g. due to depletion), then you should use a dx smaller than 1.0

[mccoys]

@YCAOsz Did you try the recommendations from @Z10Frank ?

[Z10Frank]

@YCAOsz In absence of answer I'll close the issue, please feel free to reopen it if you have news about it.