Plotting energy spectra at multiple times

[rashid611]

Having utilized VisIt for visualizing parameters like number density and electric field from the SDF files obtained through EPOCH simulations for Laser Wakefield Acceleration (LWFA), I'm interested in plotting or generating a contour representation that shows the relation between the energy of accelerated electrons and their number density. Could you please provide guidance on utilizing VisIt or matlab to create a plot or contour demonstrating the relationship between the energy of accelerated electrons and their number density from these SDF files?" Attached is an example plot....

[Status-Mirror]

I don't use VisIt myself, but I can give you some MATLAB instructions. I have slightly modified the input.deck file given in our momentum-distribution particle loading demo. We will need to output px, py, pz and the macro-particle weight:

begin:control
    nx = 500
    ny = 500
    t_end = 1.0e-15
    x_min = 0
    x_max = 25e-6
    y_min = 0
    y_max = 25e-6
    stdout_frequency = 100
end:control

begin:boundaries
    bc_x_min = open
    bc_x_max = open
    bc_y_min = open
    bc_y_max = open
end:boundaries

begin:species
    name = Electron_temp
    mass = 1.0
    charge = -1.0
    npart = 5 * nx * ny
    density = 1.0e24
    temp_ev = 1000
end:species

begin:output
    dt_snapshot = t_end
    px = always
    py = always
    pz = always
    weight = always
end:output

To plot the macro-particle energy distribution in MATLAB, copy over the files in epoch/SDF/Matlab to your working directory (or add them to the path). Once you have access to GetDataSDF(), you can run this MATLAB plotting script:

% energy_spectrum.m
% Plots the energy spectrum of all particles present in a species at a
% snapshot in time

c = 299792458;
qe = 1.60217663e-19;
me = 9.1093837e-31;

% Choose filename and species
data = GetDataSDF('0000.sdf');
species_name = "Electron_temp";
species_mass = me;

% Extract data
px = data.Particles.Px.(species_name).data;
py = data.Particles.Py.(species_name).data;
pz = data.Particles.Pz.(species_name).data;
weight = data.Particles.Weight.(species_name).data;
energy = sqrt((px.^2 + py.^2 + pz.^2)*c^2 + species_mass^2 * c^4);
energy_MeV = energy / 1.0e6 / qe;

% Bin macro-particles by energy
[~, bin_edges, ibin] = histcounts(energy_MeV);
bin_centres = 0.5*(bin_edges(2:end) + bin_edges(1:end-1));
bin_vals = 0*bin_centres;

% Number of real particles in each bin
for i = 1:length(ibin)
    bin_vals(ibin(i)) = bin_vals(ibin(i)) + weight(i);
end

% Normalise to dN/dE
bin_width = bin_edges(2) - bin_edges(1);
dN_dE = bin_vals / bin_width;

% Create plot
plot(bin_centres, dN_dE, "LineWidth", 2);
xlabel('Energy [MeV]');
ylabel('dN/dE [1/MeV]');
box on
grid on
ax = gca;
ax.FontSize = 16;
ax.YAxis.Scale = 'log';

This yields the following plot:

Hope this helps,
Stuart

[Status-Mirror]

That line suggests to me that the GetDataSDF() was unable to open the file, so data was not set-up correctly, and attempts to extract variables from data have failed.

Some things to check:

• Is GetDataSDF.m visible in your path? You'll find the function in epoch/SDF/Matlab.
• The function takes a file-name as an argument. Have you given it the file-name of a valid .sdf file?
• If your .sdf file is in a different folder, have you provided the path to the folder in the filename?

[DanRRRR]

Rashid611,
I see you are using moving window and claiming it runs OK. I still failed with moving window in my cases, EPOCH always crashes at the time of moving window start. Can you please post your input.deck which runs (Ideally in 2D for speed, or in 3D with smaller processor demands, or any simple demos which run) ? The one you have posted above has numerous typo problems by some reason.
What compiler + version do you use and what Makefile settings you have modified from their default values?

Here is example of error report

Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
Program received signal SIGSEGV: Segmentation fault - invalid memory reference.
Backtrace for this error:
#0  0x7f327692ead0 in ???
#1  0x7f327692dc35 in ???
#2  0x7f327652151f in ???
	at ./signal/../sysdeps/unix/sysv/linux/x86_64/libc_sigaction.c:0
#3  0x55578ccc45f1 in __particles_MOD_push_particles
	at src/include/triangle/e_part.inc:10
#4  0x55578cc96227 in pic
	at src/epoch2d.F90:216
#5  0x55578cc2ccbe in main
	at src/epoch2d.F90:33

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 376540 RUNNING AT dr-System-Product-Name
=   EXIT CODE: 139
=   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
YOUR APPLICATION TERMINATED WITH THE EXIT STRING: Segmentation fault (signal 11)
This typically refers to a problem with your application.
Please see the FAQ page for debugging suggestions

[rashid611]

yes , i am using moving window. here is the code

begin:constant

  lambda0        = 0.8*micron
  I              = 3.42e19  #a0=4
  waist          = 9.50*micron
  pulse_duration = 17*femto
  sigma_t        = pulse_duration/(sqrt(2*loge(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
  sim_end        = 2000*femto
  snap_freq      = 10*femto
end:constant


begin:control
  nx               = (x2-x1)*cpwlx/lambda0
  ny               = (y2-y1)*cpwly/lambda0
  nz               = (z2-z1)*cpwlz/lambda0
  npart            = part_per_cell*nx*ny*nz
  x_min            = x1
  x_max            = x2
  y_min            = y1
  y_max            = y2
  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)
  t_profile       = gauss(time,3*sigma_t,sigma_t)
end:laser

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

begin:output
  dt_snapshot     = snap_freq
  use_offset_grid = T


  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]

Thanks rashid611.
By the way do you have problems with this line t_profile = ....?
EPOCH refuses it in my case with absolutely no reason i could resolve.

[rashid611]

Hi Dear, I've checked the line " t_profile = --------- " in my EPOCH code, and I'm not encountering any issues. The code runs successfully without any errors on my end.