{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Custom reference state for the dimensional anelastic MHD formulation in a convective spherical shell based on polytropes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem Setup\n", "\n", "We define the reference state of the convective region based on a polytrope that has a form given by \n", "\n", "$\\rho=\\rho_0 z^n$,\n", "\n", "$P=P_0 z^{n+1}$, \n", "\n", "and \n", "\n", "$T=T_0 z$, \n", "\n", "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 \n", "\n", "$P=R\\rho T$,\n", "\n", "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 \n", "\n", "$R=c_p-c_v=c_p(1-\\frac{1}{\\gamma})$, \n", "\n", "with \n", "\n", "$\\gamma\\equiv \\frac{c_p}{c_v}$.\n", "\n", "For a monotomic ideal gas, we have that\n", "\n", "$\\gamma\\equiv \\frac{5}{3}$. \n", "\n", "\n", "Finally, we require that the polytrope satisfies hydrostatic balance. Namely,\n", "\n", "$\\rho\\frac{GM}{r^2}=-\\frac{\\partial P}{\\partial r}$, \n", "\n", "where $G$ is the gravitational constant, and $M$ is the mass of the star.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Polytropic Solution\n", "\n", "Substituting our relations for $P$, $\\rho$, and $T$ into the equation of hydrostatic balance, we arive at\n", "\n", "$\\frac{\\partial z}{\\partial r}= - \\frac{GM}{(n+1)RT_0 r^2} = - \\frac{2GM}{5(n+1)c_pT_0 r^2}$ .\n", "\n", "This motivates us to seek a form for $z$ of \n", "\n", "$z = a +\\frac{b}{r}$,\n", "\n", "and immediately, we see that b must be given by\n", "\n", "$b = \\frac{2GM}{5(n+1)c_pT_0}$.\n", "\n", "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", "$N_\\rho = \\frac{ \\rho_b }{ \\rho_t } = \\frac{z_b^n}{z_t^n}$, \n", "\n", "or equivalently, using $C$ to denote the exponential factor, and $\\beta \\equiv \\frac{r_b}{r_t}$, we have\n", "\n", "$C\\equiv e^{\\frac{N_\\rho}{n}} = \\frac{a+b/r_b}{a+b/r_t}$.\n", "\n", "Rearranging, we find our expression for $a$\n", "\n", "$a = \\frac{fb}{r_b}$, \n", "\n", "where \n", "\n", "$f \\equiv \\frac{\\beta C-1}{1-C}$. \n", "\n", "This yields our expression for $z$ in terms of $T_0$\n", "\n", "$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)$.\n", "\n", "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\n", "\n", "$T_0 = T_b = \\frac{2GM}{5(n+1)c_p}\\left(\\frac{1+f}{r_b}\\right)$,\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#######################\n", "import numpy\n", "import matplotlib.pyplot as plt\n", "from matplotlib.pyplot import plot, draw, show\n", "\n", "import os, sys\n", "sys.path.insert(0, os.path.abspath('../../'))\n", "\n", "import post_processing.reference_tools as rt # You will need the refernce_tools.py to run this notebook\n", "import post_processing.rayleigh_diagnostics as rdiag" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Grid Parameters\n", "# Here, we use solar values as an example\n", "\n", "nr = 512 # Number of radial points\n", "ri = 5e10 # Inner boundary of radial domain\n", "ro = 6.586e10 # Outer boundary of radial domain\n", "rcz = 5e10 # base of the CZ \n", "\n", "#apsect ratio\n", "beta=ri/ro \n", "\n", "#Polytrope Parameters\n", "ncz = 1.5 # polytropic index of convection zone\n", "nrho = 3. # number of density scaleheights across convection zone (not full domain)\n", "mass = 1.98891e33 # Mass of the star\n", "G = 6.67e-8 # gravitational constant\n", "rhoi = 0.18053428 # density at the base of convection zone\n", "cp = 3.5e8 # specific heat at constant pressure\n", "gamma = 5.0/3.0 # Ratio of specific heats for ideal gas\n", "\n", "lsun = 3.846e33 # solar luminosity for scaling the volumetric heating" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "radius = numpy.linspace(ri,ro,nr) #Radial domain of the CZ\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Compute a CZ polytrope (see reference_tools)\n", "\n", "poly1 = rt.gen_poly(radius,ncz,nrho,mass,rhoi,G,cp,rcz)\n", "gas_constant = cp*(1.0-1.0/gamma) # R\n", "\n", "temperature= poly1.temperature # temperature T\n", "density = poly1.density # density rho\n", "pressure = poly1.pressure # pressure P, this won't matter -- set it equal to something\n", "dsdr = poly1.entropy_gradient # dS/dr\n", "entropy= poly1.entropy # entropy S, this won't matter -- set it equal to something\n", "dpdr = poly1.pressure_gradient # dP/dr\n", "gravity= mass*G/radius**2 # gravity g" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Calculation of the first and second derivative of lnrho and the first derivative of lnT.\n", "# If we do not calculate those here, then they are calculated within Rayleigh!\n", "\n", "d_density_dr = numpy.gradient(density,radius, edge_order=2)\n", "dlnrho = d_density_dr/density\n", "d2lnrho = numpy.gradient(dlnrho,radius, edge_order=2)\n", "dtdr = numpy.gradient(temperature,radius, edge_order =2)\n", "dlnt = dtdr/temperature" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plots of the dimensional reference state profiles, i.e. the density profile, the temperature profile, etc.\n", "\n", "fig, ax = plt.subplots(nrows =3,ncols = 3, figsize=(15,10) )\n", "ax[0][0].plot(radius,density,'r')\n", "ax[0][0].set_xlabel('Radius')\n", "ax[0][0].set_ylabel('Density')\n", "\n", "ax[0][1].plot(radius,entropy)\n", "ax[0][1].set_xlabel('Radius')\n", "ax[0][1].set_ylabel('Entropy')\n", "\n", "ax[0][2].plot(radius,temperature)\n", "ax[0][2].set_xlabel('Radius')\n", "ax[0][2].set_ylabel('Temperature')\n", "\n", "ax[1][0].plot(radius,dsdr)\n", "ax[1][0].set_xlabel('Radius')\n", "ax[1][0].set_ylabel('Entropy Gradient')\n", "\n", "ax[1][1].plot(radius,pressure)\n", "ax[1][1].set_xlabel('Radius')\n", "ax[1][1].set_ylabel('Pressure')\n", "\n", "ax[1][2].plot(radius,gravity)\n", "ax[1][2].set_xlabel('Radius')\n", "ax[1][2].set_ylabel('Gravity')\n", "\n", "ax[2][0].plot(radius,dlnrho)\n", "ax[2][0].set_xlabel('Radius')\n", "ax[2][0].set_ylabel('dlnrho')\n", "\n", "ax[2][1].plot(radius,d2lnrho)\n", "ax[2][1].set_xlabel('Radius')\n", "ax[2][1].set_ylabel('d2lnrho')\n", "\n", "ax[2][2].plot(radius,dlnt)\n", "ax[2][2].set_xlabel('Radius')\n", "ax[2][2].set_ylabel('dlnt')\n", "\n", "plt.tight_layout()\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Calculation of the heating profile based on the pressure, if we want to have a heating function in our setup.\n", "\n", "hprofile = numpy.zeros(nr,dtype='float64')\n", "hprofile[:] = pressure[:]\n", "\n", "#################################################################\n", "# We normalize the heating function so that it integrates to 1.\n", "\n", "integrand= numpy.pi*4*radius*radius*hprofile\n", "hint = numpy.trapz(integrand,x=radius)\n", "hprofile = hprofile/hint\n", "#plt.plot(radius,hprofile)\n", "plt.plot(radius, hprofile*lsun)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "my_ref = rt.equation_coefficients(radius)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Here we define all the functions and constants that will be written in our data file and\n", "# read by Rayleigh if we choose the custom reference state (=4)\n", "\n", "unity = numpy.ones(nr,dtype='float64')\n", "buoy =density*gravity/cp # calculation of the buoyancy term used in the momentum equation\n", "\n", "my_ref.set_function(density,1) # density rho\n", "my_ref.set_function(buoy,2) # buoyancy term\n", "my_ref.set_function(unity,3) # nu(r) -- can be overwritten via nu_type in Rayleigh\n", "my_ref.set_function(temperature,4) # temperature T\n", "my_ref.set_function(unity,5) # kappa(r) -- works like nu\n", "my_ref.set_function(hprofile,6) # heating -- this is normalized to one\n", "\n", "my_ref.set_function(dlnrho,8) # dlnrho/dr\n", "my_ref.set_function(d2lnrho,9) # d^2lnrho/dr^2\n", "my_ref.set_function(dlnt,10) # dlnT/dr\n", "my_ref.set_function(unity,7) # eta -- works like nu and kappa -- magnetic diffusivity\n", "my_ref.set_function(dsdr,14) # dS/dr \n", "\n", "# The constants can all be set/overridden in the input file\n", "# NOTE that they default to ZERO, but we want\n", "# most of them to be UNITY. \n", "# These aren't very useful in the dimensional anelastic formulation but are VERY useful\n", "# for the non-dimensional version of the anelastic custom reference state!\n", "\n", "my_ref.set_constant(1.0,1) # multiplies the Coriolis force--Should be 2 Omega (complication here)\n", "my_ref.set_constant(1.0,2) # multiplies buoyancy\n", "my_ref.set_constant(1.0,3) # multiplies pressure gradient\n", "my_ref.set_constant(0.0,4) # multiplies lorentz force\n", "my_ref.set_constant(1.0,5) # multiplies viscosity\n", "my_ref.set_constant(1.0,6) # multiplies entropy diffusion (kappa)\n", "my_ref.set_constant(0.0,7) # multiplies eta in induction equation\n", "my_ref.set_constant(1.0,8) # multiplies viscous heating\n", "my_ref.set_constant(1.0,9) # multiplies ohmic heating\n", "my_ref.set_constant(lsun,10) # multiplies the heating (if normalized to 1, this should be the luminosity)\n", "my_ref.write('dimensional.dat') # Here we write our data file to be used to run our simulation with Rayleigh!\n", "print(my_ref.fset)\n", "print(my_ref.cset)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Once you've run for one time step, set have_run = True !!\n", "\n", "# Here we check if the output reference state is the same as the one we used as an input (sanity check)\n", "\n", "# NOTE: We need the output file \"equation_coefficients\" to run this, as well as the PDE_Coefficients\n", "# from rayleigh_diagnostics.py\n", "\n", "try:\n", " cref = rdiag.PDE_Coefficients() # This will give us the output reference state\n", " have_run = True\n", "except:\n", " have_run = False\n", " \n", "if (have_run):\n", "\n", " fig, ax = plt.subplots(ncols=3,nrows=4, figsize=(16,4*4))\n", " # density, derivatives of lnrho\n", " ax[0][0].plot(cref.radius,cref.density,'yo')\n", " ax[0][0].plot(radius,density)\n", " ax[0][0].set_xlabel('Radius (cm)')\n", " ax[0][0].set_title('Density')\n", " \n", " ax[0][1].plot(cref.radius, cref.dlnrho,'yo')\n", " ax[0][1].plot(radius, dlnrho)\n", " ax[0][1].set_xlabel('Radius (cm)')\n", " ax[0][1].set_title('Log density gradient')\n", " \n", " ax[0][2].plot(cref.radius,cref.d2lnrho,'yo')\n", " ax[0][2].plot(radius,d2lnrho)\n", " ax[0][2].set_xlabel('Radius (cm)')\n", " ax[0][2].set_title('d_dr{Log density gradient}')\n", " \n", " # temperature and derivative of lnT \n", " ax[1][0].plot(cref.radius,cref.temperature,'yo')\n", " ax[1][0].plot(radius,temperature)\n", " ax[1][0].set_xlabel('Radius (cm)')\n", " ax[1][0].set_title('Temperature')\n", " \n", " ax[1][1].plot(cref.radius, cref.dlnT,'yo')\n", " ax[1][1].plot(radius, dlnt)\n", " ax[1][1].set_xlabel('Radius (cm)')\n", " ax[1][1].set_title('Log temperature gradient')\n", "\n", " # entropy gradient \n", " ax[2][1].plot(cref.radius, cref.dsdr,'yo')\n", " ax[2][1].plot(radius, dsdr,'b')\n", " ax[2][1].set_xlabel('Radius (cm)')\n", " ax[2][1].set_title('Entropy gradient') \n", " \n", " # Note that you must build the buoyancy term from the functions/constants\n", " ax[3][1].plot(cref.radius, cref.functions[:,1]*cref.constants[1],'yo')\n", " ax[3][1].plot(radius, gravity*density/cp)\n", " ax[3][1].set_xlabel('Radius (cm)')\n", " ax[3][1].set_title('Buoyancy') \n", " \n", " # Note that the output heating (cref.heating) is hprofile/density/temperature*luminosity\n", " ax[3][2].plot(cref.radius, cref.heating*cref.rho*cref.T,'yo')\n", " ax[3][2].plot(radius, hprofile*lsun)\n", " ax[3][2].set_xlabel('Radius (cm)')\n", " ax[3][2].set_title('Heating') \n", " \n", " plt.tight_layout()\n", " plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 2 }