Lots of improvements (#143)

* refactor time distance plotting to be more general

* unit fixes; strand em; flexible bounds on em

* function for making clean copy of hydrad

* cleanup bkgrnd heating config; add overwrite option

* fix b-field/gravity configuration

* fix pass by reference issue

* fix heating template tests

* clean up CI config now that HYDRAD is public

* adjust time array tests

* ci debugging

* fix glob sorting

.github/workflows/test.yml

@@ -25,6 +25,17 @@ jobs:
         python -m pip install --upgrade pip
         python -m pip install -r requirements/requirements-dev.txt
         python setup.py install
+    - name: Clone HYDRAD
+      uses: actions/checkout@v2
+      with:
+        repository: rice-solar-physics/HYDRAD
+        path: HYDRAD
     - name: Test
       run: |
-        pytest --cov=pydrad --cov-report=xml
+        pytest --hydrad-dir=$GITHUB_WORKSPACE/HYDRAD --cov=pydrad --cov-report=xml
+    - name: Report coverage
+      uses: codecov/codecov-action@v1
+      with:
+        fail_ci_if_error: true
+        verbose: true
+        file: ./coverage.xml

README.md

@@ -1,7 +1,6 @@
 # pydrad
 
 ![pydrad CI status](https://github.com/rice-solar-physics/pydrad/workflows/Test/badge.svg)
-![pydrad CI status](https://github.com/rice-solar-physics/pydrad/workflows/Test%20with%20HYDRAD/badge.svg)
 [![Documentation Status](https://readthedocs.org/projects/pydrad/badge/?version=latest)](https://pydrad.readthedocs.io/en/latest/?badge=latest)
 [![codecov](https://codecov.io/gh/rice-solar-physics/pydrad/branch/master/graph/badge.svg)](https://codecov.io/gh/rice-solar-physics/pydrad)
 

pydrad/configure/configure.py

@@ -14,7 +14,7 @@
 import asdf
 
 from . import filters
-from .util import run_shell_command, on_windows
+from .util import run_shell_command, on_windows, get_equilibrium_heating_rate
 
 __all__ = ['Configure']
 
@@ -69,8 +69,12 @@ def save_config(self, filename):
         """
         asdf.AsdfFile(self.config).write_to(filename)
 
-    def setup_simulation(self, output_path, base_path,
-                         run_initial_conditions=True, **kwargs):
+    def setup_simulation(self,
+                         output_path,
+                         base_path,
+                         run_initial_conditions=True,
+                         overwrite=False,
+                         **kwargs):
         """
         Setup a HYDRAD simulation with desired outputs from a clean copy
 
@@ -79,6 +83,8 @@ def setup_simulation(self, output_path, base_path,
         base_path (`str`): Path to existing HYDRAD
         run_initial_conditions (`bool`): If True, compile and run the initial
         conditions code
+        overwrite : `bool`
+            If True and the `output_path` directory already exists, overwrite it
         """
         with tempfile.TemporaryDirectory() as tmpdir:
             # NOTE: this is all done in a temp directory and then copied over
@@ -89,6 +95,9 @@ def setup_simulation(self, output_path, base_path,
                 self.setup_initial_conditions(tmpdir, execute=execute)
             self.setup_hydrad(tmpdir)
             self.save_config(os.path.join(tmpdir, 'pydrad_config.asdf'))
+            if overwrite:
+                if os.path.exists(output_path):
+                    shutil.rmtree(output_path)
             shutil.copytree(tmpdir, output_path)
 
     def setup_initial_conditions(self, root_dir, execute=True):
@@ -146,20 +155,7 @@ def setup_initial_conditions(self, root_dir, execute=True):
                 root_dir,
             )
             if self.config['heating']['background'].get('use_initial_conditions', False):
-                self.equilibrium_heating_rate = self.get_equilibrium_heating_rate(root_dir)
-
-    def get_equilibrium_heating_rate(self, root_dir):
-        """
-        Read equilibrium heating rate from initial conditions results
-
-        # Parameters
-        root_dir (`str`): Path to HYDRAD directory
-        """
-        filename = os.path.join(root_dir,
-                                'Initial_Conditions/profiles/initial.amr.sol')
-        with open(filename, 'r') as f:
-            equilibrium_heating_rate = float(f.readline()) * u.erg / u.s / (u.cm**3)
-        return equilibrium_heating_rate
+                self.equilibrium_heating_rate = get_equilibrium_heating_rate(root_dir)
 
     def setup_hydrad(self, root_dir):
         """
@@ -292,26 +288,11 @@ def heating_cfg(self):
         values, you must run the initial conditions and set the
         `equilibrium_heating_rate` attribute first.
         """
-        if self.config['heating'].get('background', False):
-            bg = self.config['heating']['background']
-            if bg.get('use_initial_conditions', False):
-                background = {
-                    'rate': self.equilibrium_heating_rate,
-                    'location': self.config['initial_conditions']['heating_location'],
-                    'scale_height': self.config['initial_conditions']['heating_scale_height'],
-                }
-            elif all((k in bg for k in ('rate', 'location', 'scale_height'))):
-                background = self.config['heating']['background']
-            else:
-                raise ValueError(
-                    'Set use_initial_conditions to True or set parameters '
-                    'explicitly in order to use background heating.')
-        else:
-            background = {
-                'rate': 0*u.erg/(u.cm**3 * u.s),
-                'location': 0*u.cm,
-                'scale_height': 0*u.cm
-            }
+        background = copy.deepcopy(self.config['heating']['background'])
+        if background.get('use_initial_conditions', False):
+            background['rate'] = self.equilibrium_heating_rate
+            background['location'] = self.config['initial_conditions']['heating_location']
+            background['scale_height'] = self.config['initial_conditions']['heating_scale_height']
         return self.env.get_template('heating.cfg').render(
             date=self.date,
             background=background,
@@ -390,10 +371,10 @@ 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')
         return self.env.get_template('coefficients.cfg').render(
             date=self.date,
-            fit=self.config['general']['poly_fit_magnetic_field'],
-            y_unit='G',
+            **fit,
         )
 
     @property
@@ -402,12 +383,31 @@ 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')
         return self.env.get_template('coefficients.cfg').render(
             date=self.date,
-            fit=self.config['general']['poly_fit_gravity'],
-            y_unit='cm s-2',
+            **fit,
         )
 
+    def _fit_poly_domains(self, name, 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
+        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])
+            minmax.append([y[i_d].min(), y[i_d].max()])
+        fit['minmax'] = minmax
+        fit['coefficients'] = coefficients
+        return fit
+
     @property
     def minimum_cells(self):
         """

pydrad/configure/templates/coefficients.cfg

@@ -1,10 +1,9 @@
-{{ fit.order }}
-{{ fit.domains | length - 1}}
-{% for d in fit.domains %}{{ d }}    {% endfor %}
-{{ fit.x | length }}
-{%- set y_val = fit.y | units_filter(y_unit) -%}
-{% for x in fit.x %}
-{{ x }}    {{ y_val[loop.index0] }}
-{%- endfor %}
+{{ order }}
+{{ domains | length - 1}}
+{% for d in domains %}{{ d }}    {% endfor %}
+{% for c in coefficients -%}
+{% for _c in c %}{{ _c }}    {% endfor %}
+{{ minmax[loop.index0][0] }}    {{ minmax[loop.index0][1] }}
+{% endfor %}
 
 Configuration file generated by pydrad on {{ date }}

pydrad/configure/tests/test_templates.py

@@ -46,7 +46,11 @@ def test_collisions_header(configuration):
 
 def test_heating_config(configuration):
     # No events, no background
-    configuration.config['heating']['background'] = {}
+    configuration.config['heating']['background'] = {
+        'location': 0*u.cm,
+        'scale_height': 0*u.cm,
+        'rate': 0.0*u.erg/(u.cm**3)/u.s 
+    }
     heating_config = f"""0.00000000e+00 0.00000000e+00 0.0
 
 0
@@ -72,14 +76,12 @@ def test_heating_config(configuration):
     }
     with pytest.raises(AttributeError):
         configuration.heating_cfg
-    # Cannot determine background model
+    # 1 event, no background
     configuration.config['heating']['background'] = {
-        'foo': 'bar'
+        'location': 0*u.cm,
+        'scale_height': 0*u.cm,
+        'rate': 0.0*u.erg/(u.cm**3)/u.s 
     }
-    with pytest.raises(ValueError):
-        configuration.heating_cfg
-    # 1 event, no background
-    configuration.config['heating']['background'] = False
     configuration.config['heating']['events'] = [
         {'time_start': 0.*u.s,
          'rise_duration': 100.*u.s,

pydrad/configure/util.py

@@ -1,11 +1,24 @@
 """
 Utilities for HYDRAD configuration
 """
+import os
 import subprocess
 import platform
+import tempfile
+import shutil
+from distutils.dir_util import copy_tree
+
+import astropy.units as u
 
 from pydrad import log
 
+__all__ = ['MissingParameter',
+           'HYDRADError',
+           'run_shell_command',
+           'on_windows',
+           'hydrad_stripped_down',
+           'get_equilibrium_heating_rate']
+
 
 class MissingParameter(Exception):
     """
@@ -48,3 +61,53 @@ def on_windows():
     Determine whether the user's operating system is Windows
     """
     return platform.system().lower() == 'windows'
+
+
+def hydrad_stripped_down(output_path, base_path=None, from_github=False):
+    """
+    Create a clean copy of HYDRAD with only the files necessary to run the code. May be
+    useful when making many copies.
+
+    Parameters
+    ----------
+    output_path : pathlike
+        Path to the new stripped down version
+    base_path : pathlike, optional
+        Path to the original copy. This will not be modified.
+    from_github : `bool`, optional
+        If True, grab the latest copy of HYDRAD from GitHub. In this case,
+        `base_path` is ignored. Note that this requires the GitPython package.
+    """
+    # NOTE: this is all done in a temp directory and then copied over
+    # so that if something fails, all the files are cleaned up
+    with tempfile.TemporaryDirectory() as tmpdir:
+        if from_github:
+            import git
+            git.Repo.clone_from('https://github.com/rice-solar-physics/HYDRAD',
+                                tmpdir)
+        elif base_path:
+            copy_tree(base_path, tmpdir)
+        else:
+            raise ValueError('Specify local path to HYDRAD or clone from GitHub')
+        rm_dirs =['Forward_Model', 'HYDRAD_GUI', 'Visualisation', '.git']
+        rm_files = ['HYDRAD_GUI.jar', 'LICENSE', 'README.md', '.gitignore']
+        for d in rm_dirs:
+            shutil.rmtree(os.path.join(tmpdir, d))
+        for f in rm_files:
+            os.remove(os.path.join(tmpdir, f))
+        shutil.copytree(tmpdir, output_path)
+
+
+def get_equilibrium_heating_rate(root_dir):
+    """
+    Read equilibrium heating rate from initial conditions results
+
+    Parameters
+    ----------
+    root_dir : `str` or pathlike
+        Path to HYDRAD directory
+    """
+    filename = os.path.join(root_dir, 'Initial_Conditions/profiles/initial.amr.sol')
+    with open(filename, 'r') as f:
+        equilibrium_heating_rate = float(f.readline()) * u.Unit('erg cm-3 s-1')
+    return equilibrium_heating_rate