use adaptive stepping for exact Riemann integration

[zingale]

PR summary

There are 2 main changes here:

1. the riemann_support.H header is now split into 2, one for shocks and one for rarefactions
2. the rarefaction integration no longer uses a fixed stepsize, but instead now does an adaptive RK4 integration

This changes results to roundoff / tol error, as expected.

PR motivation

PR checklist

• test suite needs to be run on this PR
• this PR will change answers in the test suite to more than roundoff level
• all newly-added functions have docstrings as per the coding conventions
• the CHANGES file has been updated, if appropriate
• if appropriate, this change is described in the docs

[zingale]

Here are comparisons of density for the 4 test problems from our 2015 paper:

[zingale]

Here are some of my notes for this:
exact_riemann.pdf

[aisclark91]

Everything looks great, there is only one issue:

In riemman_star_state.H, after ustar_l and ustar_r are computed of each iteration, in line 141 we obtain :
ustar = 0.5_rt * (ustar_l + ustar_r);
This expression is not right, the actual expression is:
ustar = (Z_l * ustar_l + Z_r * ustar_r)/(Z_l + Z_r)
In accordance with Van Leer [22] eq A3. I double checked it to be sure.

[zingale]

I don't think that the expression I am using is wrong. If the iteration were successful, then u_star,l = u_star,r, so this average is just dealing with the tolerance of the solver.

[zingale]

it is actually just a different average, the Z's there are the slope at the intersection, which should be the same if we've converged the root finding.