Refactor parsing of results files (#176)

* refactor poly fit interface

* refactor fitting methods for gravity and B fields

* refactor trm parsing

* move all file parsing to external functions, return tables instead of arrays

* fix amr parsing for variable columns

* add smoke test for conserved quantity calculation

* add constants module

* respond to code review feedback

* add pkg_resources deprectation warning

* lots of docstrings

* avergage ion mass comment

.gitignore

@@ -3,6 +3,7 @@ site/
 docs/api/*
 docs/config_tables.rst
 docs/generated
+docs/sg_execution_times.rst
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

examples/configure_simulation.py

@@ -27,7 +27,7 @@
 # Mm loop lasting 5000 s heated by a single 200 s nanoflare
 # solved on an adaptive grid.
 # A complete list of configuration parameters can be found in
-# the `configuration-tables`_ page.
+# the :ref:`configuration-tables` page.
 config_dict = {
     'general': {
         'loop_length': 80*u.Mm,

examples/parse_simulation.py

@@ -103,7 +103,7 @@
 # grid. Note that each profile is a `~astropy.units.Quantity` and can
 # be easily converted to any compatible unit.
 print(p.electron_temperature)
-print(p.ion_density)
+print(p.hydrogen_density)
 print(p.velocity)
 
 #################################################################
@@ -125,7 +125,7 @@
 # single time step as a function of field-aligned coordinate.
 # The `~pydrad.parse.Profile` object provides a simple quicklook
 # method for this,
-s[0].peek()
+p.peek()
 
 #################################################################
 # Similarly, we can call this method on a `~pydrad.parse.Strand` to

pydrad/__init__.py

@@ -10,12 +10,6 @@
 
 
 def _init_log():
-    """
-    Initializes the SunPy log.
-    In most circumstances this is called automatically when importing
-    SunPy. This code is based on that provided by Astropy see
-    "licenses/ASTROPY.rst".
-    """
     orig_logger_cls = logging.getLoggerClass()
     logging.setLoggerClass(AstropyLogger)
     try:

pydrad/configure/configure.py

@@ -352,10 +352,13 @@ def poly_fit_magnetic_field(self):
         Sixth-order polynomial fit coefficients for computing flux tube
         expansion
         """
-        fit = self._fit_poly_domains('poly_fit_magnetic_field', 'G')
+        fit_results = self._fit_poly_domains(
+            **self.config['general']['poly_fit_magnetic_field'],
+            target_unit='G',
+        )
         return self.env.get_template('coefficients.cfg').render(
             date=self.date,
-            **fit,
+            **fit_results,
         )
 
     @property
@@ -364,30 +367,37 @@ def poly_fit_gravity(self):
         Sixth-order polynomial fit coefficients for computing gravitational
         acceleration
         """
-        fit = self._fit_poly_domains('poly_fit_gravity', 'cm s-2')
+        fit_results = self._fit_poly_domains(
+            **self.config['general']['poly_fit_gravity'],
+            target_unit='cm s-2'
+        )
         return self.env.get_template('coefficients.cfg').render(
             date=self.date,
-            **fit,
+            **fit_results,
         )
 
-    def _fit_poly_domains(self, name, unit):
+    def _fit_poly_domains(self, x, y, domains, order, target_unit):
         """
         Perform polynomial fit to quantity as a function of field aligned coordinate
         over multiple domains and return fitting coefficients.
         """
-        # TODO: refactor to be independent of dictionary
-        fit = copy.deepcopy(self.config['general'][name])
-        x = (fit['x'] / self.config['general']['loop_length']).decompose().to(u.dimensionless_unscaled).value
-        y = fit['y'].to(unit).value
+        x /= self.config['general']['loop_length']
+        x = x.decompose().to_value(u.dimensionless_unscaled)
+        y = y.to_value(target_unit)
         coefficients = []
         minmax = []
-        for i in range(len(fit['domains'])-1):
-            i_d = np.where(np.logical_and(x>=fit['domains'][i], x<=fit['domains'][i+1]))
-            coefficients.append(np.polyfit(x[i_d], y[i_d], fit['order'])[::-1])
+        for i in range(len(domains)-1):
+            i_d = np.where(np.logical_and(x>=domains[i], x<=domains[i+1]))
+            # NOTE: The order is reversed because HYDRAD expects the opposite order of
+            # what NumPy outputs
+            coefficients.append(np.polyfit(x[i_d], y[i_d], order)[::-1])
             minmax.append([y[i_d].min(), y[i_d].max()])
-        fit['minmax'] = minmax
-        fit['coefficients'] = coefficients
-        return fit
+        return {
+            'domains': domains,
+            'order': order,
+            'minmax': minmax,
+            'coefficients': coefficients,
+        }
 
     @property
     def minimum_cells(self):
@@ -415,7 +425,8 @@ def maximum_cells(self):
         refinement level and :math:`n_{min}` is the minimum allowed number of
         grid cells. Optionally, if the maximum number of cells is specified
         in ``config['grid']['maximum_cells']``, this value will take
-        precedence.
+        precedence. Alternatively, using a value of 30000 will be sufficiently large for
+        any amount of refinement used in HYDRAD.
         """
         if 'maximum_cells' in self.config['grid']:
             return int(self.config['grid']['maximum_cells'])