[Mesa-users] Inconsistency between Omega/Omega_crit reported by MESA and that computed with Paxton et al. (2019) equations
Larry Molnar
lmolnar at calvin.edu
Tue Jul 14 15:33:51 EDT 2020
Hi Pablo,
Sorry, I didn't realize I was looking at the arXiv version.
I modified the Fortran and reran it. Now both the Omega_crit and Omega/Omega_crit report sensible values.
Curiously, the v_rot is still incorrect. This one is the easiest to check (v_rot should be Omega*r_equatorial). It is not clear to me where the value is set in the code.
Larry
________________________________
From: Pablo Marchant <pamarca at gmail.com>
Sent: Tuesday, July 14, 2020 1:12 PM
To: Larry Molnar <lmolnar at calvin.edu>
Cc: mesa-users at lists.mesastar.org <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,
was using the equation numbering of the published version at the ApJ which differs for the appendices. In the arXiv version its B10.
Regarding the file I sent, I dumbly made a change before sending it and undid the relevant part, you should switch the last line in the omega_crit function to
omega_crit = sqrt(gamma_factor*s% cgrav(k)*s% m_grav(k)/pow3(rmid))
On Tue, Jul 14, 2020 at 6:56 PM Larry Molnar <lmolnar at calvin.edu<mailto:lmolnar at calvin.edu>> wrote:
Pablo,
Thanks for your note.
1. I am aware of the factor with the Eddington luminosity, but as you say it is negligible in this case. It is good to know the ratio can be given not including it.
2. I am aware that the relation between r_Psi and r_equatorial is a fit. (There is not an equation 56 in Paxton et al. 2019. Equations B8, B10, and B11 give three versions of this fit. I am not sure which of these three is used by the code, but they all agree to within a few tenths of a percent as advertised.)
3. We tried installing and running the star_utils.f90 that you sent, but the output was unchanged. Looking at the omega_crit function I see that ‘rmid’ is now defined in the function relative to r_equatorial. However, the final line defining omega_crit does not use this rmid. Instead it uses s%rmid(k) which is based on r_Psi.
Larry
From: Pablo Marchant <pamarca at gmail.com<mailto:pamarca at gmail.com>>
Sent: Tuesday, July 14, 2020 2:36 AM
To: Larry Molnar <lmolnar at calvin.edu<mailto:lmolnar at calvin.edu>>
Cc: mesa-users at lists.mesastar.org<mailto: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,
found the problem. This is an issue with output, so there's no impact in your structure calculation. This comes from the code that writes the value of omega_crit in the profiles using r_phi instead of r_e. You can actually show using Eq. (56) of Paxton+ (2019) that your equation omega_Eqn*(1 – 0.25*omega_Eqn^2) is exactly what should happen in that case, with additional terms in the parentheses of order omega^4 and above. The attached star_utils.f90 can fix this issue, it goes into star/private and you should recompile and export star after applying it (./mk && ./export in the star folder). Check for the function omega_crit to see the change.
Two extra things you should be aware of:
- the value used to compute the rotation corrections can be written to profiles as w_div_w_crit_roche, which is different than omega_div_omega_crit. The latter includes the contribution from the Eddington factor. For your 1 Msun model the Eddington factor is pretty puny so it does not make much of a difference.
- even with this fix you'll still see a small difference between the w_div_w_crit_roche you compute independently and that from the routine omega_crit. This is because the equation for the equatorial radius (56 in Paxton+ 2019) is a fit with a small error,
Will make sure that this fix is applied in future versions. Thanks for reporting, and let me know if there's some other lingering issues :).
Cheers
On Tue, Jul 14, 2020 at 1:50 AM Larry Molnar <lmolnar at calvin.edu<mailto:lmolnar at calvin.edu>> wrote:
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<mailto:pamarca at gmail.com>>
Sent: Monday, July 13, 2020 2:44 PM
To: Larry Molnar <lmolnar at calvin.edu<mailto:lmolnar at calvin.edu>>
Cc: mesa-users at lists.mesastar.org<mailto: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
--
Pablo Marchant Campos
M.Sc on Astrophysics, Universidad Católica de Chile
PhD on Astrophysics, Argelander-Institut für Astronomie, Universität Bonn
--
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/20200714/cda11591/attachment.htm>
More information about the Mesa-users
mailing list