# Custom reference state for the dimensional anelastic MHD formulation in a convective spherical shell based on polytropes#

## Problem Setup#

We define the reference state of the convective region based on a polytrope that has a form given by

$$\rho=\rho_0 z^n$$,

$$P=P_0 z^{n+1}$$,

and

$$T=T_0 z$$,

where $$\rho$$, $$P$$, and $$T$$ are the density, pressure, and temperature respectively, and where the polytropic variable $$z$$ is an as-of-yet undetermined function of radius. We further assume that these quantities are related by the ideal gas law, with

$$P=R\rho T$$,

where $$R$$ is the gas constant. It is related to the specific heats at constant pressure ($$c_p$$) and contant volume ($$c_v$$) through the relation

$$R=c_p-c_v=c_p(1-\frac{1}{\gamma})$$,

with

$$\gamma\equiv \frac{c_p}{c_v}$$.

For a monotomic ideal gas, we have that

$$\gamma\equiv \frac{5}{3}$$.

Finally, we require that the polytrope satisfies hydrostatic balance. Namely,

$$\rho\frac{GM}{r^2}=-\frac{\partial P}{\partial r}$$,

where $$G$$ is the gravitational constant, and $$M$$ is the mass of the star.

## Polytropic Solution#

Substituting our relations for $$P$$, $$\rho$$, and $$T$$ into the equation of hydrostatic balance, we arive at

$$\frac{\partial z}{\partial r}= - \frac{GM}{(n+1)RT_0 r^2} = - \frac{2GM}{5(n+1)c_pT_0 r^2}$$ .

This motivates us to seek a form for $$z$$ of

$$z = a +\frac{b}{r}$$,

and immediately, we see that b must be given by

$$b = \frac{2GM}{5(n+1)c_pT_0}$$.

Note that while $$T_0$$ remains undetermined, we can now compute $$\partial T/\partial r$$. In order to determine $$a$$, we need one more constraint. In our case, we will specify the number of density scaleheights, $$N_\rho$$, across the convection zone. We denote the top of the convection zone by a subscript $$t$$ and the base of the convection zone by a subscript $$b$$. We then have $$N_\rho = \frac{ \rho_b }{ \rho_t } = \frac{z_b^n}{z_t^n}$$,

or equivalently, using $$C$$ to denote the exponential factor, and $$\beta \equiv \frac{r_b}{r_t}$$, we have

$$C\equiv e^{\frac{N_\rho}{n}} = \frac{a+b/r_b}{a+b/r_t}$$.

Rearranging, we find our expression for $$a$$

$$a = \frac{fb}{r_b}$$,

where

$$f \equiv \frac{\beta C-1}{1-C}$$.

This yields our expression for $$z$$ in terms of $$T_0$$

$$z = b (\frac{1}{r} +\frac{f}{r_b} )= \frac{2GM}{5(n+1)c_pT_0}\left(\frac{1}{r} +\frac{f}{r_b} \right)$$.

Factors of $$T_0$$ cancel out, when calculating the temperature, leaving us with a complete description of its functional form. We are free to choose any value of $$T_0$$ as a result; we use $$T_b$$, the temperature at the base of the convection zone. We have that

$$T_0 = T_b = \frac{2GM}{5(n+1)c_p}\left(\frac{1+f}{r_b}\right)$$,

completing our description of z. Values for $$\rho_0$$ and $$P_0$$ can now similarly be computed by enforcing the value of $$\rho$$ and $$P$$ as a particular point. As with $$T$$, we choose the base of the convection zone in the code that follows.

[ ]:
#######################
import numpy
import matplotlib.pyplot as plt
from matplotlib.pyplot import plot, draw, show

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

import post_processing.reference_tools as rt # You will need the refernce_tools.py to run this notebook
import post_processing.rayleigh_diagnostics as rdiag
[ ]:
# Grid Parameters
# Here, we use solar values as an example

nr    = 512         # Number of radial points
ri    = 5e10        # Inner boundary of radial domain
ro    = 6.586e10    # Outer boundary of radial domain
rcz   = 5e10        # base of the CZ

#apsect ratio
beta=ri/ro

#Polytrope Parameters
ncz   = 1.5         # polytropic index of convection zone
nrho  = 3.          # number of density scaleheights across convection zone (not full domain)
mass  = 1.98891e33  # Mass of the star
G     = 6.67e-8     # gravitational constant
rhoi  = 0.18053428  # density at the base of convection zone
cp    = 3.5e8        # specific heat at constant pressure
gamma = 5.0/3.0     # Ratio of specific heats for ideal gas

lsun = 3.846e33     # solar luminosity for scaling the volumetric heating
[ ]:
radius = numpy.linspace(ri,ro,nr)           #Radial domain of the CZ

[ ]:
#Compute a  CZ polytrope (see reference_tools)

gas_constant = cp*(1.0-1.0/gamma)  # R

temperature= poly1.temperature # temperature T
density = poly1.density  # density rho
pressure = poly1.pressure # pressure P, this won't matter -- set it equal to something
dsdr = poly1.entropy_gradient # dS/dr
entropy= poly1.entropy # entropy S, this won't matter -- set it equal to something
dpdr = poly1.pressure_gradient # dP/dr
gravity= mass*G/radius**2  # gravity g
[ ]:
# Calculation of the first and second derivative of lnrho and the first derivative of lnT.
# If we do not calculate those here, then they are calculated within Rayleigh!

dlnrho = d_density_dr/density
dlnt = dtdr/temperature
[ ]:
# Plots of  the dimensional reference state profiles, i.e. the density profile, the temperature profile, etc.

fig, ax = plt.subplots(nrows =3,ncols = 3, figsize=(15,10) )
ax[0][0].set_ylabel('Density')

ax[0][1].set_ylabel('Entropy')

ax[0][2].set_ylabel('Temperature')

ax[1][1].set_ylabel('Pressure')

ax[1][2].set_ylabel('Gravity')

ax[2][0].set_ylabel('dlnrho')

ax[2][1].set_ylabel('d2lnrho')

ax[2][2].set_ylabel('dlnt')

plt.tight_layout()

plt.show()
[ ]:
# Calculation of the heating profile based on the pressure, if we want to have a heating function in our setup.

hprofile = numpy.zeros(nr,dtype='float64')
hprofile[:] = pressure[:]

#################################################################
# We normalize the heating function so that it integrates to 1.

hprofile = hprofile/hint

[ ]:
[ ]:
# Here we define all the functions and constants that will be written in our data file and
# read by Rayleigh if we choose the custom reference state (=4)

unity = numpy.ones(nr,dtype='float64')
buoy =density*gravity/cp # calculation of the buoyancy term used in the momentum equation

my_ref.set_function(density,1) # density rho
my_ref.set_function(buoy,2)    # buoyancy term
my_ref.set_function(unity,3)   # nu(r) -- can be overwritten via nu_type in Rayleigh
my_ref.set_function(temperature,4) # temperature T
my_ref.set_function(unity,5)  # kappa(r) -- works like nu
my_ref.set_function(hprofile,6)  # heating -- this is normalized to one

my_ref.set_function(dlnrho,8)  # dlnrho/dr
my_ref.set_function(d2lnrho,9) # d^2lnrho/dr^2
my_ref.set_function(dlnt,10)   # dlnT/dr
my_ref.set_function(unity,7)  # eta -- works like nu and kappa --  magnetic diffusivity
my_ref.set_function(dsdr,14)   # dS/dr

# The constants can all be set/overridden in the input file
# NOTE that they default to ZERO, but we want
# most of them to be UNITY.
# These aren't very useful in the dimensional anelastic formulation but are VERY useful
# for the non-dimensional version of the anelastic custom reference state!

my_ref.set_constant(1.0,1)  # multiplies the Coriolis force--Should be 2 Omega (complication here)
my_ref.set_constant(1.0,2)  # multiplies buoyancy
my_ref.set_constant(1.0,3)  # multiplies pressure gradient
my_ref.set_constant(0.0,4)  # multiplies lorentz force
my_ref.set_constant(1.0,5)  # multiplies viscosity
my_ref.set_constant(1.0,6)  # multiplies entropy diffusion (kappa)
my_ref.set_constant(0.0,7)  # multiplies eta in induction equation
my_ref.set_constant(1.0,8)  # multiplies viscous heating
my_ref.set_constant(1.0,9)  # multiplies ohmic heating
my_ref.set_constant(lsun,10) # multiplies the heating (if normalized to 1, this should be the luminosity)
my_ref.write('dimensional.dat')  # Here we write our data file to be used to run our simulation with Rayleigh!
print(my_ref.fset)
print(my_ref.cset)
[ ]:
# Once you've run for one time step, set have_run = True !!

# Here we check if the output reference state is the same as the one we used as an input (sanity check)

# NOTE: We need the output file "equation_coefficients" to run this, as well as the PDE_Coefficients
# from rayleigh_diagnostics.py

try:
cref = rdiag.PDE_Coefficients() # This will give us the output reference state
have_run = True
except:
have_run = False

if (have_run):

fig, ax = plt.subplots(ncols=3,nrows=4, figsize=(16,4*4))
# density, derivatives of lnrho
ax[0][0].set_title('Density')

# temperature and derivative of lnT
ax[1][0].set_title('Temperature')

# Note that you must build the buoyancy term from the functions/constants