Add Electron Transport

[melekderman]

✅ (uploaded/inserted (almost) as is; I will check these later to make them compatible with the new data refactoring)
✅ The file element.py to object_/
✅ Electron-related items in reaction.py to object_/reaction.py (append or insert, not overwrite the entire file)
✅ The folder physics/electron/ to transport/physics/
👩🏻‍💻 Electron-related items in material.py to object_/material.py (append or insert, not overwrite the entire file)

[jpmorgan98]

Hey Melek!
I just pushed in some stuff to past the tests Fetch the new status and run the lintier and you should be good.

[ilhamv]

Thanks, Melek!

As we planned, you can go ahead and drop the remaining items (material.py, etc). Don't worry about the failing tests. They are supposed to fail, and then from there, we collaboratively fix and organize them here.

Some comments on some edits I made in my recent commits:

• I see that you are already aware of the type-hinting stuff needed for Numba mode and code factory. 👍🏽
• Note that now we have DistributionBase and DataBase, described in object_/distribution.py and object_/data.py, respectively. (See the changes in reaction.py)
• Note that I added the Element class and elements list in the simulation object description in object_/simulation.py. This is needed to get the new object discoverable by the code factory and simulation.

[ilhamv]

@melekderman - Can you please move the ISOTOPIC_ABUNDANCE data to object_/util.py?
Thanks!

[melekderman]

@ilhamv

Thank you so much for your review! I refactored the reactions.py, element.py, and physics/electron/native.py. I also moved the ISOTOPIC_ABUNDANCE data to object_/util.py, as you suggested. Your feedback would be greatly appreciated. 🙂 The remaining work is to insert material.py functions and to add the EDEP tally and related functions.

[ilhamv]

Looking good, @melekderman!

Can you also drop some examples in MCDC/examples.
For now, I'll go and make some adjustments to MCDC/mcdc/transport/physics/electron.

[ilhamv]

About ELECTRON_CUTOFF_ENERGY used in collision() in mcdc/transport/physics/electron/native.py.

1. Want to make sure: So, an electron with energy below the threshold will be absorbed if it collides? Or, when the said electron is created, it will be absorbed right away (does not have the chance to move to a collision site)?
2. Should the threshold be element dependent?

===

After learning more about this, I think it is best to make the cutoff adjustable per region. But we can do that later.
But still, I think the cut-off check should not be done in the collision() function, unless we have a mechanism that reduces the electron energy while free-streaming.

[ilhamv]

The neutron_production_xs in mcdc/transport/physics/neutron is only used for neutron k-eigenvalue simulations. So, I'm not sure if electron_production_xs would be needed..?

[ilhamv]

To have a consistent sub-function organization, can you please move compute_scattering_eta and sample_small_angle_mu_coulomb to after elastic_scattering? Thanks @melekderman !

[ilhamv]

Some questions on the elastic scattering.

What is the 0.999999 in mu_cut = float(large_angle.attrs.get("mu_cut", 0.999999)) for?

I see that you do some checks here:

    # If large-angle, xs from data table (EEDL MT=525)
    xs_large = 0.0
    xs_large_ID = int(reaction["xs_large_ID"])

    if xs_large_ID >= 0:
        xs_large = elastic_large_xs(E, reaction, mcdc, data)

    # Important to check because of numerical issues
    if xs_large < 0.0:
        xs_large = 0.0
    if xs_large > xs_total:
        xs_large = xs_total

Can we just manage the xs_large in mcdc/object_/reaction.py instead, so that xs_large_ID is always >=0 and we avoid any possible numerical issue?

[ilhamv]

In inonization, when you sample the subshell, you first calculate the total xs:

    # Sample subshell
    N = int(reaction["N_subshell"])
    xs_vals = np.empty(N, dtype=np.float64)

    total = 0.0
    for i in range(N):
        xs_sub_id = int(reaction["subshell_xs_IDs"][i])
        xs_sub_table = mcdc["data"][xs_sub_id]
        xs_sub_i = evaluate_data(E, xs_sub_table, mcdc, data)
        xs_vals[i] = xs_sub_i
        total += xs_sub_i

    xi = rng.lcg(particle_container) * total
    total_acc = 0.0
    chosen = 0
    for i in range(N):
        total_acc += xs_vals[i]
        if total_acc >= xi:
            chosen = i
            break

Should the total XS be the ionization XS? If not, should we normalize the subshell XS so that they sum up to the ionization XS? If we have them consistent with each other, we don't need to recalculate the total XS from the subshell XS.

[ilhamv]

In ionization we have:

    # Deposit all energy if below binding energy
    if E <= B:
        edep = E
        particle["alive"] = False
        particle["E"] = 0.0
        return edep

Isn't it impossible for a sampled subshell to have a binding energy that is larger than the electron energy? I'm wondering, if it does occur in the simulation, does that mean something is incorrect in the subshell cross-sections used for sampling?

[ilhamv]

Have you made sure whether the scattering kinematics data for elastic scattering and ionization are provided in the laboratory or the center-of-mass frame?

[ilhamv]

In Ionization we have

    primary_alive_after = True
    if E_out <= ELECTRON_CUTOFF_ENERGY:
        edep += E_out
        particle["E"] = 0.0
        particle["alive"] = False
        primary_alive_after = False

    if T_delta <= ELECTRON_CUTOFF_ENERGY:
        edep += T_delta
        return edep

I think if the delta electron got cut off, but the primary electron is still alive, we should not return as we still need to determine the primary electron direction.

[melekderman]

Have you made sure whether the scattering kinematics data for elastic scattering and ionization are provided in the laboratory or the center-of-mass frame?

I think it is in the LAB frame. The angular distribution is given only for large-angle elastic scattering, and recoil is neglected because it is assumed that the electron energies are much smaller than Mc2 (and there’s no LCT flag in the original EEDL data).

For ionization, the angular distribution is calculated inside MCDC using the formula defined in the LAB frame.

[melekderman]

In ionization we have:

    # Deposit all energy if below binding energy
    if E <= B:
        edep = E
        particle["alive"] = False
        particle["E"] = 0.0
        return edep

Isn't it impossible for a sampled subshell to have a binding energy that is larger than the electron energy? I'm wondering, if it does occur in the simulation, does that mean something is incorrect in the subshell cross-sections used for sampling?

Sure, that is not possible. It seems like all the extra checks I added while debugging to make sure energy is conserved are still in the code. I will remove them. Thank you!

[melekderman]

About ELECTRON_CUTOFF_ENERGY used in collision() in mcdc/transport/physics/electron/native.py.

1. Want to make sure: So, an electron with energy below the threshold will be absorbed if it collides? Or, when the said electron is created, it will be absorbed right away (does not have the chance to move to a collision site)?
2. Should the threshold be element dependent?

===

After learning more about this, I think it is best to make the cutoff adjustable per region. But we can do that later. But still, I think the cut-off check should not be done in the collision() function, unless we have a mechanism that reduces the electron energy while free-streaming.

I noticed that (I think I added an extra check again during debugging), I put a cutoff check both inside the collision and in the reactions. The one inside the collision must have been for debugging. And since the particle cannot actually travel that path, since the particle is absorbed immediately after it’s created in the reaction (or after the reaction for primary electron), it cannot reach the cutoff check in the collision routine. I'll delete it. Thank you!

[melekderman]

Some questions on the elastic scattering.

What is the 0.999999 in mu_cut = float(large_angle.attrs.get("mu_cut", 0.999999)) for?

I see that you do some checks here:

    # If large-angle, xs from data table (EEDL MT=525)
    xs_large = 0.0
    xs_large_ID = int(reaction["xs_large_ID"])

    if xs_large_ID >= 0:
        xs_large = elastic_large_xs(E, reaction, mcdc, data)

    # Important to check because of numerical issues
    if xs_large < 0.0:
        xs_large = 0.0
    if xs_large > xs_total:
        xs_large = xs_total

Can we just manage the xs_large in mcdc/object_/reaction.py instead, so that xs_large_ID is always >=0 and we avoid any possible numerical issue?

In EEDL2023, 0.999999 is the default value for mu_cut used to separate small-angle from large-angle elastic scattering. In this dataset (EEDL2023), angles with mu in [0.999999, 1.0] are considered small-angle, and anything below that is considered large-angle.

I only left it like this because I was considering adding mu_cut into the h5 file in case a future dataset uses a different cutoff value. We can definitely move this into mcdc/object_/reaction.py; I had thought about it, but I was not sure where it would fit best in the current code structure, so I left it there. I will change it. 👍