Hydrostatic Equilibrium in input model

[ehsan]

Dear Rich et al.
I am comparing adiabatic oscillation frequencies from GYRE and GraCo (Moya et al. 2008), and have received fairly good agreement between the two. I will come to that later. For now, I like to discuss how well the MESA input models to GYRE satisfy the hydrostatic equilibrium. To assess that, I do a simplistic run with MESA for a 5Msun model, and store the first model (very young MS star), along with a GYRE profile.
Then, I read in the GYRE file, and follow the prescription by Pieter Degroote
http://www.ster.kuleuven.be/~pieterd/ru ... xtras.html
I have attached my MESA inlist and run_star_extras.f.

I have attached the GYRE model (tar.gz file), and the MESA profile output (profile1.txt). The following plot shows that the hydrostatic equilibrium (HSE) is perfectly fulfilled in the MESA profile.

Hydr-Equil-profile.pdf

plot for MESA profile

(252.7 KiB) Downloaded 490 times

However, this is not the case with the GYRE model.

Hydr-Equil-GYRE.pdf

plot for GYRE model

(192.38 KiB) Downloaded 366 times

I also attach a short Python script that produces a plot of the departure from HSE inside the star from the core to the surface (black line) for the GYRE model. The left hand side (LHS) and right hand side (RHS) of the equation are also overplotted with blue and red, respectively. To run the script, unzip GYRE-test.tgz, and make sure the profile1.prof.GYRE is sitting beside the .py script.
As you notice, the model is not in a perfect HSE!

For now (if I have not made silly mistakes), I guess the difference between the two is that pressure, density, temperature etc in MESA profiles are cell-averaged, while in GYRE files, they are interpolated to cell edges. I wonder if this can violate the perfect HSE reached in the profiles. For now, I have no other explanations.

As you can see in run_star_extras, I call

Code: Select all

           call star_get_gyre_info(s% id, keep_surface_point, add_atmosphere, &
              G, M_star, R_star, L_star, r, w, L, p, rho, T, &
              N2, Gamma_1, nabla_ad, delta, nabla,  &
              kappa, kappa_rho, kappa_T, &
              epsilon, epsilon_rho, epsilon_T, omega, ierr)

to get p, rho, T, etc. and this is where I suspect the culprit lies. I know that all intensive quantities are consistently interpolated to cell edges, but the concern is that HSE is maintained with this anymore? On the other hand, if you use logP, logRho etc provided in the default profile_column.list, then the model is in perfect HSE!

I would be grateful to hear back your comments on how to improve this, and if other users also have similar result in their GYRE files.

Best regards
Ehsan

[pieterd]

Hi Ehsan,

There is not much you can do about this, my guess is that this is all due to numerics. Recall that a MESA model is only in hydrostatic equilibrium according to the very strict definition you use. If you change the definition, e.g. from Langrangian to Eulerian, you're not in HSE. If you change the definition of numerical derivative (left, right, centered...), you're not in HSE. So it should come as no surprise that redefining the locations of the variables has a huge effect. You could try to recalculate e.g. mass and pressure from density to force HSE in your model. Chances are that you will end up with a model with a slightly different mass or radius though. Also you will need to make the choice as to which definition of numerical derivative you use (left, right, centered...). Preferably choose the same one as the pulsation code you use -- but this is (I think) the beauty of GYRE: it doesn't use any numerical derivatives. I take the liberty of copying a reply Rich once gave me:

This is a tricky question to answer, and I'm not sure if I myself have the full truth of it. But let me make some suggestions...

Internally, GYRE doesn't know anything about pressures, densities, etc --- instead, it works with the dimensionless structure coefficients V, U, A* and c_1 (see Unno et al 1989 if you aren't familiar with the definitions). In hydrostatic equilibrium, certain differential relationships between these variables must hold. For instance,

dlnV/dlnr = (1 - 1/Gamma_1)V - A* + U - 1

dc_1/dlnr = 3 - U

Of course, there's still ambiguity in these equations -- it's not clear how we should calculate the logarithmic gradients of V and c_1. GYRE itself never explicitly calculates these gradients --- but in deriving the pulsation equations that GYRE implements, these identity *are* used.

Luckily, from my simulations, I see that this is issue is not so important. As long as your variable locations are consistent (i.e. mass, pressure... all located at the same position in a numerical cell), you're fine. If you want hard numbers on this, I'm sure I can throw something together.

Cheers,
Pieter

[ehsan]

Thanks Pieter, for your reply,

If you want hard numbers on this, I'm sure I can throw something together.

I really cannot resist this offer; Yes!

It would be interesting to verify if you as another experienced MESA and GYRE user can independently find (with hard numbers and simple plots) that:
- the HSE is different between a MESA profile and GYRE input model for the same snapshot during the evolution,
- the overall look of the HSE across the model (like my plot in the first post).
It would be a huge loss of effort and information if MESA produces models at such a perfect level of HSE, but the perfection is lost when migrating the structure into oscillation format. At least, this is my impression.

If GYRE is (almost) derivative-independent, then what if we keep perfectionism going, and retain intensive quantities in the cell center, and the rest at the cell edges? I actually can try this and see the effect on frequencies; though seems a brute idea.

Of course, there's still ambiguity in these equations -- it's not clear how we should calculate the logarithmic gradients of V and c_1. GYRE itself never explicitly calculates these gradients

Then, how is this implicit differentiation carried out? In the input GYRE inlist, we need to specify differentiation scheme.

Best regards
Ehsan.

[pieterd]

Hi Ehsan,

I just wrote a long reply but I got timed out and lost my text. I'm not gonna input it again, so I try to keep it short.

I computed frequencies for the raw MESA model (*.data file, not *.data.GYRE file) and a version that I interpolated myself (going from 1600->5000 points). This interpolation procedure destroys the condition of HSE completely (see left and right subplots in the bottom row of comparison.png, the green zone means HSE, the red zone means non-HSE). The reasoning is that if I get the same frequencies before and after interpolation, the code is not sensitive to the exact numerical definition of HSE. I expect huge differences between the three codes I used (LOSC, ADIPLS and GYRE), just because the model is bullocks (it is defined on an interlaced grid).

In the frequency comparison plots (top row), I compare the set of frequencies of radial order modes for a 5Msol main sequence model to those computed with GYRE on the interpolated model. The difference between the left and right plots show you have sensitive the codes are to the numerical derivatives. They are not at the 0.006% level (which includes the influence of mesh density).

Conclusion: HSE is a good diagnostic to check the quality of the stellar evolution calculations, but the numerical derivatives involved are not relevant for the pulsation codes.

If you wish to see other exercises, you know where I live!

Cheers,
Pieter