Work integral

[ehsan]

Dear Rich,
I'm looking at the instabilities of few modes, and tried to reproduce the net work, W, which is an output attribute with dW_dx and x which can be found in the output mode files.
Looking at the two files /nad/gyre_mode.f90 and /nad/gyre_util.f90, I see how W is simply calculated:

Code: Select all

    W = integrate(this%x, this%dW_dx())

and

Code: Select all

  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: Select all

      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.

[rhtownsend]

Hi Ehsan --

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

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

[rhtownsend]

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

[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/di ... /dW_dx.zip

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

Cheers
Ehsan.

[rhtownsend]

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/di ... /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

[ehsan]

Hi Rich,
To start, I updated the zip file, so, you may still use the same ftp link.
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.

[rhtownsend]

Hi Ehsan --

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

Work function for Ehsan's eigenfunction

dW_dx_rhdt.png (16.46 KiB) Viewed 4409 times

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

[rhtownsend]

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: Select all

#!/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