pyro-kinetics / pyrokinetics

Python library to run and analyse gyrokinetics simulations
https://pyrokinetics.readthedocs.io/en/latest/#
GNU Lesser General Public License v3.0
25 stars 6 forks source link

Feature request: Clever Ion naming #335

Open bpatel2107 opened 6 months ago

bpatel2107 commented 6 months ago

We should mostly be able to determine which ion we are dealing with from the charge/mass alone. This get complicated with partial ionisation states but maybe we can handle this in a clever way. See conversation in #315

@arkabokshi FYI

ZedThree commented 6 months ago

Might be worth looking at PlasmaPy's Particle module, or the mendeleev package to see if there's something that already does this?

ZedThree commented 6 months ago

PlasmaPy sort of:

import plasmapy.particles as pp

particle = pp.ionic_symbol(1, mass_numb=3, Z=1)
print(particle)
# 'T 1+'

pp.element_name(particle)
# 'hydrogen'

Note that it does need to know the atomic number, and I don't see an obvious way to get the name "tritium" although it does understand that element_name("tritium") == "hydrogen".

I think without the atomic number, this can only ever be a "best guess"

arkabokshi commented 6 months ago

Thanks Bhavin. As we have been discussing on the sidelines, these are my current thoughts on the main issues:

  1. How do we distinguish between He-ash and fast-alphas when converting from one code to another? Current thought is to label the species with the higher temp as fast (and the other is then ash)?

  2. About naming the species, I think it may be possible to identify all the tokamak relevant species purely on the basis of the mass number (with may be the exception of He3 and Tritium). Consider the following code which tries to open up gaps in the periodic table based on 2 exclusionary criteria:

import numpy as np
from periodictable import elements
likely_elements = ['H','He','Li','Be','C','N','O','Ar','Ne','W']
for Z,element in enumerate(elements):
    if str(element) in likely_elements:
        for isotope in element:
            if isotope.abundance > 0.0:
                print(f'{element},{isotope.mass:.0f}')

With this the only possible candidates are below and we needn't use another criterion such as Z which could complicate things due to partial ionization etc.

H,2
He,3
He,4
Li,6
Li,7
Be,9
C,12
C,13
N,14
N,15
O,16
O,17
O,18
Ne,20
Ne,21
Ne,22
Ar,36
Ar,38
Ar,40
W,180
W,182
W,183
W,184
W,186

@fcasson Thoughts?

bpatel2107 commented 6 months ago

The choice of fast vs thermal should be made upon the data we load from (say a TRANSP run where a fast species would have a different name). Otherwise when we load from a gk file we can look at temperature ratios or in some cases the code will specify the distribution function for the species.

I've tweaked the logic above and found something that is pretty good

import numpy as np
from periodictable import elements

likely_symbols = ["H", "He", "Li", "Be", "C", "N", "O", "Ar", "Ne", "W"]

def determine_isotope(charge, mass):
    likely_match = {}
    for element in elements:
        for isotope in element:
            # Need to loop through all elements first to find perfect match
            if np.isclose(isotope.mass, mass, atol=0.1):
                if np.isclose(isotope.number, charge, rtol=0.01):
                    return f"{isotope.name}{round(isotope.mass)}"
                else:
                    likely_match[isotope.symbol] = f"{isotope.name}{round(isotope.mass)}"

    for symbol in likely_symbols:
        if symbol in likely_match.keys():
            return likely_match[symbol]

    return "ion"

charge = 6
mass = 12

isotope = determine_isotope(charge, mass)
print(f"Isotope: {isotope}")

charge = 1
mass = 3

isotope = determine_isotope(charge, mass)
print(f"Isotope: {isotope}")

charge = 41
mass = 184

isotope = determine_isotope(charge, mass)
print(f"Isotope: {isotope}")

charge = 43
mass = 184

isotope = determine_isotope(charge, mass)
print(f"Isotope: {isotope}")

This gets the following

Isotope: carbon12
Isotope: tritium3
Isotope: tungsten184
Isotope: tungsten184

First we check for a match on the charge and mass and if nothing is found then we loop through likely elements and find a matching mass.

This could be implemented in get_local_species for each gk code (or in normalise which always gets called)

fcasson commented 6 months ago

Overall I'm not sure "naming" ions a GK context is helpful, because of the ambiguities. Why do you want to do it? Mass and charge are all we should ever need.

Some thoughts

  1. Ions can be partially ionised, which complicates attempts to identify element name by charge alone, since it is only the partial charge that matters in GK context, but the atomic charge that matters in naming (periodic table) context.
  2. Atoms have different isotopes, which complicates attempts to identify by mass number alone.
  3. Often in simulation, we want to create ficticious species, like tokamakium (z=1, m=2.5), or "lumped heavy impurity", "lumped light impurity", or average over charge states, or average over isotopes.

All in all, life will be a lot simpler if you can avoid doing this.

Related: https://jira.iter.org/browse/IMAS-4900

d7919 commented 6 months ago

I guess the motivation is that GENE v3 seems to run into some challenges if one has multiple species with the same "name" input. Currently we just add a counter to "ion", which I think is probably fine. The approach proposed here might make things appear more intuitive but I think I agree that it can start to open up additional complexity and potential to be misleading.

fcasson commented 6 months ago

Here is the list of isotopes one can expect to see in a tokamak (with any charge up to their atomic number)

H,(D),(T)
(He-3),He-4,
(Li-6),Li-7

The elements below can occur with any isotope from their naturally occurring isotopic composition

Be,B,C,N,O,Ne  
Al,Ar
Ti,V,Cr,Mn,Fe,Co,Ni,Cu,Zn,Kr
Ni,Mo,Ag,Sn,Xe
Hf,W,Pt,Au,Pb
bpatel2107 commented 6 months ago

Maybe the solution is to have an optional method that tries to determine the ion but by default sticks to ion1, ion2 etc.

The code I wrote above does check all isotopes for each element and works well for fully ionised species and is pretty good for partially ionised species if the list of likely_species is relatively small. Would need to see how it does if we extend it to the list above. If it can't find a match then it could leave it be which seems safe.

fcasson commented 6 months ago

Most likely species, based on JET / STEP / ITER / MAST / CMOD / AUG

Assume full ionisation [assumption getting weaker...]

H,D,T,He-3,He,Li-7,Be,B,C,[N,O,Ne]

Likely partial ionisation

Ar,Fe,Ni,Kr,Mo,Xe,W

If you want to get really clever you can use lookup table for average ionisation state at each temp. I have examples of this for W at least .

This paper is a nice overview. https://iopscience.iop.org/article/10.1088/1741-4326/ab0384/meta

arkabokshi commented 6 months ago

Merged impurities could indeed be a problem, and if we find such an exotic species, we could default to ion1, ion2 etc.

Consider the histogram which has a bin width of 1 amu, made with the most likely species and their isotopes from JET / STEP / ITER etc above. isotope

The conclusion is that if can identify / define the atomic mass with a precision of 1 amu, only two isotopes are indistinguishable on 3 occasions. If the atomic mass is specified with a precision of 0.001 amu, all naturally occurring isotopes possible in the tokamak are distinguishable.

My conclusion is that the mass number is a more robust way of determining the isotope and we can rely less on the charge state (secondary criteria)?