[Mesa-users] Resolving TP-AGB

Wed Apr 1 12:57:25 UTC 2026

Hello Theo,

Following  your initial description "the first thermal pulses go smoothly, but the late evolution struggles and the timestep becomes small and goes into a loop of increase/decrease between -3 < log dt <  -2.5."

Have you plotted quantities from your model visually or inspected its evolution in the pgstar, perhaps look at a Kippenhahn and take note of exactly when the model begins to encounter numerical issues? I noticed you are using a large (not small) amount of convective overshooting on all the convective shells including burning shells in your model. Perhaps some shells (H-He shells or H and nonburning envelope) are merging/touching during the late thermal pulses inducing the sudden timestep dive? Thermohaline driven mixing might contribute to this behavior as well. (Perhaps you could confirm if this is the case or not).

Also, since you are using mlt++. You might be transporting an unphysical amount of energy into the surface layers, hence you might need to raise Pextra_factor = 1d0 ----> 1.5d0 to stabilize the surface of your stellar model. High Pextra_factor becomes increasingly important for stabilizing the WD atmosphere as well, although increasing beyond values of 1.5d0 can result in unphysical radii for more extended stars such as those evolving through the TP-AGB, see section 6.1 in MESA VI<https://ui.adsabs.harvard.edu/abs/2023ApJS..265...15J/abstract>.

Lastly, as opposed to directly stripping the model, boosting the mass loss rate along the AGB will generally result is less thermal pulses and a smoother evolution into a the white dwarf phase. One could try increasing the Blocker_scaling_factor to 0.5d0 or even 1d0 to advance through this phase more quickly, although this will change the mapping between your initial model and final WD mass and composition.

If you are going to use hydrodynamics and drag, consider extending min_dq_for_drag down to 0.8d0. The purpose is to damp out an kappa-epsilon driven pulsations from He/H ionization and recombination, which MESA is likely to resolve when the timestep dives to resolve the He shell driven flash. H-recombination is one of the main reasons for the supposed recombination driven runaway described in https://ui.adsabs.harvard.edu/abs/1994A%26A...290..807W/abstract. In reality, this runaway behavior is likely related to the existence of Mira long period radially pulsating variable stars along the AGB. Whether the last few layers are ejected in one episode is up for debate, however I would lean on multi-dimensional simulations like CO5-Bold models for insight in this territory, See, Freytag et al. 2017<https://ui.adsabs.harvard.edu/abs/2017A%26A...600A.137F/abstract> or
Ahmad et al. 2023<https://ui.adsabs.harvard.edu/abs/2023A%26A...669A..49A/abstract>. Falk might also have some thoughts here.

Lastly, Natalie's work<https://ui.adsabs.harvard.edu/abs/2024MNRAS.531.4033R/abstract> attempts to smooth over the dynamical behavior encountered during this phase, by evolving the model hydrostatically and injecting/removing heat energy into the right layers of the stars to prevent ionization/recombination driven instabilities entirely. So you if you employ her methods, I would suggest evolving with v_flag = .false.

-EbF
________________________________
From: Mesa-users <mesa-users-bounces at lists.mesastar.org> on behalf of Théo-Courty via Mesa-users <mesa-users at lists.mesastar.org>
Sent: Wednesday, April 1, 2026 7:49 AM
To: Frank Timmes <fxtimmes at gmail.com>
Cc: mesa-users <mesa-users at lists.mesastar.org>
Subject: Re: [Mesa-users] Resolving TP-AGB

Hi,

Thank you everyone and sorry for the late reply.

The MIST paper used an earlier version of mesa so I don’t use exactly the same parameters. I have same convection treatment (Henyey + alpha mlt = 1.8) but different overshooting scheme. I tried to play with different convection/overshooting parameters during TP-AGB, but it didn’t seem to solve the issues that much.

I think I will try what Ana suggested, do some pulses until the core has grown, and then remove the enveloppe. Unfortunately, I want to evolve multiple systems, so week/days per simulation is not something I can afford for my project.

Best,

Théo

> Le 27 mars 2026 à 22:47, Frank Timmes <fxtimmes at gmail.com> a écrit :
>
> hi théo,
>
> mesa can evolve through all the thermal pulses to a wd (for example, as done in https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fscixplorer.org%2Fabs%2F2022ApJ...935...21C%2Fabstract&data=05%7C02%7CEbraheem.Farag%40yale.edu%7Caa0643dd866d4a3868c108de8fe4e8bb%7Cdd8cbebb21394df8b4114e3e87abeb5c%7C0%7C0%7C639106410514148638%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=0kpVHLx4ZF7JNr5cal7psofVTfXCjem7dUliQXharIw%3D&reserved=0<https://scixplorer.org/abs/2022ApJ...935...21C/abstract> ), but the cost is to use very small time steps that can result in run times of order a week or longer. that your timestep is flopping around between -3 < log dt <  -2.5 is a large clue that the timesteps are too large. the best way to decrease the timestep (https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdocs.mesastar.org%2Fen%2F25.12.1%2Freference%2Fcontrols.html%23timestep-controls&data=05%7C02%7CEbraheem.Farag%40yale.edu%7Caa0643dd866d4a3868c108de8fe4e8bb%7Cdd8cbebb21394df8b4114e3e87abeb5c%7C0%7C0%7C639106410514177349%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=eAtzWDLF2yy%2BeLYzOLlKewLn3%2FYcEYpyHP01Lgq8k3A%3D&reserved=0<https://docs.mesastar.org/en/25.12.1/reference/controls.html#timestep-controls>) is to limit changes in the relevant physical quantities (not numerical quantities). for example, experiment with significantly lowering the default values for changes in the density delta_lgRho_limit, temperature delta_lgT_limit, luminosity dL_divₗ_limit, helium luminosity delta_lgL_He_limit, and the major elements of the composition dX_limit. you've found the right combination when the convergence of a single timestep takes ~3 newton iterations and the timestep is varying smoothly and slowly. with patience for longer run times, and a little determination, it can be done!
>
> fxt
>
>
>
>
>> On Mar 26, 2026, at 6:20 AM, Théo-Courty via Mesa-users <mesa-users at lists.mesastar.org> wrote:
>>
>>
>> Hi Ana,
>>
>> Thanks a lot for your answer! I’ll definitely look at that paper. I’m primarily interested in the white dwarf stage, however when using the make_co_wd test suite (with the remove_env routine) I cannot reproduce the initial-final mass relation of https://ui.adsabs.harvard.edu/abs/2016ApJ...823..102C/abstract and https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Facademic.oup.com%2Fmnras%2Farticle%2F527%2F2%2F3602%2F7420520&data=05%7C02%7CEbraheem.Farag%40yale.edu%7Caa0643dd866d4a3868c108de8fe4e8bb%7Cdd8cbebb21394df8b4114e3e87abeb5c%7C0%7C0%7C639106410514212661%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=kUnmfxMrU2YaqQ0RuINP3HkflSsYfOZ9c24vRZlNO7w%3D&reserved=0<https://academic.oup.com/mnras/article/527/2/3602/7420520>. The final mass is significantly underestimated, especially around 2-3Msun (0.5-0.6 Msun vs expected 0.6-0.7Msun). At higher masses it seems to fit better however. I went through the TP-AGB route to use a more physically accurate mass loss scheme to correct that issue. Anyway, thanks again for your reply, I’ll try to find some workarounds.
>>
>> Best,
>>
>> Théo
>>
>>> Le 26 mars 2026 à 13:17, Ana Antonini <anacarolina.srantonini at gmail.com> a écrit :
>>>
>>> Hi Theo,
>>> The late TP-AGB is a very tricky evolutionary stage, especially at higher masses.
>>> I'd recommend having a look at Rees&Izzard 2024  https://ui.adsabs.harvard.edu/abs/2024MNRAS.531.4033R/abstract , as they go through the main instabilities TP-AGB models experience, and offer possible solutions as mesa subroutines.
>>> Given  the max residuals  at the surface and the envelope mass, you might be encountering HRI .
>>>
>>> However, their subroutines might not (and at higher masses certainly won't)  solve your problem, merely delaying it to a later point in the TP-AGB.
>>> If your interest lies *in* the TP-AGB stage, the default approach is to carry calculations on for as long as possible and go through post-processing. If you're interested in the WD stage,  the only workaround is to remove the envelope before the instabilities crash the code.
>>>
>>> Good luck!
>>> Ana Antonini
>>>
>>> Théo-Courty via Mesa-users <mesa-users at lists.mesastar.org> escreveu (quinta, 26/03/2026 à(s) 08:52):
>>> Dear Mesa users,
>>>
>>> I’m trying to evolve stars with mass varying between 1-5 Msun from ZAMS to white dwarf cooling track with version r25.12.1. However, I encounter numerical difficulties during the late TP-AGB phase for a 5Msun model with Z=0.012 (remaining mass ~1.25 Msun). The first thermal pulses go smoothly, but the late evolution struggles and the timestep becomes small and goes into a loop of increase/decrease between -3 < log dt <  -2.5. The limiting factor for dt is labelled as max increase until a retry and loops like this. I tried to lower max_timestep_factor to 1.01 but it does not solve the issue. I observe a lot of retries during that part of evolution until crash due to min_timestep_limit.
>>>
>>> I am fairly new to mesa and not very experienced with debugging. Here are typical examples of error messages I get:
>>>
>>> avg resid  0.223E-06     max resid dv_dt     1  0.35367E-03 mix type xx000 avg corr  0.474E-03      max corr lnd   451 -0.21100E-01 mix type 00000  avg+max corr+resid -- give up
>>>
>>> hydro_mtx: logT too small       26764         407           3    6.7237179921904777D-01    5.3757388036508846D+00
>>>
>>> hydro_mtx: logRho too large       26738         421           5    5.6268674362592407D+01   -7.8503397605837941D+00
>>>
>>> I tried to include okay_to_reduce_gradT_excess = .true. which I found used in some inlists online, but with this parameter activated the star shows non-physical behavior (log L > 10^5-6 with few small timesteps) so I keep it to false.
>>>
>>> Please find below the enlists I use which are based on make_co_wd test_suite which I modified to include element diffusion during MS and thermal pulses to remove the enveloppe instead of remove_env routine. Thanks a lot!
>>>
>>> Best regards,
>>>
>>> Théo
>>>
>>> Common inlist for shared parameters
>>> &star_job
>>>      show_log_description_at_start = .false.
>>>
>>>
>>>      change_net = .true.
>>>      new_net_name = 'co_burn_extras.net'
>>>
>>>      ! age start at ZAMS
>>>      set_initial_age = .false.
>>>      set_initial_model_number = .false.
>>>
>>>      num_special_rate_factors = 2
>>>      reaction_for_special_factor(1) = 'r_c12_ag_o16'
>>>      special_rate_factor(1) = 1
>>>      filename_of_special_rate(1) = 'c12ag_deboer_sigma_0p0_2000_Tgrid.dat'
>>>
>>>      reaction_for_special_factor(2) = 'r_he4_he4_he4_to_c12'
>>>      special_rate_factor(2) = 1
>>>      filename_of_special_rate(2) = 'r_he4_he4_he4_to_c12_cf88.txt'
>>>
>>> / ! end of star_job namelist
>>>
>>> &eos
>>>
>>>
>>> / ! end of eos namelist
>>>
>>> &kap
>>>      Zbase = 1.2d-02
>>>
>>>      kap_file_prefix = 'gs98'
>>>      use_Type2_opacities = .true.
>>>
>>> / ! end of kap namelist
>>>
>>> &controls
>>>
>>> initial_mass = 5.00d+00
>>>
>>> ! winds
>>>      cool_wind_full_on_T = 9.99d9
>>>      hot_wind_full_on_T = 1d10
>>>      cool_wind_RGB_scheme = 'Reimers'
>>>      cool_wind_AGB_scheme = 'Blocker'
>>>      RGB_to_AGB_wind_switch = 1d-4
>>>      Reimers_scaling_factor = 0.1d0
>>>      Blocker_scaling_factor = 0.2d0
>>>
>>> ! when to stop
>>>      max_age = 1.4d+10
>>>
>>> ! mesh
>>>      max_allowed_nz = 20000 ! preventing 'fail to adjust mesh'
>>>
>>> ! diffusion
>>>      diffusion_v_max = 1d-5
>>>      do_element_diffusion = .false. ! Activated during MS and settling only
>>>      diffusion_use_full_net = .true.
>>>
>>>      ! Solver controls.
>>>      diffusion_use_cgs_solver = .true.
>>>      diffusion_use_isolve = .true.
>>>      diffusion_rtol_for_isolve = 1d-4
>>>      diffusion_atol_for_isolve = 1d-5
>>>      diffusion_min_X_hard_limit = -1d-4
>>>      diffusion_maxsteps_for_isolve = 1000
>>>      diffusion_isolve_solver = 'ros2_solver'
>>>
>>>      ! Timestep controls to prevent steps that are difficult for diffusion.
>>>      diffusion_steps_limit = 20
>>>      diffusion_steps_hard_limit = 150
>>>      diffusion_iters_limit = 50
>>>      diffusion_iters_hard_limit = 100
>>>
>>> ! convection
>>>      mlt_option = 'Henyey'
>>>      mixing_length_alpha = 1.8
>>>      use_Ledoux_criterion = .true.
>>>      thermohaline_coeff = 1d2
>>>
>>> ! timesteps
>>>      max_years_for_timestep = 1d7
>>>      time_delta_coeff = 1.0d0
>>>      use_gold2_tolerances = .true.
>>>
>>>      varcontrol_target = 1d-3
>>>      dX_nuc_drop_limit = 1d-2
>>>      delta_HR_limit = 0.1
>>>
>>> ! Atm
>>>      atm_option = 'T_tau'
>>>      atm_T_tau_relation = 'Eddington'
>>>      atm_T_tau_opacity = 'iterated'
>>>      atm_T_tau_max_iters = 100
>>>
>>> ! Overshoot
>>>      overshoot_scheme(1) = 'exponential'
>>>      overshoot_zone_type(1) = 'any'
>>>      overshoot_zone_loc(1) = 'any'
>>>      overshoot_bdy_loc(1) = 'any'
>>>      overshoot_f(1) = 0.014
>>>      overshoot_f0(1) = 0.004
>>>
>>> / ! end of controls namelist
>>>
>>> AGB inlist
>>>
>>> &star_job
>>>
>>>      load_saved_model = .true.
>>>      load_model_filename = 'end_he_core_burn.mod'
>>>      save_model_when_terminate = .true.
>>>      save_model_filename = 'end_AGB.mod'
>>>      required_termination_code_string = ''
>>>      set_initial_model_number = .true.
>>>      initial_model_number = 0
>>>      pgstar_flag = .true.
>>>
>>>
>>> / ! end of star_job namelist
>>>
>>> &eos
>>>
>>> / ! end of eos namelist
>>>
>>> &kap
>>>
>>> / ! end of kap namelist
>>>
>>> &controls
>>>      log_directory = 'LOGS_AGB'
>>>
>>> ! energy conservation
>>>      energy_eqn_option = 'eps_grav'
>>>      use_time_centered_eps_grav = .true.
>>>
>>> ! when to stop
>>>      envelope_mass_limit = 1e-3
>>>
>>> ! winds
>>>      max_wind = 1d-2
>>>
>>> ! solver
>>>      retry_hold = 0
>>>      neg_mass_fraction_hold = 3
>>>
>>> ! mesh
>>>      mesh_delta_coeff = 1.0
>>>
>>> ! convection & dragging
>>>      num_cells_for_smooth_gradL_composition_term = 0
>>>      drag_coefficient = 1d0
>>>      min_q_for_drag = 0.98d0
>>>      alpha_semiconvection = 0.1
>>>
>>> ! timesteps
>>>
>>>      varcontrol_target = 1d-4
>>>      delta_lgL_He_limit = 0.005
>>>      lgL_He_burn_min = 2.0
>>>      dX_nuc_drop_limit = 1d-3
>>>      delta_HR_limit = 0.1
>>>      max_timestep_factor = 1.01d0
>>>
>>>
>>>      delta_lgTeff_limit = 0.005
>>>      delta_lgTeff_hard_limit = 0.01
>>>      delta_lgL_limit = 0.02
>>>      delta_lgL_hard_limit = 0.05
>>>
>>>
>>> ! output
>>>      num_trace_history_values = 2
>>>      trace_history_value_name(1) = 'rel_E_err'
>>>      trace_history_value_name(2) = 'log_rel_run_E_err'
>>>
>>>      !report_solver_progress = .true. ! set true to see info about solver iterations
>>>      !report_ierr = .true. ! if true, produce terminal output when have some internal error
>>>      !stop_for_bad_nums = .true.
>>>
>>> / ! end of controls namelist
>>>
>>>
>>> _______________________________________________
>>> mesa-users at lists.mesastar.org
>>> https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.mesastar.org%2Fmailman%2Flistinfo%2Fmesa-users&data=05%7C02%7CEbraheem.Farag%40yale.edu%7Caa0643dd866d4a3868c108de8fe4e8bb%7Cdd8cbebb21394df8b4114e3e87abeb5c%7C0%7C0%7C639106410514248929%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=hyySX%2F3je9XXUCC%2BhCsQnBAmeo0zbl7k%2F7bv5lc0RAs%3D&reserved=0<https://lists.mesastar.org/mailman/listinfo/mesa-users>
>>>
>>
>> _______________________________________________
>> mesa-users at lists.mesastar.org
>> https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.mesastar.org%2Fmailman%2Flistinfo%2Fmesa-users&data=05%7C02%7CEbraheem.Farag%40yale.edu%7Caa0643dd866d4a3868c108de8fe4e8bb%7Cdd8cbebb21394df8b4114e3e87abeb5c%7C0%7C0%7C639106410514284133%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=r0M%2BGB%2B7A19tOsqvb5qR5gQNYOacRFmYKUrDr8Xp7Gs%3D&reserved=0<https://lists.mesastar.org/mailman/listinfo/mesa-users>
>>
>

_______________________________________________
mesa-users at lists.mesastar.org
https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.mesastar.org%2Fmailman%2Flistinfo%2Fmesa-users&data=05%7C02%7CEbraheem.Farag%40yale.edu%7Caa0643dd866d4a3868c108de8fe4e8bb%7Cdd8cbebb21394df8b4114e3e87abeb5c%7C0%7C0%7C639106410514323128%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=cH7IPEG3exaylkmLnnzXnsfun2uSnXEC9kip%2FCEopC8%3D&reserved=0<https://lists.mesastar.org/mailman/listinfo/mesa-users>

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