-
Notifications
You must be signed in to change notification settings - Fork 17
/
helita.txt
235 lines (180 loc) · 10.9 KB
/
helita.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
.. _helita-label:
********************
``helita`` interface
********************
The `helita <https://github.com/ITA-Solar/helita>`_ Python package contains several routines to interface with RH 1.5D. Installation instructions are `available in its website <https://helita.readthedocs.io/en/latest/installation.html>`_.
Reading and writing input files
===============================
Writing atmosphere files
------------------------
The ``rh15d`` module in ``helita.sim`` contains a function to write an input atmosphere in RH 1.5D format, assuming the user already has the required data to write at hand. Its function definition is:
.. code-block:: python
def make_xarray_atmos(outfile, T, vz, z, nH=None, x=None, y=None, Bz=None, By=None,
Bx=None, rho=None, ne=None, vx=None, vy=None, vturb=None,
desc=None, snap=None, boundary=None, append=False):
"""
Creates HDF5 input file for RH 1.5D using xarray.
Parameters
----------
outfile : string
Name of destination. If file exists it will be wiped.
T : n-D array
Temperature in K. Its shape will determine the output
dimensions. Shape is generally (nt, nx, ny, nz), but any
dimensions except nz can be omitted. Therefore the array can
be 1D, 2D, or 3D, 4D but ultimately will always be saved as 4D.
vz : n-D array
Line of sight velocity in m/s. Same shape as T.
z : n-D array
Height in m. Can have same shape as T (different height scale
for each column) or be only 1D (same height for all columns).
nH : n-D array, optional
Hydrogen populations in m^-3. Shape is (nt, nhydr, nx, ny, nz),
where nt, nx, ny can be omitted but must be consistent with
the shape of T. nhydr can be 1 (total number of protons) or
more (level populations). If nH is not given, rho must be given!
ne : n-D array, optional
Electron density in m^-3. Same shape as T.
rho : n-D array, optional
Density in kg m^-3. Same shape as T. Only used if nH is not given.
vx : n-D array, optional
x velocity in m/s. Same shape as T. Not in use by RH 1.5D.
vy : n-D array, optional
y velocity in m/s. Same shape as T. Not in use by RH 1.5D.
vturb : n-D array, optional
Turbulent velocity (Microturbulence) in km/s. Not usually needed
for MHD models, and should only be used when a depth dependent
microturbulence is needed (constant microturbulence can be added
in RH).
Bx : n-D array, optional
Magnetic field in x dimension, in Tesla. Same shape as T.
By : n-D array, optional
Magnetic field in y dimension, in Tesla. Same shape as T.
Bz : n-D array, optional
Magnetic field in z dimension, in Tesla. Same shape as T.
x : 1-D array, optional
Grid distances in m. Same shape as first index of T.
y : 1-D array, optional
Grid distances in m. Same shape as second index of T.
x : 1-D array, optional
Grid distances in m. Same shape as first index of T.
snap : array-like, optional
Snapshot number(s).
desc : string, optional
Description of file
boundary : Tuple, optional
Tuple with [bottom, top] boundary conditions. Options are:
0: Zero, 1: Thermalised, 2: Reflective.
append : boolean, optional
If True, will append to existing file (if any).
"""
Note that while in this routine the writing of the hydrogen populations is optional (they can be derived from the mass density, if available), RH 1.5D does not support this yet.
.. note::
The variables passed to ``make_xarray_atmos`` must be consistent with the height scale. The first height index must be the top of the atmosphere (closest to observer), and the height scale must be strictly decreasing.
Reading atmosphere files
------------------------
Once written, the input atmosphere files can be read in Python with ``xarray``, and do not require ``helita``.
For example:
.. code-block:: python
>>> import xarray
>>> atmos = xarray.open_dataset('my_atmos.hdf5')
>>> atmos
<xarray.Dataset>
Dimensions: (depth: 82, nhydr: 6, snapshot_number: 1, x: 5, y: 5)
Coordinates:
* x (x) int64 0 1 2 3 4
* y (y) int64 0 1 2 3 4
z (snapshot_number, depth) float32 ...
* snapshot_number (snapshot_number) int32 0
Dimensions without coordinates: depth, nhydr
Data variables:
temperature (snapshot_number, x, y, depth) float32 ...
velocity_z (snapshot_number, x, y, depth) float32 ...
electron_density (snapshot_number, x, y, depth) float64 ...
hydrogen_populations (snapshot_number, nhydr, x, y, depth) float32 ...
velocity_turbulent (snapshot_number, x, y, depth) float32 ...
Attributes:
comment: Created with make_xarray_atmos on 2018-01-25 15:28:10.4...
boundary_top: 0
boundary_bottom: 1
has_B: 0
description: FAL C model with 82 depth points replicated to 5x5 colu...
nx: 5
ny: 5
nz: 82
nt: 1
The amount of detail loaded by ``xarray`` will depend how the atmosphere was written. Older atmosphere files may not have as much verbose attributes or labeled coordinates (especially if written by plain HDF5 with no attaching of dimension scales), but they are still valid. Older netCDF atmospheres should work fine with ``xarray``.
It is also possible to modify the data with ``xarray``, and saving and updated atmosphere is done via the ``to_netcdf()`` method:
.. code-block:: python
>>> atmos.to_netcdf("newfile.hdf5", format='NETCDF4')
Be sure to use ``format='NETCDF4'`` so that the file is internally HDF5!
Writing wavelength files
------------------------
Another utility function in ``rh15d.py`` is ``make_wave_file``. This creates an RH wavelength file (to be used with the option ``WAVETABLE`` in ``keyword.input``) that contains additional wavelengths to be calculated. The function's usage is documented in its function call:
.. code-block:: python
def make_wave_file(outfile, start=None, end=None, step=None, new_wave=None,
ewave=None, air=True):
"""
Writes RH wave file (in xdr format). All wavelengths should be in nm.
Parameters
----------
start: number
Starting wavelength.
end: number
Ending wavelength (non-inclusive)
step: number
Wavelength separation
outfile: string
Name of file to write.
ewave: 1-D array, optional
Array of existing wavelengths. Program will make discard points
to make sure no step is enforced using these points too.
air: boolean, optional
If true, will at the end convert the wavelengths into vacuum
wavelengths.
"""
You can either supply an array with the wavelengths, or give a range of wavelengths and a fixed spacing, e.g.:
.. code-block:: python
>>> from helita.sim import rh15d
>>> rr = rh15d.Rh15dout()
# this will write wavelenghts from 650 to 650 nm, 0.01 nm spacing:
>>> rh15d.make_wave_file('my.wave', 650, 660, 0.01)
# this will write an existing array "my_waves", if it exists
>>> rh15d.make_wave_file('my.wave', ewave=my_waves)
Reading output files
====================
The main class to read the output is called ``Rh15dout``. It uses ``xarray`` under the hood and populates an object with all the different datasets. It can be initiated in the following way:
.. code-block:: python
>>> from helita.sim import rh15d
>>> rr = rh15d.Rh15dout()
--- Read ./output_aux.hdf5 file.
--- Read ./output_indata.hdf5 file.
--- Read ./output_ray.hdf5 file.
By default, it will look for the three files in the directory specified as main argument (defaults to current directory). Additionally, the method ``read_group(infile)`` can be used to manually load the ``output_aux.hdf5`` or ``output_indata.hdf5`` and the method and ``read_ray(infile)`` can be used to manually load the ``output_ray.hdf5`` file. The variables themselves are not read into memory, but are rather a `memmap object <http://docs.scipy.org/doc/numpy/reference/generated/numpy.memmap.html>`_ (file pointer; only read when needed) that ``xarray`` opens.
After loading the files, the ``Rh15dout`` instance loads each file as an ``xarray`` dataset with the base name of each group (e.g. ``ray``, ``atmos``, ``atom_CA``, ``mpi``).The ``ray`` attribute contains the same dataset as shown in the ``xarray`` example above.
The attributes of each file are still accessible under the attributes of each object, e.g.:
.. code-block:: python
>>> rr.ray.creation_time
'2018-01-10T16:16:42+0100'
>>> rr.atmos.nrays
5
>>> rr.mpi.nprocesses
2048
With ``xarray`` it is easy to quickly inspect and plot different quantities. For example, to plot the intensity at ``(x, y) = (0, 0)``:
>>> rr.ray.intensity[0, 0].plot()
Or the intensity at a fixed wavelength:
>>> rr.ray.intensity.sel(wavelength=279.55, method='nearest').plot()
(This only shows a 2D image if you calculated the intensity from a 3D model, otherwise an histogram or line plot is shown.)
Visualisation and notebooks
===========================
``helita`` includes a visualisation module, ``helita.sim.rh15d_vis``, with widgets that are meant to be used inside the `Jupyter notebook <https://jupyter.org/>`_. To use these, you will need to install not only ``helita`` but also the `Matplotlib Jupyter Extension <https://github.com/matplotlib/jupyter-matplotlib>`_ and the `IPython widgets for Jupyter <https://github.com/jupyter-widgets/ipywidgets>`_. If you have Anaconda, both can be installed with conda::
conda install -c conda-forge ipywidgets ipympl widgetsnbextension
jupyter nbextension enable --py widgetsnbextension
You can also install them with ``pip`` (check their pages for details).
Currently we have the following Jupyter notebooks for visualisation of RH 1.5D output:
* `Basic output <https://github.com/ITA-Solar/rh/blob/master/doc/notebooks/BasicOutput.ipynb>`_
* `Visualisation widgets <https://github.com/ITA-Solar/rh/blob/master/doc/notebooks/Visualisation.ipynb>`_
To use the above notebooks, you need to have run RH 1.5D and have the output files ready!
You can also explore the input atmosphere files with the Jupyter widget ``rh15d_vis.InputAtmosphere`` in ``helita``:
>>> from helita.sim import rh15d_vis
>>> rh15d_vis.InputAtmosphere('my_atmos.hdf5');