ShockTests/Tutorial-GRMHD-2D_ShockTests.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# 2D Balsara Shock Tests\n",
    "[comment]: <> (Author: TODO)\n",
    "## This module provides initial data options for GRMHD evolutions.\n",
    "\n",
    "## This notebook develops and implements 2D shock tests, specifically the cylindrical explosion, magnetic rotor, and magnetic loop advection tests, used for magnetohydrodynamic codes.\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). \n",
    "\n",
    "### NRPy+ Source Code for this module: [ShockTests_2D.py](../../edit/in_progress-ShockTests/ShockTests_2D.py)\n",
    "\n",
    "## Introduction:\n",
    "\n",
    "Define shocks.\n",
    "\n",
    "Shocks can occur commonly in astrophysical contexts. Among other methods, shocks may be modeled using codes that solve the equations of magnetohydrodynamics. Here we will develop 1-3D shock tests for MHD codes. Within NRPy+, we will specifically use these tests to validate the NRPy+ version of IllinoisGRMHD. We will borrow the same tests used by the Spritz code for its validation, described [here](https://iopscience.iop.org/article/10.1088/1361-6382/ab8be8).\n",
    "\n",
    "In this notebook, we document and implement 2D shock tests commonly used within the literature. Specifically, we will implement the cylindrical explosion, magnetic rotor, and magnetic loop advection tests.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents:\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#initializenrpy): Import core NRPy+ modules and set NRPy+ parameters\n",
    "1. [Step 2](#cyl_expl): Cylindrical Explosion\n",
    "1. [Step 3](#mag_rotor): Magnetic Rotor\n",
    "1. [Step 4](#loop_adv): Loop Advection\n",
    "1. [Step 5](#code_validation): Code Validation against `ShockTests_2D` NRPy+ Module\n",
    "1. [Step 6](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Import core NRPy+ modules and set NRPy+ parameters \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "Here, we will import the NRPy+ core modules, set the reference metric to Cartesian, and set commonly used NRPy+ parameters. We will also set up a parameter to determine what initial data is set up, although it won't do much yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 0: Add NRPy's directory to the path\n",
    "# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n",
    "import os,sys\n",
    "nrpy_dir_path = os.path.join(\"..\")\n",
    "if nrpy_dir_path not in sys.path:\n",
    "    sys.path.append(nrpy_dir_path)\n",
    "\n",
    "# Step 0.a: Import the NRPy+ core modules and set the reference metric to Cartesian\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "\n",
    "par.set_parval_from_str(\"grid::DIM\", 3)\n",
    "DIM = par.parval_from_str(\"grid::DIM\")\n",
    "\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n",
    "rfm.reference_metric()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### <a id='cyl_expl'></a>\n",
    "\n",
    "# Step 2: Cylindrical Explosion \\[Back to [top](#toc)\\]\n",
    "$$\\label{cyl_expl}$$\n",
    "\n",
    "This test is described in section 4.2.1 of [Cipolletta et al](https://iopscience.iop.org/article/10.1088/1361-6382/ab8be8). Briefly, the test consists of a dense, over-pressurized central region surrounded by an ambient medium, with some smoothing function connecting the two regions. The entire domain has a uniform magnetic field\n",
    "\n",
    "Similar to our implementation of the 1D Balsara tests, we will make use of `Min_Max_and_Piecewise_Expressions` NRPy+ module, to replace `if` statements with `fabs` calls. The inner region goes out to radius $r_{in} = 0.8$, and the outer region begins at $r_{out} = 1.0$. For the inner region we have $\\rho_{in} = 10^{-2}, \\ p_{in} = 1.0$. The outer region has $\\rho_{out} = 10^{-4}, \\ p_{out} = 3 \\times 10^{-5}$.\n",
    "\n",
    "The density and pressure profiles are then given by\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "\\rho\\left( r \\right) &= \\left \\{ \\begin{array}{lll} \\rho_{in} & \\mbox{if} & r \\leq r_{in} \\\\\n",
    "\\exp \\left[ \\frac{\\left( r_{out} - r \\right) \\ln\\rho_{in} + \\left( r - r_{in} \\right) \\ln\\rho_{out}}{r_{out} - r_{in}} \\right] & \\mbox{if} & r_{in} < r < r_{out} \\\\ \n",
    "\\rho_{out} & \\mbox{if} & r \\geq r_{out} \\end{array} \\right.\\\\\n",
    "p\\left( r \\right) &= \\left \\{ \\begin{array}{lll} p_{in} & \\mbox{if} & r \\leq r_{in} \\\\\n",
    "\\exp \\left[ \\frac{\\left( r_{out} - r \\right) \\ln p_{in} + \\left( r - r_{in} \\right) \\ln p_{out}}{r_{out} - r_{in}} \\right] & \\mbox{if} & r_{in} < r < r_{out} \\\\ \n",
    "p_{out} & \\mbox{if} & r \\geq r_{out} \\end{array} \\right. \\\\\n",
    "\\end{align}\n",
    "\n",
    "We also have $v^i = 0$ and $B^x = 0.1, \\ B^y = B_z = 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Min_Max_and_Piecewise_Expressions as noif\n",
    "from sympy import Rational as rl\n",
    "\n",
    "def cylindrical_explosion(r, r_in=0.8, r_out=1.0):\n",
    "    vU = ixp.zerorank1()\n",
    "    BU = ixp.zerorank1()\n",
    "\n",
    "    r_out_minus_r_in = r_out - r_in\n",
    "    r_out_minus_r = r_out - r\n",
    "    r_minus_r_in = r - r_in\n",
    "\n",
    "    rho_in = rl(1,100)\n",
    "    rho_out = rl(1,10000)\n",
    "    rho_mid = sp.exp( (r_out_minus_r*sp.ln(rho_in)  +  r_minus_r_in*sp.ln(rho_out) ) / (r_out_minus_r_in) )\n",
    "\n",
    "    rho = noif.coord_leq_bound(r, r_in)*rho_in\\\n",
    "         +noif.coord_greater_bound(r, r_in)*noif.coord_less_bound(r, r_out)*rho_mid\\\n",
    "         +noif.coord_geq_bound(r, r_out)*rho_out\n",
    "\n",
    "\n",
    "    press_in = rl(1.0)\n",
    "    press_out = rl(3,100000)\n",
    "    press_mid = sp.exp( (r_out_minus_r*sp.ln(press_in)  +  r_minus_r_in*sp.ln(press_out) ) / (r_out_minus_r_in) )\n",
    "\n",
    "    press = noif.coord_leq_bound(r, r_in)*press_in\\\n",
    "           +noif.coord_greater_bound(r, r_in)*noif.coord_less_bound(r, r_out)*press_mid\\\n",
    "           +noif.coord_geq_bound(r, r_out)*press_out\n",
    "\n",
    "    BU[0] = rl(1,10)\n",
    "\n",
    "    return rho, press, vU, BU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### <a id='mag_rotor'></a>\n",
    "\n",
    "# Step 3: Magnetic Rotor \\[Back to [top](#toc)\\]\n",
    "$$\\label{mag_rotor}$$\n",
    "\n",
    "This test, described in section 4.2.2 of [Cipolletta et al](https://iopscience.iop.org/article/10.1088/1361-6382/ab8be8), consists of a rapidly spinning dense, central region surrounded by a static ambient medium. Pressure and magnetic field are uniform throughout the entire domain. For the fluid density we have $r_{in} = 0.1, \\ \\rho_{in} = 10.0,\\ \\rho_{out} = 1.0$, giving\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "\\rho\\left( r \\right) &= \\left \\{ \\begin{array}{lll} \\rho_{in} & \\mbox{if} & r \\leq r_{in} \\\\ \n",
    "\\rho_{out} & \\mbox{if} & r > r_{in} \\end{array} \\right.\\\\\n",
    "\\end{align}\n",
    "\n",
    "The inner region is rotating with uniform angular velocity $\\Omega = 9.95$, i.e.\n",
    "\n",
    "\\begin{align}\n",
    "v^\\phi\\left( r \\right) &= \\left \\{ \\begin{array}{lll} \\Omega & \\mbox{if} & r \\leq r_{in} \\\\ \n",
    "\\rho_{out} & \\mbox{if} & r > r_{in} \\end{array} \\right.\\\\\n",
    "\\end{align}\n",
    "\n",
    "Below we basis transform $v^i = \\left(0, \\Omega, 0 \\right)$ from cylindrical coordinates to Cartesian coordinates. The entire domain has $B^i = \\left( 1.0, 0, 0 \\right), \\ p = 1.0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def magnetic_rotor(r, r_in=0.1, Omega=9.95, cartx=rfm.Cartx, carty=rfm.Carty):\n",
    "    cart_list = [cartx, carty]\n",
    "\n",
    "    vU_cyl = ixp.zerorank1()\n",
    "    BU = ixp.zerorank1()\n",
    "\n",
    "    rho_in = rl(10.0)\n",
    "    rho_out = rl(1.0)\n",
    "\n",
    "    rho = noif.coord_leq_bound(r, r_in)*rho_in\\\n",
    "         +noif.coord_greater_bound(r, r_in)*rho_out\n",
    "\n",
    "    press = rl(1.0)\n",
    "\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cylindrical\")\n",
    "    rfm.reference_metric()\n",
    "\n",
    "    vU_Cyl = ixp.zerorank1()\n",
    "    vU_Cyl[1] = noif.coord_leq_bound(r, r_in)*Omega\n",
    "\n",
    "    Jac_dUCart_dDrfmUD, Jac_dUrfm_dDCartUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Cartesian()\n",
    "    vU = rfm.basis_transform_vectorU_from_rfmbasis_to_Cartesian(Jac_dUCart_dDrfmUD, vU_Cyl)\n",
    "\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            vU[i] = vU[i].subs(rfm.xx[j], rfm.Cart_to_xx[j])\n",
    "\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n",
    "    rfm.reference_metric()\n",
    "\n",
    "    for i in range(DIM):\n",
    "        for j in range(2):\n",
    "            vU[i] = vU[i].subs(rfm.Cart[j], cart_list[j])\n",
    "\n",
    "    BU[0] = rl(1.0)\n",
    "\n",
    "    return rho, press, vU, BU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### <a id='loop_adv'></a>\n",
    "\n",
    "# Step 4: Loop Advection \\[Back to [top](#toc)\\]\n",
    "$$\\label{loop_adv}$$\n",
    "\n",
    "This test, described in section 4.2.3 of [Cipolletta et al](https://iopscience.iop.org/article/10.1088/1361-6382/ab8be8), consists of a magnetized circular loop propagating through an ambient, unmagnetized medium at constant velocity.\n",
    "\n",
    "The magnetic field is given by\n",
    "\n",
    "\\begin{align}\n",
    "B^x,B^y &= \\left \\{ \\begin{array}{lll} A_{loop}y/r, A_{loop}x/r & \\mbox{if} & r < R_{loop} \\\\ \n",
    "0 & \\mbox{if} & r \\geq R_{loop} \\end{array} \\right.\\\\\n",
    "B^z &= 0\n",
    "\\end{align}\n",
    "\n",
    "The corresponding vector potential may be written as $A_i = \\left(0,0, \\max\\left[0, A_{loop} \\left(R_{loop} - r\\right)\\right] \\right)$. [Cipolletta et al](https://iopscience.iop.org/article/10.1088/1361-6382/ab8be8) chooses $\\rho=1.0, p=3.0, A_{loop}=0.001, R_{loop}=0.3$\n",
    "\n",
    "The velocity field is given by\n",
    "\n",
    "\\begin{align}\n",
    "v^x,v^y &= \\left \\{ \\begin{array}{lll} 1/12, 1/24 & \\mbox{if} & r < R_{loop} \\\\ \n",
    "0 & \\mbox{if} & r \\geq R_{loop} \\end{array} \\right.\\\\\n",
    "\\end{align}\n",
    "\n",
    "For $v^z$, we have the option to set it as $0$ or $1/24$. Cipolletta et al show results for the latter.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loop_advection(r, R_loop=rl(3,10), A_loop=rl(1/1000), vz_equals_zero=False):\n",
    "    vU = ixp.zerorank1()\n",
    "    AD = ixp.zerorank1()\n",
    "\n",
    "    rho = rl(1.0)\n",
    "\n",
    "    press = rl(3.0)\n",
    "\n",
    "    vx_in = rl(1,12)\n",
    "    vy_in = rl(1,24)\n",
    "\n",
    "    vU[0] = noif.coord_less_bound(r, R_loop)*vx_in\n",
    "    vU[1] = noif.coord_less_bound(r, R_loop)*vy_in\n",
    "    if not vz_equals_zero:\n",
    "        vU[2] = noif.coord_less_bound(r, R_loop)*vy_in\n",
    "\n",
    "    AD[2] = noif.max_noif(0, A_loop*(R_loop - r))\n",
    "\n",
    "    return rho, press, vU, AD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 5: Code Validation against `ShockTests_2D` NRPy+ Module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for the 1D shock tests we intend to use between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [`ShockTests_2D.py`](../../edit/in_progress-ShockTests/ShockTests_2D.py) module.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we modify the check_zero function\n",
    "from UnitTesting.assert_equal import check_zero\n",
    "\n",
    "def assert_zero(exp):\n",
    "    assert check_zero(exp)==True\n",
    "\n",
    "def compare_results(functionA, functionB, x, bound=0.0):\n",
    "    rho_A, press_A, vU_A, BU_A = functionA(x, bound)\n",
    "    rho_B, press_B, vU_B, BU_B = functionB(x, bound)\n",
    "\n",
    "    assert_zero(rho_A - rho_B)\n",
    "    assert_zero(press_A - press_B)\n",
    "\n",
    "    for i in range(DIM):\n",
    "        assert_zero(vU_A[i] - vU_B[i])\n",
    "        assert_zero(BU_A[i] - BU_B[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All assertions have passed!\n"
     ]
    }
   ],
   "source": [
    "import ShockTests_2D as st_2d\n",
    "\n",
    "r = rfm.xxSph[0]\n",
    "\n",
    "compare_results(cylindrical_explosion, st_2d.cylindrical_explosion, r)\n",
    "compare_results(magnetic_rotor, st_2d.magnetic_rotor, r)\n",
    "compare_results(loop_advection, st_2d.loop_advection, r)\n",
    "\n",
    "print('All assertions have passed!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-GRMHD-2D_ShockTests.pdf](Tutorial-GRMHD-2D_ShockTests.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-GRMHD-2D_ShockTests.tex, and compiled LaTeX file to PDF\n",
      "    file Tutorial-GRMHD-2D_ShockTests.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-GRMHD-2D_ShockTests\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}