Setting the background heating to 0 with adaptive step causes crash

[jwreep]

When I set the background heating to zero while using adaptive step, or some small number (1e-50), the code crashes with the following error message:

[xerxes:~/ebtelPlusPlus] jreep% bin/ebtel++.run libc++abi.dylib: terminating with uncaught exception of type std::runtime_error: Adaptive solver exceeded maximum number of allowed failures. Abort

The code runs fine with adaptive step turned off or with a larger background heating. Example configuration:

<?xml version="1.0" ?>
<root>
  <total_time>5000.0</total_time>
  <tau>1</tau>
  <tau_max>10</tau_max>
  <loop_length>40.0e+8</loop_length>
  <saturation_limit>1.0</saturation_limit>
  <force_single_fluid>False</force_single_fluid>
  <use_c1_loss_correction>True</use_c1_loss_correction>
  <use_c1_grav_correction>True</use_c1_grav_correction>
  <use_power_law_radiative_losses>True</use_power_law_radiative_losses>
  <use_flux_limiting>False</use_flux_limiting>
  <calculate_dem>False</calculate_dem>
  <save_terms>False</save_terms>
  <use_adaptive_solver>True</use_adaptive_solver>
  <output_filename>ebtel++_results_file.txt</output_filename>
  <adaptive_solver_error>1e-6</adaptive_solver_error>
  <adaptive_solver_safety>0.5</adaptive_solver_safety>
  <c1_cond0>2.0</c1_cond0>
  <c1_rad0>0.6</c1_rad0>
  <helium_to_hydrogen_ratio>0.075</helium_to_hydrogen_ratio>
  <dem>
    <use_new_method>True</use_new_method>
    <temperature bins="451" log_min="4" log_max="8.5"/>
  </dem>
  <heating>
    <background>0.0</background>
    <partition>1.0</partition>
    <events>
      <event magnitude="0.1" rise_start="0.0" rise_end="50.0" decay_start="50.0" decay_end="100.0"/>
      <event magnitude="0.15" rise_start="200.0" rise_end="250.0" decay_start="250.0" decay_end="300.0"/>
      <event magnitude="0.05" rise_start="1000.0" rise_end="1050.0" decay_start="1050.0" decay_end="1100.0"/>
    </events>
  </heating>
</root>

[wtbarnes]

Ok I can reproduce the failure with the adaptive step. When set "use_adaptive_solver" to False, the code does not return an error, but the answers are all NaN. The reason for this is that NaNs are marked as failures in the adaptive solver, but not in the fixed-step solver. So when the adaptive solver sees a NaN, it tries to fix it (until it exceeds the max number of tries) whereas the fixed-step solver just lets it pass on through. This probably needs fixing.

The bigger issue here is the lack of background heating. The reason this leads to a NaN is that the heating at t=0 is used to configure the intial temperature and density (see the CalculateInitialConditions() function). It requires that the heating at t=0 be finite in order to calculate an equilibrium temperature and density.

If you have a heating function which is 0 at t=0 (as you do with these triangular profiles), then the code relies on the background heating to set the initial conditions. If that is 0, you get nonsense (or rather NaN!) for the initial temperature and density values. So this behavior is actually expected.

If you'd like to avoid using background heating, I'd suggest adding a heating event which is nonzero at t=0 (e.g. a square profile) that persists until the next heating event. This will allow you to start with a desired temperature and density (based on the heating rate you start with), but will not impose any background heating later on.

As an example, the heating configuration

<events>
      <event magnitude="1e-5" rise_start="0.0" rise_end="0.0" decay_start="1000.0" decay_end="1000.0"/>
      <event magnitude="0.05" rise_start="1000.0" rise_end="1050.0" decay_start="1050.0" decay_end="1100.0"/>
</events>

with a total simulation time of 9000 s produces the following electron temperature and density profiles,

Notice how at late times, the temperature and density fall below their initial values because there is no background heating to support them at equilibrium.

Let me know if you have any questions about this. I actually don't think the failure for 0 background heating is a bug per se (other than the fixed-step solver not raising an error for NaN).

[jwreep]

Thanks! It looks like this solves the issue for my purposes.

In general, do you think it would be worth adding a line like this?

double heat = heater->Get_Heating(0.0); if( heat < SMALL_THRESHOLD) heat = DEFAULT_NONZERO_VALUE;

It would have the advantage of preventing crashes due to user error, but it has the drawback of defining some default profile. My guess is that most of the time, there would be a background heating, so this might not be a problem anyway.

[wtbarnes]

I sort of prefer the idea of explicitly (and loudly) failing rather than silently adding some background heating. I think that could be potentially confusing. At least with a crash, you know there is a problem and can fix it.

That being said, there probably should be some sort of warning (or error) if the heating at t=0 is not properly configured. Or at the very least, a note in the docs.

[wtbarnes]

Thanks for reporting this @jwreep! It's brought up two other important issues.