[Mesa-users] Question about GR corrections for SMS

Rob Farmer r.j.farmer at uva.nl
Mon Jul 27 03:13:31 EDT 2020 

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> wrote:

>
> On Jul 26, 2020, at 10:01 AM, Nicholas Herrington <
> 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> 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> 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> 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> 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> 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
>> 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/bc734741/attachment.htm>

More information about the Mesa-users mailing list