{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Custom reference state for a non-dimensional Boussinesq MHD setup in a convective spherical shell " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Custom reference state for a Boussinesq convective setup. We non-dimensionalize the MHD Boussinesq equations\n", "using the rotation period such that $[t]=1/\\Omega_o$ is the timescale, the shell depth $[r]=r_o-r_i=L$ is the lengthscale (where $r_o$ is the outer radius and\n", "$r_i$ is the inner radius), $[u]=L\\Omega_o$ is the velocity scale, and $[T]=\\Delta T$\n", "is the temperature scale*.\n", "Assuming that $\\Theta$ are the temperature perturbations,\n", "the non-dimensional Boussinesq equations can be written as:\n", "\n", "\\begin{equation}\n", "\\dfrac{\\partial\\vec{u}}{\\partial t}+\\vec{u}\\cdot\\nabla\\vec{u}+2\\hat{z}\\times\\vec{u}=-\\dfrac{\\nabla P}{\\rho_m}+{\\rm{Ra^*}}\\left(\\dfrac{r_o}{r}\\right)^n\\Theta \\hat{e}_r+\\dfrac{1}{4\\pi}\\dfrac{E}{\\rm{P_m}}(\\nabla\\times\\vec{B})\\times{\\vec{B}}+E\\nabla^2\\vec{u}\n", "\\end{equation}\n", "\n", "\n", "\\begin{equation}\n", "\\nabla\\cdot\\vec{u}=0\n", "\\end{equation}\n", "\n", "\\begin{equation}\n", "\\dfrac{\\partial\\Theta}{\\partial t}+\\vec{u}\\cdot\\nabla\\Theta=\\dfrac{E}{\\rm{Pr}}\\nabla^2\\Theta\n", "\\end{equation}\n", "\n", "and\n", "\n", "\\begin{equation}\n", "\\dfrac{\\partial\\vec{B}}{\\partial t}-\\nabla\\times(\\vec{u}\\times\\vec{B})=\\dfrac{E}{\\rm{P_m}}\\nabla^2\\vec{B}.\n", "\\end{equation}\n", "\n", "We have assumed that $\\nu$, $\\kappa$, and $\\eta$ are constants, $\\rho_m$ is the mean denisty and $n$ is the gravity power (hence e.g. $n=0$ for constant gravity). We also have the modified Rayleigh number Ra$^*$ given by\n", "\\begin{equation}\n", "{\\rm{Ra^*}}=\\dfrac{\\alpha g_o \\Delta T}{L\\Omega_o^2}=\\dfrac{\\rm{Ra}}{\\rm{Pr}}E^2,\n", "\\end{equation}\n", "where $\\vec{g}(r)=g_o(r_o/r)^n$ , Pr=$\\nu/\\kappa$ is the Prandtl number, $E=\\dfrac{\\nu}{L^2\\Omega_o}$ is the Ekman number and Ra$=\\alpha g_o\\Delta T L^3/(\\kappa\\nu)$. Finally, the magnetic Prandtl number is P$_m=\\nu/\\eta$.\n", "Then the corresponding fuctions $f$ used here are:\n", "\n", "$f_1(r)\\rightarrow 1$, \n", "\n", "$f_2(r)\\rightarrow (r_o/r)^n$,\n", "\n", "$f_3(r)\\rightarrow 1$,\n", "\n", "$f_4(r)\\rightarrow 1$, \n", "\n", "$f_5(r)\\rightarrow 1$, \n", "\n", "$f_6(r)\\rightarrow 0$, \n", "\n", "$f_7(r)\\rightarrow 1$, \n", "\n", "$f_8(r)\\rightarrow 0$,\n", "\n", "$f_9(r)\\rightarrow 0$, \n", "\n", "$f_{10}(r)\\rightarrow 0$, \n", "\n", "\n", "$f_{14}(r)\\rightarrow 0$,\n", "\n", "$f_{15}(r)\\rightarrow 0$, \n", "\n", "$f_{16}(r)\\rightarrow 1$ \n", "\n", "and the constants $c$ are:\n", "\n", "$c_1\\rightarrow 2$,\n", "\n", "$c_2\\rightarrow Ra E^2/Pr$, \n", "\n", "$c_3\\rightarrow 1$,\n", "\n", "$c_4\\rightarrow E/(4\\pi{\\rm{P_m})}$,\n", "\n", "$c_5\\rightarrow E$,\n", "\n", "$c_6\\rightarrow E/{\\rm{Pr}}$, \n", "\n", "$c_7\\rightarrow E/{\\rm{P_m}}$, \n", "\n", "$c_8\\rightarrow 0$,\n", "\n", "$c_9\\rightarrow 0$, \n", "\n", "$c_{10}\\rightarrow 0$.\n", "\n", "*This is the relevant temperature scale for isothermal BCs -- for fixed flux, fixed temperature BCs, the temperature scale should be something like $[T]=L dT_o/dr$." ] }, { "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" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Grid Parameters\n", "nr = 500 # Number of radial points\n", "ri = 0.7e0 # Inner boundary of radial domain\n", "ro = 1.0e0 # Outer boundary of radial domain\n", "\n", "# radial grid\n", "r=numpy.linspace(ri,ro,nr)\n", "\n", "#aspect ratio\n", "beta=ri/ro \n", "\n", "# shell depth depending on the non-dimensionalization\n", "d=1.e0 \n", "\n", "#non-dimensional r_i\n", "ri_nd=beta*d/(1-beta)\n", "\n", "#non-dimensional r_o\n", "ro_nd=d/(1-beta)\n", "\n", "#non-dimensional radial grid\n", "radius1=numpy.linspace(ri_nd,ro_nd,nr)\n", "print(ri_nd,ro_nd)\n", "#print(radius1[0],radius1[nr-1])\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ones = numpy.ones(nr,dtype='float64')\n", "zeros = numpy.zeros(nr,dtype='float64')\n", "\n", "# Here we define the reference state i.e. density, temperature,etc.\n", "# For a classic RBC setup, this is the reference state to be used.\n", "\n", "density = ones # density rho\n", "dlnrho= zeros # dlnrho/dr\n", "d2lnrho= zeros # d^2lnrho/dr^2\n", "temperature=ones # temperature T\n", "dlnt=zeros # dlnT/dr\n", "pressure=ones # pressure P\n", "entropy = zeros # entropy S, not used in Boussinesq -- set it = 0\n", "gravity=zeros # gravity -- it is part of the buoyancy term in the non-dimensional momentum equation (see notes above)\n", "hprofile=zeros # heating function (if we want one)\n", "dsdr=zeros # dS/dr -- not useful in Boussinesq -- set it = 0" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "my_ref = rt.equation_coefficients(radius1)" ] }, { "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). For more info, check notes above!\n", "## Also, chech main_input_Boussinesq to see how to run a simulation with Rayleigh and this custom\n", "## reference state ( \"Boussinesq.dat\" input file generated below!)\n", "\n", "unity = numpy.ones(nr,dtype='float64')\n", "gravity_power=0.0\n", "buoy = (radius1[nr-1]/radius1)**gravity_power # buoyancy term calculation\n", "my_ref.set_function(density,1) # denisty 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 function \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\n", "my_ref.set_function(dsdr,14) # This is not used in Boussinesq -- set it = 0\n", "\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. These constants will explicitly depend on the non-dimensionalization chosen.\n", "\n", "# The comments corresponding to each one of the constants are generic but here, we also specify \n", "# what they are exactly in our example which is based on the non-dimensionalization used in this notebook!\n", "\n", "## This is a non-magnetic example with Pr=1, E=0.001, and Ra=10^6, such that Ra*=1 !\n", "my_ref.set_constant(2.0,1) # multiplies the Coriolis term, here it is: 2\n", "my_ref.set_constant(1.0,2) # multiplies the buoyancy, here it is Ra*=Ra.E^2/Pr as defined above\n", "my_ref.set_constant(1.0,3) # multiplies the pressure gradient\n", "my_ref.set_constant(0.0 , 4) # multiplies the lorentz force, here it is: E/(4*pi*Pm)\n", "my_ref.set_constant(0.001e0,5) # multiplies the viscosity, here it is: E\n", "my_ref.set_constant(0.001e0,6) # multiplies the entropy diffusion (kappa), here it is: E/Pr\n", "my_ref.set_constant(0.0,7) # multiplies eta in induction equation, here it is E/Pm\n", "my_ref.set_constant(0.0,8) # multiplies viscous heating, here it is always 0, since we assume the Boussinesq approximation\n", "my_ref.set_constant(1.0,9) # multiplies ohmic heating, here it is always 0, since we assume the Boussinesq approximation\n", "my_ref.set_constant(1.0,10) # multiplies the heating, here it is 0, since we have assumed that there is no heating function!\n", "my_ref.write('Boussinesq.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)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "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.5" } }, "nbformat": 4, "nbformat_minor": 2 }