[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
-2.7729163044273741D-01
s% gradT(k) 1    5.4141579620067472-225
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
>
>     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