[Mesa-users] Inconsistency between Omega/Omega_crit reported by MESA and that computed with Paxton et al. (2019) equations
Pablo Marchant
pamarca at gmail.com
Tue Jul 14 02:36:29 EDT 2020
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> 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>
> *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> 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
>
> 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
> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20200714/ce9129d5/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: star_utils.f90
Type: text/x-fortran
Size: 120462 bytes
Desc: not available
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20200714/ce9129d5/attachment.bin>
More information about the Mesa-users
mailing list