{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boussinesq dynamo: nondimensionalized using viscous timescale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "\n",
    "This code provides an example for using a custom non-dimensionalization\n",
    "of Rayleigh in Boussinesq dynamo mode.\n",
    "\n",
    "The non-dimensionalization used here is based on the Rayleigh default\n",
    "viscous diffusion scaling for convection, thus providing a check\n",
    "on the custom reference state.\n",
    "\n",
    "The parameters correspond to a case listed in Table 2 of \n",
    "Soderlund et al.: \"The influence of magnetic fields in planetary\n",
    "dynamo models\", Earth Planet. Sci. Lett., v.333-334, p.9-20 (2012)\n",
    "\n",
    "The numbers referenced below for the various functions and constants \n",
    "refer to equation (5) in \"Rayleigh_Output_Variables.pdf\". Please refer to \n",
    "that document for further details.\n",
    "\n",
    "Requirements: (1) \"rayleigh_diagnostics.py\" ; and (2) \"reference_tools.py\"\n",
    "\n",
    "'''\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "import os, sys\n",
    "sys.path.insert(0, os.path.abspath('../../'))\n",
    "\n",
    "import post_processing.reference_tools as rt\n",
    "\n",
    "# name of output file containing custom reference data\n",
    "filename = 'custom_ref_viscous.dat'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Non-dimensional input parameters\n",
    "\n",
    "Ra = 1.12e5             # Rayleigh number\n",
    "Pr = 1.0                # Prandtl number\n",
    "Ek = 2.0e-3             # Ekman number\n",
    "Pm = 5.0                # Magnetic Prandtl number\n",
    "beta = 0.4              # Aspect ratio = r_inner/r_outer\n",
    "gravity_power = 1.0     # power law variation of gravitational acceleration\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create radial grid\n",
    "\n",
    "# Numer of radial grid points for radial functions (f_1, f_2, etc.)\n",
    "# Make large enough for accurate interpolation onto Chebyshev grid\n",
    "nr = 2000 \n",
    "\n",
    "# non-dimensional r_inner\n",
    "ri = beta/(1-beta)\n",
    "\n",
    "# non-dimensional r_outer\n",
    "ro=1.0/(1-beta)\n",
    "\n",
    "# non-dimensional radial grid\n",
    "radius=np.linspace(ri,ro,nr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the reference state functions and constants\n",
    "\n",
    "ones = np.ones(nr,dtype='float64')\n",
    "zeros = np.zeros(nr,dtype='float64')\n",
    "\n",
    "# the function list below is default for Boussinesq\n",
    "f_1 = ones\n",
    "f_2 = (radius/radius[nr-1])**gravity_power \n",
    "f_3 = ones\n",
    "f_4 = ones\n",
    "f_5 = ones\n",
    "f_6 = zeros\n",
    "f_7 = ones\n",
    "f_8 = zeros\n",
    "f_9 = zeros\n",
    "f_10 = zeros\n",
    "f_11 = zeros\n",
    "f_12 = zeros\n",
    "f_13 = zeros\n",
    "\n",
    "c_1 = 2.0/Ek                # Coriolis force\n",
    "c_2 = Ra/Pr                 # Buoyancy force\n",
    "c_3 = 1.0/Ek                # Pressure gradient \n",
    "c_4 = 1.0/(Ek*Pm)           # Lorentz force\n",
    "c_5 = 1.0                   # Viscous force\n",
    "c_6 = 1.0/Pr                # Thermal diffusion\n",
    "c_7 = 1.0/Pm                # Ohmic diffusion\n",
    "c_8 = 0.0                   # Viscous heating\n",
    "c_9 = 0.0                   # Ohmic heating\n",
    "c_10 = 0.0                  # Internal heating\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set all of the functions and constants\n",
    "\n",
    "my_ref = rt.equation_coefficients(radius)\n",
    "\n",
    "# Set functions here\n",
    "my_ref.set_function(f_1, 1)\n",
    "my_ref.set_function(f_2, 2)\n",
    "my_ref.set_function(f_3, 3)\n",
    "my_ref.set_function(f_4, 4)\n",
    "my_ref.set_function(f_5, 5)\n",
    "my_ref.set_function(f_6, 6)\n",
    "my_ref.set_function(f_7, 7)\n",
    "my_ref.set_function(f_8, 8)\n",
    "my_ref.set_function(f_9, 9)\n",
    "my_ref.set_function(f_10, 10)\n",
    "my_ref.set_function(f_11, 11)\n",
    "my_ref.set_function(f_12, 12)\n",
    "my_ref.set_function(f_13, 13)\n",
    "\n",
    "# Set constants here\n",
    "my_ref.set_constant(c_1, 1)\n",
    "my_ref.set_constant(c_2, 2)  \n",
    "my_ref.set_constant(c_3, 3) \n",
    "my_ref.set_constant(c_4, 4) \n",
    "my_ref.set_constant(c_5, 5) \n",
    "my_ref.set_constant(c_6, 6) \n",
    "my_ref.set_constant(c_7, 7) \n",
    "my_ref.set_constant(c_8, 8) \n",
    "my_ref.set_constant(c_9, 9) \n",
    "my_ref.set_constant(c_10, 10) \n",
    "\n",
    "my_ref.write(filename)\n",
    "\n",
    "print('Custom reference file', filename, 'was written successfully.')\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
}