search

Warwick-Plasma / epoch

Particle-in-cell code for plasma physics simulations
https://epochpic.github.io
GNU General Public License v3.0
176 stars 56 forks source link

Ionisation in Input deck example doesn't work #537

Closed erny123 closed 1 year ago

erny123 commented 1 year ago

I have an input deck where I basically just copy the ionisation example from: https://epochpic.github.io/documentation/input_deck/input_deck_species.html

However I get the following errors: 


        d########P  d########b        .######b          d#######  d##P      d##P
       d########P  d###########    d###########     .##########  d##P      d##P 
      ----        ----     ----  -----     ----   -----         ----      -- P  
     d########P  d####,,,####P ####.      .#### d###P          d############P   
    d########P  d#########P   ####       .###P ####.          d############P    
   d##P        d##P           ####     d####   ####.         d##P      d##P     
  d########P  d##P            ###########P     ##########P  d##P      d##P      
 d########P  d##P              d######P          #######P  d##P      d##P       

[32m Welcome to EPOCH1D version 4.19.2   (commit v4.19.2-9-g724471f4-clean)

[39m0m *************************************************************
 The code was compiled with no compile time options
 *************************************************************
 Code is running on 40 processing elements

 Specify output directory
 Name of species 1 is Carbon
 Name of species 2 is Hydrogen
 Name of species 3 is Electron
 Name of species 4 is Carbon1
 Name of species 5 is Carbon2
 Name of species 6 is Carbon3
 Name of species 7 is Hydrogen1

 *** ERROR ***
 No mass specified for particle species "Carbon1"

 *** ERROR ***
 No mass specified for particle species "Carbon2"

 *** ERROR ***
 No mass specified for particle species "Carbon3"

 *** ERROR ***
 No mass specified for particle species "Hydrogen1"

 *** ERROR ***
 Not all required elements of input deck specified.
 Please fix input deck and rerun code
 *** WARNING ***
 The collision routine now uses the Nanbu-Perez scheme by default, rather than
 Sentoku-Kemp. This method is faster and does not appear to suffer from some
 unusual behaviour exhibited by Sentoku-Kemp under certain conditions.
 To revert to Sentoku-Kemp, specify "use_nanbu = F" in the collisions block.
 To remove this warning message, specify "use_nanbu = T".

I'm going to try adding each ionised species one by one. But honestly, ya'll should have better documentation on this.

And here is the input deck. 

# =================================================================================
begin:constant

 # --------------------------------------------------------------------------------
 # CONTROLLED CONSTANTS -----------------------------------------------------------

 # Plasma -------------------------------------------------------------------------
 Te_eV = 1.0                                # Electron temperature (eV)

 # Laser --------------------------------------------------------------------------
 l0 = 2.00 * micron                             # Wavelength (m)
 I0 = 10.0^14                               # Max single fiber laser intensity (W/cm^2)
 nc_ratio = 1.111                           # Ratio of critical den. to plasma den.
 n_fibers = 1.0 
 b = 1.0                            # Number of contributing laser fibers

 # Pulse --------------------------------------------------------------------------
 a  = 0.003                             # Envelope amp. defined as pulse edge

 # Run ----------------------------------------------------------------------------
 x_max_coeff = 100                      # Length of X domain in units of la
 n_output_files = 200                       # Total number of output files
 npart_cell_e = 40                      # Electrons per cell
 npart_cell_i = 40                          # Ions per cell

 # --------------------------------------------------------------------------------
 # Derived Constants --------------------------------------------------------------

 # Laser --------------------------------------------------------------------------
 w0 = 2.0 * pi * c / l0                     # Frequency (rad / s)
 #E0 = sqrt(2.0*I0*10^4/(epsilon0*c))       # Single-fiber amplitude (V / m) 
 #E0 = me*w0*c /qe * a
 #a0 = a #qe*E0/(me*w0*c)                       # Normalized single-fiber laser vector potential
 vg = c * sqrt(1.0 - 1.0/nc_ratio)          # Laser pulse group velocity (m / s)

 # Plasma -------------------------------------------------------------------------
 wpe = w0 / sqrt(nc_ratio)                  # Plasma freqeuncy (1 / s) 
 ne0 = me * epsilon0 / qe^2 * wpe^2         # Plasma density (1 / m^3)
 Te = Te_eV * ev                            # Electron temperature (J)
 tpe = 2.0 * pi / wpe                       # Plasma period (s)
 lp = 2.0 * pi * c / wpe                    # Wakefield wavelength (m)
 la = 2 * w0^2 * c / wpe^3                  # Acceleration saturation length (m)
 lD = sqrt(epsilon0 * Te / (ne0 * qe^2))    # Debye length (m)
 Ek_max = (2 * nc_ratio - 1) * me * c^2     # Theoretical maximum kinetic energy
 bg_max = sqrt((2.0 * nc_ratio)^2.0 - 1.0)  # Theoretical maximum beta-gamma

 # Second Laser --------------------------------------------- 
 wpump = w0 + wpe
 apump = a                     #* w0/wpump
 lpump = 2.0*pi *c/wpump
 Epump = me*wpump*c /qe * apump
 a0 = apump * l0/lpump*0.316227766  # 0.65465367        #seed 
 E0 = me*w0*c /qe * a0          #seed  amplitude

 # Pulse --------------------------------------------------------------------------
 Lt = pi * c / wpe                          # Laser pulse size (m)
 #Tt = Lt / c
 Tt = 1.0*pico                              # Laser pulse time (s)
 dt_pulse = Tt / (2.0 * sqrt(loge(2)))      # Pulse width (fwhm) (s)
 ETD = me * wpe * c / qe                    # Tajima-Dawson Field (V / m)
 #t0 = 0.5*wt_coeff*tpe*(loge(f0)/loge(a))^(1/(2*b))    # Laser start time (s)
 t0 = Tt*0.9
 #wt = wt_coeff*tpe/(2*(-loge(a))^(1/(2*b)))    # Laser pulse width (s)
 #wt = tpe / (4.0 * loge(1.0/a)^(1/(2*b)))  # Laser pulse width (s)
 #wt = wt_coeff*tpe/(2*(loge(1/a))^(1/(2*b)))# Laser pulse width (s)
 #wx = lp / (4.0 * (loge(1.0/a))^(1/(2*b))) # Laser pulse width (m)
 #x0 = x0_coeff * lp                            # Starting center position of pulse (m)

 # Run ----------------------------------------------------------------------------
 x_max0 = 200*micro #x_max_coeff * lp                   # Length of X domain (m)
 nx0 = ceil(x_max0 / lD/100)                    # Number of cells in X domain
 #nx0 = ceil(x_max_coeff * 100)         # Number of cells in X domain
 t_end0 = 2.5*pico  #1.0 * x_max0 / vg                  # Time limit (s)
 dt_snapshot0 = t_end0 / n_output_files     # Output interval (s)
end:constant
# =================================================================================

begin:control
 nx = nx0
 x_min = 0.0
 x_max = x_max0
 t_end = t_end0
 dt_multiplier = 0.95
 print_constants = T
 field_ionisation = T
 use_multiphoton = T
 physics_table_location = /data/homezvol0/ernestob/epoch/epoch1d/src/physics_packages/TABLES
 #restart_snapshot = 0001.sdf
end:control

begin:boundaries
 bc_x_min_field = simple_laser
 bc_x_min_particle = open
 bc_x_max_field = simple_outflow
 bc_x_max_particle = open
end:boundaries

# ----------------------------------------------------------------

begin:laser
 boundary = x_min
 amp = E0 
 #lambda = l0
 omega = w0
 t_profile = supergauss(time, t0, dt_pulse, 2*b)
end:laser

begin:laser
 boundary = x_min
 amp = Epump 
 #lambda = l0
 omega = wpump
 t_profile = supergauss(time, t0, dt_pulse, 2*b)
end:laser

# ----------------------------------------------------------------

begin:collisions
  use_collisions = T   
  collisional_ionisation = T 
   ci_n_step = 3
   coulomb_log = auto
   collide = all
end:collisions

begin:species
  name = Carbon
  charge = 0
  atomic_no = 6
  ionise = T
  ionise_limit = 3
  ionisation_electron_species = Electron
 density = if((x gt 20.0*micro) ,ne0*1.0/5.0,0)*if((x lt 160.0*micro) ,1,0)
 temp_ev = Te_eV

end:species   

begin:species
  name = Hydrogen
  charge = 0
  atomic_no = 1
  ionise = T
  ionise_limit = 1
  ionisation_electron_species = Electron
 density = if((x gt 20.0*micro) ,ne0*4.0/5.0,0)*if((x lt 160.0*micro) ,1,0)
 temp_ev = Te_eV

end:species  

begin:species
  name = Electron
  identify:electron
end:species

#begin:species
# name = Electrons
# charge = -1.0
# mass = 1.0
# density = if((x gt 10.0*micro) ,ne0,0)*if((x lt 20.0*micro) ,1,0)
# temp_ev = Te_eV
# npart_per_cell = npart_cell_e
#end:species

#begin:species
# name = Protons
# charge = 1.0
# mass = 12.0*1836.153
# density = if((x gt 10.0*micro) ,ne0,0)*if((x lt 20.0*micro) ,1,0)
# temp_ev = Te_eV
# npart_per_cell = npart_cell_i
# immobile = F
#end:species

begin:dist_fn
 name = Energy
 ndims = 1
 direction1 = dir_en
 range1 = (0,0)
 resolution1 = 100000
 include_species: Electron
end:dist_fn

begin:dist_fn
 name = x_px
 ndims = 1
 direction1 = dir_px
 range1 = (0,0)
 resolution1 = 100000
 include_species: Electron
end:dist_fn

begin:probe
 name = electron_probe
 include_species: Electron
 point = 195.0e-6
 normal = 1.0
 ek_min = 0.0
 ek_max = -1.0
end:probe

begin:probe
 name = hydrogen_probe
 include_species: Hydrogen1
 point = 195.0e-6
 normal = 1.0
 ek_min = 0.0
 ek_max = -1.0
end:probe

begin:probe
 name = carbon_probe
 include_species: Carbon1
 point = 195.0e-6
 normal = 1.0
 ek_min = 0.0
 ek_max = -1.0
end:probe

begin:probe
 name = carbon2_probe
 include_species: Carbon2
 point = 195.0e-6
 normal = 1.0
 ek_min = 0.0
 ek_max = -1.0
end:probe

begin:probe
 name = carbon3_probe
 include_species: Carbon3
 point = 195.0e-6
 normal = 1.0
 ek_min = 0.0
 ek_max = -1.0
end:probe

#begin:probe
# name = electron_accel_probe_min
# include_species: Electrons
# point = x_min
# normal = 0.0
# ek_min = 0.0
# ek_max = -1.0
#end:probe

begin:output
 name = rstart
 restartable = T
 dump_cycle = 1
 dt_snapshot = dt_snapshot0*5
 dump_first = F
 dump_last = T
end:output

begin:output
 name = dist
 file_prefix = dist
 dt_snapshot = dt_snapshot0
 distribution_functions = always + no_sum
 absorption = always + no_sum
 # distinguish energy in laser and particle
 total_energy_sum = always + no_sum
end:output

begin:output
 name = Probe
 file_prefix = probe
 dt_snapshot = dt_snapshot0
 particle_probes = always
end:output

begin:output
 name = Energy
 file_prefix = energy
 dt_snapshot = dt_snapshot0
 grid = always
 absorption = always
 total_energy_sum = always
 ekflux = always + species
 poynt_flux = always
end:output

begin:output
 name = Fields
 file_prefix = flds
 dt_snapshot = 1.0*femto
 grid=always
 ex = always + single
 ey = always + single
 #ez = always + single
 #bx = always + single
 #by = always + single
 #bz = always + single
end:output

#begin:output
# name = Density
# file_prefix = dens
# dt_snapshot = dt_snapshot0
# grid=always
# number_density=always + no_sum + single + species
# #charge_density = always + single
#end:output

begin:output
 name = Electron_momentum
 file_prefix = phas
 dt_snapshot = dt_snapshot0
 #particle_weight = single
 px = single
 ex = always + single
 ey = always + single
 number_density=always + no_sum + single + species
 #py = SubElec + single
 #pz = SubElec + single
 #particle_energy =  single
end:output

begin:output
 name = Electron_KEnergy
 file_prefix = KE
 dt_snapshot = dt_snapshot0
 particle_energy = single
 px =  single
 #py = SubElec + single

end:output
Status-Mirror commented 1 year ago

Hey @erny123,

I think the code you're referring to is in the Ionisation section of the species block documentation. I also think there's been some confusion here, as the input deck lines here only show where ionisation-relevant keys go. The above text reads: A basic example of using both ionisation mechanisms is given below, where non-relevant lines have been omitted. Perhaps we should include an ionisation example in our "basic examples" section to avoid future confusion.

There are lines which are "non-relevant" for ionisation, but still required for a standard EPOCH input deck. These include things like x_min, x_max, t_end in the control block, but there are also things needed by the species block.

In the attached input deck, you should include mass in your Carbon and Hydrogen blocks, such that they read:

begin:species
  name = Carbon
  charge = 0
  mass = 1836.2 * 12
  atomic_no = 6
  ionise = T
  ionise_limit = 3
  ionisation_electron_species = Electron
 density = if((x gt 20.0*micro) ,ne0*1.0/5.0,0)*if((x lt 160.0*micro) ,1,0)
 temp_ev = Te_eV

end:species   

begin:species
  name = Hydrogen
  charge = 0
  mass = 1836.2 * 1
  atomic_no = 1
  ionise = T
  ionise_limit = 1
  ionisation_electron_species = Electron
 density = if((x gt 20.0*micro) ,ne0*4.0/5.0,0)*if((x lt 160.0*micro) ,1,0)
 temp_ev = Te_eV
end:species

Note that mass must be specified for a species (unless using the identify alias like your electron species), and is given in units of the electron mass. The mass, charge and atomic number of ionised species will be generated automatically.

If you have any further problems, let me know Stuart