Moving Window crash with DEFINES += $(D)PER_PARTICLE_CHARGE_MASS
DanRRRR opened this issue · 12 comments
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
By the way here is how it looks before it crash. Further history is impossible.
That plotting was100% done in Fortran
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
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 ?
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.?
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
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
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?
Stuart,
Can you please point me at the offending place with missing code?
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.
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
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?
This is done with Clearwin+ of Silverfrost Fortran and their built-in OpenGL Fortran library
This is done with Clearwin+ of Silverfrost Fortran and their built-in OpenGL Fortran library
OK, thank you.