SmileiPIC / Smilei

Particle-in-cell code for plasma simulation
https://smileipic.github.io/Smilei
344 stars 120 forks source link

electron species name #468

Closed iplasma closed 2 years ago

iplasma commented 3 years ago

Hi,

I encountered some issue with crashing for the diagnostic at t=0 with both ionization and preionized species, and it seemed to be related to the name of the electron of the preionized species. It has to be called "elec", otherwise it will crashed. Just wondering if this is intended.

Thanks!

mccoys commented 3 years ago

No this is not intended. Please provide an example input file

iplasma commented 3 years ago

Here is the input file to reproduce the problem. It crashed when the electron name is "elec_C" instead of "elec". It actually has nothing to do with ionization so I remove those species in the input deck.

Thanks!

import numpy as np
twopi = 2*np.pi
twopi2 = twopi**2.0

box = [600,60]
res = [60,16]   # resolution in space: cell number per unit length
number_of_patches = [256,128]
tmax = 15.00  # max sim. time

dr =[ 1./i for i in res ]  
dt_CFL = 1./np.sqrt(sum(i**2 for i in res))  # Courant condition
res_t = (1./dt_CFL)/0.98  # resolution in t
dt =1./res_t
Nt = int(res_t) # steps needed to finish one time unit

for i in range(len(box)):  
        fac = float(box[i])*res[i]/number_of_patches[i]
        if fac!=fac//1: fac = int(fac)+1
        box[i] = float(fac)*number_of_patches[i]/res[i]
cells = list( np.array(box)*np.array(res) ) # total number of cells

lambda0=1000.0 #laser wavelegnth, nm
plasma_start = 100.0

a0=np.sqrt(1e19/1.1e18); tau=3e-12/(lambda0*1e-9/3e8);
focus=[plasma_start, box[1]/2.0]; waist=10.0

def time_envelope(t):
    if 0<t<tau: return( (np.sin(np.pi*t/tau))**2 )
    else: return (0)

def time_envelope_planewave(t):
    global a0
    if 0<t<tau: return( a0*(np.sin(np.pi*t/tau))**2 )
    else: return (0)

dump_freq = 5*Nt
point = [box[0]-100.0, box[1]/2.0]
vector=[1.0, 0.0]

Main(
    geometry = "2Dcartesian",
    interpolation_order = 2,
    grid_length  = box,
    cell_length = dr,
    number_of_patches = number_of_patches, # each must be power of 2, total be greater than MPI ranks
    timestep = dt,
    simulation_time = tmax,
    EM_boundary_conditions = [
        ['silver-muller'],
        ['periodic'], 
    ],
        solve_poisson = False,  # False to have neutralized background
    reference_angular_frequency_SI = 3.e8/1.0e-6, #speed of light / laser wavelegth, so length is normalized to laser wavelength, time to laser cycle, field to mcw0/e/(2pi), density to n_c/(2pi^2)                                                                                                                                                  
        print_every = 10,
)

LaserGaussian2D(box_side="xmin", a0=a0, omega=twopi, polarization_phi=0.0, ellipticity=0, focus=focus, waist=waist, incidence_angle=0.0, time_envelope=time_envelope,)

Species(name = 'C6', position_initialization = 'regular',momentum_initialization = 'cold', particles_per_cell = 1, mass = 12*1836.0, charge = 6.0, number_density = 1.0, boundary_conditions = [ ["remove"], ['periodic'],], ) 

Species(name = 'elec_C',position_initialization = 'regular',momentum_initialization = 'cold', particles_per_cell = 1, mass = 1.0, charge = -1.0, number_density = 1.0, temperature = [1.0e-10], boundary_conditions = [ ["remove"], ["periodic"]], )

DiagFields(
    every = [0*Nt, tmax*Nt, int(0.5*Nt)],
    flush_every = 50*Nt, 
        fields = ['Rho_elec_C', ]
)
mccoys commented 3 years ago

Thanks for the report. I will push a fix soon

mccoys commented 2 years ago

This is fixed already.