Questions about arbitrary distribution functions (species block)

[wnguyen1312]

How can one create a piece wise distribution in the species block and is the distribution supposed to be normalized (integrated to 1)? For example, I want to implement a "semi-quadratic" distribution in energy (kinetic energy) (I want the slope to be steeper on right side when E > 10 MeV). Is there a way to do this in EPOCH with the IF function? An example would be very helpful, thank you!

[Status-Mirror]

Hey @wnguyen1312,

This can be done using a similar method to that outlined in our momentum distribution demo. In my attached input.deck, I create a kinetic-energy variable part_ke, and I use this to construct two quardratic distribution functions, both of which peak at 1 for a kinetic energy KE_edge. The first function rises from 0 at 0 KE, and the second drops back to 0 at KE_max.

This example isn't perfect, as I haven't accounted for something important. EPOCH randomly draws px, py and pz values, compares it to the distribution function (assumed between 0 and 1), and draws a random number to see if we keep this px, py and pz set or not. However, this would only work for a KE distribution if all kinetic energies are equally likely to be drawn, but this is not the case! This is why our loaded distribution isn't the KE distribution we expected.

To do this properly, you would need to divide the distribution function by the probability of drawing a kinetic energy with a certain value from 3 uniformly sampled momenta.

Hopefully this is enough to get you started - if you're still having problems I might be able to spare some time later on to look over how to do this properly.

Cheers,
Stuart

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:constant
    part_p2 = px^2 + py^2 + pz^2
    part_ke = sqrt(part_p2*c^2 + me^2*c^4) - me*c^2
    KE_edge = 10.0e6 * qe
    KE_max = 15.0e6 * qe
end:constant

begin:constant
    p_max = sqrt((KE_max + me*c^2)^2 - me^2*c^4)/c
    b0 = 2 / KE_edge
    a0 = -b0 / (2 * KE_edge)
    c0 = 0
    b1 = 1/(0.5*KE_edge + KE_max^2/(2*KE_edge) - KE_max)
    a1 = - b1 / (2*KE_edge)
    c1 = b1 * (KE_max^2/(2*KE_edge) - KE_max)
    func_0_to_edge = (a0*part_ke^2 + b0 * part_ke + c0)
    func_edge_to_max = (a1*part_ke^2 + b1 * part_ke + c1)
end:constant

begin:species
    name = Electron_user
    mass = 1.0
    charge = -1.0
    npart = 5 * nx * ny
    density = 1.0e24
    dist_fn = if (part_ke lt KE_edge, func_0_to_edge, func_edge_to_max)
    dist_fn_px_range = (-p_max, p_max)
    dist_fn_py_range = (-p_max, p_max)
    dist_fn_pz_range = (-p_max, p_max)
end:species

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

[wnguyen1312]

Thank you, this works well for my purpose! I also managed to reproduce a similar distribution with the injector block, which seems to be more flexible