Understand differences between IDL and C++ codes

[wtbarnes]

Now that #39 has been merged, there are some tests for comparing against the IDL implementation by using the single-fluid option in ebtel++. However, this test is currently marked as an expected failure because there small differences that I do not quite understand.

These differences occur mostly within 10 s of the start of the heating event and are very short-lived though there are also some minor differences in the late cooling stage as well. Except for the very early heating-stage differences, which can be up to ~20%, all differences are a order ~1% or less.

These differences become larger when using a non-zero He abundance.

• [ ] Understand differences in single-event case
• [ ] Add multi-event test against IDL code
• [ ] Test against non-zero He abundance

[wtbarnes]

Since merging #43, differences seem to be limited to only the first 10 s, with very good agreement (<0.1%) everywhere else.

[wtbarnes]

After the modifications for area expansion in the IDL code, I've found that the differences are once again significant, even with no expansion and only a coronal loop length. As an example, here is a plot of ebtel++ (dashed), HYDRAD (solid), and EBTEL IDL (dotted). Note the differences between the IDL code and ebtel++,

If instead I set $c_1=2$ at all times, I get the following,

As such, the differences here must be down to how $c_1$ is calculated. I need to go back and double check both the IDL code and the C++ code to make sure there is not a typo in how this is calculated.

[wtbarnes]

Another update: After the fix to $c_1$ implemented in #81, the comparison is now:

where the last panel represents the relative difference between the two. It is still unclear to me why the two diverge so much at both a)the initial rise phase and b) the late cooling and draining phase.

It's possible that this is down to the solver that's being used. In the IDL code, it's just a simple first order Euler method. In the C++ case, I'm using a Cash-Karp method which in principle should be more accurate.

[wtbarnes]

Note that for doing exact comparisons, should probably be using the same solver. In Boost, you can use a 1st order Euler method,

typedef boost::numeric::odeint::euler< state_type > euler_stepper_type;

// Constant timestep integration
num_steps = boost::numeric::odeint::integrate_const(euler_stepper_type(),
                                                    loop->CalculateDerivatives,
                                                    state,
                                                    loop->parameters.tau,
                                                    loop->parameters.total_time,
                                                    loop->parameters.tau,
                                                    obs->Observe);

But even for small timesteps, this still doesn't fully solve the problem,

In fact, things look worse at the end.

[wtbarnes]

Even still using the Cash-Karp method, I find better agreement as I lower the timestep (here at 0.125 s),