Convergence issues with mixing with a new eos

We are trying to contribute a new equation of state (eos) to MESA for modeling giant planets, but we are experiencing convergence issues when there is mixing. More details are below. We’re happy to answer questions or provide more details if needed.

### Implementation of the new eos
For the modeling of planets which often have large mass fractions of heavy elements, we implemented a new eos into MESA: Our eos combines the SCvH or CMS hydrogen-helium eos with the quotidian eos for water, using the linear mixing approximation. This is implemented outside of MESA as a python module, which we use to calculate eos tables in the MESA format.

To use our tables in MESA, we use the other_eos_frac and other_eos_component hooks. The code that loads and interpolates the tables is in large parts simply copied from the private eos routines (since our tables are in the same format). The main modification we made is to port the way X tables are interpolated to the Z tables (instead of only linear) and made the choice of interpolation method an inlist option.

For homogeneously mixed planets, or if we suppress mixing by setting mix_factor = 0, the module runs just fine and is consistent with MESA’s results. For example, we re-implemented MESA’s SCvH tables and got the same result with and without our custom eos module. The main two differences we encountered during testing were that 1.) create_initial_model takes a lot longer and 2.) evolving a planet takes a little bit longer.

### The error
The error occurs when mixing a gradient in Z (or a steep gradient in Y) for a gaseous planet. Over many inlist setups, the simulation breaks from one of two errors:

After mixing some material, the solver’s residual of dlnE_dt (or sometimes equ_o16; we use o16 as our heavy element) increases drastically, accompanied by a large max correction in lnd (sometimes also max corr o16). Typical values for the residuals are about 1E+4 and for the max corrections about 5. This is followed by many (typically adjacent) cells throwing a hydro_mtx: logT too large (or small) error. As the error indicates the cells jump to completely unreasonable values. For example, a cell with logT = 4.52 would change to logT = -2.2 or logT = 8.5. These values are then rejected by the solver, throwing a retry: logT > hydro_mtx_max_allowed_logT error, and reducing the timestep. This pattern continues until the time step is too small and the simulation aborts.

Sometimes, the solver doesn’t terminate because of the hydro error, but because the residuals and max correction become extremely large and the solver simply rejects the solution. A typical example here would be:

> 27 3 coeff 1.0000
> avg resid 0.108E+01
> max resid dlnE_dt 1889 0.93888E+04
> mix type 00000
> avg corr 0.149E+03
> max corr lnd 1886 0.62059E+06
> mix type 00000

This also triggers a reduction in the timestep, which ultimately leads to the termination of the simulation. This error pattern occurs for both, our eos as well as MESA’s SCvH and CMS eos. We attached two files of the solver output for a model with a Z gradient and a Y gradient respectively.

### Inlist Setup
As a default case, we use a gas giant of 1 Jupiter mass and an initial specific entropy of s = 10 kb/bary, but we also tested other entropies and found the same behavior. The inlist options we use are motivated by other test suite examples. We attached a typical inlist to evolve the model, too. Without going through every single inlist option we tested, we can say that this convergence problem persists over a wide range of options (e.g., different convergence tolerances or energy equation options). We tried smoothing the eos using different techniques, but this also didn’t help much. While we are able to run a Y gradient model with the right set of inlist options using our smoothed version of the CMS eos, we were only able to get one Z model to work: A model were we reduce the mixing factor by 1e-7 and don’t use convective premixing (or predictive mixing). However, all these “makeshift” solutions don’t address the underlying issue that the code struggles to mix even small amounts of Z material.

Any help is much appreciated, thanks!

[output_Y_gradient_run.txt](https://github.com/MESAHub/mesa/files/13219282/output_Y_gradient_run.txt)
[output_Z_gradient_run.txt](https://github.com/MESAHub/mesa/files/13219284/output_Z_gradient_run.txt)
[inlist_evolve.txt](https://github.com/MESAHub/mesa/files/13219286/inlist_evolve.txt)

Start by reproducing the failure with inlist_evolve.txt and comparing output_Y_gradient_run.txt with output_Z_gradient_run.txt. Review the other_eos_frac and other_eos_component integration and the reported solver residuals, using the attached runs as diagnostics. Done means identifying the convergence cause and demonstrating stable evolution with composition gradients and mixing.

Written by the indexing model from the issue text.