Improve Derivatives for Cubic Interpolation in Composition for Opacity
evbauer opened this issue · 5 comments
Right now, our default is to do linear interpolation in composition (X,Z) for opacity tables. Even though cubic interpolation returns more realistic opacities, it returns much worse derivatives. Turning on cubic interpolation by default results in a number of test suite failures: https://testhub.mesastar.org/test_kap_cubic_splines/commits
It seems a big part of the issue is that we interpolate the derivatives from the tables, rather than taking derivatives of the interpolants. This may require a fairly significant amount of effort to refactor and overhaul how the code approaches composition interpolation.
mesa/kap/private/kap_eval_fixed.f90
Lines 518 to 534 in c020320
To get a sense of how bad the derivatives are, here are some dfridr plots from @Debraheem at X = 0.625, Z = 0 when using cubic interpolation.
Incorporating autodiff for the composition interpolation helps substantially (see f373537).
The dfridr plot at X=0.625, Z=0 now looks like this:
If we add in some metallicity that is between tabulation points, the plot for X=0.625, Z=5e-4 looks somewhat worse:
So we still want to investigate where these residual bad stripes are coming from, but even this is substantially better than what things looked like previously.
Zooming in on the area with the stripes shows that these bad derivatives occur right along lines of constant logT and logR corresponding to grid points in the underlying opacity tabulations.
I think the dfridr badness is actually mostly an artifact of the fact that the underlying derivatives have discontinuities across grid points, which perhaps shouldn't be shocking. Here's a zoom in on the structure of the logT derivative when doing cubic composition interpolation (with autodiff):
And for comparison here's the derivative when doing linear composition interpolation:
This might improve if we allowed the composition interpolation to use something other than the "piecewise monotonic" cubic interpolation routines. Could we relax that so that we maybe get less discontinuities in the derivatives of the interpolants?
Or should we declare victory and say that the quality of these derivatives is good enough?
@fxt44 do you have any thoughts?
as long as this effort is in experimental mode, maybe do the experiment of trying a different cubic. even if it ends up doing little in this opacity case, it might be informative for other cubic interpolation use cases.
Using the monotonicity preserving interpolation routines e.g. interp_mp -> interp_m3a, interp_m3b, interp_m3q in place of the peace-wise monotonic spline used to generate Evan's plots (interp_pm) from $MESA_DIR/interp_1d/public/interp_lib.f90
where a, q, b stand for the (average, quartic, and superbee limiters) MESA uses citing:
Huynh, H.T., "Accurate Monotone Cubic Interpolation", SIAM J Numer. Anal. (30) 1993, 57-100.
Suresh, A, and H.T. Huynh, "Accurate Monotonicity-Preserving Schemes with Runge-Kutta Time Stepping", JCP (136) 1997, 83-99.
I remade the same two (dfridr and dlogkapdlogT) plots for each.
m3a
m3b
m3q
These monotinicity presrvering routines don't seem to do much better in terms reducing discontinuities in the derivative of the interpolant.