An example for custom reference states from MESA#

This script will take a MESA stellar evolution profile and convert it into a format that can be read in as a custom reference state in Rayleigh. You will need the rayleigh_diagnostics.py, reference_tools.py, and mesa.py files. You will also need a suitable MESA profile file, such as profile_mesa.data.

[ ]:

import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as spi
import scipy.integrate as spint
import scipy.signal as spsig
import sys, os

sys.path.insert(0, os.path.abspath('../../'))

import post_processing.rayleigh_diagnostics as rd
import post_processing.reference_tools as rt
import mesa
%matplotlib inline

[ ]:

def interp(r, v):
    prad = p.rmid[::-1] * mesa.rsol
    #You can also use 10**p.logR[::-1] or p.radius[::-1] instead of rmid[::-1], but rmid is the most accurate choice
    return np.interp(r, prad, v[::-1])

Set the work_dir variable to the location of the Python files listed above and MESA profile you would like to use.

[ ]:

work_dir = './'
sys.path.append(work_dir)

[ ]:

p = mesa.profile('profile_mesa.data')

Choose a suitable number of radial grid points. They do not need to be regularly spaced. You should err on the side of high resolution since Rayleigh’s Chebyshev domains have very fine grid spacing at the top and bottom of the domain.

[ ]:

nr = 5000
r0 = 5.1e10 # in cm
r1 = 6.8e10 # in cm
radius = np.linspace(r0, r1, nr)

From the MESA model, Rayleigh will need the density, buoyancy function \(\rho g / C_P\), temperature, vicoscity, thermal diffusion, electrical resistivity (for magnetic cases), heating profile (for cases with \(Q \ne 0\)), entropy gradient (for cases with reference state advection). Note that MESA radial indicies start at the bottom, while Rayleigh radial indicies start at the top.

[ ]:

r_MESA = p.rmid*mesa.rsol
density = interp(radius, 10**p.logRho)
temperature = interp(radius, 10**p.logT)
grav = interp(radius, p.grav)
cp = interp(radius, p.cp)
buoy = density * grav / cp
nu = 1e14 * np.ones_like(radius)
kappa = 1e14 * np.ones_like(radius)
eta = 1e14 * np.ones_like(radius)
hprofile = np.zeros_like(radius)
dsdr = np.zeros_like(radius)

WARNING You should be very careful how you think about the entropy gradient when moving from MESA to Rayleigh due to the differing equations of state. Here we have chosen to simply set the reference state entropy gradient to zero and let the convection establish it’s own entropy gradient. This may or may not be satisfactory for your application.

ANOTHER WARNING You should be very careful with radiative luminosity and/or nuclear energy generation. There are a number of ways to compute the heating functions you need. For this example, we have chosen to simply compute the luminosity profile needed for flux balance if the convective transport matches the values from MESA.

[ ]:

q_rad = -np.gradient((p.luminosity - p.conv_L_div_L*p.luminosity)*mesa.solarlum, r_MESA)/(4.0*np.pi*r_MESA**2)
heatingp = interp(radius, q_rad)
luminosity = np.trapz(4.0*np.pi*radius**2*heatingp, radius)
hprofile = heatingp/luminosity

Plot the density from MESA and the newly interpolated density that will be fed into Rayleigh to make sure they are consistent. We also plot the heating profile as a sanity check.

[ ]:

plt.plot(radius, density,'or')
plt.plot(p.rmid[::-1] * mesa.rsol, 10**p.logRho[::-1], '-b')
plt.xlabel('Radius (cm)')
plt.ylabel(r'Desnity (g/cm$^3$)')
plt.xlim([r0,r1])
plt.ylim([np.min(density), np.max(density)])

[ ]:

plt.plot(radius, luminosity*hprofile, '-r')
plt.xlabel('Radius (cm)')
plt.ylabel(r'$Q$ (erg/(cm$^3$ s))')

Now create the data structure that will be written to a file that Rayleigh can read, and then load in the needed radial functions.

[ ]:

my_ref = rt.equation_coefficients(radius)

my_ref.set_function(density,'density')
my_ref.set_function(buoy,'buoy')
my_ref.set_constant(1.0, 'buoy_fact')
my_ref.set_function(nu,'nu')
my_ref.set_constant(1.0, 'visc_fact')
my_ref.set_function(temperature,'temperature')
my_ref.set_function(kappa,'kappa')
my_ref.set_constant(1.0, 'diff_fact')
my_ref.set_constant(1.0, 'p_fact')
my_ref.set_function(hprofile,'heating')
my_ref.set_constant(1.0, 'luminosity')
my_ref.set_function(eta,'eta')
my_ref.set_constant(1.0, 'resist_fact')
my_ref.set_function(dsdr,'ds_dr')

my_ref.set_constant(luminosity,'luminosity')
print(my_ref.constants)

[ ]:

file_write='cref_from_MESA.dat'
my_ref.write(file_write)

Now you can use this file to run a Rayleigh model. Once your Rayleigh model has run you can use the equation_coefficients file to check how your specified reference state looks when transfered into Rayleigh.

[ ]:

#Once you're run for one time step, set have_run = True
radius1 = radius
gravity = grav
try:
    cref = rd.PDE_Coefficients()
    have_run = True
except:
    have_run = False

if (have_run):

    fig, ax = plt.subplots(ncols=3,nrows=3, figsize=(9,3*3))
    # Density variables
    ax[0][0].plot(cref.radius,cref.density,'yo')
    ax[0][0].plot(radius1,density)
    ax[0][0].set_xlabel('Radius')
    ax[0][0].set_title('Density')

    ax[0][1].plot(cref.radius, cref.nu,'yo')
    ax[0][1].plot(radius1, nu)
    ax[0][1].set_xlabel('Radius')
    ax[0][1].set_title(r'$\nu$')

    ax[0][2].plot(cref.radius,cref.kappa,'yo')
    ax[0][2].plot(radius1,kappa)
    ax[0][2].set_xlabel('Radius')
    ax[0][2].set_title(r'$\kappa$')

    ax[1][1].plot(cref.radius,cref.temperature,'yo')
    ax[1][1].plot(radius1,temperature)
    ax[1][1].set_xlabel('Radius')
    ax[1][1].set_title('Temperature')

    ''''
    #Activate this if your case is magnetic
    ax[1][0].plot(cref.radius, cref.eta,'yo')
    ax[1][0].plot(radius1, eta)
    ax[1][0].set_xlabel('Radius')
    ax[1][0].set_title(r'$\eta$')
    '''

    ax[2][1].plot(cref.radius, cref.dsdr,'yo')
    ax[2][1].plot(radius1, dsdr)
    ax[2][1].set_xlabel('Radius')
    ax[2][1].set_title('Log entropy gradient')

    ax[1][2].plot(cref.radius, cref.functions[:,1]*cref.constants[1],'yo')
    ax[1][2].plot(radius1, gravity*density/cp)
    ax[1][2].set_xlabel('Radius')
    ax[1][2].set_title('Gravity')

    ax[2][0].plot(cref.radius, cref.heating*cref.rho*cref.T,'yo')
    ax[2][0].plot(radius1, hprofile*luminosity)
    ax[2][0].set_xlabel('Radius')
    ax[2][0].set_title('Heating')

    plt.tight_layout()
    plt.show()

[ ]:

Rayleigh documentation

An example for custom reference states from MESA

An example for custom reference states from MESA#