Moving Window crash with DEFINES += $(D)PER_PARTICLE_CHARGE_MASS

[DanRRRR]

Stuart, It took me more than a month to find why moving window crashed my code. Looks like this was not caused by any my additions.

If in Makefile i use the option DEFINES += $(D)PER_PARTICLE_CHARGE_MASS then the EPOCH crashes at the time moving window starts with no single changes in EPOCH source code (latest version).

I tried several different input files and here is the input file i took from our users. This works OK if

DEFINES += $(D)PER_PARTICLE_CHARGE_MASS is commented out. But it crashes at the time "window_start_time = 5.e-14" if i remove the

This demo takes few minutes on a single core:

begin:constant

lambda0 = 0.8*micron
I = 3.42e19 #a0=4
waist = 9.50*micron
pulse_duration = 17*femto
sigma_t = pulse_duration/(sqrt(2loge(2)))

r1 = 0e-5
r2 = 50e-6
den_peak = 5.0e24

x1 = 0 * micron
x2 = 65*micron
y1 = -38*micron
y2 = 38*micron
z1 = y1
z2 = y2
cpwlx = 25
cpwly = 6
cpwlz = cpwly
part_per_cell = 4
z = 0
r = sqrt(y^2+z^2)

sim_end = 2000*femto
snap_freq = 10*femto
end:constant

begin:control
nx = 2032 # (x2-x1)*cpwlx/lambda0
ny = 570  # (y2-y1)*cpwly/lambda0
# nz = (z2-z1)*cpwlz/lambda0
npart = part_per_cell*nx*ny               #   *nz
x_min =  0. *micron                       
x_max =  65.*micron                   
y_min = -38.*micron
y_max =  38.*micron
#z_min = z1
#z_max = z2
t_end = sim_end
stdout_frequency = 100
field_order = 6

end:control

begin:boundaries
bc_x_min = simple_laser
bc_x_max = simple_outflow
bc_y_min = simple_outflow
bc_y_max = simple_outflow
bc_z_min = simple_outflow
bc_z_max = simple_outflow

end:boundaries

begin:species
name = electron
charge = -1.0
mass = 1.0
density = if(x lt r2,0.0,1.0)
density = if(x lt r2 and x gt r1,(x-r1)/(r2-r1),density(electron))
density = density(electron)*den_peak
frac = 1.0
end:species

begin:laser
boundary = x_min
intensity_w_cm2 = I
lambda = lambda0
phase     = 0.0
pol_angle = 0
#profile = supergauss(sqrt(y*y+z*z), 0., waist, 8.)
profile = supergauss(r, 0., waist, 8.)
t_profile = gauss(time, 3.*sigma_t, sigma_t)
end:laser

begin:window
move_window = T
window_v_x = 299233811
window_start_time = 5.e-14  #  2.05e-13
bc_x_min_after_move = simple_outflow
bc_x_max_after_move = simple_outflow
end:window

begin:output
dt_snapshot = 10*femto
use_offset_grid = T

file_prefix = z

grid = always
#ex = always + single
ey = always + single
#ez = always + single
#bx = always + single
#by = always + single
#bz = always + single

particles = always
#px = always
#!py = always
#pz = always
#id = always
#particle_weight = always + single
charge = never
mass = never
#number_density = always + single
end:output

[DanRRRR]

By the way here is how it looks before it crash. Further history is impossible. That plotting was100% done in Fortran

[Status-Mirror]

Good spot, and nice plot!

When the window moves, new cells are filled with particles based on the density and temperature specified by the input.deck. This is done in the insert_particles subroutine in housekeeping/window.F90. There are 3 lines which set the momentum in this subroutine which read species%mass, but this is not set when -DPER_PARTICLE_CHARGE_MASS is used, and whatever nonsense mass the code uses instead eventually leads to your crash.

This problem may go away if we use a subroutine like setup_particle_temperature instead. I'll add this to my fix list.

Cheers, Stuart

[DanRRRR]

Many thanks, Stuart. Hope to see this finally fixed ! I thought all these months that the errors was due to my additions in the subr. diagnostics until it really got on my nerves and finally this conflict was found :). By the way, i also use, or plan to use DEFINES += $(D)PARTICLE_ID. Is it OK with moving frame ?

[DanRRRR]

Stuart, Because I still not quite understand how exactly this moving frame functionality was implemented because never succeeded with it so far, so my naive question would be : shouldn't the added cells in front of laser just take exactly the same parameters as the last existing boundary cell before it? These last older cells were created self-consistently during all previous time history of the plasma. Hence may be things to add new cells will be potentially easier and the setup_particle_temperature would be not needed.?

[Status-Mirror]

Hey Dan,

Just to make sure we're on the same page, let me describe a simple 1D toy system

Say we have a window with 100 cells in $x$. In the first time-step, let this describe cells with indices 1-100. Now say we have a moving window, so that in the next timestep, these 100 cells should represent cells 2-101. It sounds like you're saying we should take all the particles in cell 100, and duplicate them in cell 101, to make sure the density and temperature are identical.

However, EPOCH has a greater functionality than this. You can define a density profile beyond the simulation window, and EPOCH will remember this and use it to set densities and temperatures in the new cells. For example, you could have a linear density ramp from cell 1-500, in a simulation window only 100 cells long. As the window moves towards cell 500, the densities in the window cells will keep increasing according to your distribution.

Hope this helps, Stuart

[DanRRRR]

Stuart, Definitely that's cool if it has this functionality! I know experimental situations when this is needed. Looking forward to try it.

Any progress and timeline with fixing the bug?

[DanRRRR]

Stuart, Can you please point me at the offending place with missing code?

[Status-Mirror]

Yes, sorry for the delay. You can fix the problem immediately by switching off the -DPER_PARTICLE_CHARGE_MASS flag. I would recommend this, as this option slows down the code - more particle variables mean slower MPI transfers between ranks. I'm not sure why EPOCH includes this compiler option at all.

But if you're determined to run with this flag, then the error is in the subroutine insert_particles in src/housekeeping/window.F90. When setting the momenta of loaded macro-particles, the code uses the function momentum_from_temperature, momentum_from_temperature_relativistic or sample_from_deck_expression. In all 3 cases, the code passes species%mass, but this is undefined for -DPER_PARTICLE_CHARGE_MASS.

Instead, the code should pass current%mass if running with -DPER_PARTICLE_CHARGE_MASS, or species%mass otherwise. It should be a simple fix, but it may be more complicated when you start to look at it.

[DanRRRR]

Stuart, I'd be glad not use this option if there were other way to output these 9 major variables of every PIC code: 3 particle coordinates, 3 momenta, weight, charge state and ID. I am surprised myself why it is needed when I do that in diagnostics subroutine because the code should operate same way with or without DPER_PARTICLE_CHARGE_MASS. I also need the mass of particle but that is just info for each given species not related to the code flow...All that is done just doe info purposes, to save the output on harddrives, not for any code necessity.

Essentially only mass and Z are missing. You probably much better know how to make a workaround respect to this, specifically ionization charge.

If i do not use this compilation switch DPER_PARTICLE_CHARGE_MASS i get the error at compile time at these two variables

particleMass = current%mas charge111 = current%charge

in this entire code i added to diagnostics sub respect to output of particles:

DO ispecies = 1, n_species

  if( species_list(ispecies)%attached_list%count.le.0) CYCLE

#ifdef PARTICLE_ID 
      CALL generate_particle_ids(species_list(ispecies)%attached_list) ! Stuart's suggestion
#endif
  current => species_list(ispecies)%attached_list%head
  particleMass = current%mass
  mc0_111      = current%mass * 3e8
  mi_to_me     = current%mass / 9.1094e-31

  DO WHILE(ASSOCIATED(current))

    x111 = current%part_pos(1)*100.
    y111 = current%part_pos(2)*100.
    z111 = current%part_pos(2)*100.

    px111 = current%part_p(1) /mc0_111
    py111 = current%part_p(2) /mc0_111
    pz111 = current%part_p(2) /mc0_111

    part_weight = current%weight
    charge111    = current%charge

    p_tot2 = px111**2+py111**2+pz111**2
    gamma_p = sqrt (1+p_tot2)
    energy111 = 0.511 *p_tot2/(1+gamma_p) * mi_to_me

    WRITE(iUnit, format111) x111, y111, z111....

  current =>current%next

  ENDDO 
ENDIF

[qiqi77346]

By the way here is how it looks before it crash. Further history is impossible. That plotting was100% done in Fortran

Hello, your picture is beautiful. What software was used to draw this picture?

[DanRRRR]

This is done with Clearwin+ of Silverfrost Fortran and their built-in OpenGL Fortran library

[qiqi77346]

This is done with Clearwin+ of Silverfrost Fortran and their built-in OpenGL Fortran library

OK, thank you.