Transfer classes for eigendata, NEI simulations, and simulation results to this repo #20
Conversation
Co-authored-by: Chengcai Shen <[email protected]>
These tests will not work since they are in files that were directly copied over from the NEI-modeling/NEI repository.
|
Hello @namurphy! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:
Comment last updated at 2020-06-03 21:55:45 UTC |
The class is adapted from the class EigenData2 which is in the eigenvaluetable.py file that will be deleted soon. The tests are adapted from the test_eigenvaluetable.py class that will also be deleted soon.
These were copied over from the NEI-modeling/NEI repository.
Co-authored-by: Chengcai Shen <[email protected]> Co-authored-by: Marcus Dupont <[email protected]>
|
I still need to parse the errors more, but it looks like the HDF5 file is not being read in on the Azure tests. Tomorrow... |
|
Looking into the test failures now... |
namurphy
changed the title
Avoid using Optional[u.s] as this seems to cause problems.
Codecov Report
@@ Coverage Diff @@
## master #20 +/- ##
============================================
- Coverage 100.00% 80.99% -19.01%
============================================
Files 2 4 +2
Lines 4 721 +717
============================================
+ Hits 4 584 +580
- Misses 0 137 +137
Continue to review full report at Codecov.
|
This is to try to pin down some doc build errors for which error messages are not being useful. I largely removed references to external packages (e.g., changing `~astropy.units.Quantity` to ``astropy.units.Quantity``). These should be added back in later.
In order to get links to different packages in the docs.
|
|
||
| Parameters | ||
| ---------- | ||
| """ |
There was a problem hiding this comment.
- add parameter descriptions
- more explicit variable names, or add inline comments for clarification
- add some background/cites to the docstrings that explains the calculation
There was a problem hiding this comment.
I agree with all of these. I'm going to slightly address these now, and raise an issue or two so that the refactoring will be done in a follow-up pull request. I also think this would be more readable if it were broken up into smaller functions, which may reduce the need for inline clarification comments.
plasmapy_nei/eigen/eigenclass.py
Outdated
| """ | ||
| Returns the inverses of the eigenvectors for the ionization and | ||
| recombination rates for the temperature specified in the class. | ||
| """ |
There was a problem hiding this comment.
add documentation for parameters
plasmapy_nei/eigen/eigenclass.py
Outdated
| """ | ||
| Return the equilibrium charge state distribution for the | ||
| temperature specified in the class. | ||
| """ |
There was a problem hiding this comment.
add documentation for parameters
Co-authored-by: Erik Everson <[email protected]>
| # TODO: Allow this to keep track of velocity and position too, and | ||
| # eventually to have density and temperature be able to be functions of | ||
| # position. (and more complicated expressions for density and | ||
| # temperature too) | ||
|
|
||
| # TODO: Expand Simulation docstring | ||
|
|
||
|
|
||
| # TODO: Include the methods in the original Visualize class which is a | ||
| # subclass of NEI in the NEI-modeling/NEI repo. These were deleted | ||
| # temporarily to make it possible to get the NEI class itself | ||
| # adapted into this package. | ||
|
|
||
|
|
||
| # TODO: In this file and test_nei.py, there are a few places with | ||
| # initial.ionic_fractions.keys(), where initial is an instance | ||
| # of IonizationStates. This workaround exists because I forgot | ||
| # to put in an `elements` attribute in IonizationStates, and | ||
| # should be corrected. |
There was a problem hiding this comment.
These really should be issues. You've already stated that - just a reminder. :)
plasmapy_nei/nei/nei.py
Outdated
| import plasmapy as pl | ||
| from scipy import interpolate, optimize | ||
| from plasmapy_nei.eigen import EigenData | ||
| from plasmapy.atomic import IonizationStates |
There was a problem hiding this comment.
The attached PR wraps this in a try/except. Why do you use this ancient version, anyway? :p
| self, | ||
| initial: IonizationStates, | ||
| n_init: u.Quantity, | ||
| T_e_init: u.Quantity, |
There was a problem hiding this comment.
I'm missing a quantity validator somewhere around here...
| The new ionization fractions for this time step. The keys | ||
| of this `dict` are the atomic symbols of the elements being | ||
| tracked, and with the corresponding value being an | ||
| ``numpy.ndarray`` representing the ionic fractions. Each |
There was a problem hiding this comment.
Looks like a nested list would probably do. "2D iterable?" ndarray is fine.
| @property | ||
| def last_step(self) -> int: | ||
| """The time index of the last step.""" | ||
| return self._last_step | ||
|
|
||
| @property | ||
| def nstates(self) -> Dict[str, int]: | ||
| """ | ||
| Return the dictionary containing atomic symbols as keys and the | ||
| number of ionic species for the corresponding element as the | ||
| value. | ||
| """ | ||
| return self._nstates |
There was a problem hiding this comment.
For this kind of field, if there's no setter anywhere... is this just Sphinx-driven development? :)
There was a problem hiding this comment.
I tend to do this when either (1) I want sphinx documentation to be more readable and explicit, (2) I want the code to be more readable, and/or (3) I want to control overwriting of the attribute. In this case, nstates would be be protect but its items() would not be. So, I'm guessing # 1 or # 2 is the rationale here.
| T_e: astropy.units.Quantity or callable | ||
| The electron temperature, which may be a constant, an array of | ||
| temperatures corresponding to the times in `time_input`, or a | ||
| function that yields the temperature as a function of time. |
There was a problem hiding this comment.
Okay, that's pretty cool!
There was a problem hiding this comment.
This was my favorite part of this whole NEI class. With all that said, that reminds me that it will be necessary to implement a different formulation...
- Create issue on allowing heat to be added to the parcel of plasma over the course of the simulation
- Create issue on enabling cooling due to adiabatic expansion
- Create issue on enabling radiative cooling based on the actual ionization states (this one is going to be a tough one)
| time step. Setting ``verbose`` to `True` is useful for testing. | ||
| Defaults to `False`. | ||
|
|
||
| abundances: dict |
There was a problem hiding this comment.
Undocumented parameter again.
| # print(f"time_sta = ", time) | ||
| # print(f"INI: ", f0) | ||
| # print(f"Sum(f0) = ", np.sum(f0)) |
| adapt_dt: `bool` | ||
| If `True`, change the time step based on the characteristic | ||
| ionization and recombination time scales and change in | ||
| temperature. Not yet implemented. |
There was a problem hiding this comment.
- Create an issue on implementing the adaptive time step
| The ionization and recombination rates are from Chianti version | ||
| 8.7. These rates include radiative and dielectronic recombination. | ||
| Photoionization is not included. | ||
| """ |
There was a problem hiding this comment.
- Create issue: Put in a reference to Shen et al. (2015) Astronomy & Computing, 12, 1, "A Lagrangian scheme for time-dependent ionization simulations of astrophysical plasmas"...in particular eqs 1 & 2
There was a problem hiding this comment.
Wait, how is it Lagrangian...? This looked rather, er, space-free...
There was a problem hiding this comment.
In practice, this technique would be used following a parcel of plasma, for example as it moves away from the Sun or passes through a supernova remnant shock wave. Another example would be taking the output of a numerical simulation and tracing the temperature and density of a parcel of plasma as a function of time.
In any case, it'll be necessary to clarify this in the forthcoming narrative documentation and be careful about language. I'll probably avoid the terms "Lagrangian" and "Eulerian" because it always takes me 17 seconds to remember which is which.
|
|
||
| The keys of this dictionary are atomic symbols. The values are | ||
| 2D arrays where the first index refers to the time step and the | ||
| second index refers to the integer charge. |
There was a problem hiding this comment.
Add an example of how to use this attribute
| if time is None: | ||
| self._time_max = ( | ||
| self.time_input[-1] if self.time_input is not None else np.inf * u.s | ||
| ) | ||
| elif isinstance(time, u.Quantity): | ||
| if not time.isscalar: | ||
| raise ValueError("time_max must be a scalar") | ||
| try: | ||
| time = time.to(u.s) | ||
| except u.UnitConversionError: | ||
| raise u.UnitsError("time_max must have units of seconds") from None | ||
| if ( | ||
| hasattr(self, "_time_start") | ||
| and self._time_start is not None | ||
| and self._time_start >= time | ||
| ): | ||
| raise ValueError("time_max must be greater than time_start") |
There was a problem hiding this comment.
The two times have a codependence on each other, but I still wonder if it might not make sense to extract this code out?
| def hydrogen_number_density(self, time: u.Quantity) -> u.Quantity: | ||
| try: | ||
| time = time.to(u.s) | ||
| except (AttributeError, u.UnitsError): | ||
| raise NEIError("Invalid time in hydrogen_density") |
There was a problem hiding this comment.
We could probably use another validator here. And a docstring, too.
There was a problem hiding this comment.
- Create issue on applying validators in
NEIclass
| # Is there a way to use the inspect package or something similar | ||
| # to only return self.results if it is in an expression where |
There was a problem hiding this comment.
What did you have in mind?
There was a problem hiding this comment.
I wouldn't have simulate() return results at all, since that encourages users to separate the results from the originating simulation. I would encourage the use of...
>>> sim.simulate()
>>> sim.results.plot_something()
Or, if the result was a xarray.Dataset then one could do...
>>> sim.simulate()
>>> sim.results.T_e.plot()
| # we can use the old timestep as a reasonable guess. If we are | ||
| # simulation, then we can either use the inputted timestep or |
There was a problem hiding this comment.
If we are simulation, I am the :=.
| """Advance the simulation by one time step.""" | ||
| # TODO: Expand docstring and include equations! | ||
|
|
||
| # TODO: Fully implement units into this. |
There was a problem hiding this comment.
If you want to take that performance hit, sure; but I'd try to do repetitive calculations without them for now, and keep units to interfaces. I think advancements in numpy dtypes may help, but I'm not entirely sure.
There was a problem hiding this comment.
Yeah...performance is going to be important here.
There was a problem hiding this comment.
- Create issue to optimize the time advance that includes this comment
| Save the `~plasmapy_nei.nei.NEI` instance to an HDF5 file. Not | ||
| implemented. | ||
| """ | ||
| raise NotImplementedError |
There was a problem hiding this comment.
Would have been simple with an xarray! ;)
Co-authored-by: Dominik Stańczak <[email protected]>
Co-authored-by: Erik Everson <[email protected]>
| for j in range(self.nstates): | ||
| self._eigenvalues[temperature_index, j] = eigenvalues[j] | ||
| self._equilibrium_states[temperature_index, j] = eqi[j] | ||
| for i in range(self.nstates): | ||
| self._eigenvectors[temperature_index, i, j] = eigenvectors[i, j] | ||
| self._eigenvector_inverses[temperature_index, i, j] = v_inverse[ | ||
| i, j | ||
| ] |
There was a problem hiding this comment.
| for j in range(self.nstates): | |
| self._eigenvalues[temperature_index, j] = eigenvalues[j] | |
| self._equilibrium_states[temperature_index, j] = eqi[j] | |
| for i in range(self.nstates): | |
| self._eigenvectors[temperature_index, i, j] = eigenvectors[i, j] | |
| self._eigenvector_inverses[temperature_index, i, j] = v_inverse[ | |
| i, j | |
| ] | |
| self._eigenvalues[temperature_index, ...] = eigenvalues | |
| self._equilibrium_states[temperature_index, ...] = eqi | |
| self._eigenvectors[temperature_index, ...] = eigenvectors | |
| self._eigenvector_inverses[temperature_index, ...] = v_inverse |
More for readability than speed.
| def temperature(self): | ||
| return self._temperature | ||
|
|
||
| @property |
There was a problem hiding this comment.
| @property | |
| @temperature.setter |
There was a problem hiding this comment.
Review of nei.py.
As @StanczakDominik has mentioned several times, I too am curious what this would look like if all the data was stored in an internal xarray.Dataset. It could also provide quick access to plotting the data. If implemented, then the implementation should probably be a separate PR.
A comment for nei.py and eigenclass.py. There is not much unit wrangling or discussion of expected units in the docstrings. This does not need to be something handled in this PR, but you should open an issue for it. This will probably include the use of plasmapy.utils.decorators.validate_quantities. The nei.py module is better about unit wrangling, but there still is not much about units in the docstrings.
| import numpy as np | ||
| from typing import Union, Optional, List, Dict, Callable | ||
| import astropy.units as u | ||
| from scipy import interpolate, optimize | ||
| from plasmapy_nei.eigen import EigenData |
| for elem in elements: | ||
| self._ionic_fractions[elem][index, :] = new_ionfracs[elem][:] | ||
|
|
||
| # Calculate elemental and ionic number densities | ||
| n_elem = {elem: new_n * self.abundances[elem] for elem in elements} | ||
| number_densities = { | ||
| elem: n_elem[elem] * new_ionfracs[elem] for elem in elements | ||
| } | ||
|
|
||
| # Calculate the electron number density | ||
| n_e = 0.0 * u.cm ** -3 | ||
| for elem in elements: | ||
| integer_charges = np.linspace( | ||
| 0, self.nstates[elem] - 1, self.nstates[elem] | ||
| ) | ||
| n_e += np.sum(number_densities[elem] * integer_charges) | ||
|
|
||
| # Assign densities | ||
| self._n_e[index] = n_e | ||
| for elem in elements: | ||
| self._n_elem[elem][index] = n_elem[elem] | ||
| self._number_densities[elem][index, :] = number_densities[elem] |
There was a problem hiding this comment.
I count this for-loop, for elem in elements, is done at least 5 time here. Should be easy enough to pack into one for-loop.
| @property | ||
| def last_step(self) -> int: | ||
| """The time index of the last step.""" | ||
| return self._last_step | ||
|
|
||
| @property | ||
| def nstates(self) -> Dict[str, int]: | ||
| """ | ||
| Return the dictionary containing atomic symbols as keys and the | ||
| number of ionic species for the corresponding element as the | ||
| value. | ||
| """ | ||
| return self._nstates |
There was a problem hiding this comment.
I tend to do this when either (1) I want sphinx documentation to be more readable and explicit, (2) I want the code to be more readable, and/or (3) I want to control overwriting of the attribute. In this case, nstates would be be protect but its items() would not be. So, I'm guessing # 1 or # 2 is the rationale here.
| The number density multiplicative factor. The number density of | ||
| each element will be ``n`` times the abundance given in | ||
| ``abundances``. For example, if ``abundance['H'] = 1``, then this | ||
| will correspond to the number density of hydrogen (including | ||
| neutral hydrogen and protons). This factor may be a constant, | ||
| an array of number densities over time, or a function that | ||
| yields a number density as a function of time. |
There was a problem hiding this comment.
Personally I don't really follow this. I'd suggest adding a non-unity example and/or an example with two ion species.
| After having inputted all of the necessary information, we can run | ||
| the simulation. | ||
|
|
||
| >>> results = sim.simulate() |
There was a problem hiding this comment.
Why not have the simulation generate a results attribute? For example,
>>> sim.results
No simulate results generated
>>> sim.simulate()
>>> sim.results
We got data!!!Then you wouldn't have a detach of the simulation parameters and the results, unless the user explicitly does results = sum.results.
| T_e(self.time_start).to(u.K) | ||
| T_e(self.time_max).to(u.K) |
| def initial(self, initial_states: IonizationStates): | ||
| if isinstance(initial_states, IonizationStates): | ||
| self._initial = initial_states | ||
| self._elements = ( |
There was a problem hiding this comment.
| self._elements = ( | |
| self._elements = list( |
| # Is there a way to use the inspect package or something similar | ||
| # to only return self.results if it is in an expression where |
There was a problem hiding this comment.
I wouldn't have simulate() return results at all, since that encourages users to separate the results from the originating simulation. I would encourage the use of...
>>> sim.simulate()
>>> sim.results.plot_something()
Or, if the result was a xarray.Dataset then one could do...
>>> sim.simulate()
>>> sim.results.T_e.plot()
| tol=1e-6, | ||
| ) | ||
|
|
||
| if not np.isclose(self.time_max / u.s, self.results.time[-1] / u.s): |
There was a problem hiding this comment.
Since you're policing units, why not just do...
| if not np.isclose(self.time_max / u.s, self.results.time[-1] / u.s): | |
| if not np.isclose(self.time_max.value, self.results.time[-1].vlaue): |
Co-authored-by: Erik Everson <[email protected]>
|
I did not yet address all of the issues made during code review, but rather created issues #25, #26, #27, #28, #29, #30, #31, and #32 in order to record remaining tasks. Each of these issues warrants its own pull request (and sometimes even multiple pull requests), so I'm going to merge the current pull request now so that it does not become more intractable. There are still valuable suggestions from code review above, so the above comments will be really helpful for when the different parts of the package are refactored and revised:
My sincere thanks to everyone who contributed to the SunNEI and NEI packages, and to everyone who provided suggestions and feedback on this PR! |
This PR adapts the
EigenData2class from the NEI-modeling/NEI repository for use in this package atEigenData. This class provides access to the ionization and recombination rates.The PR also includes the
NEIclass (and companion classes) as adapted from the same repository.This should be a merge commit without squashing to preserve the original files from the other repository here.