In other words, my old line was:

- Code:
`yvals = np.sign(dic['dW_dx']) * np.log10(np.abs(dic['dW_dx']))`

but now it looks more and less like:

- Code:
`yvals = np.sign(copy.deepcopy(dic['dW_dx'])) * np.log10(np.abs(copy.deepcopy(dic['dW_dx'])))`

and problem is solved. Now my integrations are unaffected by the former lines. Please consider the attached plot.

Sorry about bothering you with my mistakes. Now, at least I'm sure I can externally reproduce the net cumulative work.

Best regards

Ehsan.

Statistics: Posted by ehsan — Sat Oct 04, 2014 10:38 am

]]>

ehsan wrote:

Thanks Rich for your attempt.

But, it is not fair to compare your plot with mine (top panel) for now. There, I show " sign(dW_dx) . log(abs(dW_dx))". The only reason I do so is just proper visualization, and not be dominated by the huge factor of ~10^7. Our difference is not a Python problem. If you change your plotting convention, then you will certainly reproduce my top panel.

However, the issue is the net work, which I cannot reproduce and compare with the GYRE internal calculations.

Cheers

Ehsan.

Thanks Rich for your attempt.

But, it is not fair to compare your plot with mine (top panel) for now. There, I show " sign(dW_dx) . log(abs(dW_dx))". The only reason I do so is just proper visualization, and not be dominated by the huge factor of ~10^7. Our difference is not a Python problem. If you change your plotting convention, then you will certainly reproduce my top panel.

However, the issue is the net work, which I cannot reproduce and compare with the GYRE internal calculations.

Cheers

Ehsan.

OK, I've reproduced your plot -- so you're right, the problem isn't there. But for the eigenfunction you supply in the .zip file (nad-Dmix-00001.h5), I find I get identical values (=-963786.305756) for W from GYRE and from the sample Python code you pasted above. Here's the script I'm running:

- Code:
`#!/usr/bin/env python`

import gyre as gy

import numpy as np

d, r = gy.read_output('nad-Dmix-00001.h5')

nz = len(r['x'])

dx = r['x'][1 : ] - r['x'][ : -1]

integrand= (r['dW_dx'][1 : ] + r['dW_dx'][ : -1]) / 2.

W_cumul = np.zeros(nz-1)

W_net = 0.0

for i in range(nz-1):

W_cumul[i] = np.sum(integrand[0:i+1] * dx[0:i+1])

W_net += integrand[i] * dx[i]

print "GYRE gives : %f" % d['W']

print "Python gives : %f" % W_net

Can you try running this code to see what happens?

cheers,

Rich

Statistics: Posted by rhtownsend — Fri Oct 03, 2014 11:28 am

]]>

However, the issue is the net work, which I cannot reproduce and compare with the GYRE internal calculations.

Cheers

Ehsan.

Statistics: Posted by ehsan — Wed Oct 01, 2014 2:49 pm

]]>

Looks like a problem with Python. Attached is a plot of the eigenfunction file you sent me; this looks essentially identical to what I get when I run GYRE myself.

cheers,

Rich

dW_dx_rhdt.png

Statistics: Posted by rhtownsend — Wed Oct 01, 2014 1:01 pm

]]>

This is for the mode with n_pg = -38, the first mode calculated in the series.

The attached plot also shows the work derivative (top), and the cumulative work (bottom).

For this specific mode, the net work from the gyre output HDF5 file gives: dic['W'] = -963786.305756, but what I calculate (using the simple python block in the former message) is

W_net = -2.78840123095.

Moreover, all my W_net estimates are negative, whereas dic['W'] changes sign for some modes (which should be).

Thus, I am starting to speculate my routine, though I really do not touch the columns and attributes read from eigenfunction files using gyre.read_output().

Best

Ehsan.

Statistics: Posted by ehsan — Wed Oct 01, 2014 10:24 am

]]>

ehsan wrote:

Thanks Rich for looking into this.

Please find the inlist and a sample model.gyre as a zip file in the following ftp address; I cannot attach to this message

ftp://anonymous@ftp.ster.kuleuven.be/dist/ehsan/dW_dx.zip

I use GYRE v.3.0. Sorry, a bit old fashioned

Cheers

Ehsan.

Thanks Rich for looking into this.

Please find the inlist and a sample model.gyre as a zip file in the following ftp address; I cannot attach to this message

ftp://anonymous@ftp.ster.kuleuven.be/dist/ehsan/dW_dx.zip

I use GYRE v.3.0. Sorry, a bit old fashioned

Cheers

Ehsan.

Hi Ehsan --

I'm still unable to reproduce your problem, whether with 3.0 or the latest development version. Could you post a link to one of the problem output files (nad-NNNNN.h5)?

cheers,

Rich

Statistics: Posted by rhtownsend — Tue Sep 30, 2014 9:03 pm

]]>

ftp://anonymous@ftp.ster.kuleuven.be/dist/ehsan/dW_dx.zip

I use GYRE v.3.0. Sorry, a bit old fashioned

Cheers

Ehsan.

Statistics: Posted by ehsan — Tue Sep 30, 2014 12:21 pm

]]>

ehsan wrote:

Hi Rich,

Thanks for your reaction.

Indeed I use and import the "gyre" module to interact with the output HDF5 files.

That's also strange that my routine works for you, but not for me

How shall we compare stuff?

Ehsan.

Hi Rich,

Thanks for your reaction.

Indeed I use and import the "gyre" module to interact with the output HDF5 files.

That's also strange that my routine works for you, but not for me

How shall we compare stuff?

Ehsan.

To start with, could you post the inlist + model.

cheers,

R

Statistics: Posted by rhtownsend — Tue Sep 30, 2014 8:21 am

]]>

Indeed I use and import the "gyre" module to interact with the output HDF5 files.

That's also strange that my routine works for you, but not for me

How shall we compare stuff?

Ehsan.

Statistics: Posted by ehsan — Tue Sep 30, 2014 5:38 am

]]>

I can't see anything wring with your Python code -- it works for me when I apply it to one of the non-adiabatic test cases.

How are you reading the GYRE data into Python? Are you using the supplied python module, or something else?

cheers,

Rich

Statistics: Posted by rhtownsend — Mon Sep 29, 2014 4:09 pm

]]>

Looking at the two files /nad/gyre_mode.f90 and /nad/gyre_util.f90, I see how W is simply calculated:

- Code:
`W = integrate(this%x, this%dW_dx())`

and

- Code:
`function integrate_r_ (x, y) result (int_y)`

real(WP), intent(in) :: x(:)

real(WP), intent(in) :: y(:)

real(WP) :: int_y

integer :: n

! Integrate y(x) using trapezoidal quadrature

n = SIZE(x)

int_y = SUM(0.5_WP*(y(2:) + y(:n-1))*(x(2:) - x(:n-1)))

! Finish

return

end function integrate_r_

Then, I try to reproduce the same net work using a simple Python script, and just cannot get there!

- Code:
`dx = dic['x'][1 : ] - dic['x'][ : -1]`

integrand= (dic['dW_dx'][1 : ] + dic['dW_dx'][ : -1]) / 2.

W_cumul = np.zeros(nz-1)

W_net = 0.0

for i in range(nz-1):

W_cumul[i] = np.sum(integrand[0:i+1] * dx[0:i+1])

W_net += integrand[i] * dx[i]

As an example, for a specific mode that GYRE gives W=1351.65750141, I get W_net=-1.93474738483.

Am I making a mistake here?

dW_dx from gyre_nad is still a real array?

Do I need to provide any file? inlist? mode info? input model?

Thanks in advance.

Ehsan.

Statistics: Posted by ehsan — Mon Sep 29, 2014 3:53 am

]]>

This is a bug-fix release, primarily focused on addressing some issues with the MESA interface.

As usual, the source code can be downloaded from

https://bitbucket.org/rhdtownsend/gyre/downloads/

Rich

Statistics: Posted by rhtownsend — Tue Sep 02, 2014 10:35 am

]]>

Apologies for the delay in answering you, I've been trying to come up with a better answer than the one I give here: the problem is with the model, and not with GYRE.

Experience over the past few years has shown that small departures from hydrostatic equilibrium, and other equilibria, in the input stellar model can cause distortions to the mode eigenfunctions, which in turn lead to mode mis-classifications. Figuring out how to suppress these departures is an area of ongoing research; but as a possible workaround, you might want to try increasing the number of grid points MESA is using for the model, to ensure that it's properly resolving the stellar structure.

Best wishes,

Rich

Statistics: Posted by rhtownsend — Sun Jul 06, 2014 10:47 am

]]>

I am having some trouble with the eigenfunction calculations for a subgiant (MESA-) model using GYRE. During the calculation of l=9 modes I noticed some improper values for n_pg. I tried to fix this issue by playing around with the different grid options as proposed on the GYRE website 'Understanding grids', but in this case it did not help. Even worse, when keeping the model grid the same (with RESAMP_DISPERSION enabled) and just using a finer frequency scan grid (n_freq=1500 -> n_freq=3500) I got the same frequencies as a result, but for some of them different eigenfunctions.

Even for some frequencies where the frequency values are exactly the same for both cases of n_freq, and the n_pg values are also the same and seem to be correct, there are differences in the eigenfunctions depending on the parameter n_freq!

I then had a closer look at the differing eigenfunctions, and it turns out that for the same frequency, they show in principle the same oscillatory behavior, but the p-mode part is shifted somehow (see attached figure "xir_mode17.png"), so the oscillations are not centered around zero anymore (which actually would explain the problem with the n_pg values in some of the cases, when the shift gets so large that the whole function is shifted above/below zero). In the file "diff_xir_mode17.png" the difference of the two eigenfunctions is plotted and you can see, that in the inner part, the eigenfunctions are nearly identical, but the difference gets larger in the outer part of the model. When zooming into the g-mode region, close to the core, there is actually also a difference present, but it is very small (~1e-10, see "diff_xir_mode17_core.png").

I find it quite unsettling that this whole problem occurs only for some of the modes but not all of them: for 10 of the 17 modes that were calculated, the eigenfunctions are exactly identical for both settings of n_freq (as I would expect).

Do you know what might cause this behavior of the eigenfunctions? Or how to solve this problem? I would be grateful for any help. Attached you can find the MESA model file and GYRE inlist I used.

Best wishes,

Wiebke

Statistics: Posted by wiebke — Tue Jun 24, 2014 9:38 am

]]>

Sorry for taking so long to check back here. As a feature, don't worry about trying to make sense of stuff above the cutoff. I realise that the adiabatic calculations stop making complete sense out there, and I was more curious about how the codes were treating this situation differently. Like I said, the main reason I asked is because I was trying to fit a particular star that has some very high reported frequencies, near the cutoff. When MESA compared the observed and modelled frequencies, I think the GYRE super-cutoff results were better at steering MESA in the right direction.

The problem, it turned out, was that I wasn't including the atmosphere in the pulsation calculation. Including it increased the acoustic cutoff by enough that the fit to start behaving better again. That is, the seismically-derived model parameters are consistent with the spectroscopic parameters.

Anyway, thanks as ever for the work and support!

W

Statistics: Posted by warrick — Wed May 21, 2014 2:02 am

]]>