Brunt B from eos partials #823
Conversation
|
At the moment, the eos partials from skye, PC, and ideal eos are 0. MESA gets around this using finite differences for these terms in the hydro solver, so i'm considering employing that approach here, while there is no alternative. At the moment, this just sets the brunt B term to 0 when using these three EOS'. If gfortan had good support for derived types for arrays, then this would become a non-issue for Skye as it could then return the eos partials using the current machinery and a dynamic array size set by the number of isotopes in the network. Still waiting for a day when this becomes a reality, but looking at the latest gfort 15.1, this compiler support issue still persists. |
|
What do you mean by “derived types for arrays”?
… On May 22, 2025, at 1:03 PM, Ebraheem Farag ***@***.***> wrote:
Debraheem left a comment (MESAHub/mesa#823)
At the moment, the eos partials from skye, PC, and ideal eos are 0. MESA gets around this using finite differences for these terms in the hydro solver, so i'm considering employing that approach here, while there is no alternative.
At the moment, this just sets the brunt B term to 0 when using these three EOS'.
If gfortan had good support for derived types for arrays, then this would become a non-issue for Skye as it could then return the eos partials using the current machinery and a dynamic array size set by the number of isotopes in the network. Still waiting for a day when this becomes a reality, but looking at the latest gfort 15.1, this compiler support issue still persists.
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because your review was requested.Message ID: ***@***.***>
--
Rich Townsend • Professor of Astronomy
Astronomy Department • University of Wisconsin-Madison
Phone: 608-262-1752 • E-mail: ***@***.***
|
|
Sorry i think it is called a parameterized derived type? I probably explained that incorrectly, and I'm a compiler newb. You might already be familiar with this. Instead of having fixed typing for example: In theory you could have an auto_diff type that can handle an arbitrary sized array of variables and the compiler would handle it at run time, something like:
In SKYE, Adam has and he sets all the eos derivatives to 0, but if you could do from https://ui.adsabs.harvard.edu/abs/2021ApJ...913...72J/abstract
As far as I'm aware, this is still an issue. But there could be a way to hack around it... |
|
Ah, got it. Yes, gfortran has been painfully slow to implement parameterized derived types. Been waiting for years.
… On May 22, 2025, at 1:33 PM, Ebraheem Farag ***@***.***> wrote:
Debraheem left a comment (MESAHub/mesa#823)
Sorry i think it is called a parameterized derived type? I probably explained that incorrectly, and I'm a compiler newb. You might already be familiar with this.
Instead of having fixed typing for example:
auto_diff_real_star_order1 (which is fixed to 33 star variables) is set with a fixed size.
type :: auto_diff_real_star_order1
real(dp) :: val
real(dp) :: d1Array(33)
end type auto_diff_real_star_order1
In theory you could have an auto_diff type that can handle an arbitrary sized array of variables and the compiler would handle it at run time, something like:
type(auto_diff_real_star_order1(nVar))
In SKYE, Adam has
auto_diff_real_2var_order3 (fixed to 2 variables for rho-T, what SKYE uses)
and he sets all the eos derivatives to 0, but if you could do auto_diff_real_order3(2+num_isos), then SKYE would compute these partials for your for free-ish?
from https://ui.adsabs.harvard.edu/abs/2021ApJ...913...72J/abstract
Screen.Shot.2025-05-07.at.10.22.53.PM.png (view on web)
As far as I'm aware, this is still an issue. But there could be a way to hack around it...
—
Reply to this email directly, view it on GitHub, or unsubscribe.
You are receiving this because your review was requested.Message ID: ***@***.***>
--
Rich Townsend • Professor of Astronomy
Astronomy Department • University of Wisconsin-Madison
Phone: 608-262-1752 • E-mail: ***@***.***
|
|
IIRC this GCC bug report is the relevant one. |
This pr adds an option to directly compute the Brunt B term from eos partials.
The “B” term appears in the Brunt-Vaisala frequency and the Ledoux criterion (see section 3.3 of Paxton+2013, aka MESA II). Currently, we evaluate B on cell faces using finite difference directional derivative between the compositions of neighboring cells (equations 7 & 8 of MESA II).
This pr attempts to compute the "B" term directly using its formal definition following equation 6.
B = −1/χT ∑i=1N−1 (∂ ln P / ∂ ln Xi)ρ,T,Xj≠i · (d ln Xi / d ln P).
In practice, I actually sum over all the isotopes, and don't enforce the (\sum_{i=1}^N X_i = 1) condition, but that might be necessary if this is eventually turned into an implicit formulation.
I am still left wondering if it is worth computing this term using a compensator or quad precision (optionally?), as this term involves a summation over all terms that appear in the network, and hence is prone to loss of precision from consecutive subtraction/addition of large and small numbers.
Tagging asteroseismology folks who probablt know better than me for a review... I'm curious what if any effect this will have...
(still missing a changelog entry, but we will get there if this even works.)