Pablo,
Thanks for the reply and the suggestion. We did another run with M = 1 Msun, omega_initial = 0.4, but this time with a maximum time step of 1E5 years. The result is much the same as before with the larger timestamp. I copy below the first few zones of profile500.data (model_number 37155, age 3.68E9 y). For zone 1, omega is reported to be 0.43529 when it should be 0.45893, so the reported value was 5.15% low. This is close to my previously found trend: 0.25*0.45893^2 = 0.0527, predicting 5.27% short. So the time step is not the issue.
Larry
Profile500.data:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
model_number num_zones initial_mass initial_z star_age time_step Teff photosphere_L photosphere_r center_eta center_h1 center_he3 center_he4 center_c12 center_n14 center_o16 center_ne20 star_mass star_mdot star_mass_h1 star_mass_he3 star_mass_he4 star_mass_c12 star_mass_n14 star_mass_o16 star_mass_ne20 he_core_mass c_core_mass o_core_mass si_core_mass fe_core_mass neutron_rich_core_mass tau10_mass tau10_radius tau100_mass tau100_radius dynamic_time kh_timescale nuc_timescale power_nuc_burn power_h_burn power_he_burn power_neu burn_min1 burn_min2 time_seconds version_number compiler build MESA_SDK_version math_backend date
37155 770 1.0000000000000000E+000 2.0000000000000000E-002 3.6776696980139370E+009 1.0000000000000000E+005 5.5064198195436538E+003 8.2480411522123964E-001 9.9929935376573187E-001 -1.5577637218453766E+000 4.4884867064318057E-001 1.7431050803908474E-005 5.3058379483133655E-001 1.7661882834157228E-005 5.1466402311758520E-003 9.1960164866795834E-003 2.0994604783345795E-003 1.0000000000000000E+000 0.0000000000000000E+000 6.7357678919379660E-001 7.6976808273400042E-004 3.0549348367536516E-001 2.4801365851414444E-003 2.1318769117749047E-003 9.3581606771990823E-003 2.0994604783341284E-003 0.0000000000000000E+000 0.0000000000000000E+000 0.0000000000000000E+000 0.0000000000000000E+000 0.0000000000000000E+000 0.0000000000000000E+000 9.9999999986903432E-001 9.9906003390403830E-001 9.9999999977714160E-001 9.9893083391412185E-001 1.0001772717994632E+004 2.8479549376492925E+007 1.2124090817997034E+010 8.2468244158947501E-001 8.2468244158947501E-001 2.2023040521293154E-045 1.8673900154586386E-002 5.0000000000000000E+001 1.0000000000000000E+003 1.1606045192153581E+017 "12778" "gfortran" "9.2.0" "x86_64-linux-20.3.2" "CRMATH" "20200713"
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
zone mass logR logT logRho logP x_mass_fraction_H y_mass_fraction_He z_mass_fraction_metals x luminosity conv_L_div_L omega log_omega log_J_inside i_rot j_rot v_rot r_polar r_equatorial r_e_div_r_p omega_crit omega_div_omega_crit
1 1.0000000000000000E+000 -3.0439344182504120E-004 3.7410950107712910E+000 -6.7400221415136272E+000 4.8053562336963447E+000 6.9999569506667336E-001 2.8000430184701625E-001 2.0000003086310336E-002 6.9999569506667336E-001 8.2480411522123964E-001 0.0000000000000000E+000 2.7344741784363536E-004 -3.5631261733006983E+000 5.0305441558992349E+001 3.6132762893748577E+021 9.8804107128518707E+017 1.9018059112971798E+002 9.3637241738738453E-001 1.0349058902269257E+000 1.1052289356348877E+000 6.2820147949070896E-004 4.3528617294139882E-001
2 9.9999999999949996E-001 -3.0495496642397277E-004 3.7416163217212626E+000 -6.7370369840172959E+000 4.8088642728522153E+000 6.9999569506667336E-001 2.8000430184701625E-001 2.0000003086310336E-002 6.9999569506667336E-001 8.2480411522123964E-001 0.0000000000000000E+000 2.7344741784363536E-004 -3.5631261733006983E+000 5.0305441558990232E+001 3.6132654001746047E+021 9.8803809366149542E+017 1.9018034523432775E+002 9.3637141594964823E-001 1.0349044147091435E+000 1.1052285418810710E+000 6.2820269784908712E-004 4.3528532873210540E-001
From: Pablo Marchant <pamarca at gmail.com>
Sent: Monday, July 13, 2020 2:44 PM
To: Larry Molnar <lmolnar at calvin.edu>
Cc: mesa-users at lists.mesastar.org
Subject: Re: [Mesa-users] Inconsistency between Omega/Omega_crit reported by MESA and that computed with Paxton et al. (2019) equations
Hi Larry,
thanks for reporting this, I'll give it a check in the following days. One thing you have to be aware of is that omega (by which I mean the little omega used in Paxton+ 2019 and not the angular frequency) is computed explicitly at the beginning of the step, we are not implicitly solving its value in a consistent way inside the newton solver as the structure of the star changes. The other detail to be aware of is that angular momentum mixing is done separate from the structure of the star. I believe a combination of these two things should be related to the variation you see in reported values.
A simple test I'd ask you to do to confirm this is to see if the large variation holds if you use a very small timestep, which you can do with max_timestep. In this case the explicit calculation of omega should not matter, and you should find near consistent values. If you keep finding a large discrepancy in omega from step to step while the entire structure of the star and its rotation remain static, then there's likely a deeper issue.
In any case, these subtleties should be better documented in the future.
Cheers
On Mon, Jul 13, 2020 at 8:21 PM Larry Molnar <lmolnar at calvin.edu<mailto:lmolnar at calvin.edu>> wrote:
I am running the latest MESA based on the revisions described in Paxton et al. (2019): MESA R12778, SDK 20.3.2, on a Linux operating system.
I am interested in rotation of main sequence stars with intermediate mass, so I ran a case with 1.0 Msun and default (solar) abundances with several initial values of the frequency ratio omega := Omega_Psi/Omega_crit [p. 28 of Paxton et al. (2019)]. I attach the inlist used and copy the first few lines of output from profile39.data (at which age = 6.79 Gy).
The run completed fine and the results largely look sensible. The left side of Equation 39 of section 4 of Paxton et al. (2019) contains j_rot, m_Psi, and r_Psi, the quantities needed to determine omega for a given zone (and hence the quantities r_e and i_rot which are found on the right side of the equation). As a consistency check I computed both sides separately for a number of zones and found significant disagreement (up to several percent, depending on the zone). So I used Equation B30 to implicitly find the value of omega that is consistent (and the values for other things like r_e, r_p, etc.). I found a slightly larger value of omega that was consistent with everything. Using this revised omega, I checked separately and found the profile values of Omega_Psi, r_e, r_p, and i_rot are all consistent with the left side of Eqn. 39 while omega, v_rot, and Omega_crit are not. I repeated this for different zones and for two initial values of omega (0.3 and 0.4). I find the omega reported by the profile is approximately equal to the following function of the omega from the equation: omega_Eqn*(1 – 0.25*omega_Eqn^2).
For the example below, omega of zone 1 is reported to be 0.32961 when it should be 0.33906, a difference of 2.8%.
In summary, the correct omega ratio seems to be computed and used for some things (like r_e, r_p, i_rot, Omega_Psi), but something else is being reported (and being used to compute v_rot and Omega_crit).
Prof. Larry Molnar
Department of Physics and Astronomy
Calvin University
1734 Knollcrest Circle SE
Grand Rapids, MI 49546-4403 USA
E-mail: lmolnar at calvin.edu<mailto:lmolnar at calvin.edu>
Phone: 616-526-6341
Profile39.data: first three zones plus header
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
model_number num_zones initial_mass initial_z star_age time_step Teff photosphere_L photosphere_r center_eta center_h1 center_he3 center_he4 center_c12 center_n14 center_o16 center_ne20 star_mass star_mdot star_mass_h1 star_mass_he3 star_mass_he4 star_mass_c12 star_mass_n14 star_mass_o16 star_mass_ne20 he_core_mass c_core_mass o_core_mass si_core_mass fe_core_mass neutron_rich_core_mass tau10_mass tau10_radius tau100_mass tau100_radius dynamic_time kh_timescale nuc_timescale power_nuc_burn power_h_burn power_he_burn power_neu burn_min1 burn_min2 time_seconds version_number compiler build MESA_SDK_version math_backend date
840 788 1.0000000000000000E+000 2.0000000000000000E-002 6.7862319465169029E+009 1.0974467458049601E+008 5.6913137780955512E+003 1.1106458896908029E+000 1.0854793927473201E+000 -1.4234108223730109E+000 1.7892625342288088E-001 1.2875481967828730E-006 8.0076368174195711E-001 3.0901329101104673E-005 6.8050315038511392E-003 7.2830595800253044E-003 2.0994604783340942E-003 1.0000000000000000E+000 0.0000000000000000E+000 6.4251076451589306E-001 9.1978547547992594E-004 3.3638630678708026E-001 2.3204305753635767E-003 2.3422350540312545E-003 9.3306927181358093E-003 2.0994604783340300E-003 0.0000000000000000E+000 0.0000000000000000E+000 0.0000000000000000E+000 0.0000000000000000E+000 0.0000000000000000E+000 0.0000000000000000E+000 9.9999999986641774E-001 1.0852251060789109E+000 9.9999999977220166E-001 1.0850847080353851E+000 1.1323115727543849E+004 1.9470734833037574E+007 9.0037698719471607E+009 1.1106109430736193E+000 1.1106109430736197E+000 1.9379726463585424E-041 3.1117358618271040E-002 5.0000000000000000E+001 1.0000000000000000E+003 2.1416092559439270E+017 "12778" "gfortran" "9.2.0" "x86_64-linux-20.3.2" "CRMATH" "20200707"
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
zone mass logR logT logRho logP x_mass_fraction_H y_mass_fraction_He z_mass_fraction_metals x luminosity conv_L_div_L omega log_omega log_J_inside i_rot j_rot v_rot r_polar r_equatorial r_e_div_r_p omega_crit omega_div_omega_crit
1 1.0000000000000000E+000 3.5621583017451837E-002 3.7553873899760544E+000 -6.8226384760895904E+000 4.7364521258960668E+000 6.9837283051018251E-001 2.8162185819852292E-001 2.0005311291294570E-002 6.9837283051018251E-001 1.1106458896908029E+000 0.0000000000000000E+000 1.8289678787990470E-004 -3.7377939217421856E+000 5.0141558734190070E+001 4.0452322195411635E+021 7.3985997918237632E+017 1.3817339257946193E+002 1.0466378035378214E+000 1.1061445120438380E+000 1.0568551110086730E+000 5.5489221113985861E-004 3.2960777644400241E-001
2 9.9999999999949996E-001 3.5621053101314769E-002 3.7558294668053462E+000 -6.8198510069919314E+000 4.7396826558142946E+000 6.9837283051018251E-001 2.8162185819852292E-001 2.0005311291294570E-002 6.9837283051018251E-001 1.1106458896908031E+000 0.0000000000000000E+000 1.8289678787990470E-004 -3.7377939217421856E+000 5.0141558734187761E+001 4.0452210785591907E+021 7.3985794153255962E+017 1.3817322398356947E+002 1.0466367126293932E+000 1.1061430561848453E+000 1.0568548215798377E+000 5.5489322674071766E-004 3.2960717317489624E-001
3 9.9999999999862488E-001 3.5620131676504144E-002 3.7563968080601273E+000 -6.8163424224394999E+000 4.7437598235062506E+000 6.9837283051018240E-001 2.8162185819852292E-001 2.0005311291294681E-002 6.9837283051018240E-001 1.1106458896908031E+000 0.0000000000000000E+000 1.8289678787990470E-004 -3.7377939217421856E+000 5.0141558734183725E+001 4.0452026442495694E+021 7.3985456995654310E+017 1.3817293082741162E+002 1.0466346781989255E+000 1.1061406031617111E+000 1.0568545321517389E+000 5.5489490448713092E-004 3.2960617659473646E-001
_______________________________________________
mesa-users at lists.mesastar.org<mailto:mesa-users at lists.mesastar.org>
https://lists.mesastar.org/mailman/listinfo/mesa-users<https://urldefense.proofpoint.com/v2/url?u=https-3A__lists.mesastar.org_mailman_listinfo_mesa-2Dusers&d=DwMGaQ&c=4rZ6NPIETe-LE5i2KBR4rw&r=rohBEfdOjnrgwB12cgkmE-lSVHI5uTtNMjZ7eHCNBFw&m=E57QJO6co1xKKB_2kGtLjamqEByBy9pqdV6tO3ue9Xk&s=wHon_ob5c77cGITVX-tvzFee-7FKFNZvJaN9AckNBoY&e=>
--
Pablo Marchant Campos
M.Sc on Astrophysics, Universidad Católica de Chile
PhD on Astrophysics, Argelander-Institut für Astronomie, Universität Bonn
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20200713/0fc1f13c/attachment.htm>