Page 1 of 2
Rotational splittings in summary files
Posted: Tue May 01, 2018 6:32 am
by warrick
I thought it'd be neat if GYRE could compute rotational splittings directly in the summary output. Given a GYRE model, all the data is there and, though it isn't hard to do in post-processing, it can be done on the fly.
I took a stab at coding up the relevant components but I quickly got quite confused about what one would need to add. The key thing is adding the calculation to `src/mode/gyre_mode.fpp`. I've just tried adding this, which at least compiles:
Code: Select all
function splitting (this)
class(mode_t), intent(in) :: this
real(WP) :: splitting
integer :: k
real(WP) :: dbeta_dx(this%n_k)
real(WP) :: Omega_rot(this%n_k)
! Evaluate rotational splitting. Based on eqns 3.354-3.357 of
! [Aer2010]
if (this%l /= 0) then
do k = 1, this%n_k
dbeta_dx(k) = this%dbeta_dx(k)
associate ( &
ml => this%cx%ml, &
pt => this%gr%pt(k))
Omega_rot(k) = ml%coeff(I_OMEGA_ROT, pt)
end associate
end do
splitting = integrate(this%gr%pt%x, dbeta_dx*Omega_rot)
else
splitting = 0._WP
end if
! Finish
return
end function splitting
I'm not much good with object-oriented Fortran so I'm not sure if I'm getting `Omega_rot` in a sensible way. I've also done a little bit of basic algebra to convince myself that `dbeta_dx` is basically the rotation kernel, possibly up to some constants that I'll figure out.
Anyway, I'll now start having a look at how to make this callable in the summary output.
Re: Rotational splittings in summary files
Posted: Tue May 01, 2018 7:31 am
by warrick
OK, here's the diff to add the splitting to the summary.
Code: Select all
--- (old)/src/output/gyre_output.fpp 2017-09-23 21:41:31.000000000 +0000
+++ (new)/src/output/gyre_output.fpp 2018-05-01 12:20:29.459089049 +0000
@@ -180,6 +180,7 @@
$OUTPUT_MODES(r,W,W())
$OUTPUT_MODES(r,W_eps,W_eps())
$OUTPUT_MODES(r,beta,beta())
+ $OUTPUT_MODES(r,splitting,splitting())
$OUTPUT_MODES(r,x_ref,gr%pt(md(i_md)%k_ref)%x)
$OUTPUT_MODES(c,xi_r_ref,xi_r(md(i_md)%k_ref))
$OUTPUT_MODES(c,xi_h_ref,xi_h(md(i_md)%k_ref))
@@ -432,6 +433,10 @@
call wr%write('beta', md%beta())
+ case ('splitting')
+
+ call wr%write('splitting', md%splitting())
+
case ('x')
$if ($GFORTRAN_PR_49636)
Also, I missed a public declaration in the first block, so here's the full diff for completeness.
Code: Select all
--- (old)/src/mode/gyre_mode.fpp 2017-09-23 21:41:31.000000000 +0000
+++ (new)/src/mode/gyre_mode.fpp 2018-05-01 12:26:01.476054947 +0000
@@ -111,6 +111,7 @@
procedure, public :: W
procedure, public :: W_eps
procedure, public :: beta
+ procedure, public :: splitting
procedure, public :: omega_int
procedure, public :: eta
end type mode_t
@@ -2089,4 +2090,42 @@
end function eta
+ !****
+
+ function splitting (this)
+
+ class(mode_t), intent(in) :: this
+ real(WP) :: splitting
+
+ integer :: k
+ real(WP) :: dbeta_dx(this%n_k)
+ real(WP) :: Omega_rot(this%n_k)
+
+ ! Evaluate rotational splitting. Based on eqns 3.354-3.357 of
+ ! [Aer2010]
+
+ if (this%l /= 0) then
+
+ do k = 1, this%n_k
+ dbeta_dx(k) = this%dbeta_dx(k)
+
+ associate ( &
+ ml => this%cx%ml, &
+ pt => this%gr%pt(k))
+ Omega_rot(k) = ml%coeff(I_OMEGA_ROT, pt)
+ end associate
+
+ end do
+
+ splitting = integrate(this%gr%pt%x, dbeta_dx*Omega_rot)
+ else
+ splitting = 0._WP
+ end if
+
+ ! Finish
+
+ return
+
+ end function splitting
+
end module gyre_mode
So I'm now writing some number to the summary file, and it's called `splitting`, but I'll need to tweak it. Also, it's very slow to save the data to disk so I'll fix that too.
Re: Rotational splittings in summary files
Posted: Tue May 01, 2018 9:07 am
by rhtownsend
Hi Warrick --
My gut reaction to your post(s) was "Well, we already can write out beta, which gives the splitting for uniform rotation. And for differential rotation, a little post-processing with dbeta_dx can achieve the desired result. So what's the point?".
However, on thinking on it I agree that having the capability to write out splittings there and then, without the need for post-processing, would in fact be pretty useful. So, I'm all for this -- I'll add your changes in when I get a chance.
Many thanks,
Rich
Re: Rotational splittings in summary files
Posted: Tue May 01, 2018 9:35 am
by rhtownsend
OK, just committed the changes to the development version of GYRE (the splitting is labeled as domega_rot).
Thanks again for your input!
cheers,
Rich
Re: Rotational splittings in summary files
Posted: Tue May 01, 2018 10:40 am
by warrick
Hang on a second: I'm sure my formula needs normalisation. I just inserted Omega = 1.0 throughout a stellar model and I'm getting a range of results in the range 0.5-1.0 * 10^6 where I expect 1.0.
I'm happy to try my hand at opening this as a pull request, which will give me more of a chance to actually write something that works properly...
The stellar model (MESA/GYRE format) is attached. Just for reference, I'm playing with the MESA r10398 bundled version (5.1).
Re: Rotational splittings in summary files
Posted: Tue May 01, 2018 11:19 am
by warrick
OK, part of the problem is simply that I'm confusing myself by using a low-luminosity red giant, in which case all the mixed modes throw me off. I've just run a Sun-like model (attached) so at least the rotational splitting is roughly constant. Now just to figure out the units...
Re: Rotational splittings in summary files
Posted: Tue May 01, 2018 11:44 am
by warrick
OK, last post for today. It looks like the splittings are in units of 2*dynamical frequency. Should I try to use the function `freq_from_omega` to get the same units as the mode frequency? Like the mode frequency itself, the user could indicate this with `dfreq_rot` rather than `domega_rot`.
Re: Rotational splittings in summary files
Posted: Tue May 01, 2018 11:56 am
by rhtownsend
warrick wrote:OK, last post for today. It looks like the splittings are in units of 2*dynamical frequency. Should I try to use the function `freq_from_omega` to get the same units as the mode frequency? Like the mode frequency itself, the user could indicate this with `dfreq_rot` rather than `domega_rot`.
The splittings are in dimensionless frequency units, same as omega:
omega = sigma / sqrt(GM/R^3)
...where sigma is the angular eigenfrequency in the inertial frame. (The 'dynamical frequency' is not a well defined quantity, so I prefer to give an explicit definition like this).
You should indeed be able to use freq_from_omega to calculate dfreq_rot, but make sure the frame is set to 'INERTIAL'.
I've never dealt with a pull request, but would be happy to learn -- so, go ahead and give it a shot. Note that you should be working against the development version of GYRE, rather than 5.1; to set up the dev version:
Code: Select all
hg clone https://bitbucket.org/rhdtownsend/gyre
mkdir gyre/src/extern
cd gyre/src/extern
hg clone https://bitbucket.org/madstar/core
cheers,
Rich
Re: Rotational splittings in summary files
Posted: Wed May 02, 2018 8:24 am
by warrick
I've worked through the algebra and am reasonably confident that the rotational splitting needs to be divided by 4PI. This is corroborated by a model in which I inserted a rotation rate of 2PI rad/s. I'm readying the pull request with this change and a dimensioned frequency splitting but here's a PDF with my working in the meantime.
Re: Rotational splittings in summary files
Posted: Wed May 02, 2018 8:24 pm
by rhtownsend
warrick wrote:I've worked through the algebra and am reasonably confident that the rotational splitting needs to be divided by 4PI. This is corroborated by a model in which I inserted a rotation rate of 2PI rad/s. I'm readying the pull request with this change and a dimensioned frequency splitting but here's a PDF with my working in the meantime.
Thanks for spotting this. I think the problem crept in because when I evaluate dbeta_dx I divide by the mode inertia E (because that's essentially what forms the denominator of eqn. 3.357 in Aerts+2010). However, at some point in the recent past I altered my definition of mode amplitudes, and that entailed scaling E by a factor 4pi so that it now agrees with the E given in eqn. 3.139 (ibid.). Clearly, I didn't propagate this change into the evaluation of dbeta_dx and/or beta.
I'll figure how to go about reviewing your pull request and integrating it into the development branch. Something we should also think about: right now, we evaluate E (involving an integral over the full star) *each time* we call dbeta_dx. This is pretty expensive; we might want to think about caching E in the wave_t derived type.
cheers,
Rich