freeze particles

[HLQzZ]

Dear developer! I want to make the plasma move normally while freezing the electron beam, but it seems that the particle_tstart can only freeze all particles at the same time, is there any solution? Thanks! HLQ

[Status-Mirror]

Hi HLQ,

At present, there is no feature which would allow you to do this. However, you could achieve the same effect by inserting a line into the source-code.

Let's say you're running epoch2d. If you navigate to epoch2d/src/particles.F90, you'll find a subroutine called push_particles, and in this subroutine, you'll find the line:

      IF (species_list(ispecies)%immobile) CYCLE

This line is responsible for skipping the particle push for any species which has the immobile = T flag set in the epoch input deck. If you have a species named "electron_beam", and you only want it to move at 100 fs, then you could add a new line just below this:

      IF (species_list(ispecies)%immobile) CYCLE
      IF (TRIM(species_list(ispecies)%name) == 'electron_beam' .AND. time < 100.0e-15_num) CYCLE

Hence, when the code goes to move particles in the electron_beam species, and the simulation time is before 100e-15 seconds, the particle push would be skipped for that species as if it were an immobile species.

Hope this helps! If this doesn't work, let me know, Stuart

[HLQzZ]

Dear developer: This method works well in epoch, meeting my requirements. However, under the same compilation and conditions in epoch, it appears to have no effect. HLQ

[Status-Mirror]

Hi HLQ,

I'm not sure what you mean here, these statements seem contradictory.

Cheers, Stuart

[HLQzZ]

Sorry! I mean this method works well in epoch2d, meeting my requirements. However, under the same compilation and conditions in epoch3d, it appears to have no effect.

[Status-Mirror]

Ahh I see. EPOCH has a separate set of source codes for epoch1d, epoch2d and epoch3d. To change the 3D version, you should make the same edit to epoch3d/src/particles.F90.

Hope this helps, Stuart

[HLQzZ]

I changed the files in the 3D directory and recompiled, but it doesn't seem to work.

[Status-Mirror]

Could you send me your 3D input deck, and copy/paste the new line you added to epoch3d/src/particles.F90 here? It would also be useful if you could tell me the line-number for your new line.

Cheers, Stuart

[HLQzZ]

      IF (species_list(ispecies)%immobile) CYCLE
      IF (TRIM(species_list(ispecies)%name) == 'positron1' .AND. time < 400.0e-15_num) CYCLE
      IF (TRIM(species_list(ispecies)%name) == 'positron2' .AND. time < 400.0e-15_num) CYCLE
      IF (species_list(ispecies)%species_type == c_species_id_photon) THEN

begin:constant
  las_lambda       = 0.8 * micron
  las_omega        = 2.0 * pi * c / las_lambda
  laser_period     = las_lambda/c
  x_m              = 2*micron
  laser_k          = 2*pi/las_lambda 

  w0               = 2.37* micron   # spot size 2.5 *micron
  rayleigh_length  = pi*w0^2/las_lambda
  wz               = w0*sqrt(1+(x_m/rayleigh_length)^2)
  radius_curv      = x_m*(1.0+(rayleigh_length/x_m)^2)

  a0               = 40.0    
  Intensity        = (1.37*10^18)*a0^2/((las_lambda^2)*10^12)   #8.56e20[W/cm^2]
  n_crit           = 1.1*10^27/((las_lambda^2)*10^12)

  las_t_fwhm       = 8 * laser_period    # 12 periods

  r = sqrt(y^2 + z^2)
  E = 100 * mev  # 350MeV
  p = E/c   # E^2>>me^2*c^4
  E0=2*mev
  p0 = E0/c
  Ed=4.6884*mev
  pd = Ed/c

end:constant

begin:control
  nx = 4000
  ny = 240
  nz = 240

  # final time of simulation
  t_end            = 0.5 * pico

  # size of domain
  y_min            = -12 * micron
  y_max            =  12 * micron

  z_min            = -12 * micron
  z_max            =  12 * micron

  x_min            =  0.0 * micron
  x_max            =  200.0 * micron

  stdout_frequency =10

end:control

begin:boundaries
  bc_x_min_field     = simple_laser
  bc_x_max_field     = simple_laser
  bc_y_min_field     = simple_laser
  bc_y_max_field     = simple_outflow
  bc_z_min_field     = simple_laser
  bc_z_max_field     = simple_outflow
  bc_x_min_particle  = open
  bc_x_max_particle  = open
  bc_y_min_particle  = open
  bc_y_max_particle  = open
  bc_z_min_particle  = open
  bc_z_max_particle  = open
end:boundaries

begin:species
  name = electron
  charge = -1.0
  mass = 1.0
  npart_per_cell = 20
  tracer = F
  density = if((abs(r) lt 30 * micron) and (x gt 2.0 * micron) and (x lt 198.0 * micron), 0.099 * n_crit, 0.0 ) 
  temp_ev = 10
  identify:electron
end:species

begin:species
  name = positron1
  charge = 1.0
  mass = 1.0
  npart_per_cell =9
  tracer = F
  dist_fn= exp (-((px - p0)/pd)^2)
  dist_fn_px_range = (0,p)
  density = if((abs(r) lt 1.0 * micron) and (x gt 95 * micron) and (x lt 105 * micron),  1e10, 0.0) 
  temp_ev = 10

end:species

begin:species
  name = positron2
  charge = 1.0
  mass = 1.0
  npart_per_cell =9 
  tracer = F
  dist_fn=exp (-((-px - p0)/pd)^2)
  dist_fn_px_range = (-p,0)
  density = if((abs(r) lt 1.0 * micron) and (x gt 95 * micron) and (x lt 105 * micron),  1e10, 0.0) 
  temp_ev = 10

end:species

begin:species
  name = C
  charge = 6.0
  mass = 1836.0*12 
  npart_per_cell =9
  tracer = F
  density = if((abs(r) lt 30.0 * micron) and (x gt 2.0 * micron) and (x lt 198.0 * micron),  0.0165 * n_crit , 0.0) 
  temp_ev = 10
#  identify:proton
end:species

begin:species
  name             = Photon  
  charge           = 0.0
  npart            = 0 
  mass             = 0.0
  dumpmask         = always
  identify:photon
end:species

begin:collisions
  use_collisions   =    F
  coulomb_log      = auto
  collide          = all
end:collisions

begin:qed
   use_qed                = T 
   qed_start_time         = 10.0 * femto       #0.1 * pico
   produce_photons        = T
   use_radiation_reaction = T
   photon_energy_min      = 10.0 * kev
   produce_pairs          = F 
   photon_dynamics        = T  
end:qed

begin:laser
   boundary          = x_min
   intensity_w_cm2   = Intensity * w0/wz
   lambda            = las_lambda
   pol_angle         = 0.5*pi    #0.5*pi
   phase             = 0.0
   t_profile         = sin(0.5*pi*time/las_t_fwhm) * sin(0.5*pi*time/las_t_fwhm)
   t_start           = 0.0
   t_end             = 2*las_t_fwhm
   profile           = gauss(r,0.0, w0)
end:laser

begin:laser
   boundary          = x_max
   intensity_w_cm2   = Intensity * w0/wz 
   lambda            = las_lambda
   pol_angle         = 0.5*pi     #0.5*pi
   phase             = 0.0
   t_profile         = sin(0.5*pi*time/las_t_fwhm) * sin(0.5*pi*time/las_t_fwhm)
   t_start           = 0.0
   t_end             = 2*las_t_fwhm
   profile           = gauss(r,0.0, w0)
end:laser

begin:window   
   move_window = F
   window_v_x = c  
   window_start_time =  (x_max - x_min) / c 
   bc_x_min_after_move = simple_outflow  
   bc_x_max_after_move = simple_outflow
end:window

begin:subset
   name             = selected
   random_fraction  = 1.0
   include_species  = Photon
end:subset

begin:subset
   name             = selected1
   random_fraction  = 1.0
   include_species  = electron
   gamma_min        = 21.0
end:subset

begin:subset
  name             = selected2
  random_fraction  = 1.0
  include_species  = positron1
end:subset

begin:subset
  name             = selected3
  random_fraction  = 1.0
  include_species  = positron2
end:subset

begin:output
  # number of timesteps between output dumps
   dt_snapshot         = 10.0 * femto
   dt_average          = 1.0*laser_period

  # Properties at particle positions
    particles          =  selected2+ selected3
    vx                 =  never
    vy                 =  never
    vz                 =  never
    gamma              =  never #+ selected2+ selected3 
    ek                 =  never
    px                 =   selected2+ selected3
    py                 =   selected2+ selected3
    pz                 =   selected2+ selected3
    charge             =  never
    mass               =  never
    particle_weight    =  never
    id                 =  never 

  # Properties on grid
    grid               = always
    ex                 = never
    ey                 = always
    ez                 = always
    bx                 = never
    by                 = always
    bz                 = always
    jx                 = never
    jy                 = never
    jz                 = never
    ekbar              = never + species + no_sum
    charge_density     = never
    number_density     = always + species + no_sum
    temperature        = never + species
    total_energy_sum   = never
    distribution_functions = never
end:output

[Status-Mirror]

So are you finding that particles in the positron1 and positron2 species are moving before 400 fs?

I'm not sure what could be going wrong here. Maybe the TRIM(species_list(ispecies)%name) == 'positron1' test is failing? If you added a new line:

      IF (species_list(ispecies)%immobile) CYCLE
      IF (TRIM(species_list(ispecies)%name) == 'positron1' .AND. time < 400.0e-15_num) PRINT*, 'Positron1 frozen'
      IF (TRIM(species_list(ispecies)%name) == 'positron1' .AND. time < 400.0e-15_num) CYCLE
      IF (TRIM(species_list(ispecies)%name) == 'positron2' .AND. time < 400.0e-15_num) CYCLE

that would test whether the CYCLE was being triggered. I'm wondering if maybe some whitespace around species_list(ispecies)%name could make the name check fail.

[Status-Mirror]

I'm going to mark this issue as closed for now, as I believe the code I sent should perform as expected. Feel free to continue the thread if you have any further comments!

Cheers, Stuart