How to restrict particle output vx in the dist_fn block?

[feelingxxxxxxx]

Dear Professor My goal is to obtain the velocity of particles in the x direction within the x range, as the number of particles is too large, I limit the range of vx. But the code I wrote reported an error, can you help me take a look? I am looking forward to your answer!

[Status-Mirror]

Hi @feelingxxxxxxx,

I've adapted the uniform load demo to show you how you can run a dist_fn. I'll include the input deck at the end.

I've made some changes to your block, these are:

• Swapped the order of your speeds in range2 - you need to provide the lower value first.
• Distribution blocks don't support dir_vx, so I've switched to relativistic dir_px.
• I've included a constant block, so you can enter speeds, and get momenta for the dist_fn block
• Because carbon and electrons have different mass, we can't use the same $p_x$ range for both, so I've separated them into two dist_fn blocks.

Our new input deck can run, and we get the following distribution function:

Here is the deck I used:

begin:control
    nx = 500
    ny = 500
    t_end = 1.0e-15
    x_min = -10.0e-6
    x_max = 15e-6
    y_min = 0
    y_max = 25e-6
    stdout_frequency = 100
    npart = 50 * nx * ny
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
    mass = 1.0
    charge = -1.0
    frac = 0.8
    density = 1.0e24
    temp_ev = 1e7
end:species

begin:species
    name = C6
    mass = 22033.0
    charge = 6.0
    frac = 0.2
    density = density(electron) / 6
    temp_ev = 1e7
end:species

begin:output
    dt_snapshot = t_end
    distribution_functions = always
end:output

begin:constant
  # Set speed range
  vx_low = 2.25e8
  vx_high = 0.999*c

  # Calculate Lorentz gamma for these speeds
  gamma_low = 1.0 / sqrt(1.0 - (vx_low/c)^2)
  gamma_high = 1.0 / sqrt(1.0 - (vx_high/c)^2)

  # Electron px range (p = gamma * mass * v)
  px_e_low = gamma_low * me * vx_low
  px_e_high = gamma_high * me * vx_high

  # Carbon px range
  px_c_low = gamma_low * 22033*me * vx_low
  px_c_high = gamma_high * 22033*me * vx_high
end:constant

begin:dist_fn
  name = electron_dist
  ndims = 2
  dumpmask = always
  direction1 = dir_x
  direction2 = dir_px
  range1 = (-0.8e-6, 0.3e-6)
  range2 = (-px_e_high, -px_e_low)
  resolution1 = 200
  resolution2 = 200
  include_species:electron
end:dist_fn

begin:dist_fn
  name = carbon_dist
  ndims = 2
  dumpmask = always
  direction1 = dir_x
  direction2 = dir_px
  range1 = (-0.8e-6, 0.3e-6)
  range2 = (-px_c_high, -px_c_low)
  resolution1 = 200
  resolution2 = 200
  include_species:C6
end:dist_fn

Hope this helps, Stuart

[bear-maker]

@Status-Mirror Dear Stuart, I am the one who asked this question and thank you very much for the answer that solved most of my questions. But right now I have one more question. So far we know that there are macroparticles in a cell, because the initial setup condition has a density gradient, which means that the macroparticles represent different densities. When I draw the density map, should I also output the weight file of each macro particle? If you need to output, can you please do a demonstration code? I look forward to your reply. Thank you very much

[Status-Mirror]

Hi @bear-maker,

You're right about the weights - in the default mode of EPOCH, each cell loads the same number of macro-particles, and density gradients are represented by differences in the macro-particle weights. If you want to plot the number density, EPOCH provides a number_density key in the output block. This will sum the macro-particle weights in each cell and divide by the cell volume, to give you an array showing the number density of real particles in each cell.

I've provided an example input deck to show this. It's the same as before, but I've removed the dist_fn and added a number_density output. By including the + species key, we can output the densities from specific species, along with the total density. I've also changed the densities in the species blocks to give an exponentially increasing number density. Running this deck allows me to plot:

begin:control
    nx = 500
    ny = 500
    t_end = 1.0e-15
    x_min = -10.0e-6
    x_max = 15e-6
    y_min = 0
    y_max = 25e-6
    stdout_frequency = 100
    npart = 50 * nx * ny
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
    mass = 1.0
    charge = -1.0
    frac = 0.8
    density = 1.0e24 * exp(x / 10.0e-6)
    temp_ev = 1e7
end:species

begin:species
    name = C6
    mass = 22033.0
    charge = 6.0
    frac = 0.2
    density = density(electron) / 6
    temp_ev = 1e7
end:species

begin:output
    dt_snapshot = t_end
    number_density = always + species
end:output

You can include both distribution_functions and number_density in the output block - I've just taken out the dist_fn to simplify the deck here.

Hope this helps, Stuart

[bear-maker]

@Status-Mirror Dear Stuart, Thank you very much for your answer, but currently I only want to draw the density of particles within a certain momentum range as in my first question, and do not need to draw the density of all particles. To solve this problem, my idea is to use the dist_fn module to restrict the output particles (momentum and position), and then sum up these restricted particles using the weights represented by macro particles. Another question is, how does epoch calculate particle density through macro particle weights in 2D or 1D? 2D divided by area? Divide 1D by 1? Looking forward to your answer. Your help is of great significance to me, and I wish you a happy life.

[Status-Mirror]

Ah I understand, you want to get number density information from a specific momentum range. The distribution functions used in our discussion will bin particles by their $x$ position and $p_x$ momentum, then sum the macro-particle weight in each bin.

Lets assume you've set up the momentum range of your distribution function to match the momentum range of interest, so you're only interested in the number density of particles plotted by the dist_fn. Since the dist_fn gives you the sum of macro-particle weights in each bin, if you sum all the bins which share the same $x$, then you'll get the total weight of particles within your momentum range in each cell. That looks like this in my first example:

The sum of macro-particle weights tells you the number of real particles in each cell. To get a number density, we divide by the cell volume. In the lower-dimension epoch codes, weights are assigned assuming the cell-size is 1m in the omitted dimensions (see discussion in our uniform load demo). In our example, the cell-sizes in $(x,y,z)$ are $(50\text{nm}, 50\text{nm}, 1\text{m})$ which gives a cell volume of $2.5\times 10^{-15} \text{m}^3$. But remember that the dist_fn only bins by $x$ and $p_x$, which means that a particle in a given bin can have any $y$ position. This means that a particle can be assigned to a dist_fn bin from a region with sides $(50\text{nm}, 25\mu\text{m}, 1\text{m})$, or volume $1.25\times 10^{-12} \text{m}^3$. When dividing the above plot by this volume, we get the number density:

If you were interested in a 2D density map, in both $x$ and $y$, you would need a 3D distribution function with axes $x,y,p_x$.

Hope this helps, Stuart