[Mesa-users] Resolving TP-AGB

Wed Apr 1 11:49:01 UTC 2026

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://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://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://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://lists.mesastar.org/mailman/listinfo/mesa-users
>>> 
>> 
>> _______________________________________________
>> mesa-users at lists.mesastar.org
>> https://lists.mesastar.org/mailman/listinfo/mesa-users
>> 
> 