SmileiPIC / Smilei

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

Segmentation fault with time_frozen #654

Open iplasma opened 1 year ago

iplasma commented 1 year ago

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

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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

iplasma commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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

iplasma commented 1 year ago

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 commented 11 months ago

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