MESAHub/mesa

Port Adrian's Python Implementation of Fraser 2023 Thermohaline into MESA

Closed this issue · 3 comments

afraser3 commented

(Recall that in BGS13, HG19, and this model, once you get w_f it's a simple formula to get Nu_mu)

For the not-so-great model in section 5.1 of that paper (the dotted curve in fig 5), w_f is calculated via the w_f function here. This relies upon the following functions in kolmogorov_EVP.py:

  • Lmat sets up the matrix we want the eigenvalues of
  • gamfromL just passes L to np.linalg.eig to get the largest eigenvalue, "gamma" or "gam" in the code here, confusingly called "sigma" in the paper
  • gamma_over_k and gammax_kscan runs the above but for a range of kz's to maximize growth rate over kz
  • gammax_minus_lambda evaluates Eq 28 in my paper (but in the form LHS - RHS = 0) so that, in w_f, scipy.optimize.root_scalar can use it to find the w_f that satisfies Eq 28

Currently working on updating this so that the L matrix is the "with T and C" one that Sam and Pascale worked out (the better one). Otherwise the machinery will all be more or less the same.

afraser3 commented

I've implemented the "better" model in my repo (the "with T, C" one in FRG23, described in sec 5.2; as opposed to the worse model, "no T, C", described in sec 5.1). Enough changes that in several instances I elected to implement it as separate functions, rather than passing flags to the existing functions. But the structure remains the same. Below I try to list out the changes we need to note for implementation in MESA. It looks tedious, but very little of the structure/logic is changing, you just end up needing to redo a lot of things because you need slightly different parameters. If you want to take it for a spin yourself, this script recreates fig 8 from the paper; sadly it takes tens of minutes at least to run.

Note that there are several places where I'm using the kwarg Sam=True, this is because I'm copy+pasting Sam Reifenstein's matrix while I try to figure out what's wrong with my own matrix. For now, we want the Sam=True stuff. Also, wherever you see withTC=True, that's a flag to use the functions corresponding to the better model.

I'll walk through this top-down, so starting with w_f and ending with Lmat.
Starting from parasite_model.py:

  • Instead of w_f, the new model uses w_f_withTC. Most of the input arguments are the same, but I got rid of the delta argument (which should always be set to 0 anyway), and the model parameter CH (see C_H in HG19 Eq 31) has been renamed to C2 (see C_2 in FRG23 Eq 28) -- note these parameters are related according to C_2 = 1/C_H. Also, there's now the argument Sam which, as mentioned above, should be set to True until I fix my own matrix.
  • Other changes in this file are only relevant to my own purposes, not for the MESA implementation.

Next, kolmogorov_EVP.py (AKA kh_instability.f90):

  • Previously, w_f used scipy's root-finding combined with the function gammax_minus_lambda. Now, w_f_withTC does the same thing but calls gammax_minus_lambda_withTC. This is essentially the same, except the same CH vs C2 and Sam arguments mentioned above are also passed along here. The old version calculated the Reynolds number, magnetic Reynolds number, and "M2", but those are no longer necessary and have been dropped; also, the old model was calculating sigma (the KH growth rate) in a different non-dimensionalization than the new model is, and that's why the old code multiplied sigma by w*lhat but this version does not. Note that this function returns sigma*C2 - lambda (where sigma is the KH growth rate, lambda is the thermohaline growth rate) while the old function returned sigma - lambda*CH.
  • The above function calls gammax_kscan_withTC, which replaces the gammax_kscan called by the old function. The same Sam argument continues to be passed along here. The only other notable change is that there's some new logic on whether or not to throw an error/warning in the event that gammax == gammas[0], which I'm happy to explain in greater detail.
  • The above function calls gamma_over_k_withTC, which replaces the gamma_over_k called by the old function. Only difference is that the Sam parameter is passed along, and we call the underlying sigma_from_fingering_params with the withTC=True argument (I could have similarly made a new sigma_from_fingering_params_withTC function)
  • sigma_from_fingering_params now includes new code that is executed if withTC is True. The first thing to note is it does a unit conversion with kz = k_star * lhat -- each wavenumber k in the ks array that you first pass to w_f is assumed to be measured in units of \hat{l}_f (the wavenumber you get from gaml2max, the function Rich has been improving); this line converts it to the units described by Eq 10 of FRG23, because this model was set up with slightly different units. Next, for Sam=True (which we care about), it sets up a few of the parameters that Sam's matrix needs as inputs, then constructs the matrix and calculates its eigenvalues (the old version shipped the matrix off to a different function to have it calculate the eigenvalues).
  • Finally, the previous Lmat is replaced by what I'm currently calling Sams_Lmat here

Note that Rich has implemented his own version of this Matrix, and has done so in a way that has all-real coefficients. We'll probably want to use his matrix, which will likely mean minor changes to the arguments that sigma_from_fingering_params passes along to Lmat.

  • 👍2
evbauer commented

MESA can now produce this using turb/plotter:
iTerm2 tpR3WZ turb_plotter

This looks pretty much identical to the plot produced by paper_plot_fig8.py that @afraser3 linked above. As far as porting the python code goes, I think that's now done.

Next we will need to hook this into MESA to be available for stellar evolution calculations, and optimize it to be efficient enough to call on the fly.

  • 🎉3