Custom reference state for non-dimensional anelastic MHD formulation in a convective spherical shell
Custom reference state for non-dimensional anelastic MHD formulation in a convective spherical shell#
We non-dimensionalize the dimensional anelastic MHD equations (see Rayleigh documentation for the full set of equations, and the definitions of the variables and parameters used here) using the width of the convection zone \([l]=L=r_o-r_{cz}\) as the lengthscale (where \(r_o\) is the outer radius and \(r_{cz}\) is the radius at the bottom of the convection zone (CZ)), the viscous timescale \([t]=L^2/\nu\), and the velocity scale \([u]=\nu/L\). For the density and temperature scales, we choose \([\rho]=\tilde{\rho}\), and \([T]=\tilde{T}\), respectively, where the tildes correspond to the volume average of the respective quantity such that \(\tilde{q}=\dfrac{1}{V}\int_V q dV\), where \(V\) is the volume of the convective region of the spherical shell (and \(q\rightarrow \rho, T\)). For the magnetic field, we use \([B]=\sqrt{\tilde{\rho}\mu\eta\Omega_o}\) and finally for the entropy \(S\) we use a scaling related to the thermal energy flux \(\tilde{F}\) such that \([S]=\dfrac{L\tilde{F}}{\tilde{\rho}\tilde{T}\kappa}\). Then, we can write the non-dimensional anelastic MHD equations as:
\begin{equation} \dfrac{\partial{\vec{u}}}{\partial{t}} +\vec{u}\cdot\nabla\vec{u}+\dfrac{1}{\rm{E}}2\hat{z}\times\vec{u}=\dfrac{\rm{Ra}}{\rm{Pr}}\dfrac{g(r)}{\tilde{g}}S\hat{r}-\nabla(p/\bar{\rho})+\dfrac{1}{\bar{\rho}}\nabla\cdot \vec{D}+ \dfrac{1}{4\pi\bar{\rho}}\dfrac{1}{\rm{P_m E}}(\nabla\times\vec{B})\times\vec{B}. \end{equation}
\begin{equation} \nabla\cdot(\bar{\rho}\vec{u})=0, \end{equation} \begin{equation} \nabla\cdot\vec{B}=0. \end{equation}
\begin{equation} \bar{\rho}\bar{T}\left( \dfrac{dS}{dt}+\vec{u}\cdot\nabla S \right )=\dfrac{1}{\rm{Pr}}\nabla\cdot(\bar{\rho}\bar{T}\nabla S)+\dfrac{1}{\rm{Pr}}Q_{nd}+\dfrac{\rm{Di Pr}}{\rm{Ra}}\Phi+\dfrac{1}{4\pi}\dfrac{\rm{Di Pr}}{\rm{E P_m^2 Ra}}|\nabla\times\vec{B}|^2, \end{equation} where \(Q_{nd}=\dfrac{L}{\tilde{F}}Q\), and \begin{equation} \dfrac{\partial{\vec{B}}}{\partial{t}}-\nabla\times(\vec{u}\times\vec{B})=-\nabla\times\left(\dfrac{1}{\rm{P_m}}\nabla\times\vec{B}\right)\left\{=\dfrac{1}{\rm{P_m}}\nabla^2\vec{B}\right\}, \end{equation} where we have assumed that \(\kappa\), \(\nu\) and \(\eta\) are constants. We now have five non-dimensional numbers: the flux Rayleigh number Ra, the Prandtl number Pr, the magnetic Prandtl number \(\rm{P_m}\), the Ekman number E and the dissipation number Di, which are defined respectively as: \({\rm{Ra}}=\dfrac{\tilde{g}\tilde{F}L^4}{c_p\tilde{\rho}\tilde{T}\kappa^2\nu},\quad {\rm{Pr}}=\dfrac{\nu}{\kappa}, \quad{\rm{P_m}}=\dfrac{\nu}{\eta},\quad {\rm{E}}=\dfrac{\nu}{\Omega_oL^2},\quad\text{and}\quad {\rm{Di}}=\dfrac{\tilde{g}L}{c_p\tilde{T}}\). Then the functions \(f\) are:
\(f_1(r)\rightarrow \bar{\rho}\),
\(f_2(r)\rightarrow \bar{\rho}g(r)/\tilde{g}\),
\(f_3(r)\rightarrow 1\),
\(f_4(r)\rightarrow \bar{T}\),
\(f_5(r)\rightarrow 1\),
\(f_6(r)\rightarrow Q_{nd}(r)\),
\(f_7(r)\rightarrow 1\),
and the constants are
\(c_1\rightarrow 2/{\rm{E}}\),
\(c_2\rightarrow {\rm{Ra/Pr}}\),
\(c_3\rightarrow 1\),
\(c_4\rightarrow \dfrac{1}{4\pi}\dfrac{1}{\rm{P_m E}}\),
\(c_5\rightarrow 1\),
\(c_6\rightarrow 1/{\rm{Pr}}\),
\(c_7\rightarrow 1/{\rm{P_m}}\),
\(c_8\rightarrow \rm{Di Pr/Ra}\),
\(c_9\rightarrow\dfrac{1}{4\pi}\dfrac{\rm{Di Pr}}{\rm{E P_m^2 Ra}}\), and
\(c_{10}\rightarrow \dfrac{1}{\rm{Pr}}\).
The non-dimensional reference density and temperature profiles are equal to their dimensional profile divided by its volume average in the CZ, i.e. \(\bar{\rho}=\bar{\rho}_{dim}/\tilde{\rho}\), and \(\bar{T}=\bar{T}_{dim}/\tilde{T}\).
Note that the reference state is based on a polytrope, similarly to the dimensional case. For more info on that, check the dimensional anelastic notebook version!
[ ]:
#######################
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
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
# shell depth according to non-dimensionalization
d=1.e0
## In the main_input file, we should use the following radial boundaries or else set the right aspect ratio and shell depth
ri_nd=beta*d/(1-beta) # non-dimensional inner boundary
ro_nd=d/(1-beta) # non-dimensional outer boundary
print(ri_nd,ro_nd)
# Polytrope Parameters
# In this example, we use solar values
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
[ ]:
# Define the Radial Grid
radius = numpy.linspace(ri,ro,nr) # Full domain radial grid for reference state
[ ]:
#Compute CZ polytrope -- see reference_tools.py script
poly1 = rt.gen_poly(radius,ncz,nrho,mass,rhoi,G,cp,rcz)
gas_constant = cp*(1.0-1.0/gamma) # R
#Dimensional profiles for reference state functions
temperature1 = poly1.temperature # T
density1 = poly1.density # rho
pressure1 = poly1.pressure # P
dsdr1 = poly1.entropy_gradient # dS/dr
s1= poly1.entropy # S
dpdr1 = poly1.pressure_gradient # dP/dr
gravity1= mass*G/radius**2 # g
# We use the volume average of rho,T and g in the CZ to non-dimensionalize our reference state functions
def vol_av(f,radius):
int1=radius**2.
int2=f*radius**2.
vol = numpy.trapz(int1,x=radius)
f_av1 = numpy.trapz(int2,x=radius)
f_av=f_av1/vol
return f_av
temperature_av=vol_av(temperature1,radius) # CZ volume average of T
density_av=vol_av(density1,radius) # CZ volume average of rho
gravity_av=vol_av(gravity1,radius) # CZ volume average of g
print(temperature_av,density_av,gravity_av)
# Here we define the non-dimensional T, rho, g, S, r, dS/dr and P
# The non-dimensional temperature: T=T_dimensional/T_av
temperature=temperature1/temperature_av
# The non-dimensional density: rho=rho_dimensional/rho_av
density=density1/density_av
# The non-dimensional gravity: g=g_dimensional/g_av
gravity=gravity1/gravity_av
# Here we define our non-dimensional radius
radius1=radius/(ro-ri) # OR: radius1=numpy.linspace(ri_nd,ro_nd,nr)
#The pressure profile won't matter, we set it equal to rho*T as a reference
pressure=temperature*density
# For a purely convective region, we use dS/dr=0
dsdr=dsdr1
# Entropy won't matter, set it to something-- here I use the non-dimensional S i.e. S_nondim=[S_dim*(L^3*g_av)]/Ra*cp*kappa*nu]
entropy=numpy.ones(len(radius))*s1*(1.586e10**3.*gravity_av)/(13303.43109666924*cp*(8e12)**2)
[ ]:
# Here we calculate the derivatives of lnrho and lnT based on the non-dimensional radius1, since
# we want their non-dimensional profiles
d_density_dr = numpy.gradient(density,radius1, edge_order=2)
dlnrho = d_density_dr/density
d2lnrho = numpy.gradient(dlnrho,radius1, edge_order=2)
dtdr = numpy.gradient(temperature,radius1, edge_order =2)
dlnt = dtdr/temperature
### Here, we plot our non-dimensional profiles for our reference state functions, i.e. for rho, T, etc.
fig, ax = plt.subplots(nrows =3,ncols = 3, figsize=(15,10) )
ax[0][0].plot(radius1,density,'r')
ax[0][0].set_xlabel('Radius')
ax[0][0].set_ylabel('Density')
ax[0][1].plot(radius1,entropy)
ax[0][1].set_xlabel('Radius')
ax[0][1].set_ylabel('Entropy')
ax[0][2].plot(radius1,temperature)
ax[0][2].set_xlabel('Radius')
ax[0][2].set_ylabel('Temperature')
ax[1][0].plot(radius1,dsdr)
ax[1][0].set_xlabel('Radius')
ax[1][0].set_ylabel('Entropy Gradient')
ax[1][1].plot(radius1,pressure)
ax[1][1].set_xlabel('Radius')
ax[1][1].set_ylabel('Pressure')
ax[1][2].plot(radius1,gravity)
ax[1][2].set_xlabel('Radius')
ax[1][2].set_ylabel('Gravity')
ax[2][0].plot(radius1,dlnrho)
ax[2][0].set_xlabel('Radius')
ax[2][0].set_ylabel('dlnrho')
ax[2][1].plot(radius1,d2lnrho)
ax[2][1].set_xlabel('Radius')
ax[2][1].set_ylabel('d2lnrho')
ax[2][2].plot(radius1,dlnt)
ax[2][2].set_xlabel('Radius')
ax[2][2].set_ylabel('dlnt')
plt.tight_layout()
plt.show()
print(density[0],temperature[0])
[ ]:
## This is where we define a heating function, if we want one in our model.
# Units are energy / volume / time such that {rho_hat T_hat dS/dt} = hprofile(r)
hprofile = numpy.zeros(nr,dtype='float64')
hprofile[:] = pressure1[:]
# Next, we need to integrate the heating profile
# We normalize it such that its integral over the volume is 1
# This way, we can set the luminosity via a constant in the input file for the dimensional case!
integrand= numpy.pi*4*radius*radius*hprofile
hint = numpy.trapz(integrand,x=radius)
hprofile = hprofile/hint
###########################################################
# Plot the integrated luminosity as a function of radius
# (should integrate to 1 at r = r_top)
###########################################################
# We then need to calculate the non-dimensional heating function
lq1 = numpy.zeros(nr)
lq1[0]=0
lq2 = numpy.zeros(nr)
lq2[0]=0
lq3 = numpy.zeros(nr)
lq3[0]=0
lsun= 3.846e33 # solar luminosity
integrand1= lsun*radius*radius*hprofile # Luminosity*r^2*heating
for i in range(1,nr):
lq1[i] =(1/(radius[i]**2.))*numpy.trapz(integrand1[0:i+1],x=radius[0:i+1])
integrand2=4*numpy.pi*radius**2.
lq2= numpy.trapz(integrand2,x=radius)
integrand3=lq1*4*numpy.pi*radius**2.
lq3=numpy.trapz(integrand3,x=radius)/lq2 # That is the volume Flux F_tilde
# The non-dimensional heating profile is hprofile_nd = lsun*hprofile*L/F_tilde=Q_dim*(r_o-r_cz)/F_tilde
hprofile_nd=lsun*hprofile*(ro-ri)/lq3
#print((ro-ri)/lq3) # This comes out of the non-dimensionalization used (L/F_tilde)
nu=8.e12
kappa=8.e12
Ra=gravity_av*lq3*(ro-ri)**4./(cp*density_av*temperature_av*nu*kappa**2.) # Ra is defined earlier in the notes
Di=gravity_av*(ro-ri)/(cp*temperature_av) # Dissipation number is defined in the notes above
print(Ra,Di/Ra)
fig, ax = plt.subplots(ncols=2,figsize=(12,4))
ax[0].plot(radius1,hprofile_nd,'ob')
ax[0].set_title('Non-Dimensional Heating Profile')
ax[0].set_xlabel('Radius')
ax[1].plot(radius,hprofile,'.b')
ax[1].set_xlabel('Radius (cm)')
ax[1].set_title('Dimensional Heating Profile')
plt.show()
[ ]:
# Have everything in terms of the non-dimensional radius
my_ref = rt.equation_coefficients(radius1)
[ ]:
# 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 # buoyancy term calculation
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_nd,6) # heating profile
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
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 constants will explicitly depend on the non-dimensionalization chosen.
# The comments corresponding to each one of the constants are generic but we also specify
# what they are exactly in our example here, assuming the non-dimensionalization we used.
# This is a non-magnetic, non-rotating example that reproduces the results from Featherstone & Hindman (2016)
# for the case with Nrho=3,and kappa=nu=8e12
my_ref.set_constant(1.0,1) # multiplies the Coriolis force, here it is: 2/E
my_ref.set_constant(13303.43109666924,2) # multiplies buoyancy -- Ra/Pr here
my_ref.set_constant(1.0,3) # multiplies the pressure gradient
my_ref.set_constant(0.0 , 4) # multiplies the lorentz force, here it is: (1/(4*pi)*1/(Pm*E))
my_ref.set_constant(1.0,5) # multiplies viscosity
my_ref.set_constant(1.0,6) # multiplies entropy diffusion , here it is 1/Pr
my_ref.set_constant(0.0,7) # multiplies magnetic diffusion in induction equation, here: 1/Pm
my_ref.set_constant(0.00012954929449971041,8) # multiplies viscous heating, here it is: Di*Pr/Ra
my_ref.set_constant(1.0,9) # multiplies ohmic heating, here it is: (1/(4*pi)*Di*Pr/(E*Pm^2*Ra))
my_ref.set_constant(1.0,10) # multiplies the heating function -- here it is 1/Pr (if normalized to 1, this is the luminosity)
my_ref.write('CZtest.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()
have_run = True
except:
have_run = False
if (have_run):
fig, ax = plt.subplots(ncols=3,nrows=4, figsize=(16,4*4))
# Density and derivatives of lnrho
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.dlnrho,'yo')
ax[0][1].plot(radius1, dlnrho)
ax[0][1].set_xlabel('Radius')
ax[0][1].set_title('Log density gradient')
ax[0][2].plot(cref.radius,cref.d2lnrho,'yo')
ax[0][2].plot(radius1,d2lnrho)
ax[0][2].set_xlabel('Radius')
ax[0][2].set_title('d_dr{Log density gradient}')
# Temperature and derivative of lnT
ax[1][0].plot(cref.radius,cref.temperature,'yo')
ax[1][0].plot(radius1,temperature)
ax[1][0].set_xlabel('Radius')
ax[1][0].set_title('Temperature')
ax[1][1].plot(cref.radius, cref.dlnT,'yo')
ax[1][1].plot(radius1, dlnt)
ax[1][1].set_xlabel('Radius')
ax[1][1].set_title('Log temperature gradient')
# dS/dr
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('Entropy gradient')
# Buoyancy, Heating
# Note that you must build the buoyancy term from the functions/constants
ax[3][1].plot(cref.radius, cref.functions[:,1]*cref.constants[1],'yo')
ax[3][1].plot(radius1, gravity*density*Ra)
ax[3][1].set_xlabel('Radius')
ax[3][1].set_title('Buoyancy')
# Note that the output heating (cref.heating) is hprofile/density/temperature
ax[3][2].plot(cref.radius, cref.heating,'yo')
ax[3][2].plot(radius1, hprofile_nd/density/temperature)
ax[3][2].set_xlabel('Radius')
ax[3][2].set_title('Heating')
plt.tight_layout()
plt.show()
[ ]:
[ ]: