[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 15:42:12 EDT 2020
That one has the same issue now that I see. The values for the profiles are
computed in
$MESA_DIR/star/private/profile_getval.f90
and you can see there that v_rot is simply computed as omega*r_phi. One
extra thing to fix.
On Tue, Jul 14, 2020 at 9:33 PM Larry Molnar <lmolnar at calvin.edu> wrote:
> 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> 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>
> *Sent:* Tuesday, July 14, 2020 2:36 AM
> *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,
>
>
>
> 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
>
>
>
> --
> 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/a7733e98/attachment.htm>
More information about the Mesa-users
mailing list