MESAHub / mesa

Modules for Experiments in Stellar Astrophysics
http://mesastar.org
GNU Lesser General Public License v2.1
138 stars 38 forks source link

Improve Derivatives for Cubic Interpolation in Composition for Opacity #593

Open evbauer opened 11 months ago

evbauer commented 11 months ago

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 commented 9 months ago

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 commented 9 months ago

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 commented 9 months ago

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?

fxt44 commented 9 months ago

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.

Debraheem commented 9 months ago

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.