[Mesa-users] Question about GR corrections for SMS
Nicholas Herrington
Nicholas.Herrington at myport.ac.uk
Mon Jul 27 10:41:49 EDT 2020
Hi,
Thanks for the corrections Rob, I understand the problems now.
I have done some quick tests with the new subroutine and sometimes when
I re-run the same star, at 10^2 or 10^3 Msun (depending on the accretion
rate) I will get,
f 1 Infinity
s% brunt_B(k) 1 0.0000000000000000D+00
s% gradT_sub_grada(k) 1
-2.7729163044273741D-01
s% gradT(k) 1 5.4141579620067472-225
s% grada_face(k) 1
2.7729163044273741D-01
STOP brunt
DATE: 2020-07-27
TIME: 13:25:01
But then if I run it again (same controls) I get normal/expected
evolution up to 10^4-10^5 Msun. I will keep testing and figure out whats
happening, but mostly this correction to my subroutine is producing good
final masses so thank you again.
Bill, I will talk to some people and look through some papers on the
post-Newtonian corrections to G and figure out why there is a dependence
on P. For the graphs shall I use the run_star_extras to output Grel
(log(Grel-1)) and plot this against mass?
Thanks again guys for your help
Nick
On 27/07/2020 08:13, Rob Farmer wrote:
> Hi,
> There is an issue with the code you have:
>
> do k = 1, s% nz !s% nz,1,-1
>
> TwoG = (2*standard_cgrav*s% m(k))/(s% r(k)*clight*clight)
> Tone = (s% P(k)/(s% rho(k)*clight*clight))
> Ttwo = (4*pi*s% P(k)*s% r(k)*s% r(k)*s% r(k))/(s%
> m(k)*clight*clight)
> Grel = standard_cgrav*((1+Tone+Ttwo)/(1-TwoG))
>
> end do
> s% cgrav(:) = Grel
>
>
> What you are doing here is computing Grel in each zone, but at the end
> you are only assigning the final value of Grel (from zone s%nz) to all
> zones in the star.
>
> If you want a different cgrav based on grel in each zone then you need
> to do:
>
> do k = 1, s% nz !s% nz,1,-1
>
> TwoG = (2*standard_cgrav*s% m(k))/(s% r(k)*clight*clight)
> Tone = (s% P(k)/(s% rho(k)*clight*clight))
> Ttwo = (4*pi*s% P(k)*s% r(k)*s% r(k)*s% r(k))/(s%
> m(k)*clight*clight)
> Grel = standard_cgrav*((1+Tone+Ttwo)/(1-TwoG))
> s% cgrav(k) = Grel
> end do
>
> To assign a different value in each zone.
>
>
> Also, with mesa if you do want to use array operations dont use (:)
> use (1:s%nz), i.e instead of:
>
> s%a(:) = s%b(:) * c
>
> do
>
> s%a(1:s%nz) = s%b(1:s%nz) * c
>
> As mesa's arrays may be bigger than the star currently is due to how
> we do memory management and thus array elements greater than s%nz may
> be filled with garbage from a previous step.
>
> Rob
>
>
>
> On Sun, 26 Jul 2020 at 19:47, Bill Paxton via Mesa-users
> <mesa-users at lists.mesastar.org <mailto:mesa-users at lists.mesastar.org>>
> wrote:
>
>
>> On Jul 26, 2020, at 10:01 AM, Nicholas Herrington
>> <Nicholas.Herrington at myport.ac.uk
>> <mailto:Nicholas.Herrington at myport.ac.uk>> wrote:
>>
>> Hi Bill
>>
>> I am trying to evolve a rapidly accreting supermassive primordial
>> star. The star starts at 10 Msun and accretes matter at 1 Msun /
>> yr until the star reaches ~450,000 Msun (no gr factors), 120,000
>> Msun (gr factors on) and ~ 220,000-320,000 Msun with my other
>> cgrav routine.
>
> Wow! ok, that’s certainly pushing mesa into new territory. so
> i can understand that you might require more accurate
> approximations of the GR effects.
>
>>
>> I am not familiar with how the corrections are produced. I’ve
>> been told they can be linearised in different ways, but from
>> previous work in the field with other codes such as Kepler and
>> Geneva, they implemented gr corrections to G in the same way I am
>> trying.
>
> i’m no expert, but it seems there are various options for
> approximating GR corrections - e.g., see
> https://www.researchgate.net/publication/318836552_Einstein's_postulate_as_a_correction_to_Newton's_law_of_gravity
> <https://eur04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.researchgate.net%2Fpublication%2F318836552_Einstein%27s_postulate_as_a_correction_to_Newton%27s_law_of_gravity&data=02%7C01%7Cr.j.farmer%40uva.nl%7C68785f4aa2394a822d3508d8318c060e%7Ca0f1cacd618c4403b94576fb3d6874e5%7C1%7C0%7C637313824758872278&sdata=rVW0Y%2BqOEYZerosxgSLr6HFo61Oobtelu80yQ4X4874%3D&reserved=0>)
>
> which mentions that for the Schwarzschild solution, the correction
> is what we have currently in mesa
> 1/sqrt(1 - 2*G*m/(r*c^2))
> for weak fields, this can be simplified to
> 1/(1 - G*m/(r*c^2))
> the author derives yet another form that includes the equivalent
> of a field mass,
> 1/(1 - G*m/(2*r*c^2))
>
> note that none of these correspond to what you are using.
>
>> Also we think the terms in Grel do not significantly effect G
>> until the stars surpass 10^4 Msun and that the between 1-3 x 10^5
>> Msun these corrections lead to strong collapse features before
>> the code stops with a min timestep error.
>>
>> Initially I saw gr factors in Mesa and wasn’t sure what they were
>> for but for the case of supermassive stars they actually produce
>> final masses in the 10^5 Msun range. I was hoping to try
>> implement the corrections from Kepler and Geneva within
>> other_cgrav to compare to gr factors in Mesa. What was the
>> intention for gr factors?
>
> the current GR factors are expected to be adequate for neutron
> star surface calcutions. the difference in your application is
> that you want to model deep interior of very massive stars,
> whereas we’ve only done surface conditions for an envelope with
> core removed.
>
>>
>> I will look further into where the gr corrections came from for
>> supermassive stars.
>
> great. and please help me understand the need for a more
> complicated form of the correction including P by making some
> plots vs mass of the values of log(correction - 1) for various
> alternatives.
>
> the code currently assumes that timesteps will be small enough
> that the GR factor for a particular mass coordinate will not
> change significantly during a single step. that allows us to
> compute the correction at the start of the step and hold it
> constant during the solver iterations. but if you are looking at
> problems where the GR factor at a mass location can change by a
> non-negligible amount during the step, it may be necessary to
> recalculate the factor at each solver iteration and use partial
> derivatives of the factor in the Jacobian matrix that guides the
> solver changes. that would mean adding a routine that returns
> partials of the factor wrt the variables it depends on (lnR for
> the simple form, plus lnd and lnT for the form including
> pressure). if we find this is necessary to enable the solver
> to find solutions, then we’ll need to add a hook for
> “other_cgrav_implicit” in addition to the current “other_cgrav”.
>
> once you provide numbers indictating showing whether or not we
> need pressure dependence, i’ll take on the task adding the new
> hook to mesa/star so you can do your experiments without needing
> to change the private code.
>
> thanks,
> bill
>
>
>>
>> Thanks
>> Nick
>>
>> On Sun, 26 Jul 2020 at 17:01, Bill Paxton <paxton at kitp.ucsb.edu
>> <mailto:paxton at kitp.ucsb.edu>> wrote:
>>
>> Hi Nick,
>>
>> As you know, mesa currently supports a simpler form of GR
>> corrections:
>>
>> !### use_gr_factors
>>
>> ! GR corrections. Currently just for pressure equation.
>> ! GR gravity factor = 1/sqrt(1-2Gm/(rc^2))
>>
>> use_gr_factors = .false.
>>
>> You have changed of this term of the correction to 1/(1 -
>> 2Gm/(rc^2)). You’ve approximated away the sqrt, but didn’t
>> drop the factor of 2 as I would have expected. Is that a
>> bug? Or do I have a bug in mesa?
>>
>> Also, I’ve been told that the other corrections you’ve
>> included are not significant for situations expected in mesa
>> (at least until we try to do neutron star interiors!). Could
>> you please provide numbers to show that in fact those extra
>> terms are important for a plausible application in mesa (such
>> as neutron star surface as opposed to interior) and that
>> limiting the correction to the term GR gravity factor
>> currently implemented in mesa would be significantly wrong.
>>
>> thanks,
>> Bill
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>> On Jul 26, 2020, at 8:26 AM, Nicholas Herrington via
>>> Mesa-users <mesa-users at lists.mesastar.org
>>> <mailto:mesa-users at lists.mesastar.org>> wrote:
>>>
>>> Hi
>>>
>>> I have transferred the corrections from the source code to
>>> the other_cgrav subroutine in the run_star_extras.f, as you
>>> suggested.
>>>
>>> When implementing the TOV corrections to G in the source
>>> code I added a few extra terms onto G = s% cgrav(k) to form
>>> G = s% cgrav(k)*Grel, where Grel is a function of the
>>> pressure, density, mass and radius. In my run_star extras.f
>>> I initially copied my Grel(P(k),rho(k),M(k),R(k)) into the
>>> subroutine other_cgrav, under s% cgrav(:) =
>>> standard_cgrav*Grel (where previously s% cgrav(:) =
>>> standard_cgrav). This did not work, so I presume I needed to
>>> iterate this throughout all the zones of the star. I have
>>> attached what I've got so far and from my tests and tweaks
>>> the results of this subroutine are close to what has been
>>> done before with this TOV implementation in other stellar
>>> evolution codes (final masses around 10^5 Msun).
>>>
>>> I am still new to FORTRAN and my understanding of the loop I
>>> have made so far is that after each step of the code the
>>> loop evaluates G from 1 up to the current number of zones of
>>> the star, appends this to cgrav(:) as a list or array, which
>>> is then fed into the stellar equations? However I am not
>>> sure this is the case as each time I run the same star with
>>> same controls inlist, the final mass of the star changes
>>> significantly. With my source code implementation the star
>>> would consistently stop at ~110,000 Msun (accreting at 1
>>> Msun/yr) when the timesteps drop significantly below the dt
>>> limit, which is expected for a supermassive star at this
>>> accretion rate. Additionally depending on the number of
>>> isotopes in the nuclear network the final mass changes
>>> significantly as well which wasn't previously the case with
>>> the source code changes. Other work with the TOV corrections
>>> to G in the KEPLER code show consistent final masses of
>>> ~200,000 Msun, which also leads me to believe I am doing
>>> something wrong in my subroutine.
>>>
>>> I would really appreciate any comments and feedback on
>>> whether the subroutine is computing Grel in the right way.
>>> Below is the TOV approximations to G I am looking to
>>> implement, <aepgpaknhdgjjhfh.png>
>>>
>>> Here is the subroutine I have created so far to try and
>>> compute this,<jflkhccofibpbple.png>
>>>
>>>
>>> ('real (dp) :: Mass, Rad, Pr, rho' were previously in the
>>> subroutine, but have been removed now please ignore them).
>>>
>>>
>>> I have attached the run_star_extras.f, inlist_1Mdot for the
>>> controls and images of Grel and the other_cgrav subroutine
>>> named 'TOV_cgrav'.
>>>
>>> Thanks
>>> Nick
>>>
>>> On 23/06/2020 16:06, Rob Farmer wrote:
>>>> Hi
>>>> Please keep all replies on mesa-users so others can learn
>>>> and/or help.
>>>>
>>>> You should never need to alter the f90 source files (except
>>>> your run_star_extras.f90). Mesa has the concept of "hooks"
>>>> which can change the internal physics without needing to
>>>> alter the mesa source code.
>>>>
>>>> In this case look at the routine in
>>>> star/other/other_cgrav.f90. If you follow the instructions
>>>> you can alter G however you want. This has the advantage
>>>> that it ensures that all occurrences of G have your
>>>> corrections applied.
>>>>
>>>> You can then call your routine again from your
>>>> run_star_extras.f90 in the data_for_extra_profile_columns
>>>> function to save the values to the profile file.
>>>>
>>>> Rob
>>>>
>>>>
>>>> On Tue, 23 Jun 2020 at 16:55, Nicholas Herrington
>>>> <Nicholas.Herrington at myport.ac.uk
>>>> <mailto:Nicholas.Herrington at myport.ac.uk>> wrote:
>>>>
>>>> Hi Rob
>>>>
>>>> Sorry for the very late reply, I haven't been able to
>>>> access my work machine in a while. After coming back to
>>>> it, I have not been able to replicate the 10^6 Msun run
>>>> in version 11701 (only 500-800x10^5). I spoke to others
>>>> in the field and they have said that without the GR
>>>> corrections these stars don't necessarily have to reach
>>>> 10^6 Msun.
>>>>
>>>> Since the 12778 doesn't have post-Newtonian corrections
>>>> to the G, I have added the corrections to G within the
>>>> subroutine associated to 'use_gr_factors' (within the
>>>> hydro_reconstruct.f90). Instead of 'G =
>>>> cgrav*gr_factors', I have added my own 'G =
>>>> cgrav*Grel'. How can I output Grel and G in my profiles
>>>> or history.data? Also is there a better way for me to
>>>> correct G rather than changing the f90 files?
>>>>
>>>> Thanks
>>>> Nick
>>>>
>>>> On Fri, 15 May 2020 at 18:35, Rob Farmer
>>>> <r.j.farmer at uva.nl <mailto:r.j.farmer at uva.nl>> wrote:
>>>>
>>>> Hi,
>>>>
>>>> >Is there anything along the lines of the TOV
>>>> approximation? Was use_gr_factors intended for SMS?
>>>>
>>>> No only the changes to s% grav are implemented.
>>>>
>>>> > which is why I wonder if the new version has this.
>>>>
>>>> It is still there in 12778. But lots of other
>>>> things will have changed between versions, which
>>>> may be stopping the star from getting to 10^6msun
>>>> (which is what i interpret your saying is the
>>>> problem?). Or maybe you've stumbled on a bug, it
>>>> would be helpful if you could make a small test
>>>> case that shows this problem.
>>>>
>>>> Rob
>>>>
>>>>
>>>>
>>>> On Fri, 15 May 2020 at 18:59, Nicholas Herrington
>>>> via Mesa-users <mesa-users at lists.mesastar.org
>>>> <mailto:mesa-users at lists.mesastar.org>> wrote:
>>>>
>>>> Hi,
>>>>
>>>> I am evolving some supermassive primordial
>>>> stars (SMS) with accretion rates up to 1
>>>> Msun/yr. In the newest version of MESA is there
>>>> general relativistic correction to G for SMS. I
>>>> have found the control use_gr_factors with
>>>> corrections to cgrav along the lines of
>>>> (1-2GM/rc^2)^-0.5. Is there anything along the
>>>> lines of the TOV approximation? Was
>>>> use_gr_factors intended for SMS?
>>>>
>>>> I was running a SMS with version r11701 and
>>>> with use_gr_factors off I get results exceeding
>>>> 10^6 Msun which is expected, but with version
>>>> r12778 and the same inlist (use_gr_factors =
>>>> .false.), the star reaches 10^5. This is close
>>>> to the final mass of the same star from other
>>>> stellar evolution codes with post-newtonian
>>>> corrections to G, which is why I wonder if the
>>>> new version has this.
>>>>
>>>> Thanks
>>>> Nick
>>>>
>>> <run_star_extras.f><Grel.png><inlist_1Mdot.txt><runstarGrel.png>_______________________________________________
>>> mesa-users at lists.mesastar.org
>>> <mailto:mesa-users at lists.mesastar.org>
>>> https://lists.mesastar.org/mailman/listinfo/mesa-users
>>> <https://eur04.safelinks.protection.outlook.com/?url=https%3A%2F%2Flists.mesastar.org%2Fmailman%2Flistinfo%2Fmesa-users&data=02%7C01%7Cr.j.farmer%40uva.nl%7C68785f4aa2394a822d3508d8318c060e%7Ca0f1cacd618c4403b94576fb3d6874e5%7C1%7C0%7C637313824758872278&sdata=UnZ3O0xU3JoRhxRX84OQPNqz5ADqas8dV3h93feDPGM%3D&reserved=0>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.mesastar.org/pipermail/mesa-users/attachments/20200727/0711bb65/attachment.htm>
More information about the Mesa-users
mailing list