[Mesa-users] Implementation of the Canuto-Mazzitelli Model

Wed Mar 18 22:04:53 UTC 2026

Hello Vishal,

Without digging too deep into your progress, perhaps you could investigate how the scale height is limited to the width of the convection zone in the overshooting routines. This might give you some insight into how to repurpose this logic into the convection model you're developing. https://docs.mesastar.org/en/25.12.1/reference/controls.html#limit-overshoot-hp-using-size-of-convection-zone.

If you're looking to add inlist controls for testing, you can always call the optional star controls like "s% x_ctrl(:)" inside the MESA, and modify it from your inlist. But if you want to add controls to the defaults, it is a little bit more involved, adding in the control in a few places across the code. In general you need to add the control into $MESA_DIR/star_data as well as in some locations inside $MESA_DIR/star/private, and $MESA_DIR/star/defaults, but I recommend picking a working mesa control and grepping around to see how it's initialized.

A good way to test, is to pipe your convection model results into the unit tests inside $MESA_DIR/turb/test and compare directly to the results from mlt. You can compare thermodynamic gradients and convective velocities etc. there.

If you desire more detailed code specific assistance, perhaps this discussion can continue in a MESA issue thread or via an external pull request, where you can link your forked mesa-version so those who are interested in helping can see more directly your code, and implementation. https://github.com/MESAHub/mesa/issues

-EbF
________________________________
From: Vishal Vetrivel <vishalsvetrivel at gmail.com>
Sent: Wednesday, March 18, 2026 9:50 AM
To: Farag, Ebraheem <ebraheem.farag at yale.edu>
Cc: MESA Users List <mesa-users at lists.mesastar.org>
Subject: Re: [Mesa-users] Implementation of the Canuto-Mazzitelli Model

Hello all,
I've been working on this for a bit now and I believe I've made some progress.

First off I changed the root finding method to a bisection method. This was initially motivated by the realization that the roots that were being computed were not exactly "good" roots (the function value at the root was not close to zero). I thought this was due to the safe_root routine but I realized that the functions, that I needed to find roots of for CM and CGM, were practically vertical near the roots for large values of delta (The function can be found in Eq (36) of [1]). Due to this, near the root the derivatives were incredibly large and even evaluating the function accurately was difficult. To solve this I have divided the whole expression by delta to keep the roots in the same location but ensure the function is easy to evaluate at all values (I don't do this when abs(delta) < 1 however so the function doesn't get large). This combined with a reasonable upper bound for the bisection made the roots much better. I am still using the analytic derivatives since I find the roots with bisection now and I didn't see the point in running a Newton-Raphson iteration just to generate the derivatives when I had analytic expressions for them. Granted I am not entirely sure if the way I've set the derivatives ensures all the required derivatives are set.

I had a few doubts regarding some other things that I wanted to implement. One of the things I wanted to do in the beginning was to set the mixing length to the distance from the top of the convective zone as suggested by Canuto in his first paper (Section 2.8 of [1]). In my first mail I had mentioned that there was a way to limit the mixing length to this value but I haven't been able to find it anywhere again and I think it would make more sense to implement it properly. To this end I took a look at the code and I'm not sure how I'm supposed to calculate the distance to the top of the convective zone. I assume the k variable passed to the MLT routines allows me to find the current zone and I did see that there were variables storing the locations of the convective zones but I wanted to know if there was a recommended way of calculating distance to the top of the current convective zone.

There is also a parameter (Kolmogorov Constant K0) in the Canuto-Goldman-Mazzitelli Model that should be adjustable by the user. In Section 5 of [2] they discuss multiple values that could be used for this and I think exposing this as an inlist setting makes sense. Is there any particular way that this should be done?

I also wanted to know if there's any particular way I should be testing the code I'm writing. This is the first time I'm working on code larger than a few files and I know I'm supposed to be writing unit tests but I'm unsure of what is required for MESA. I've tested that the model works at least with solar calibration as that seems to give frequencies that match the observed values better than MLT. If there's any other tests I should be running please do let me know.

[1] https://articles.adsabs.harvard.edu/pdf/1991ApJ...370..295C
[2] https://articles.adsabs.harvard.edu/pdf/1996ApJ...473..550C

On Wed, Nov 19, 2025 at 9:59 PM Farag, Ebraheem <ebraheem.farag at yale.edu<mailto:ebraheem.farag at yale.edu>> wrote:
Hello Vishal, (please keep mesa-users cc'd)

I would suggest turning turbulent pressure off during the pre-ms relaxation, and then reenabling it when the model begins evolving.

If you're keen enough to try, (now that you have some experience modifying mesa to fix your mlt gradT bug), you're welcome to try adding said function into the auto_diff variable types in $MESA_DIR/star/auto_diff, to see if this alleviates your issue with adopting autodiff. I can say for certain that the numerical stability of your routines will likely be affected by your choice to include or not include auto_diff, or to return analytic derivatives. The mlt_vc_ad and gradT_ad returned by mlt to the star newton solver typically carry with them derivatives that aid in convergence.

-EbF
________________________________
From: Vishal Vetrivel <vishalsvetrivel at gmail.com<mailto:vishalsvetrivel at gmail.com>>
Sent: Wednesday, November 19, 2025 5:43 AM
To: Farag, Ebraheem <ebraheem.farag at yale.edu<mailto:ebraheem.farag at yale.edu>>
Subject: Re: [Mesa-users] Implementation of the Canuto-Mazzitelli Model

Thank you for your replies.
I read the paper on the auto_diff variables and I get how it works now. I am not sure if I can even implement it the way it is written now since I am using the safe_root for the Gamma value (which seems to accept only a real value). Unless there is a similar utility for auto_diff values I think I might just stick to what I have right now.
As for the turbulent pressure I tried it out but I wasn't able to get a model to get beyond the pre main sequence model creating for C values from Table 1 in the 1991 CM paper (~0.2). The solver seems to give up in the pre main sequence model generation. To check I tried lowering the value of C and it seems to work up to around 0.1. I am unsure exactly what the issue is but the final errors state that it stopped due to dt < min_timestep_limit. I can provide more information from the output but I'm not sure what would be relevant so I have included the last few lines of the MESA output I got when I ran a max_Pturb_factor above 0.1 (even 0.12 doesn't seem to work).

stopping because of problems dt < min_timestep_limit
 failed in relax num steps
 failed in do_relax_num_steps
 star_create_pre_ms_model ierr          -1
 do_load1_star ierr          -1
 do_before_evolve_loop ierr          -1
 do_before_evolve_loop ierr          -1

On Mon, Nov 17, 2025 at 9:36 PM Farag, Ebraheem <ebraheem.farag at yale.edu<mailto:ebraheem.farag at yale.edu>> wrote:
Hello Vishal,

Is the question, where to implement the turbulent pressure in other structure equations such as the energy and momentum equation?

Following https://adsabs.harvard.edu/pdf/1991ApJ...370..295C, I see (eq40) Pt = 2*C*(1/2*rho*v_turb^2). Where C is given in Table 1.

I believe in your implementation  v_turb -> conv_vel -> (in mesa star) mlt_vc.

In that case, the turbulent pressure will take similar form to the mlt_turbulent_pressure (see https://docs.mesastar.org/en/24.08.1/reference/controls.html#mlt-pturb-factor), the documentation here is actually slightly outdated, as we mass weight on to faces instead of simple averaging (see hydro_energy):

where Pt_mlt = 1/3 * rho * v_turb^2 = 1/3 * rho * mlt_vc^2.

To implement turbulent pressure in the momentum equation, all you need to do is change 'mlt_Pturb_factor' to something that matches the CM model.

mlt_Pturb_factor = 3*C.

-EbF
________________________________
From: Mesa-users <mesa-users-bounces at lists.mesastar.org<mailto:mesa-users-bounces at lists.mesastar.org>> on behalf of Philip Mocz via Mesa-users <mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>>
Sent: Monday, November 17, 2025 9:08 AM
To: Vishal Vetrivel <vishalsvetrivel at gmail.com<mailto:vishalsvetrivel at gmail.com>>
Cc: mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org> <mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>>
Subject: Re: [Mesa-users] Implementation of the Canuto-Mazzitelli Model

Hi Vishal,

I'm excited you are trying to implement this into MESA!

Section 2 of the MESA VI paper shows how to use the autodiff variables to automatically get derivatives of a function with respect to it's input variables. https://arxiv.org/abs/2208.03651  . It's no problem to avoid using autodiff and use analytic derivatives instead if you have those expressions.

I'll have to dig more to find where to connect/update turbulent pressure. Does anyone know off the top of their head?

Best,
Philip

On Sun, Nov 16, 2025 at 5:42 AM Vishal Vetrivel via Mesa-users <mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>> wrote:
Hello all,
I have been attempting to implement the Canuto-Mazzitelli and Canuto-Goldman-Mazzitelli model into MESA for the past few months. I managed to find an old MESA user list thread from Dr. Warrick Ball (https://sourceforge.net/p/mesa/mailman/message/33552129/<https://urldefense.com/v3/__https://sourceforge.net/p/mesa/mailman/message/33552129/__;!!DSb-azq1wVFtOg!S9dz0uQsHWzxET8X62meLlV0zU_BkQ1c-CzGSBs9bFVEgpiQ5_VJ91WGHYR0ekhNw7xjYbzDdCXwQW3qRKb55XXBe2eGUWET$>) detailing a similar implementation. However, since this was a while ago a lot of the code has since changed. I initially hoped to implement it in the same way by using a `run_stars_extra.f90` file. However, a bug in the MESA prevented me from doing so as the routine get_gradT does not call the custom mlt hook. Instead I have been implementing this into a local copy of the MESA source code (I am very grateful that you can run multiple versions of MESA so easily).

In this process I have managed to implement all the modifications made in the original thread. However, I am facing some issues figuring out where to implement the turbulent pressure. I have managed to replace the convective velocity with what I believe should be the equivalent expressions for CM and CGM. I have not yet implemented setting the mixing length to the distance to the convective surface yet (I believe there is already an option to limit the mixing length to the depth so this might not be necessary). A major issue I had was I am not sure how the new auto_diff variables are supposed to work so I have manually assigned the derivatives using the old reference implementation. This is a pretty bad way to do it so I'd be grateful if anyone could help me with how these variables are meant to be used.

I have attached my modified mlt.f90 file from the turb module which I believe should be most if not all the modifications I have made to source code. I am using version 24.08.1 of MESA. Thank you in advance for your help.

--
Regards
Vishal Sai Vetrivel
5th Year Integrated Masters
School of Physical Sciences
UM DAE CEBS, Mumbai

_______________________________________________
mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>
https://urldefense.com/v3/__https://lists.mesastar.org/mailman/listinfo/mesa-users__;!!DSb-azq1wVFtOg!S9dz0uQsHWzxET8X62meLlV0zU_BkQ1c-CzGSBs9bFVEgpiQ5_VJ91WGHYR0ekhNw7xjYbzDdCXwQW3qRKb55XXBe8VerUaX$

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20260318/d2e7d4a5/attachment.htm>