Segmentation fault with time_frozen

[iplasma]

Hi,

With the attached test input file and latest master branch of Smilei, the simulation will crash with message as follows. This appears to be related to the time_frozen variable, for example, if it is set to a value between 0 and 6, the code seems to run fine, but it crashes for time_frozen = 7.0, which is very strange. Any idea of what could be wrong?

Segmentation fault (Address not mapped to object [0x801878c88])
Stack trace (most recent call last):
#13 Object "[0xffffffffffffffff]", at 0xffffffffffffffff, in
#12 Object "/global/Work/smilei_github/Smilei/smilei", at 0x40c959, in _start
#11 Object "/lib64/libc.so.6", at 0x7f62fb9ec24c, in __libc_start_main
#10 Object "/global/Work/smilei_github/Smilei/smilei", at 0x8e3e3a, in main
#9 Object "/opt/AMD/aocc-compiler-3.2.0/bin/../lib/libomp.so", at 0x7f62fbbf7340, in __kmpc_fork_call
#8 Object "/opt/AMD/aocc-compiler-3.2.0/bin/../lib/libomp.so", at 0x7f62fbc09711, in __kmp_fork_call
#7 Object "/opt/AMD/aocc-compiler-3.2.0/bin/../lib/libomp.so", at 0x7f62fbc00ace, in
#6 Object "/opt/AMD/aocc-compiler-3.2.0/bin/../lib/libomp.so", at 0x7f62fbc6a3b2, in __kmp_invoke_microtask
#5 Object "/global/Work/smilei_github/Smilei/smilei", at 0x8e5fb0, in
#4 Object "/global/Work/smilei_github/Smilei/smilei", at 0x7cb761, in VectorPatch::dynamics(Params&, SmileiMPI*, SimWindow*, RadiationTables&, MultiphotonBreitWheelerTables&, double, Timers&, int)
#3 Object "/global/Work/smilei_github/Smilei/smilei", at 0x7cbb76, in VectorPatch::dynamicsWithoutTasks(Params&, SmileiMPI*, SimWindow*, RadiationTables&, MultiphotonBreitWheelerTables&, double, Timers&, int)
Stack trace (most recent call last) in thread 869352:
#2 Object "/global/Work/smilei_github/Smilei/smilei", at 0x91f2fa, in Species::dynamics(double, unsigned int, ElectroMagn*, Params&, bool, PartWalls*, Patch*, SmileiMPI*, RadiationTables&, MultiphotonBreitWheelerTables&)
#1 Object "/global/Work/smilei_github/Smilei/smilei", at 0x8905d0, in Projector1D2Order::currentsAndDensityWrapper(ElectroMagn*, Particles&, SmileiMPI*, int, int, int, bool, bool, int, int, int)
#0 Object "/global/Work/smilei_github/Smilei/smilei", at 0x891448, in Projector1D2Order::currents(double*, double*, double*, Particles&, unsigned int, double, int*, double*, int)

Parameters

• Smilei version: 4.7-1037-g81d858dd9-HEAD

• C++ compiler : both AOCC and gnu

• Input file

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

# ---------- pre-processing -------------------
## global info
box = [800]
res = [200]   # resolution in space: cell number per unit length
number_of_patches = [2048]
tmax = 20000  # max sim. time
## processing global info
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

time_frozen= 18.0 #plasma_start-2.0

## adjust box size (cells) to be divided by number_of_pathces
##   ** do not change **
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


# ------------ global settings --------------------
Main(
    geometry = "1Dcartesian",
    interpolation_order = 2,
    grid_length  = box,
    cell_length = dr,
    number_of_patches = number_of_patches, 
    timestep = dt,
    simulation_time = tmax,
    EM_boundary_conditions = [
        ['silver-muller'],
    ],
        solve_poisson = False,  
        reference_angular_frequency_SI = 3.e8/1.0e-6, 
        print_every = 10,
)


#--------------------- plasma -------------
ionize_model = "tunnel"
Species(name = 'carbon', ionization_model=ionize_model, ionization_electrons="elec", maximum_charge_state=6.0, position_initialization = 'regular', momentum_initialization = 'maxwell-juettner', temperature = [0.01/511]*3, particles_per_cell = 100, atomic_number = 6.0, mass = 12*1836.0, charge = 2.0, number_density = 1.0, boundary_conditions = [ ["remove", 'remove'],], time_frozen=time_frozen,)

Species(name = 'hydrogen', ionization_model=ionize_model, ionization_electrons="elec", maximum_charge_state=1.0, position_initialization = 'regular', momentum_initialization = 'maxwell-juettner', temperature = [0.01/511]*3, particles_per_cell = 100, atomic_number=1.0, mass = 1.0*1836.0, charge = 1.0, number_density = 1.0, boundary_conditions = [ ["remove", 'remove'],], time_frozen=time_frozen,)

Species(name = 'elec',position_initialization = 'regular',momentum_initialization = 'maxwell-juettner', temperature = [0.01/511]*3, particles_per_cell = 100, mass = 1.0, charge = -1.0, charge_density = 3.0, boundary_conditions = [ ["remove", 'remove'],], time_frozen=time_frozen,)

[mccoys]

Time freezing can be tricky to use. If particles should have been moving, the plasma response would be wrong, and letting them move suddenly could make discontinuities. I am not sure how stable the code would be in that case

[iplasma]

In my case, time_frozen is used to save some simulation time before a laser propagates to the plasma. For this test input file, the plasma should be neutral and the field should be small, so I don't understand why time_frozen would matter. I just tried reducing the plasma density by 100x, the simulation crashes with the same error.

[beck-llr]

Would it be possible to have the first ionization process occuring on a frozen species at time t=7. ?

[iplasma]

It should be related to the ionization process. The simulation with ionize_model="none" seems to run fine with time_frozen. BTW, smilei_test also crashes with this test input file.

[beck-llr]

Then unfreezing a frozen particle which has been ionized is probably the cause of the problem. You probably need to unfreeze the species before the first ionization occurs. Or keep it frozen the whole time.
smilei_test runs fine on my system with your input file. Which error do you get ?

[iplasma]

I may misunderstand the meaning of time_frozen here. Is it only applied to the pusher? I was thinking it should be applied to the ionization processes too, which would be convenient to skip those along with particle pushes until a specific time. If the ionization/collisions calculations are active since the beginning of the simulation, unfreezing those species perhaps would be problematic. So is there another way to control when to turn on the ionization/collision?

The smilei_test problem is probably specific to the aocc compiler, with gnu compiler it works fine. Below is the error message with aocc compiler.

     ----- TEST MODE WITH INTENDED PARTITION 1x1-----
                    _            _
  ___           _  | |        _  \ \   Version : 4.7-1037-g81d858dd9-HEAD
 / __|  _ __   (_) | |  ___  (_)  | |   
 \__ \ | '  \   _  | | / -_)  _   | |
 |___/ |_|_|_| |_| |_| \___| |_|  | |  
                                 /_/    
 
 

 Reading the simulation parameters
 --------------------------------------------------------------------------------
 HDF5 version 1.12.2
 Python version 3.9.16
         Parsing pyinit.py
         Parsing 4.7-1037-g81d858dd9-HEAD
         Parsing pyprofiles.py
         Parsing inp.master.py
         Parsing pycontrol.py
         Check for function preprocess()
         python preprocess function does not exist
         Calling python _smilei_check
         Calling python _keep_python_running() :
CAREFUL: Patches distribution: hilbertian


[WARNING](0) src/Params/Params.cpp:1123 (compute) simulation_time has been redefined from 20000.000000 to 19999.996800 to match timestep.
 

 Geometry: 1Dcartesian
 --------------------------------------------------------------------------------
         Interpolation order : 2
         Maxwell solver : Yee
         simulation duration = 19999.996800,   total number of iterations = 4081632
         timestep = 0.004900 = 0.980000 x CFL,   time resolution = 204.081633
         Grid length: 808.96
         Cell length: 0.005, 0, 0
         Number of cells: 161792
         Spatial resolution: 200
 

 Electromagnetic boundary conditions
 --------------------------------------------------------------------------------
         xmin silver-muller, absorbing vector [1]
         xmax silver-muller, absorbing vector [-1]
 

 Vectorization: 
 --------------------------------------------------------------------------------
         Mode: off
 

 Initializing MPI
 --------------------------------------------------------------------------------
         MPI_THREAD_MULTIPLE enabled
         Number of MPI processes: 1
         Number of threads per MPI process : 1
         OpenMP task parallelization not activated
 
         Number of patches: 2048
         Number of cells in one patch: 79
         Dynamic load balancing: never
 

 Initializing the restart environment
 --------------------------------------------------------------------------------
 

 Initializing species
 --------------------------------------------------------------------------------
         
         Creating Species #0: carbon
                 > Pusher: boris
                CAREFUL: For species 'carbon' possible conflict between time-frozen & not cold initialization

                 > Species frozen until time: 18.000000
                 > Boundary conditions: remove remove
                 > Density profile: 1D built-in profile `constant` (value: 0.010000)
         
         Creating Species #1: hydrogen
                 > Pusher: boris
                CAREFUL: For species 'hydrogen' possible conflict between time-frozen & not cold initialization

                 > Species frozen until time: 18.000000
                 > Boundary conditions: remove remove
                 > Density profile: 1D built-in profile `constant` (value: 0.010000)
         
         Creating Species #2: elec
                 > Pusher: boris
                CAREFUL: For species 'elec' possible conflict between time-frozen & not cold initialization

                 > Species frozen until time: 18.000000
                 > Boundary conditions: remove remove
                 > Density profile: 1D built-in profile `constant` (value: 0.030000)
 

 Initializing Patches
 --------------------------------------------------------------------------------
         First patch created
         All patches created
Stack trace (most recent call last):
#7    Object "[0xffffffffffffffff]", at 0xffffffffffffffff, in 
#6    Object "/global/Work/smilei_github/Smilei/smilei_test", at 0x40c959, in _start
#5    Object "/lib64/libc.so.6", at 0x7fcca65e424c, in __libc_start_main
#4    Object "/global/Work/smilei_github/Smilei/smilei_test", at 0x8e1585, in main
#3    Object "/global/Work/smilei_github/Smilei/smilei_test", at 0x8e56fd, in executeTestMode(VectorPatch&, Region&, SmileiMPI*, SimWindow*, Params&, Checkpoint&, OpenPMDparams&, RadiationTables*)
#2    Object "/global/Work/smilei_github/Smilei/smilei_test", at 0x8e6e3b, in PatchesFactory::createVector(VectorPatch&, Params&, SmileiMPI*, OpenPMDparams&, RadiationTables*, unsigned int, unsigned int)
#1    Object "/global/Work/smilei_github/Smilei/smilei_test", at 0x7e6829, in VectorPatch::updateFieldList(SmileiMPI*)
#0    Object "/global/Work/smilei_github/Smilei/smilei_test", at 0x91022b, in AsyncMPIbuffers::defineTags(Patch*, SmileiMPI*, int)
Segmentation fault (Address not mapped to object [0x168])

[mccoys]

For a reason I cannot remember, the ionization process is not frozen. @beck-llr is there a good reason, for example in the case of LWFA?

[beck-llr]

Yes there are simulations where you want a frozen ion background that can still be ionized and generate (non frozen) electrons.

As I see it frozen species skip pusher and projector because the don't move therefore do not deposit current. Interpolation is still done if ionization is turned on otherwise it is skipped too.

[beck-llr]

Something to try: turn off the ionization and dump a checkpoint at the time you want to start ionization. Do a restart with ionization turned on ? that might work ...

[iplasma]

Sounds like a good idea. I will try and report back.

[iplasma]

Yes, turning on ionization in a restart run from a previous frozen run without ionization works fine. Thanks for this idea. Still, it requires manual changes and would be nice to have a variable to turn on/off the ionization just like time_frozen. .

[mccoys]

@beck-llr I am not sure I understand why your solution works. What is the problem about unfreezing ionized atoms?