mesa icon indicating copy to clipboard operation
mesa copied to clipboard

Improve Derivatives for Cubic Interpolation in Composition for Opacity

Open evbauer opened this issue 1 year ago • 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. https://github.com/MESAHub/mesa/blob/c0203203faadbbaeb218d5c9bf5a37d9aff29886/kap/private/kap_eval_fixed.f90#L518-L534

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.

image image image image

evbauer avatar Oct 19 '23 16:10 evbauer

Incorporating autodiff for the composition interpolation helps substantially (see f373537d75d79950f165132ce67a3f64638375c8).

The dfridr plot at X=0.625, Z=0 now looks like this: iTerm2 gskssA kap_plotter

If we add in some metallicity that is between tabulation points, the plot for X=0.625, Z=5e-4 looks somewhat worse: iTerm2 yv3rvE cubic

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.

evbauer avatar Jan 03 '24 23:01 evbauer

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. iTerm2 MPzYyv kap_plotter

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): iTerm2 ewoC2l dlogT_cubic

And for comparison here's the derivative when doing linear composition interpolation: iTerm2 e8vWIS dlogT_linear

evbauer avatar Jan 04 '24 18:01 evbauer

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?

evbauer avatar Jan 04 '24 18:01 evbauer

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.

fxt44 avatar Jan 04 '24 19:01 fxt44

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 dfridr_m3a

derivative_m3a dfridr_m3a

m3b dfridr_m3b derivative_m3b

m3q

dfridr_m3q

derivative_m3q

These monotinicity presrvering routines don't seem to do much better in terms reducing discontinuities in the derivative of the interpolant.

Debraheem avatar Jan 04 '24 22:01 Debraheem