[Mesa-users] Reproducing Paxton et al. 2013 fig. 1

[Emily Sandford]

Hello,

I'm a new MESA user, attempting to use r24.03.1 to model brown dwarf and 
giant planet evolution. As a first step I want to reproduce the models 
plotted in Paxton et al. 2013 figure 1, which are initiated with 
create_initial_model with radius equal to 5 R_J and masses ranging from 
0.2 to 20 M_J, and then evolved for 20 Gyr.

The only controls enumerated in the paper are:
&star_job
    create_initial_model = .true.
    mass_in_gm_for_create_initial_model = 3.7972d31  ! 20 m_jupiter;  
m_jupiter = 1.8986d30
    radius_in_cm_for_create_initial_model = 3.596d10 ! 5 r_jupiter; 
r_jupiter = 7.192d9
/ ! end of star_job namelist

&eos
    eosDT_file_prefix = 'mesa' ! replaces deprecated eos_file_prefix = 
'mesa', as I understand it
/ ! end of eos namelist

&kap
    Zbase = 0.02d0             ! not mentioned explicitly in Paxton et 
al. 2013, but now seems to be necessary
    kap_file_prefix = 'gs98'
    kap_lowT_prefix = 'lowT_Freedman11'
/ ! end of kap namelist

&atm
    which_atm_option = 'simple_photosphere'
/ ! end of atm namelist

&controls
    initial_y = 0.27d0
    initial_z = 0.02d0
    mixing_length_alpha = 2.0d0
    max_age = 2.0d10
/ ! end of controls namelist

This works fine for initial mass 20 M_J; however, it fails almost 
immediately for initial mass <= 14 M_J (the initial model is created, 
but then the evolution stalls at age ~a few years). I looked at the 
make_brown_dwarf test suite example for inspiration, and found that 
adding

convergence_ignore_equL_residuals = .true.

allows models with initial masses down to 6.9 M_J to evolve for the full 
20 Gyr without problems. However, for initial masses <= 4.8 M_J, the 
evolution stalls once again (at late ages--for the 4.8 M_J model, it's 
at 12.4 Gyr), and I don't know how to fix this. I've tried turning on 
the other controls used in make_brown_dwarf (mlt_make_surface_no_mixing 
= .false. , make_gradr_sticky_in_solver_iters = .true. , 
energy_eqn_option = 'dedt'), but none of them have any effect. I've also 
played around with the make_planets test suite example, turning off the 
step where an inert core is added, but with no success.

The one thing that worked for the 4.8 M_J model was to adjust the 
initial radius downward to 2 R_J, but then that choice of radius also 
fails for initial masses <= 3.4 M_J (and I'm hoping there's an easier 
way forward than to keep hand-tuning the initial radius for every new 
choice of initial mass).

Does anyone have advice for how to make progress from here, ideally in a 
way that allows me to run a suite of models at different initial masses 
but with the same initial radius, as with the models plotted in 
Paxton+2013? I'm happy to provide more detailed debugging information if 
it's helpful, but given that these models have successfully been run 
before, I'm hoping it's just a matter of choosing the right inlist 
settings.

Thanks!
Emily Sandford
sandford at strw.leidenuniv.nl