Improve Derivatives for Cubic Interpolation in Composition for Opacity

[evbauer]

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

	call interp1(logKs, logK, ierr)
	if (ierr /= 0) then
	call mesa_error(__FILE__,__LINE__,'failed in interp1 for logK')
	return
	end if
	call interp1(dlogKs_dlogRho, dlogK_dlogRho, ierr)
	if (ierr /= 0) then
	call mesa_error(__FILE__,__LINE__,'failed in interp1 for dlogK_dlogRho')
	return
	end if
	call interp1(dlogKs_dlogT, dlogK_dlogT, ierr)
	if (ierr /= 0) then
	call mesa_error(__FILE__,__LINE__,'failed in interp1 for dlogK_dlogT')
	return
	end if

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.

[evbauer]

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.

[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.

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:

[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?

[fxt44]

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]

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.