Residual warnings -- can't evolve a 1 Mjup planet past 1 Gyr.

[Rob685]

Hello,

I'm very new to MESA and Fortran, so this issue is definitely something that I'm not doing right.

Basically, I'm trying to evolve a 1 Mjup planet usingcreate_initial_model, and following the settings of Paxton et al. (2013). After about 1 Gyr, I get these warnings:

residual_norm > tol_residual_norm        
max_residual > tol_max_residual         
correction_norm > tol_correction_norm*coeff          
max_abs_correction > tol_max_correction*coeff         
retry: avg+max corr+resid -- give up in solver

and the model just doesn't converge after that.

I've set the following lines under inlist_project:

create_initial_model = .true.
mass_in_gm_for_create_initial_model = 1.89861E30 ! 1 M_Jup
radius_in_cm_for_create_initial_model = 14E9 ! 2 R_Jup

and haven't changed much besides that, except for use_CMS = .true., and setting the freedman opacities (like Section 2 of Paxton et al. (2013)).

Can anyone point me in the right direction regarding this issue?

Thank you,

Rob

[fxt44]

hi rob,

there is not enough information here for anyone to help. is this a modification of a test suite example, a mesa directory from a zenodo archive, or something else? its also useful to know which version of mesa is being used. plots can also help you assess the model.

fxt

[Rob685]

Apologies, this is a modification of the star directory in the work folder in MESA version r.22.05.1 (downloaded from zenodo archive). I've plotted up quantities of interest (Teff, radius, age, etc. ) but nothing seems out of the ordinary, except that the models do not go past 1 Gyr.

[warrickball]

Could you attach your entire inlist, please?

[Rob685]

Hm, I can't upload it but I can paste the inlist_project, which is what I've been editing:

! For the sake of future readers of this file (yourself included),
! ONLY include the controls you are actually using.  DO NOT include
! all of the other controls that simply have their default values.

&star_job
  ! see star/defaults/star_job.defaults

  ! begin with a pre-main sequence model
    create_pre_main_sequence_model = .false.
    create_initial_model = .true.
    mass_in_gm_for_create_initial_model = 1.89861E30 ! 1 M_Jup
    radius_in_cm_for_create_initial_model = 14E9 ! 2 R_Jup
    !rhoc = 15.0 !gcm^-3

  ! save a model at the end of the run
    save_model_when_terminate = .true.
    save_model_filename = '1Mj_planet.mod'

  ! display on-screen plots
    pgstar_flag = .true.

/ ! end of star_job namelist

&eos
  ! eos options
  ! see eos/defaults/eos.defaults
  ! Rob: these are all controls in the documentation
  use_CMS = .true

/ ! end of eos namelist

&kap
  ! kap options
  ! see kap/defaults/kap.defaults
  use_Type2_opacities = .false.
  kap_lowT_prefix = 'lowT_Freedman11'
  kap_file_prefix = 'gs98'
  Zbase = 0.02

/ ! end of kap namelist

&controls
  ! see star/defaults/controls.defaults

  ! starting specifications
    initial_mass = 9.45E-4 ! in Msun units
    initial_z = 0.02
    initial_y = 0.27

  ! when to stop

    ! stop when the star nears ZAMS (Lnuc/L > 0.99)
    ! Lnuc_div_L_zams_limit = 0.99d0
    ! stop_near_zams = .true.

    ! stop when the center mass fraction of h1 drops below this limit
    !xa_central_lower_limit_species(1) = 'h1'
    !xa_central_lower_limit(1) = 1d-3

    !Teff_lower_limit = 167
    ! max_age = 1E10

  ! wind

  ! atmosphere
    atm_option = 'T_tau'
    atm_T_tau_relation = 'Eddington'
    atm_T_tau_opacity = 'fixed'

  ! rotation

  ! element diffusion

  ! mlt

  ! mixing

  ! timesteps

    !dt_years_for_steps_before_max_age = 1E7

  ! mesh

  ! solver
     ! options for energy conservation (see MESA V, Section 3)
     energy_eqn_option = 'dedt'
     use_gold_tolerances = .true.

  ! output

/ ! end of controls namelist

[warrickball]

What happens if you don't use CMS, i.e., leave use_CMS at its default value of .false.?

[Rob685]

I get the same residual warnings and the planet model still can't get past 1 Gyr.

All quantities remain the same after the time step where I get the warning, and I also get:

retry: max correction jumped -- give up in solver.

[rjfarmer]

Whats the surface temperature of your planet at 1Gyr?

[warrickball]

I just ran this case and can reproduce the result. The evolution suddenly starts struggling at an effective temperature of 158 K, which is basically 10^2.2, so I think you've found the edge of an EoS table.

[warrickball]

I don't think this is something we can solve without somehow extending the EoS into this very low temperature regime. Unless one of the EoS experts (@fxt44? @evbauer?) tells me otherwise or this is a regression (i.e., did this ever work before?), I'll close this soon.

[earlbellinger]

This seems to me quite similar to the make-planets test case and so it may be a regression. https://docs.mesastar.org/en/release-r22.11.1/test_suite/make_planets.html#make-planets

[warrickball]

I think that's kept hot by being irradiated: https://github.com/MESAHub/mesa/blob/3c9bfba29f2aaa70f407f2c716e40e64666eb0c1/star/test_suite/make_planets/inlist_evolve#L50-L52

[evbauer]

Perhaps we might have lost a blend that makes that boundary smoother between SCVH and what's below it in logT. It looks like there's a hard edge there right now.

Previously this would have fallen back to HELM. Now it falls back to an even simpler ideal gas. But is either meaningful here?

[evbauer]

Perhaps there is an argument for at least adding a blend there to make the numerics more friendly, while letting the user decide how much they want/need to believe the physics they're getting out.

[fxt44]

But is either meaningful here?

no. maybe ideal gas if T is below ionization but above molecules and the density is not too large.

[warrickball]

I don't mind the code basically crashing when it goes outside the bounds of what it's capable of. Some more helpful information about what's gone wrong (e.g. "You've gone off the EoS table to the fallback") would be good but from a practical point of view, unless we're going to do something about this issue imminently, I'm inclined to close it.

[fxt44]

i don't see myself imminently working on this.

[evbauer]

Part of the reason we developed ideal gas as a fallback is that the solver will very occasionally wander into unexpected territory during the first few Newton iterations, and we want to have something to gently point it back in the right direction. Previously we had some hard stops for very bad EOS calls, and this would cause MESA models to crash when instead they should really just keep iterating, or at worst retry. So for fundamental inputs like EOS and kap that get called many times during the solver iterations, I would caution against crashing even for pretty crazy inputs (e.g. logT = 15, logRho = -20).

But anyway, I'm not opposed to closing this for now, and keeping in mind that maybe we should add a small blend over to ideal in some of those spots where there's a hard switch.

[warrickball]

Part of the reason we developed ideal gas as a fallback is that the solver will very occasionally wander into unexpected territory during the first few Newton iterations, and we want to have something to gently point it back in the right direction.

This was my understanding too but this looks to me like a case where the solver is either doomed or, even if it isn't, we wouldn't believe the converged model anyway.

Let's see if this comes up again.