Reading in the Data¶
webnucleo data files are in either in XML or HDF5 format. wnutils routines can read either format.
In the following tutorials, you will enter Python commands. In your terminal, type python (or python3, or, perhaps, python2.x or python3.x, depending on your version). You will see something like:
Python 3.6.5 (default, Mar 29 2018, 15:38:28)
[GCC 4.2.1 Compatible Apple LLVM 7.3.0 (clang-703.0.31)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>>
The >>> is the interactive python prompt (it will typically show up in the tutorials in red). You type your commands at this prompt. To exit python, type:
>>> exit()
and hit enter.
XML¶
The format of webnucleo XML files is described in the libnucnet technical
report XML Input to libnucnet, available at the
libnucnet Home page.
The class wnutils.xml.Xml
has methods that read these XML files.
To begin, import the namespace by typing:
>>> import wnutils.xml as wx
To illustrate the use of wnutils.xml routines, use the files my_output1.xml and my_output2.xml, which you should have downloaded according to the data tutorial. For each file, create an Xml object. For example, type:
>>> my_xml = wx.Xml('my_output1.xml')
Read the nuclide data.¶
You can retrieve the nuclide data in the webnucleo XML file by typing:
>>> nuclides = my_xml.get_nuclide_data()
This returns a dictionary of data with the key being the nuclide name. You may print out all the data for a specific nuclide, say, o16, by typing:
>>> print(nuclides['o16'])
Or, to get specific data, try typing:
>>> print('The mass excess in MeV of o16 is', nuclides['o16']['mass excess'])
It is possible to use an XPath expression to select out only certain nuclides. For example, to get the data for nitrogen isotopes only, type:
>>> n = my_xml.get_nuclide_data(nuc_xpath='[z = 7]')
To confirm that you only retrieved nitrogen data, type:
>>> for isotope in n:
... print(isotope, ':', 'Z =', n[isotope]['z'], 'A =', n[isotope]['a'])
...
Partition function data for the nuclei
are stored in two numpy.array
objects. The first
array, with key t9, gives the temperature points (in billions of k) at which
the partition function is evaluated. The second array, with key
partf, gives the partition function G evaluated at each of the
temperature points.
To see how this works, try printing out the partition function
for one of the iron isotopes, say, fe56. Begin
by extracting the data for the iron isotopes by typing:
>>> fe = my_xml.get_nuclide_data(nuc_xpath='[z = 26]')
Then print out the partition function G as a function of t9 by typing:
>>> sp = 'fe56'
>>> for i in range(len(fe[sp]['t9'])):
... print('t9 = ', fe[sp]['t9'][i], 'G(t9) = ', fe[sp]['partf'][i])
...
Read the network limits.¶
It is often useful to know the limits of the network that comprises the nuclei in the nuclear data collection. To get this information, type:
>>> lim = my_xml.get_network_limits()
This returns a dict
of numpy.array
objects. The array
retrieved with key z gives the atomic numbers. The array retrieved with
key n_min gives the smallest neutron number present for the corresponding
atomic number, while the array retrieved with key n_max gives the largest
neutron number present for the corresponding atomic number. You can print
out the retrieved data by typing:
>>> for z in range(len(lim['z'])):
... print('Z =', z, ': N_min =', lim['n_min'][z], ', N_max =', lim['n_max'][z])
...
You can retrieve a subnetwork with an XPath expression. For example, you can type:
>>> lim = my_xml.get_network_limits(nuc_xpath = '[z <= 5 or z >= 25]')
Now print out the data:
>>> for i in range(len(lim['z'])):
... print('Z =', lim['z'][i], ': N_min =', lim['n_min'][i], ', N_max =', lim['n_max'][i])
...
Read the reaction data.¶
You can retrieve the reaction data in the webnucleo XML file by typing:
>>> reactions = my_xml.get_reaction_data()
This returns a dictionary with the key being the reaction string and each
value being a Reaction
. To see a list of the reactions, type:
>>> for r in reactions:
... print(r)
...
You can use an XPath expression to select the reactions. For example, you can type:
>>> reactions = my_xml.get_reaction_data('[count(non_smoker_fit) = 1]')
Since the reaction data include the reaction type, you can confirm your request by typing:
>>> for r in reactions:
... data = reactions[r].get_data()
... print(r, ': type is', data['type'])
...
You may choose a particular reaction from the dictionary by typing, for example:
>>> reac = reactions['n + he4 + he4 -> be9 + gamma']
It is then possible to retrieve the reactants, nuclide_reactants, products, nuclide_products, the reaction string, and code giving the source by typing:
>>> print(reac.reactants)
>>> print(reac.nuclide_reactants)
>>> print(reac.products)
>>> print(reac.nuclide_products)
>>> print(reac.get_string())
>>> print(reac.source)
You can also compute the rate for the reaction (among interacting multiplets and assuming one of the standard rate forms single_rate, rate_table, or non_smoker_fit) at a variety of temperatures by typing:
>>> import numpy as np
>>> t9s = np.power(10., np.linspace(-2,1))
>>> for t9 in t9s:
... print(t9, reac.compute_rate(t9))
...
To compute the rate for user-defined rate functions, each defined with a user_rate key, first write a python routine for each rate function, then bind any data to that function (which must still take t9 as an argument), and then create a dictionary of the functions associated with each key. Pass that dictionary into the compute_rate method with the keyword user_funcs.
Read all properties in a zone.¶
In a webnucleo XML file,
a zone is a collection of the mutable quantities during a network
calculation. For a single-zone network calculation, a zone is often a
time step in the calculation. The zone will contain mass fractions of
the network species at the time step to which the zone corresponds and
properties, which can be any quantity, such as time, temperature, or
density. The properties themselves have a name and up to two tags,
called tag1 and tag2. If the property only has a name, it can
be retrieved by a str
. If the property has tags, the identifier
for the property is a tuple
of up to three strings, namely,
the name, tag1, and tag2.
To retrieve all the properties of a given zone, say, the 10th zone, type:
>>> props = my_xml.get_all_properties_for_zone('[position() = 10]')
Now you can print out the properties and their values in this zone by typing:
>>> for prop in props:
... print(str(prop).rjust(25), ':', props[prop])
...
Notice the conversion to str
to print out the
(‘exposure’, ‘n’) tuple correctly.
Read properties in all zones.¶
You can retrieve selected properties in all zones. For the present example, you retrieve the time, t9 (temperature in billions of Kelvins), and rho (mass density in g/cc) by typing:
>>> props = my_xml.get_properties( ['time','t9','rho'] )
The properties are returned in the dictionary props. Each dictionary element is a list of strings giving the properties in the zones. To see this, type:
>>> print(props['time'])
This prints all the times. Print the first time entry by typing:
>>> print(props['time'][0])
To see the types, print:
>>> type(props)
which shows that it is a hash (dict
). Next, type:
>>> type(props['time'])
which shows that each dictionary entry is a list
. Next, type:
>>> type(props['time'][0])
which shows each list entry is a str
.
To retrieve properties with tags, you need to enter the appropriate tuple. For example, type:
>>> props = my_xml.get_properties(['time', ('exposure', 'n')])
To print out the exposures, type:
>>> for i in range(len(props[('exposure', 'n')])):
... print('time:', props['time'][i], 'exposure:', props[('exposure', 'n')][i])
...
Read properties of selected zones.¶
You can select out the zones whose properties you wish to read by using an XPath expression. For example, you can retrieve the time, t9, and rho properties, as in the above example, but only for the last 10 zones. Type:
>>> props = my_xml.get_properties(
... ['time','t9','rho'], zone_xpath='[position() > last() - 10]'
... )
You can print the zone properties, for example, by typing:
>>> print(props['t9'])
Confirm that there are only the properties for 10 zones by typing:
>>> print(len(props['t9']))
Read zone properties as floats.¶
Properties are by default strings. When you wish to manipulate them
(for example, to plot them), you want
them to be float
objects. You can retrieve them as floats by typing:
>>> props = my_xml.get_properties_as_floats( ['time','t9','rho'] )
The returned hash has entries that are numpy.array
, which you confirm
with:
>>> type(props['rho'])
You can confirm that the array entries are floats by typing:
>>> type(props['rho'][0])
You can print out the entries by typing:
>>> for i in range(len(props['time'])):
... print(
... 'Zone = {0:d} time(s) = {1:.2e} t9 = {2:.2f} rho(g/cc) = {3:.2e}'.format(
... i, props['time'][i], props['t9'][i], props['rho'][i]
... )
... )
...
This will output the time, temperature (in billions of K), and mass density (in g/cc) in all zones (time steps).
Read mass fractions in zones.¶
You can retrieve the mass fractions in zones. For example, to get the mass fractions of o16, si28, and s36, type:
>>> x = my_xml.get_mass_fractions(['o16','si28','s36'])
The method returns a dict
of numpy.array
. Each array element
is a float
. You can print the mass fraction of silicon-28 in all
zones by typing:
>>> print(x['si28'])
The method also accepts the zone_xpath keyword to select specific zones. For example, to retrieve the mass fraction in the first 10 zones, type:
>>> x = my_xml.get_mass_fractions(
... ['o16','si28','s36'], zone_xpath='[position() <= 10]'
... )
Read all abundances in zones.¶
You can retrieve abundances in the zones as a function of atomic and
neutron number. The retrieved data are stored in a three-dimensional
numpy.array
. The first index gives the zone, the second
gives the atomic number, and the third gives the neutron number. The
array value is the abundance (per nucleon). Zones can be selected by
XPath.
To see how this works, retrieve the abundances in all zones by typing:
>>> abunds = my_xml.get_all_abundances_in_zones()
Now print out the abundances in the 50th zone (remember the zero-indexing) by typing:
>>> for z in range(abunds.shape[1]):
... for n in range(abunds.shape[2]):
... print('Z =', z, ', N =', n, ', Y(Z,N) =', abunds[49,z,n])
...
You could do the same by typing:
>>> abunds = my_xml.get_all_abundances_in_zones(zone_xpath='[position() = 50]')
>>> for z in range(abunds.shape[1]):
... for n in range(abunds.shape[2]):
... print('Z =', z, ', N =', n, ', Y(Z,N) =', abunds[0,z,n])
...
This is because the XPath selects only one zone, which will have index 0 in the retrieved data.
Retrieve abundances summed over nucleon number in zones.¶
It is often convenient to retrieve the abundances of the nuclei in a network file summed over proton number (z), neutron number (n), or mass number (a). To do so, type:
>>> y = my_xml.get_abundances_vs_nucleon_number()
This returns a two-dimensional numpy.array
in which the first
index gives the zone and the second the mass number a. To print out
the abundances versus mass number in the eighth zone, type:
>>> for i in range(y.shape[1]):
... print('A:', i, 'Y(A):', y[7,i])
...
To retrieve the abundances summed over atomic (proton) number (z), use the keyword nucleon:
>>> y = my_xml.get_abundances_vs_nucleon_number(nucleon='z')
To retrieve the abundances in particular zones, for example, in the last 10 zones, use an XPath expression:
>>> y = my_xml.get_abundances_vs_nucleon_number(nucleon='n', zone_xpath='[position() > last() - 10]')
Retrieve abundances for a chain of species.¶
To retrieve the abundances for a set of isotopes or isotones, use the method to get chain abundances. For example, to retrieve the isotopic abundances for Z = 30 for all timesteps, type:
>>> n, y = my_xml.get_chain_abundances(('z', 30))
The method returns a tuple
with the first element being an array of
neutron numbers for the isotopes and the second element being a two dimensional
numpy.array
with the abundances for each step. To print the isotopic
abundances in the final step, type:
>>> step = y.shape[0] - 1
>>> for i in range(y.shape[1]):
... print('N =', n[i], ', Y[N] =', y[step, i])
...
To return the isotonic abundances for N = 25 in the first thirty timesteps, type:
>>> z, y = my_xml.get_chain_abundances(('n', 25), zone_xpath="[position() <= 30]")
To return the same isotonic abundances, but as a function of the mass number, set the keyword variable vs_A to True:
>>> a, y = my_xml.get_chain_abundances(('n', 25), zone_xpath="[position() <= 30]", vs_A=True)
To print these abundances in the thirtieth step, type:
>>> step = y.shape[0] - 1
>>> for i in range(y.shape[1]):
... print('A =', a[i], ', Y[A] =', y[step, i])
...
Multi_XML¶
The wnutils.multi_xml.Multi_Xml
class allows you to access and plot data
from multiple webnucleo XML files. First import the namespace by typing:
>>> import wnutils.multi_xml as mx
Then create a class instance from a list
of XML files.
For this tutorial, type
>>> my_multi_xml = mx.Multi_Xml(['my_output1.xml','my_output2.xml'])
Methods allow you to access or plot data from the files.
Read data from the individual XML instances.¶
To retrieve the individual XML instances from a Multi_Xml instance, type:
>>> xmls = my_multi_xml.get_xml()
To retrieve the original file names, type:
>>> files = my_multi_xml.get_files()
Of course the number of XML instances must equal the number of files. To confirm, type:
>>> print(len(xmls) == len(files))
Use the methods on the individual instances. For example, type:
>>> for i in range(len(xmls)):
... props = xmls[i].get_properties(['time'])
... print(files[i],'has',len(props['time']),'zones.')
...
H5¶
Methods that read webnucleo HDF5 files are in the namespace
wnutils.h5. The class that contains these methods is
wnutils.h5.H5
. Begin by importing the namespace by typing:
>>> import wnutils.h5 as w5
Then create an object for your file my_output1.h5 (which you already downloaded according to the instructions in the data tutorial) by typing:
>>> my_h5 = w5.H5('my_output1.h5')
Read the nuclide data.¶
The nuclide data are in a group of their own in the file. To retrieve the
data (as a dict
of dict
with the nuclide names as the top-level
dictionary keys), type:
>>> nuclides = my_h5.get_nuclide_data()
Print out the data for, say, o16, by typing:
>>> print(nuclides['o16'])
Print out the mass excess and spin for all species by typing:
>>> for nuclide in nuclides:
... print(nuclide, nuclides[nuclide]['mass excess'], nuclides[nuclide]['spin'])
...
Read the names of the iterable groups.¶
Iterable groups are the groups in the HDF5 file that typically represent
timesteps (that is, the groups that are not the nuclide data group).
To retrieve their names (as a list
of str
), type:
>>> groups = my_h5.get_iterable_groups()
Print them out by typing:
>>> for group in groups:
... print(group)
...
Read the zone labels for a group.¶
In a webnucleo HDF5 file, a zone is contained in a group and typically
represents a spatial region. Zones are specified by three labels, which
we denote by a tuple
. To retrieve and print out the labels for a given
group, say, Step 00010, type:
>>> labels = my_h5.get_zone_labels_for_group('Step 00010')
>>> for i in range(len(labels)):
... print('Zone',i,'has label',labels[i])
...
Read all properties in a zone for a group.¶
To retrieve all the properties from a zone in a group, type, for example:
>>> zone = ('2','0','0')
>>> props = my_h5.get_group_zone_properties('Step 00010', zone)
You can print those properties out by typing:
>>> for prop in props:
... print(str(prop).rjust(25), ':', props[prop])
...
Read properties in all zones for a group.¶
It is possible to retrieve the properties in all zones for a group as
as dict
of list
. Each list entry is a str
. For example,
to retrieve and print the properties time, t9, and rho
in all zones for a given group, say, Step 00024, type:
>>> p = ['time','t9','rho']
>>> props = my_h5.get_group_properties_in_zones('Step 00024',p)
>>> labels = my_h5.get_zone_labels_for_group('Step 00024')
>>> for i in range(len(labels)):
... print('In',labels[i],'time=',props['time'][i],'t9=',props['t9'][i],'rho=',props['rho'][i])
...
Read properties in all zones for a group as floats.¶
It is often desirable to retrieve the properties in zones for a group as floats. For example, one may again retrieve time, t9, and rho from Step 00024 but, this time, as floats by typing:
>>> p = ['time','t9','rho']
>>> props = my_h5.get_group_properties_in_zones_as_floats('Step 00024',p)
>>> type(props['time'])
>>> type(props['time'][0])
Read mass fractions in all zones for a group.¶
You can read all the mass fractions in all the zones for a given group. For a group Step 00021, type:
>>> x = my_h5.get_group_mass_fractions('Step 00021')
The array x is a 2d HDF5 h5py:Dataset
. The first index gives the zone
and the second the species. To print out the mass fraction of ne20 in all
the zones, type:
>>> i_ne20 = (my_h5.get_nuclide_data())['ne20']['index']
>>> labels = my_h5.get_zone_labels_for_group('Step 00021')
>>> for i in range(x.shape[0]):
... print('Zone',labels[i],'has X(ne20) =',x[i,i_ne20])
...
Read properties of a zone in the groups.¶
It is possible to retrieve properties from a given zone in all groups. To retrieve the properties time, t9, and rho in all groups for the zone with labels 1, 0, 0, type:
>>> zone = ('1','0','0')
>>> props = my_h5.get_zone_properties_in_groups(zone, ['time','t9','rho'])
This returns a dict
of list
of str
.
To print the properties out in the groups, type:
>>> groups = my_h5.get_iterable_groups()
>>> for i in range(len(groups)):
... print(
... groups[i], ': ', props['time'][i], props['t9'][i], props['rho'][i]
... )
...
Read properties of a zone in the groups as floats.¶
One often wants the properties of a zone in the groups as floats. To retrieve the properties time, t9, and rho in all group for the zone with labels 1, 0, 0, type:
>>> zone = ('1','0','0')
>>> props = my_h5.get_zone_properties_in_groups_as_floats(zone, ['time','t9','rho'])
This returns a dict
of numpy.array
. Each array entry is a
float
. To print the properties out in the groups, type:
>>> groups = my_h5.get_iterable_groups()
>>> for i in range(len(groups)):
... print(
... '{0:s}: time(s) = {1:.2e} t9 = {2:.2f} rho(g/cc) = {3:.2e}'.format(
... groups[i], props['time'][i], props['t9'][i], props['rho'][i]
... )
... )
...
Read mass fractions in a zone in the groups.¶
You can retrieve the mass fractions of specific species for a given zone in all the iterable groups. For example, to retrieve o16, o17, and o18 in the zone with labels 1, 0, 0, type:
>>> species = ['o16','o17','o18']
>>> zone = ('1','0','0')
>>> x = my_h5.get_zone_mass_fractions_in_groups( zone, species )
To print them out, you can now type:
>>> groups = my_h5.get_iterable_groups()
>>> for i in range(len(groups)):
... print(groups[i],':','X(o16)=',x['o16'][i],'X(o17)=',x['o17'][i],'X(o18)=',x['o18'][i])
...
Multi_H5¶
The wnutils.multi_h5.Multi_H5
class allows you to access and plot data
from multiple webnucleo HDF5 files. First import the namespace by typing:
>>> import wnutils.multi_h5 as m5
Then create a class instance from a list
of HDF5 files.
For this tutorial, type
>>> my_multi_h5 = m5.Multi_H5(['my_output1.h5','my_output2.h5'])
Methods allow you to access or plot data from the files.
Read data from the individual HDF5 instances.¶
To retrieve the individual HDF5 instances from a Multi_H5 instance, type:
>>> h5s = my_multi_h5.get_h5()
To retrieve the original file names, type:
>>> files = my_multi_h5.get_files()
Of course the number of HDF5 instances must equal the number of files. To confirm, type:
>>> print(len(h5s) == len(files))
Use the methods on the individual instances. For example, type:
>>> for i in range(len(h5s)):
... props = h5s[i].get_zone_properties_in_groups(('0','0','0'), ['time'])
... print(files[i],'has',len(props['time']),'groups.')
...