Polytropic model with excised core

[morganmacleod]

Dear Rich & Gyre-Users,

I've been simulating the gas dynamics of stellar envelopes in binary systems and would like to use gyre to compute the oscillation frequencies of my (artificial) model stars. I'd like to compare these frequencies to the modes excited in the simulation.

I have polytropic stars, with non-zero central mass and non-zero central radius (simulations are in spherical coordinates centered on the star). EOS is a gamma law: P=rho*e(\gamma-1). What is the best way to map these to a format that Gyre understands?

The approach I had in mind was to compute the relevant values in one of the typical formats, for example GSM or FGONG, like the poly_to_fgong.fpp example routine included in the code (leaving some entries zero since they aren't relevant for a polytropic model).

Is there a way or work-around to accommodate a model with non-zero central mass, radius? (Image attached for context).

Thanks for your advice and help.

Morgan

[rhtownsend]

Hi Morgan --

Thanks for your questions and your plot!

GYRE can readily handle pulsation calculations where the core is excised -- simply set x_i in the &grid namelist group to the fractional radial coordinate of the desired inner boundary. One important wrinkle is that you must also specify an appropriate inner mechanical boundary condition. Usually (i.e., in cases where the core is not excised) GYRE's inner mechanical boundary condition is just a regularity condition, which prevents solutions blowing up at the origin. However, when the core is excised you have to decide what physical constraint to place on the solutions at the inner boundary.

GYRE currently provides two options in this respect, specified by the inner_bound parameter in the &osc namelist group:

• 'ZERO_R' sets the radial displacement perturbation to zero (a no-penetration boundary condition),
• 'ZERO_H' sets the horizontal displacement perturbation to zero (a no-slip boundary condition).

You'll have to decide which one is appropriate for your stars. As an example of using the inner_bound parameter, have a look at the test/ad/poly/1.5-no-ctr test case.

On the subject of polytropes, GYRE has its own internal 'POLY' file format for representing composite polytropes -- a sequence of contiguous regions with fixed polytropic index, separated by discrete density jumps. Unfortunately, this format isn't properly documented in the current GYRE release. But as you've already figured out, GYRE does include some tools for manipulating POLY-format files.

By default, only poly_to_fgong is compiled during the build process; however, if you install the Mad SDK as your Software Development Kit, instead of the usual MESA SDK, then the build_poly executable will also be created. build_poly can be used to create files in the POLY format; see test/ad/poly/build_poly.py (a Python wrapper around build_poly) for a demonstration of how to use it.

Of course, if you want to avoid using the POLY format, then a polytropic model can also be represented in the usual FGONG format -- as you suggest, just set the undefined values (e.g., the opacity) to zero.

Hope this answers most of your questions -- please do post if you have more.

cheers,

Rich

[morganmacleod]

Dear Rich,

*Many* thanks for the detailed guidance! I especially appreciate the insight on inner BCs and links to examples in the test suite.

I will follow up on these and be in touch after I make some progress.

Morgan

[morganmacleod]

Dear Rich & Gyre Users,

Following up with a few details for the record. I've so far tested a gamma=1.5, n=1.5 simple polytrope (and have compared l=0 frequencies with those in test/poly).

I had numpy arrays: rad, rho, m, p that represent radius, density, enclosed mass, and pressure for a given model. Then the following python snippet produces an FGONG format file that GYRE can read. Posting here in case it's helpful to the next reader. (No guarantees, but in extremely limited testing this seems correct).

Thanks again for the help and very grateful to have this tool available!

Morgan



To generate the FGONG format arrays:

Code: Select all

# based on gyre's poly_to_fgong.fpp
#http://www.astro.up.pt/corot/ntools/docs/CoRoT_ESTA_Files.pdf
# assumes that gamma, poly_gamma, Ggrav, rad(nk), rho(nk),  m(nk), p(nk) are all defined. 

# Params 
IVERS = 300
ICONST = 15
IVAR = 40

# Variables 

nk = len(rad)

# Dimensionless structure vars
# 1/ gamma dlnP/dlnr - dln\rho/dlnr
#As = - G m rho/ r P (1/gamma - 1/gamma_polytrope)
As = - Ggrav*m*rho/(rad*p) * (1/gamma - 1/poly_gamma)


# Physical structure vars
M_r = m.copy()
M_r[m==0] = 1.e-38


# global 
glob = np.zeros(ICONST)
glob[0] = m[-1]
glob[1] = rad[-1]


# vars
var = np.zeros((IVAR,nk))

# weird indexing to correspond to fortran naming in PDF
var[1 -1,:] = rad
var[2 -1,:] = np.log(M_r/M_r[-1])
var[4 -1,:] = p
var[5 -1,:] = rho
var[10 -1,:] = gamma
var[15 -1,:] = As

# reverse order
var = var[:,nk::-1].copy()

Then to save this array (using the package fortranformat):

Code: Select all

import fortranformat as ff

filename = "mypolytrope.fgong"
f = open(filename,'w')

header_line = ff.FortranRecordWriter('(15A)')

for i in range(4):
    print>>f, header_line.write("comment")

index_line = ff.FortranRecordWriter('4I10')
print>>f, index_line.write([nk,ICONST,IVAR,IVERS]) 

line = ff.FortranRecordWriter('1P5E16.9')
print>>f,  line.write(glob) 

for k in range(nk):
    print>>f, line.write(var[:,k]) 
    
f.close()

[rhtownsend]

Hi Morgan --

Glad things are working for you now, and very many thanks for posting the Python script!

cheers,

Rich