FGONG output is incomplete

[earlbellinger]

MESA's FGONG output is missing variable 20 (local gravitational luminosity L_g) and variables 26-28 (partial derivatives of Gamma1 with respect to rho, P, and Y from the EOS).

FGONG spec: https://phys.au.dk/~jcd/solar_models/file-format.pdf (see page 3 for the relevant list of missing point variables) MESA code: https://github.com/MESAHub/mesa/blob/cbc645113d23c619ac44499c311f328309094c04/star/private/pulse_fgong.f90#L193 The functions store_point_data_atm, store_point_data_env, and store_point_data_ctr would need to be updated accordingly.

The EOS derivatives were previously difficult (as far as I know) but now with autodiff might be easier? (@adamjermyn what do you think? And if it's indeed straightforward and you have the time, could you perhaps please show me how to go about implementing something this?)

Also, these EOS derivatives aren't in the GYRE file format, but they would be useful (at least to me) if they were, or if it were possible to otherwise output them.

Also, MESA outputs Teff and standard_cgrav for global variables 14 and 15. I personally have no objection to this -- they are obviously very useful -- but I don't see them being listed as such in the spec. (Some sort of reference/link to the spec used might also be valuable.)

[fxt44]

helmeos has always been able to return these third derivatives of the helmholtz free energy; other eos's not so much. presumably skye can also return them.

[adamjermyn]

Derivatives with respect to rho are returned by all EOS’s.

Derivatives wrt P are returned by none of them but can be constructed from other return quantities with a bit of algebra. Happy to look into that.

Derivatives wrt Y are returned by some but not all. We could fill it in with finite differences if we had to…

-Adam

On Tue, Dec 7 2021 at 7:48 AM, fxt @.***> wrote:

helmeos has always been able to return these third derivatives of the helmholtz free energy; other eos's not so much. presumably skye can also return them.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/MESAHub/mesa/issues/352#issuecomment-988050519, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABPR5H7RFZYSUXJ2T3RPJO3UPYUENANCNFSM5JRLCSLA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[fxt44]

what is being held constant for P derivatives?

dP = ∂P/∂T dT + ∂P/∂(rho) d(rho)

for constant pressure dP = 0, so

d(rho)/dT = - ∂P/∂T / ∂P/∂(rho)

for constant temperature, dT = 0, so

dX/dP = ∂(rho)/∂P_T * dX/drho

for constant density, d(rho) = 0, so

dX/dP = ∂T/∂P|_rho * dX/dT

[warrickball]

Also, MESA outputs Teff and standard_cgrav for global variables 14 and 15. I personally have no objection to this -- they are obviously very useful -- but I don't see them being listed as such in the spec. (Some sort of reference/link to the spec used might also be valuable.)

The reference is buried in the docs but seems to have ended up somewhat displaced from the FGONG controls. I'll try to tidy this up a bit later.

https://github.com/MESAHub/mesa/blob/0bf360fd4834b71efd2e35be81d491e8ebfd5973/star/defaults/controls.defaults#L637-L641

I thought the addition of G in the header was a newer part of the FGONG spec (like the wider floating point format) but I can't find it documented. I'll have another look for it.

[warrickball]

I think this defines the updated "standard", including the wide format and Teff and G as header items 14 and 15. I'll update the docs now.

[warrickball]

Docs are somewhat improved in cb3d46cdb92d8bac4c1880439ce20596420a4398, though I don't have references for GR1D and SAIO formats. The latter might be Saio & Cox (1980).

But this is also all a tangent from the actual issue of implementing those partial derivatives.

[fxt44]

the gr1d file format is in tex and pdf at

https://github.com/evanoconnor/GR1D/tree/master/docs

if i'm understanding correctly what you seek. curiously i note these docs recommend using the mesasdk.

[warrickball]

I've just had a swing at this but could (embarassingly) use some EoS wizards' help on getting the partial derivatives right.

I noticed that the EoS subroutines return derivative information, so had a first go using those. Here are some extra profile columns I've implemented:

      subroutine data_for_extra_profile_columns(id, n, nz, names, vals, ierr)
         use eos_def, only: i_Gamma1, i_lnPgas
         use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results
         use eos_lib, only: eosDT_get, radiation_pressure

         integer, intent(in) :: id, n, nz
         character (len=maxlen_profile_column_name) :: names(n)
         real(dp) :: vals(nz,n)
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         integer :: k

         real(dp), dimension(num_eos_basic_results) :: &
              res, dres_dlnRho, dres_dlnT
         real(dp), allocatable :: dres_dxa(:,:)

         include 'formats'

         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return

         ! note: do NOT add the extra names to profile_columns.list
         ! the profile_columns.list is only for the built-in profile column options.
         ! it must not include the new column names you are adding here.

         names(1) = 'Gamma1'
         names(2) = 'Gamma1_Rho'
         names(3) = 'Gamma1_T'
         names(4) = 'Gamma1_Y'
         names(5) = 'Gamma1_P'

         allocate(dres_dxa(num_eos_basic_results, s% species))

         do k = 1, nz
            call eosDT_get(s% eos_handle, s% species, s% chem_id, s% net_iso, &
               s% xa(:,k), s% Rho(k), s% lnd(k)/ln10, s% T(k), s% lnT(k)/ln10, &
               res, dres_dlnRho, dres_dlnT, dres_dxa, ierr)

            vals(k,1) = res(i_Gamma1)
            vals(k,2) = dres_dlnRho(i_Gamma1)/res(i_Gamma1)
            vals(k,3) = dres_dlnT(i_Gamma1)/res(i_Gamma1)
            vals(k,4) = dres_dxa(i_Gamma1,1) + dres_dlnT(i_Gamma1)*dres_dxa(i_lnPgas,1)/dres_dlnT(i_lnPgas)
            vals(k,4) = -vals(k,4)/res(i_Gamma1) ! d/dY ~ -d/dX?
            vals(k,5) = dres_dlnT(i_Gamma1)/(dres_dlnT(i_lnPgas))/res(i_Gamma1)
         end do

         deallocate(dres_dxa)

      end subroutine data_for_extra_profile_columns

The partial with respect to P is the only one that I think is right (sort of, more on that in a moment) because @fxt44 spelled out the right partial derivative in a previous post. So I compared this with output that @earlbellinger provided from a separate, external program gamder that he uses to get these derivatives. This is a solar mass, initial (X,Z)=(0.7,0.02) model around 4.6 Gyr. This is sort of promising. gamder uses finite differences, whereas I think MESA is returning the derivative of the interpolating polynomial, so it looks a bit "ringy". So there's mileage here but a few things need to be ironed out.

1. Are the interpolator's derivatives (if that's what we're getting) good enough for research purposes or do we need to do finite differences too? I'm happy to try the finite differences now that I have a better idea of what the EoS calls look like.
2. Is dres_dxa(:,1) actually a derivative with respect to hydrogen abundance X? It looked like it when I followed the EoS call sequence down to OPAL/SCVH but I'm not 100% sure.
3. The other derivatives are probably wrong because they're computed at fixed (T,Y) or (ρ,T), not (P,Y) and (ρ,P). I'm hoping the EoS wizards can post how to do that faster than I can dust off thermodynamics I.

[fxt44]

which eos's are mainly involved, opal & freeeos?

vals(k,2) = dres_dlnRho(i_Gamma1)/res(i_Gamma1)

looks correct for ∂(Gamma_1)/∂(rho). how does it compare?

[warrickball]

which eos's are mainly involved, opal & freeeos?

It's all FreeEOS down to R≈0.09, after which there's a blend to Skye, which takes over below R≈0.05. I actually thought we'd moved the EOS blend to higher temperature so it wouldn't turn up in Sun-like models so I'll try to check why we're getting Skye here. This is part of the discussion I remember, where I thought we'd moved the logT blend to 7.3. But it's showing up around logT=7.15 in this model.

vals(k,2) = dres_dlnRho(i_Gamma1)/res(i_Gamma1)

looks correct for ∂(Gamma_1)/∂(rho). how does it compare?

I'm not sure it is. dres_dlnRho is presumably at constant T, rather than at constant P, which is what we need. It doesn't compare well to the other program: (Excuse the titles: I'm not sure if this is ∂Γ₁ or ∂lnΓ₁ but either way I'm pretty sure outputs are the same things.)

[fxt44]

ah, it wasn't immediately obvious to me quantities at constant pressure were desired. in this case one picks up an extra term, but the base approach is the same:

dP = ∂P/∂T dT + ∂P/∂(rho) d(rho) + ∂P/∂(mu) d(mu)

for constant pressure dP = 0 and constant composition

∂P/∂T dT + ∂P/∂(rho) d(rho) = 0

or

d(rho)/dT = - ∂P/∂T / ∂P/∂(rho)

then for some quantity f, which could be Gamma_1, at constant composition

df = ∂f/∂T dT + ∂f/∂(rho) d(rho)

dividing by dT

df/dT = ∂f/∂T + ∂f/∂(rho) d(rho)/dT

at constant pressure, substituting for dT/d(rho)

df/dT = ∂f/∂T - ∂f/∂(rho) * (∂P/∂T / ∂P/∂(rho))

similarily,

df/d(rho) = ∂f/∂(rho) - ∂f/∂(T) * (∂P/∂(rho) / ∂P/∂T)

try these expressions.

[warrickball]

That did it! Here's the state of play with (∂lnΓ₁/∂Y)_{P,ρ}:

I'll have a go at fixing the partial to be at constant (P,ρ) instead of (ρ,T) but note that the derivative goes to zero below about R~0.09, which coincides with the blend to Skye. So Skye is returning something different, and possibly just no derivative at all. I'll look into that once the partial looks closer to gamder for R>0.1.

[fxt44]

minor - both plots have titles of ∂/∂(ln(rho)). is one supposed to be ∂/∂(ln(T))? also, is it ∂(Gamma_1) or ∂(ln(Gamma_1)) being plotted. Y is the helium mass fraction?

[warrickball]

First pair are both ∂lnΓ₁/∂lnρ. The second is a zoom of the first to better show the behaviour away from the surface.

Final plot should be ∂lnΓ₁/∂Y, where Y is the helium mass fraction. But I think I've introduced a bug into the indexing of names and vals such that the columns aren't what I think they are. Let me check that...

[warrickball]

I assume the logic for ∂lnΓ₁/∂Y is the same for what you did for ∂lnΓ₁/∂lnρ above, in which case I had a sign error. I presume I'm actually getting ∂Γ₁/∂X from the EoS, where X is the hydrogen abundance, and that dX=-dY because Z is fixed and X+Y+Z=1. So now I've got

vals(k,4) = dres_dxa(i_Gamma1,1) - dres_dlnT(i_Gamma1)*dres_dxa(i_lnPgas,1)/dres_dlnT(i_lnPgas) ! ∂Γ₁/∂X
vals(k,4) = -vals(k,4)/res(i_Gamma1) ! ∂lnΓ₁/∂Y=-(∂Γ₁/∂X)/Γ₁

That looks okay in 0.1<R<0.9 but differs substantially near the center (because I don't think Skye is returning the same derivative) and near the surface (for reasons I'm not sure about). Plot titles should now be accurate. Again, right is vertical zoom of left, so the smaller differences in 0.1<R<0.9 are clearer. I think that's where I'll leave this for now.

[fxt44]

thanks. how does ∂lnΓ₁/∂lnT look?

[warrickball]

Here you go. Nothing to compare to.

Edit (again): dres_dxa is only not zero for lnE and lnPgas. In eosDT_get:

! only return 1st two d_dxa results (lnE and lnPgas) to star
d_dxa(1:num_eos_d_dxa_results,1:species) = d_dxa_eos(1:num_eos_d_dxa_results, 1:species)

I'll try taking this limit off and see what carnage ensues.

[warrickball]

If we don't zero out all those derivatives, something reasonable comes out for the derivative with respect to Y. Again, the agreement isn't fantastic but it's close enough that I believe I'm computing the right thing. So we're down to two questions. The first is as before:

1. Are the interpolator's derivatives (if that's what we're getting) good enough for research purposes or do we need to do finite differences too? I'm happy to try the finite differences now that I have a better idea of what the EoS calls look like.

I think we need @earlbellinger to answer that.

Though Skye doesn't seem to be returning composition derivatives, I'm sure it can be rigged to do so if we think this analytic route is the way forward. I'm going to give this a break while we decide on that.

PS: Apparently GitHub markdown now supports LaTeX snippets! $$\left(\frac{\partial\Gamma1}{\partial Y}\right){P,\rho}$$

[earlbellinger]

Marvelous! I will be happy to take a look at this soon. I'm on vacation this and next week though.

[adamjermyn]

Just chiming in to say Skye produces no composition derivatives because gfortran doesn't support the machinery we need for autodiff to enable that.

It is possible to add a fixed number of partials with respect to a a fixed number of species if desired. The species in question could be determined cell-by-cell and iteration-by-iteration if need be, but the number (and hence the overhead) would be fixed. The overhead is roughly "multiply Skye runtime by # of species / 2" (/2 comes from the fact that we only do third total order, meaning composition derivatives only go out to second order in (rho,T), and there are ~half as many (0,1,2)-order mixed partials in (rho,T) as there are (0,1,2,3)-order ones.)

This is probably not something I can take on, but I'd be happy to provide guidance on implementation...

Alternatively, Josiah gifted us machinery for doing finite differences to extract EOS composition partials, and that currently gets used to fill in Skye partials for the newton solver. Maybe that machinery could be reused here?

-Adam

On Sat, May 21, 2022 at 11:03 AM, Earl Patrick Bellinger < @.***> wrote:

Marvelous! I will be happy to take a look at this soon. I'm on vacation this and next week though.

— Reply to this email directly, view it on GitHub https://github.com/MESAHub/mesa/issues/352#issuecomment-1133648767, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABPR5H2ALSCF2TE7JVHIU7DVLD3NJANCNFSM5JRLCSLA . You are receiving this because you were mentioned.Message ID: @.***>

[warrickball]

I sunk a few more hours into this and think we're more or less there.

Here's the current code in my extra profile columns.

```f90 subroutine data_for_extra_profile_columns(id, n, nz, names, vals, ierr) use eos_def, only: i_Gamma1, i_lnPgas, i_chiRho, i_chiT use eos_def, only: num_eos_basic_results, num_eos_d_dxa_results use eos_lib, only: eosDT_get integer, intent(in) :: id, n, nz character (len=maxlen_profile_column_name) :: names(n) real(dp) :: vals(nz,n) integer, intent(out) :: ierr type (star_info), pointer :: s integer :: k real(dp), dimension(num_eos_basic_results) :: & res, res_loX, res_hiX, dres_dlnRho, dres_dlnT real(dp), allocatable :: dres_dxa(:,:) real(dp), allocatable :: xa(:) real(dp), parameter :: delta = 1d-3 include 'formats' ierr = 0 call star_ptr(id, s, ierr) if (ierr /= 0) return names(1) = 'Gamma1' names(2) = 'Gamma1_Rho' names(3) = 'Gamma1_P' names(4) = 'Gamma1_Y' allocate(dres_dxa(num_eos_basic_results, s% species)) allocate(xa(s% species)) do k = 1, nz xa = s% xa(:,k) ! evaluate EoS at low X xa(1) = s% xa(1,k) - delta xa(3) = s% xa(3,k) + delta call eosDT_get(s% eos_handle, s% species, s% chem_id, s% net_iso, & xa(:), s% Rho(k), s% lnd(k)/ln10, s% T(k), s% lnT(k)/ln10, & res_loX, dres_dlnRho, dres_dlnT, dres_dxa, ierr) ! evaluate EoS at high X xa(1) = s% xa(1,k) + delta xa(3) = s% xa(3,k) - delta call eosDT_get(s% eos_handle, s% species, s% chem_id, s% net_iso, & xa(:), s% Rho(k), s% lnd(k)/ln10, s% T(k), s% lnT(k)/ln10, & res_hiX, dres_dlnRho, dres_dlnT, dres_dxa, ierr) ! evaluate actual X call eosDT_get(s% eos_handle, s% species, s% chem_id, s% net_iso, & s% xa(:,k), s% Rho(k), s% lnd(k)/ln10, s% T(k), s% lnT(k)/ln10, & res, dres_dlnRho, dres_dlnT, dres_dxa, ierr) ! use dres_dxa(:,1) to store finite difference, ! which is evaluated and constant (ρ,T) dres_dxa(:,1) = (res_hiX-res_loX)/2/delta vals(k,1) = res(i_Gamma1) vals(k,2) = dres_dlnRho(i_Gamma1) - dres_dlnT(i_Gamma1)*res(i_chiRho)/res(i_chiT) vals(k,2) = vals(k,2)/res(i_Gamma1) vals(k,3) = dres_dlnT(i_Gamma1)/res(i_chiT)/res(i_Gamma1) vals(k,4) = dres_dxa(i_Gamma1,1) - dres_dlnT(i_Gamma1)*dres_dxa(i_lnPgas,1)/res(i_chiT) vals(k,4) = -vals(k,4)/res(i_Gamma1) ! d/dY ~ -d/dX end do deallocate(dres_dxa) deallocate(xa) end subroutine data_for_extra_profile_columns ```

If you don't want the details, I'm using the interpolated gradients $\chi_\rho$ and $\chi_T$ where possible and the analytic value from the interpolants for partial derivatives of $\Gamma_1$ w.r.t. thermodynamic variables. I use finite differences for the gradient w.r.t. X at constant $\rho$, $T$ and $Z$, assuming that species 1 is ¹H and species 3 is ⁴He. I then switched off FreeEOS in the inlist so that I actually use OPAL, for a fairer comparison with the gamder routine Earl was using. Here are plots of the differences in the three partials, one full view and one zoomed in for detail.

My hunch is that the "wiggliness" in the MESA results is because the MESA OPAL tables are oversampled relative to the raw OPAL and MESA uses a higher-order interpolation (cubic vs quadratic). As far as I can tell, the OPAL tables are provided with a spacing of about ~0.04-0.06 dex vs MESA's 0.02. It's not entirely clear how gamder interpolates. The "docstring" says

c..... eos(i) are obtained from a quadradic interpolation at
c      fixed T6 at three values of Rho; followed by quadratic
c      interpolation along T6. Results smoothed by mixing
c      overlapping quadratics.

I haven't pored over the F77 closely enough to figure out exactly what this "smoothing" means.

My conclusion is that I'm computing the correct thing, at least. I don't know offhand how to make the partials any smoother, or even if that should be MESA's responsibility. If this is already good enough, I can open a PR sometime in the next few weeks

Alternatively, Josiah gifted us machinery for doing finite differences to extract EOS composition partials, and that currently gets used to fill in Skye partials for the newton solver. Maybe that machinery could be reused here?

Where should I look for this in the source?

[fxt44]

it is the eos's responsibility to return all quantities needed to assess, at minimum, thermodynamic consistency. opal is incomplete in this regard. i guarantee that interpolating over multiple cells produces a thermodynamic inconsistent eos. smoothing will not make it any better, and probably worse.