Rotational splittings in summary files

[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.

[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.

[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

[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

[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).

[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...

[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`.

[rhtownsend]

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

[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.

[rhtownsend]

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