Throw error for insufficient initial heating

examples/util.py

@@ -46,31 +46,35 @@ def run_ebtel(config, ebtel_dir):
         )
         if cmd.stderr:
             raise EbtelPlusPlusError(f"{cmd.stderr.decode('utf-8')}")
-        data = np.loadtxt(results_filename)
-
-        results = {
-            'time': data[:, 0]*u.s,
-            'electron_temperature': data[:, 1]*u.K,
-            'ion_temperature': data[:, 2]*u.K,
-            'density': data[:, 3]*u.cm**(-3),
-            'electron_pressure': data[:, 4]*u.dyne/(u.cm**2),
-            'ion_pressure': data[:, 5]*u.dyne/(u.cm**2),
-            'velocity': data[:, 6]*u.cm/u.s,
-            'heat': data[:, 7]*u.erg/(u.cm**3*u.s),
-        }
-
-        results_dem = {}
-        if config['calculate_dem']:
-            results_dem['dem_tr'] = np.loadtxt(
-                config['output_filename'] + '.dem_tr')
-            results_dem['dem_corona'] = np.loadtxt(
-                config['output_filename'] + '.dem_corona')
-            # The first row of both is the temperature bins
-            results_dem['dem_temperature'] = results_dem['dem_tr'][0, :]*u.K
-            results_dem['dem_tr'] = results_dem['dem_tr'][1:, :]*u.Unit('cm-5 K-1')
-            results_dem['dem_corona'] = results_dem['dem_corona'][1:, :]*u.Unit('cm-5 K-1')
-
-    return {**results, **results_dem}
+
+        try:
+            data = np.loadtxt(results_filename)
+        except:
+            raise EbtelPlusPlusError("Failed to output a data file.\n")
+        else:
+            results = {
+                'time': data[:, 0]*u.s,
+                'electron_temperature': data[:, 1]*u.K,
+                'ion_temperature': data[:, 2]*u.K,
+                'density': data[:, 3]*u.cm**(-3),
+                'electron_pressure': data[:, 4]*u.dyne/(u.cm**2),
+                'ion_pressure': data[:, 5]*u.dyne/(u.cm**2),
+                'velocity': data[:, 6]*u.cm/u.s,
+                'heat': data[:, 7]*u.erg/(u.cm**3*u.s),
+            }
+
+            results_dem = {}
+            if config['calculate_dem']:
+                results_dem['dem_tr'] = np.loadtxt(
+                    config['output_filename'] + '.dem_tr')
+                results_dem['dem_corona'] = np.loadtxt(
+                    config['output_filename'] + '.dem_corona')
+                # The first row of both is the temperature bins
+                results_dem['dem_temperature'] = results_dem['dem_tr'][0, :]*u.K
+                results_dem['dem_tr'] = results_dem['dem_tr'][1:, :]*u.Unit('cm-5 K-1')
+                results_dem['dem_corona'] = results_dem['dem_corona'][1:, :]*u.Unit('cm-5 K-1')
+
+            return {**results, **results_dem}
 
 
 def read_xml(input_filename,):

source/loop.cpp

@@ -137,6 +137,20 @@ state_type Loop::CalculateInitialConditions(void)
   double c1 = 2.0;
   double c2 = CalculateC2();
   double heat = heater->Get_Heating(0.0);
+
+  /* Check that the initial heating value is high enough to calculate the equilibrium conditions.  At low
+   * initial heating rates, the code would output NaNs because it cannot find an equilibrium.  This minimum 
+   * threshold is found from the RTV scaling laws for coronal loops of an equilibrium temperature of 0.01 MK,
+   * (combining the laws in the form found in Reale 2014), H = (3 / L^2) * (T / 1.4)^(7/2), 
+   * which corresponds to a heating rate of approximately 9.24e-8 erg/s/cm^3 for a 10 Mm loop, falling 
+   * quadratically with length.  This is slightly higher than where the code actually fails, but puts the
+   * equilibrium conditions into a questionable temperature regime, regardless.  */
+  double minimum_heat = 9.24e10 / (parameters.loop_length * parameters.loop_length);
+  if( heat < minimum_heat )
+  {
+      std::string error_message = "Insufficient initial heating to calculate the equilibrium conditions.\nIncrease the heating at time 0.";
+      throw std::runtime_error(error_message);
+  }
   state_type state;
 
   if( parameters.use_lookup_table_losses )

tests/test_solver.py

@@ -7,7 +7,7 @@
 import numpy as np
 
 from .helpers import run_ebtelplusplus
-
+from util import EbtelPlusPlusError
 
 @pytest.fixture(scope='module')
 def base_config():
@@ -58,7 +58,7 @@ def static_results(base_config):
     config = base_config.copy()
     config['use_adaptive_solver'] = False
     return run_ebtelplusplus(config)
-
+    
 @pytest.mark.parametrize(['name', 'atol'], [
     ('electron_temperature', 1e4),
     ('ion_temperature', 1e4),
@@ -75,3 +75,10 @@ def test_quantities_equal_adaptive_static(adaptive_results, static_results, name
     # an error > 10%; due to static case not rising fast enough
     assert np.allclose(adapt_interp[5:], static_results[name][5:].to_value(adaptive_results[name].unit),
                        rtol=1e-2, atol=atol)
+
+@pytest.mark.parametrize('value', [-1e-5, 0, 1e-15])
+def test_minimum_heating(base_config, value):
+    config = base_config.copy()
+    config['background'] = value
+    with pytest.raises(EbtelPlusPlusError):
+        run_ebtelplusplus(config)