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

Thu Aug 8 11:46:02 UTC 2024

On 2024-07-12 23:05, Francis Timmes wrote:
> hi emily,
> 
> yes, this is a known issue. see this thread
> 
> https://lists.mesastar.org/pipermail/mesa-users/2024-March/015070.html
> 
> and this potential workaround
> 
> https://lists.mesastar.org/pipermail/mesa-users/2024-April/015133.html
> 
> some preliminary mesa-dev investigations suggested the phenomena may
> be due to changes since 2013 how the low temperature eos is handled.
> if you find the problem persists, perhaps file an issue on mesa github
> https://github.com/MESAHub/mesa .
> 
> fxt
> 
> 
> 
> 
> 
>> On Jul 10, 2024, at 8:24 AM, Emily Sandford via Mesa-users 
>> <mesa-users at lists.mesastar.org> wrote:
>> 
>> 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
>> _______________________________________________
>> mesa-users at lists.mesastar.org
>> https://lists.mesastar.org/mailman/listinfo/mesa-users
>> 

Dear Frank,

Thanks for this, it was very helpful--I've done some more investigating, 
and I will file a github issue, but I wanted to post my findings here 
too in case someone else runs into the same problem in the future. My 
results are consistent with the hypothesis that a change to the low-EOS 
behavior post-r12115 (released October 2019) is responsible for the 
issues I was having.

(Note that, in the below, I am only testing at the specific initial mass 
values plotted in Paxton et al. 2013 fig 1: 20.0, 14.0, 9.8, 6.9, 4.8, 
3.4, 2.3, 1.6, 1.1, 0.82, 0.57, 0.4, 0.28, and 0.2 M_J.)

In brief, I have found that:
-Adapting the r24.03.01 make_brown_dwarf test suite case works fine for 
initial masses >= 1.6 M_J. At the lower end of that mass range, it is 
necessary to tune the initial radius downward, in line with the 
suggestion in this thread: 
https://lists.mesastar.org/pipermail/mesa-users/2024-April/015133.html . 
(Unfortunately, the other suggestion proffered in that thread, to change 
kap_lowT_prefix from 'lowT_Freedman11' to 'lowT_fa05_gs98', hasn't 
worked at all for me; I get instant 'failed in kap_get 
get_kap_from_rhoT; temp too low in integration' errors, regardless of 
initial mass or radius.)
-At initial masses <= 1.1 M_J but > 0.2 M_J, the only way to make 
progress is to turn on irradiation to heat the outer layers of the 
planet and prevent temperatures from dropping below ~150K. (Thanks to 
Henrik Knierim for this suggestion!)
-At initial mass = 0.2 M_J, the create_initial_model step fails, 
regardless of choice of initial radius or irradiation parameters. 
Nothing I have tried in r24.03.01 works at this mass.
-However, reverting to MESA r12115 allows objects with initial masses 
 >=0.2 M_J to evolve to 20Gyr without irradiation. It is necessary to set use_eosDT2 = .false. and use_eosELM = .true., as in the r12115 make_brown_dwarf test suite example.

I have attached the inlists I used for r24.03.01 (split into "create" 
and "evolve" steps in order to try to isolate the create_initial_model 
problem at extremely low masses, but without success) and for r12115 in 
case this is helpful to anyone.

Best,
Emily
-------------- next part --------------
! inlist_make_brown_dwarf

&star_job

   show_log_description_at_start = .false.

   save_model_when_terminate = .true.
   save_model_filename = "planet_create_20.0_MJ_5.0_RJ.mod"
   required_termination_code_string = 'max_age'

   ! new adiabatic, contracting initial model
   create_initial_model = .true.
   mass_in_gm_for_create_initial_model = 3.79625d31 
   radius_in_cm_for_create_initial_model = 3.5746d10
   initial_model_relax_num_steps = 10
   max_tries_for_create_initial_model = 200

   set_initial_age = .true.
   initial_age = 0.d0 ! in years

/ ! end of star_job namelist

&eos
   eosDT_file_prefix = 'mesa'
/ ! end of eos namelist

&kap
   Zbase = 0.02d0
   kap_file_prefix = 'gs98'
   kap_lowT_prefix = 'lowT_Freedman11'     
/ ! end of kap namelist

&atm
   atm_option = 'T_tau'
   atm_T_tau_relation = 'Eddington'
   atm_T_tau_opacity = 'fixed'
/ ! end of atm namelist

&controls

   initial_y = 0.27d0
   initial_z = 0.02d0
   mixing_length_alpha = 2.0d0
   max_age = 5.0d1

   column_depth_for_irradiation = 300d0
   irradiation_flux = 5.0d4 ! erg/cm^2/s

   use_gold2_tolerances = .false. ! wait a few steps
   steps_before_use_gold2_tolerances = 31

   warning_limit_for_max_residual = 1d99 ! off until evolve
   warn_when_large_virial_thm_rel_err = 1d99 ! off until evolve

   limit_for_rel_error_in_energy_conservation = 1d-5
   hard_limit_for_rel_error_in_energy_conservation = 1d-3

   convergence_ignore_equL_residuals = .true.

   max_resid_jump_limit = 1d18
   max_corr_jump_limit = 1d18

/ ! end of controls namelist

&pgstar
/ ! end of pgstar namelist
-------------- next part --------------
! inlist_evolve

&star_job

   show_log_description_at_start = .false.

   load_saved_model = .true.
   load_model_filename = "planet_create_20.0_MJ_5.0_RJ.mod"   

   save_model_when_terminate = .true.
   save_model_filename = "planet_evolve_20.0_MJ_5.0_RJ.mod"   
   required_termination_code_string = 'max_age'

   set_initial_cumulative_energy_error = .true.
   new_cumulative_energy_error = 0d0

/ ! end of star_job namelist

&eos
   eosDT_file_prefix = 'mesa'
/ ! end of eos namelist

&kap
   Zbase = 0.02d0
   kap_file_prefix = 'gs98'
   kap_lowT_prefix = 'lowT_Freedman11'
/ ! end of kap namelist

&atm
   atm_option = 'T_tau'
   atm_T_tau_relation = 'Eddington'
   atm_T_tau_opacity = 'fixed'
/ ! end of atm namelist

&controls

   initial_y = 0.27d0
   initial_z = 0.02d0
   mixing_length_alpha = 2.0d0
   max_age = 2.0d10

   column_depth_for_irradiation = 300d0
   irradiation_flux = 5.0d4 ! erg/cm^2/s

   convergence_ignore_equL_residuals = .true.
   make_gradr_sticky_in_solver_iters = .true.
   mlt_make_surface_no_mixing = .false.   
   energy_eqn_option = 'dedt'
   use_gold2_tolerances = .true.

/ ! end of controls namelist

&pgstar
/ ! end of pgstar namelist
-------------- next part --------------
! inlist_create+evolve
&star_job

   show_log_description_at_start = .false. 
   create_initial_model = .true.

   save_model_when_terminate = .true.
   save_model_filename = 'planet_20.0_Mj_5.0_Rj.mod'

   mass_in_gm_for_create_initial_model = 3.79625d31
   radius_in_cm_for_create_initial_model = 3.5746d10

   kappa_lowT_prefix = 'lowT_Freedman11'

/ ! end of star_job namelist

&controls

   initial_y = 0.27d0
   initial_z = 0.02d0
   mixing_length_alpha=2.0d0
   max_age = 2.0d10

   !check_elm_helm_agreement = .true.

   use_eosDT2 = .false. !OPAL/SCVH
   use_eosELM = .true.  !HELM

   report_eos_settings_at_start_of_run = .true.

   use_gold_tolerances = .false.
   warn_when_stop_checking_residuals = .false.
   warn_when_large_rel_run_E_err = 0.05d0

   num_trace_history_values = 2
   trace_history_value_name(1) = 'rel_E_err'
   trace_history_value_name(2) = 'log_rel_run_E_err'

   max_resid_jump_limit = 1d12
   max_corr_jump_limit = 1d12

/ ! end of controls namelist

&pgstar
/ ! end of pgstar namelist