{
"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",
"# Time Evolution of Maxwell's Equations in Flat Spacetime and Curvilinear Coordinates\n",
"\n",
"## Authors: Terrence Pierre Jacques, Zachariah Etienne and Ian Ruchlin\n",
"\n",
"## This module constructs the evolution equations for Maxwell's equations as symbolic (SymPy) expressions, for an electromagnetic field in vacuum, as defined in [Tutorial-VacuumMaxwell_formulation_Curvilinear](Tutorial-VacuumMaxwell_formulation_Curvilinear.ipynb).\n",
"\n",
"**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
"\n",
"**Validation Notes:** All expressions generated in this here have been validated against the [VacuumMaxwell_Flat_Evol_Curvilinear_rescaled](../edit/Maxwell/VacuumMaxwell_Flat_Evol_Curvilinear_rescaled.py) module, as well as the [Maxwell/VacuumMaxwell_Flat_Evol_Cartesian](../edit/Maxwell/VacuumMaxwell_Flat_Evol_Cartesian.py) module when setting the coordinate system to Cartesian.\n",
"\n",
"### NRPy+ Source Code for this module: [VacuumMaxwell_Flat_Evol_Curvilinear_rescaled](../edit/Maxwell/VacuumMaxwell_Flat_Evol_Curvilinear_rescaled.py)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='top'></a>\n",
"$$\\label{top}$$\n",
"\n",
"# Table of Contents: \n",
"\n",
"1. [Step 1](#step1): Set core NRPy+ parameters for numerical grids and reference metric\n",
"1. [Step 2](#step2): System II in curvilinear coordinates, using the rescaled quantities $a^i$ and $e^i$\n",
"1. [Step 3](#cart_transform): Convert $A^i$ and $E^i$ to the Cartesian basis\n",
"1. [Step 4](#step4): Code Validation\n",
"1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='step1'></a>\n",
"# Step 1: Preliminaries \\[Back to [top](#top)\\]\n",
"$$\\label{step1}$$\n",
"\n",
"Set up the needed NRPy+ infrastructure, such the number of dimensions and finite differencing order."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:53.472644Z",
"iopub.status.busy": "2021-03-07T17:18:53.471634Z",
"iopub.status.idle": "2021-03-07T17:18:54.312583Z",
"shell.execute_reply": "2021-03-07T17:18:54.311983Z"
}
},
"outputs": [],
"source": [
"# Import needed Python modules\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",
"import grid as gri # NRPy+: Functions having to do with numerical grids\n",
"import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
"\n",
"# Set the spatial dimension parameter to 3.\n",
"par.set_parval_from_str(\"grid::DIM\", 3)\n",
"DIM = par.parval_from_str(\"grid::DIM\")\n",
"\n",
"# Set coordinate system\n",
"# Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,\n",
"# SymTP, SinhSymTP\n",
"\n",
"CoordSystem = \"Spherical\"\n",
"par.set_parval_from_str(\"reference_metric::CoordSystem\",CoordSystem)\n",
"# Set reference metric related quantities\n",
"rfm.reference_metric()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='step2'></a>\n",
"\n",
"# Step 2: System II in Curvilinear Coordinates, using the rescaled quantities $a^i$ and $e^i$ \\[Back to [top](#top)\\]\n",
"$$\\label{step2}$$\n",
"(Following discussion reproduced from [Tutorial-VacuumMaxwell_formulation_Curvilinear](Tutorial-VacuumMaxwell_formulation_Curvilinear.ipynb))\n",
"\n",
"Consider an arbitrary vector $\\Lambda^i$, with smooth (continous) Cartesian components $\\Lambda^x$, $\\Lambda^y$, and $\\Lambda^z$. Transforming $\\Lambda^i$ to, e.g. spherical coordinates, introduces terms that spoil the smoothness of $\\Lambda^i$;\n",
"\n",
"$$\n",
"\\Lambda^\\phi = \\frac{1}{r \\sin \\theta} \\times \\left[ \\text{smooth part} \\right].\n",
"$$\n",
"\n",
"Evolving $\\Lambda^\\phi$ will introduce instabilities along the $z$-axis. To avoid this, we instead evolve the _rescaled_ quantity $\\lambda^i$, defined by \n",
"\n",
"$$\n",
"\\bar{\\Lambda}^i = \\frac{\\lambda^i}{\\text{scalefactor}[i]}.\n",
"$$\n",
"\n",
"where we use the [Hadamard product](https://en.wikipedia.org/w/index.php?title=Hadamard_product_(matrices)&oldid=852272177) of matrices, and no sums are implied by the repeated indices.\n",
"\n",
"Thus, we evolve the smoothed variable $\\lambda^i$, via \n",
"\n",
"$$\n",
"\\lambda^i = \\bar{\\Lambda}^i \\text{scalefactor}[i].\n",
"$$\n",
"\n",
"Within Nrpy+, ReU[i] = 1/scalefactor[i], giving \n",
"\n",
"$$\n",
"\\lambda^i = \\frac{\\bar{\\Lambda}^i}{\\text{ReU}[i]}.\n",
"$$\n",
"\n",
"We now define the rescaled quantities $a^i$ and $e^i$ and rewrite our formulation of Maxwell's equations in curvilinear coordinates;\n",
"\n",
"\\begin{align}\n",
"a^i &= \\frac{A^i}{\\text{ReU}[i]},\\\\ \\\\\n",
"e^i &= \\frac{E^i}{\\text{ReU}[i]},\n",
"\\end{align}\n",
"\n",
"Taking a time derivative on both sides,\n",
"\n",
"\\begin{align}\n",
"\\partial_t a^i &= \\frac{\\partial_t A^i}{\\text{ReU}[i]} = \\frac{ -E^i - \\hat{g}^{ij}\\partial_j \\varphi}{\\text{ReU}[i]} = -e^i - \\frac{\\hat{g}^{ij}\\partial_j \\varphi}{\\text{ReU}[i]},\\\\ \\\\\n",
"\\partial_t e^i &= \\frac{\\partial_t E^i}{\\text{ReU}[i]} = \\frac{\\hat{g}^{ij}\\partial_j \\Gamma - \\hat{g}^{jk} \\hat{\\nabla}_{j} \\left(\\hat{\\nabla}_{k} A^{i}\\right)}{\\text{ReU}[i]} = \\frac{\\hat{g}^{ij}\\partial_j \\Gamma}{\\text{ReU}[i]} - \\frac{\\hat{g}^{jk} \\hat{\\nabla}_{j} \\left(\\hat{\\nabla}_{k} a^{i} \\text{ReU}[i] \\right)}{\\text{ReU}[i]}.\n",
"\\end{align}\n",
"\n",
"Given that\n",
"\n",
"$$\n",
"\\partial_t E^i = {\\underbrace {\\textstyle \\hat{g}^{ij}\\partial_j \\Gamma}_{\\text{Term 1}}} - \\hat{\\gamma}^{jk} \\left({\\underbrace {\\textstyle A^i_{,kj}}_{\\text{Term 2}}} + {\\underbrace {\\textstyle \\hat{\\Gamma}^i_{mk,j} A^m + \\hat{\\Gamma}^i_{mk} A^m_{,j} + \\hat{\\Gamma}^i_{dj} A^d_{,k} - \\hat{\\Gamma}^d_{kj} A^i_{,d}}_{\\text{Term 3}}} + {\\underbrace {\\textstyle \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} A^m - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} A^m}_{\\text{Term 4}}}\\right),\n",
"$$\n",
"\n",
"we can make the following replacements within the above, in terms of NRPy+ code;\n",
"\n",
"\\begin{align}\n",
"A^i = \\text{AU[i]} &\\to \\text{aU[i] * rfm.ReU[i]} \\\\\n",
"\\partial_j A^i = \\text{AUdD[i][j]} &\\to \\text{aU_dD[i][j] * rfm.ReU[i]} +\\text{aU[i] * rfm.ReUdD[i][j]} \\\\\n",
"\\partial_k \\partial_j A^i = \\text{AUdDD[i][j][k]} &\\to \n",
"\\text{aU_dDD[i][j][k] * rfm.ReU[i]} + \n",
"\\text{aU_dDD[i][j] * rfm.ReUdD[i][k]} \\\\\n",
"&+ \\text{aU_dD[i][k] * rfm.ReUdD[i][j]} +\n",
"\\text{aU[i] * rfm.ReUdDD[i][j][k]}\n",
"\\end{align}\n",
"\n",
"The remainder of Maxwell's equations are unchanged;\n",
"\n",
"$$\n",
"\\partial_t \\Gamma = -\\hat{g}^{ij} \\left( \\partial_i \\partial_j \\varphi - \\hat{\\Gamma}^k_{ji} \\partial_k \\varphi \\right),\n",
"$$\n",
"\n",
"$$\n",
"\\partial_t \\varphi = -\\Gamma,\n",
"$$\n",
"\n",
"subject to constraints\n",
"\n",
"\\begin{align}\n",
"\\mathcal{G} &\\equiv \\Gamma - \\partial_i A^i + \\hat{\\Gamma}^i_{ji} A^j &= 0,\\\\\n",
"\\mathcal{C} &\\equiv \\partial_i E^i + \\hat{\\Gamma}^i_{ji} E^j &= 0.\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.324036Z",
"iopub.status.busy": "2021-03-07T17:18:54.323038Z",
"iopub.status.idle": "2021-03-07T17:18:54.325333Z",
"shell.execute_reply": "2021-03-07T17:18:54.325842Z"
}
},
"outputs": [],
"source": [
"# Register gridfunctions that are needed as input.\n",
"\n",
"# Declare the rank-1 indexed expressions e^{i}, e^{i},\n",
"# and \\partial^{i} \\psi, that are to be evolved in time.\n",
"# Derivative variables like these must have an underscore\n",
"# in them, so the finite difference module can parse\n",
"# the variable name properly.\n",
"\n",
"# e^i\n",
"eU = ixp.register_gridfunctions_for_single_rank1(\"EVOL\", \"eU\")\n",
"\n",
"# \\partial_k ( E^i ) --> rank two tensor\n",
"eU_dD = ixp.declarerank2(\"eU_dD\", \"nosym\")\n",
"\n",
"# a^i\n",
"aU = ixp.register_gridfunctions_for_single_rank1(\"EVOL\", \"aU\")\n",
"\n",
"# \\partial_k ( a^i ) --> rank two tensor\n",
"aU_dD = ixp.declarerank2(\"aU_dD\", \"nosym\")\n",
"\n",
"# \\partial_k partial_m ( a^i ) --> rank three tensor\n",
"aU_dDD = ixp.declarerank3(\"aU_dDD\", \"sym12\")\n",
"\n",
"# \\psi is a scalar function that is time evolved\n",
"psi = gri.register_gridfunctions(\"EVOL\", [\"psi\"])\n",
"\n",
"# \\Gamma is a scalar function that is time evolved\n",
"Gamma = gri.register_gridfunctions(\"EVOL\", [\"Gamma\"])\n",
"\n",
"# \\partial_i \\psi\n",
"psi_dD = ixp.declarerank1(\"psi_dD\")\n",
"\n",
"# \\partial_i \\Gamma\n",
"Gamma_dD = ixp.declarerank1(\"Gamma_dD\")\n",
"\n",
"# partial_i \\partial_j \\psi\n",
"psi_dDD = ixp.declarerank2(\"psi_dDD\", \"sym01\")\n",
"\n",
"ghatUU = rfm.ghatUU\n",
"GammahatUDD = rfm.GammahatUDD\n",
"GammahatUDDdD = rfm.GammahatUDDdD\n",
"ReU = rfm.ReU\n",
"ReUdD = rfm.ReUdD\n",
"ReUdDD = rfm.ReUdDD\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\partial_t a^i = -e^i - \\frac{\\hat{g}^{ij}\\partial_j \\varphi}{\\text{ReU}[i]},\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.337840Z",
"iopub.status.busy": "2021-03-07T17:18:54.336800Z",
"iopub.status.idle": "2021-03-07T17:18:54.340280Z",
"shell.execute_reply": "2021-03-07T17:18:54.339434Z"
}
},
"outputs": [],
"source": [
"# \\partial_t a^i = -e^i - \\frac{\\hat{g}^{ij}\\partial_j \\varphi}{\\text{ReU}[i]}\n",
"arhsU = ixp.zerorank1()\n",
"for i in range(DIM):\n",
" arhsU[i] -= eU[i]\n",
" for j in range(DIM):\n",
" arhsU[i] -= (ghatUU[i][j]*psi_dD[j])/ReU[i]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\partial_t e^i = \\frac{\\hat{g}^{ij}\\partial_j \\Gamma - \\hat{g}^{jk} \\hat{\\nabla}_{j} \\left(\\hat{\\nabla}_{k} A^{i}\\right)}{\\text{ReU}[i]} = \\frac{\\hat{g}^{ij}\\partial_j \\Gamma}{\\text{ReU}[i]} - \\frac{\\hat{g}^{jk} \\hat{\\nabla}_{j} \\left(\\hat{\\nabla}_{k} a^{i} \\text{ReU}[i] \\right)}{\\text{ReU}[i]}.\n",
"$$\n",
"\n",
"Given that\n",
"\n",
"$$\n",
"\\partial_t E^i = {\\underbrace {\\textstyle \\hat{g}^{ij}\\partial_j \\Gamma}_{\\text{Term 1}}} - \\hat{\\gamma}^{jk} \\left({\\underbrace {\\textstyle A^i_{,kj}}_{\\text{Term 2}}} + {\\underbrace {\\textstyle \\hat{\\Gamma}^i_{mk,j} A^m + \\hat{\\Gamma}^i_{mk} A^m_{,j} + \\hat{\\Gamma}^i_{dj} A^d_{,k} - \\hat{\\Gamma}^d_{kj} A^i_{,d}}_{\\text{Term 3}}} + {\\underbrace {\\textstyle \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} A^m - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} A^m}_{\\text{Term 4}}}\\right),\n",
"$$\n",
"\n",
"we can make the following replacements within the above, in terms of NRPy+ code;\n",
"\n",
"\\begin{align}\n",
"A^i = \\text{AU[i]} &\\to \\text{aU[i] * rfm.ReU[i]} \\\\\n",
"\\partial_j A^i = \\text{AUdD[i][j]} &\\to \\text{aU_dD[i][j] * rfm.ReU[i]} +\\text{aU[i] * rfm.ReUdD[i][j]} \\\\\n",
"\\partial_k \\partial_j A^i = \\text{AUdDD[i][j][k]} &\\to \n",
"\\text{aU_dDD[i][j][k] * rfm.ReU[i]} + \n",
"\\text{aU_dD[i][j] * rfm.ReUdD[i][k]} \\\\\n",
"&+ \\text{aU_dD[i][k] * rfm.ReUdD[i][j]} +\n",
"\\text{aU[i] * rfm.ReUdDD[i][j][k]}\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.358895Z",
"iopub.status.busy": "2021-03-07T17:18:54.358219Z",
"iopub.status.idle": "2021-03-07T17:18:54.360489Z",
"shell.execute_reply": "2021-03-07T17:18:54.361050Z"
}
},
"outputs": [],
"source": [
"# A^i\n",
"AU = ixp.zerorank1()\n",
"\n",
"# \\partial_k ( A^i ) --> rank two tensor\n",
"AU_dD = ixp.zerorank2()\n",
"\n",
"# \\partial_k partial_m ( A^i ) --> rank three tensor\n",
"AU_dDD = ixp.zerorank3()\n",
"\n",
"for i in range(DIM):\n",
" AU[i] = aU[i]*ReU[i]\n",
" for j in range(DIM):\n",
" AU_dD[i][j] = aU_dD[i][j]*ReU[i] + aU[i]*ReUdD[i][j]\n",
" for k in range(DIM):\n",
" AU_dDD[i][j][k] = aU_dDD[i][j][k]*ReU[i] + aU_dD[i][j]*ReUdD[i][k] +\\\n",
" aU_dD[i][k]*ReUdD[i][j] + aU[i]*ReUdDD[i][j][k]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\text{Term 1} = \\hat{g}^{ij}\\partial_j \\Gamma\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.367238Z",
"iopub.status.busy": "2021-03-07T17:18:54.366561Z",
"iopub.status.idle": "2021-03-07T17:18:54.369165Z",
"shell.execute_reply": "2021-03-07T17:18:54.368639Z"
}
},
"outputs": [],
"source": [
"# Term 1 = \\hat{g}^{ij}\\partial_j \\Gamma\n",
"Term1U = ixp.zerorank1()\n",
"for i in range(DIM):\n",
" for j in range(DIM):\n",
" Term1U[i] += ghatUU[i][j]*Gamma_dD[j]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\text{Term 2} = A^i_{,kj}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.374802Z",
"iopub.status.busy": "2021-03-07T17:18:54.374176Z",
"iopub.status.idle": "2021-03-07T17:18:54.376323Z",
"shell.execute_reply": "2021-03-07T17:18:54.376881Z"
}
},
"outputs": [],
"source": [
"# Term 2: A^i_{,kj}\n",
"Term2UDD = ixp.zerorank3()\n",
"for i in range(DIM):\n",
" for j in range(DIM):\n",
" for k in range(DIM):\n",
" Term2UDD[i][j][k] += AU_dDD[i][k][j]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\text{Term 3} = \\hat{\\Gamma}^i_{mk,j} A^m + \\hat{\\Gamma}^i_{mk} A^m_{,j} + \\hat{\\Gamma}^i_{dj}A^d_{,k} - \\hat{\\Gamma}^d_{kj} A^i_{,d}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.401689Z",
"iopub.status.busy": "2021-03-07T17:18:54.401081Z",
"iopub.status.idle": "2021-03-07T17:18:54.403838Z",
"shell.execute_reply": "2021-03-07T17:18:54.403249Z"
}
},
"outputs": [],
"source": [
"# Term 3: \\hat{\\Gamma}^i_{mk,j} A^m + \\hat{\\Gamma}^i_{mk} A^m_{,j}\n",
"# + \\hat{\\Gamma}^i_{dj}A^d_{,k} - \\hat{\\Gamma}^d_{kj} A^i_{,d}\n",
"Term3UDD = ixp.zerorank3()\n",
"for i in range(DIM):\n",
" for j in range(DIM):\n",
" for k in range(DIM):\n",
" for m in range(DIM):\n",
" Term3UDD[i][j][k] += GammahatUDDdD[i][m][k][j]*AU[m] \\\n",
" + GammahatUDD[i][m][k]*AU_dD[m][j] \\\n",
" + GammahatUDD[i][m][j]*AU_dD[m][k] \\\n",
" - GammahatUDD[m][k][j]*AU_dD[i][m]\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\text{Term 4} = \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} A^m - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} A^m\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.427336Z",
"iopub.status.busy": "2021-03-07T17:18:54.422187Z",
"iopub.status.idle": "2021-03-07T17:18:54.431105Z",
"shell.execute_reply": "2021-03-07T17:18:54.430567Z"
}
},
"outputs": [],
"source": [
"# Term 4: \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} A^m -\n",
"# \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} A^m\n",
"Term4UDD = ixp.zerorank3()\n",
"for i in range(DIM):\n",
" for j in range(DIM):\n",
" for k in range(DIM):\n",
" for m in range(DIM):\n",
" for d in range(DIM):\n",
" Term4UDD[i][j][k] += ( GammahatUDD[i][d][j]*GammahatUDD[d][m][k] \\\n",
" -GammahatUDD[d][k][j]*GammahatUDD[i][m][d])*AU[m]\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we build up the RHS of $E^i$,\n",
"\n",
"$$\n",
"\\partial_t E^i = {\\underbrace {\\textstyle \\hat{g}^{ij}\\partial_j \\Gamma}_{\\text{Term 1}}} - \\hat{\\gamma}^{jk} \\left({\\underbrace {\\textstyle A^i_{,kj}}_{\\text{Term 2}}} + {\\underbrace {\\textstyle \\hat{\\Gamma}^i_{mk,j} A^m + \\hat{\\Gamma}^i_{mk} A^m_{,j} + \\hat{\\Gamma}^i_{dj} A^d_{,k} - \\hat{\\Gamma}^d_{kj} A^i_{,d}}_{\\text{Term 3}}} + {\\underbrace {\\textstyle \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} A^m - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} A^m}_{\\text{Term 4}}}\\right),\n",
"$$\n",
"\n",
"and divide through by ReU[i] to get $e^i$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.484338Z",
"iopub.status.busy": "2021-03-07T17:18:54.464866Z",
"iopub.status.idle": "2021-03-07T17:18:54.486972Z",
"shell.execute_reply": "2021-03-07T17:18:54.486406Z"
}
},
"outputs": [],
"source": [
"# \\partial_t E^i = \\hat{g}^{ij}\\partial_j \\Gamma - \\hat{\\gamma}^{jk}*\n",
"# (A^i_{,kj}\n",
"# + \\hat{\\Gamma}^i_{mk,j} A^m + \\hat{\\Gamma}^i_{mk} A^m_{,j}\n",
"# + \\hat{\\Gamma}^i_{dj} A^d_{,k} - \\hat{\\Gamma}^d_{kj} A^i_{,d}\n",
"# + \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} A^m\n",
"# - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} A^m)\n",
"\n",
"ErhsU = ixp.zerorank1()\n",
"for i in range(DIM):\n",
" ErhsU[i] += Term1U[i]\n",
" for j in range(DIM):\n",
" for k in range(DIM):\n",
" ErhsU[i] -= ghatUU[j][k]*(Term2UDD[i][j][k] + Term3UDD[i][j][k] + Term4UDD[i][j][k])\n",
"\n",
"erhsU = ixp.zerorank1()\n",
"for i in range(DIM):\n",
" erhsU[i] = ErhsU[i]/ReU[i]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\partial_t \\Gamma = -\\hat{g}^{ij} \\left( \\partial_i \\partial_j \\varphi - \\hat{\\Gamma}^k_{ji} \\partial_k \\varphi \\right)\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.495298Z",
"iopub.status.busy": "2021-03-07T17:18:54.494592Z",
"iopub.status.idle": "2021-03-07T17:18:54.496790Z",
"shell.execute_reply": "2021-03-07T17:18:54.497472Z"
}
},
"outputs": [],
"source": [
"# \\partial_t \\Gamma = -\\hat{g}^{ij} (\\partial_i \\partial_j \\varphi -\n",
"# \\hat{\\Gamma}^k_{ji} \\partial_k \\varphi)\n",
"Gamma_rhs = sp.sympify(0)\n",
"for i in range(DIM):\n",
" for j in range(DIM):\n",
" Gamma_rhs -= ghatUU[i][j]*psi_dDD[i][j]\n",
" for k in range(DIM):\n",
" Gamma_rhs += ghatUU[i][j]*GammahatUDD[k][j][i]*psi_dD[k]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\partial_t \\varphi = -\\Gamma\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.501992Z",
"iopub.status.busy": "2021-03-07T17:18:54.501049Z",
"iopub.status.idle": "2021-03-07T17:18:54.503243Z",
"shell.execute_reply": "2021-03-07T17:18:54.503773Z"
}
},
"outputs": [],
"source": [
"# \\partial_t \\varphi = -\\Gamma\n",
"psi_rhs = -Gamma"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Constraints:\n",
"\n",
"\\begin{align}\n",
"\\mathcal{G} &\\equiv \\Gamma - \\partial_i A^i + \\hat{\\Gamma}^i_{ji} A^j, \\\\\n",
"\\mathcal{C} &\\equiv \\partial_i E^i + \\hat{\\Gamma}^i_{ji} E^j.\n",
"\\end{align}\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.517925Z",
"iopub.status.busy": "2021-03-07T17:18:54.516841Z",
"iopub.status.idle": "2021-03-07T17:18:54.520491Z",
"shell.execute_reply": "2021-03-07T17:18:54.521084Z"
}
},
"outputs": [],
"source": [
"# \\mathcal{G} \\equiv \\Gamma - \\partial_i A^i + \\hat{\\Gamma}^i_{ji} A^j\n",
"G = Gamma\n",
"for i in range(DIM):\n",
" G -= AU_dD[i][i]\n",
" for j in range(DIM):\n",
" G += GammahatUDD[i][j][i]*AU[j]\n",
"\n",
"# E^i\n",
"EU = ixp.zerorank1()\n",
"\n",
"# \\partial_k ( A^i ) --> rank two tensor\n",
"EU_dD = ixp.zerorank2()\n",
"for i in range(DIM):\n",
" EU[i] = eU[i]*ReU[i]\n",
" for j in range(DIM):\n",
" EU_dD[i][j] = eU_dD[i][j]*ReU[i] + eU[i]*ReUdD[i][j]\n",
"\n",
"C = sp.sympify(0)\n",
"for i in range(DIM):\n",
" C += EU_dD[i][i]\n",
" for j in range(DIM):\n",
" C += GammahatUDD[i][j][i]*EU[j]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='cart_transform'></a> \n",
"# Step 3: Convert $A^i$ and $E^i$ to the Cartesian basis \\[Back to [top](#top)\\]\n",
"$$\\label{cart_transform}$$\n",
"\n",
"Here we convert $A^i$ and $E^i$ to the Cartesian basis, to make convergence tests within [Tutorial-Start_to_Finish-Solving_Maxwells_Equations_in_Vacuum-Curvilinear](Tutorial-Start_to_Finish-Solving_Maxwells_Equations_in_Vacuum-Curvilinear.ipynb) easier. Specifically, we will use the coordinate transformation definitions provided by [reference_metric.py](../edit/reference_metric.py) to build the Jacobian:\n",
"\n",
"\\begin{align} \n",
"\\frac{\\partial x_{\\rm Cart}^i}{\\partial x_{\\rm Orig}^j},\n",
"\\end{align}\n",
"\n",
"where $x_{\\rm Cart}^i \\in \\{x,y,z\\}$. We then apply it to $A^i$ and $E^i$ to transform into Cartesian coordinates, via\n",
"\n",
"\\begin{align}\n",
"A^i_{\\rm Cart} = \\frac{\\partial x_{\\rm Cart}^i}{\\partial x_{\\rm Orig}^j} A^j_{\\rm Orig}.\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.580266Z",
"iopub.status.busy": "2021-03-07T17:18:54.560097Z",
"iopub.status.idle": "2021-03-07T17:18:54.668974Z",
"shell.execute_reply": "2021-03-07T17:18:54.668430Z"
}
},
"outputs": [],
"source": [
"def Convert_to_Cartesian_basis(VU):\n",
" # Coordinate transformation from original basis to Cartesian\n",
" rfm.reference_metric()\n",
"\n",
" VU_Cart = ixp.zerorank1()\n",
" Jac_dxCartU_dxOrigD = ixp.zerorank2()\n",
" for i in range(DIM):\n",
" for j in range(DIM):\n",
" Jac_dxCartU_dxOrigD[i][j] = sp.diff(rfm.xx_to_Cart[i], rfm.xx[j])\n",
"\n",
" for i in range(DIM):\n",
" for j in range(DIM):\n",
" VU_Cart[i] += Jac_dxCartU_dxOrigD[i][j]*VU[j]\n",
" return VU_Cart\n",
"\n",
"AU_Cart = Convert_to_Cartesian_basis(AU)\n",
"EU_Cart = Convert_to_Cartesian_basis(EU)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='step4'></a> \n",
"# Step 4: NRPy+ Module Code Validation \\[Back to [top](#top)\\]\n",
"$$\\label{step4}$$\n",
"\n",
"Here, as a code validation check, we verify agreement in the SymPy expressions for the RHSs of Maxwell's equations between\n",
"1. this tutorial and \n",
"2. the NRPy+ [VacuumMaxwell_Flat_Evol_Curvilinear_rescaled](../edit/Maxwell/VacuumMaxwell_Flat_Evol_Curvilinear_rescaled.py) module.\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.677903Z",
"iopub.status.busy": "2021-03-07T17:18:54.677188Z",
"iopub.status.idle": "2021-03-07T17:18:54.870590Z",
"shell.execute_reply": "2021-03-07T17:18:54.871110Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Consistency check between Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling tutorial and NRPy+ module: ALL SHOULD BE ZERO.\n",
"C - mwevol.C = 0\n",
"G - mwevol.G = 0\n",
"psi_rhs - mwevol.psi_rhs = 0\n",
"Gamma_rhs - mwevol.Gamma_rhs = 0\n",
"arhsU[0] - mwevol.arhsU[0] = 0\n",
"erhsU[0] - mwevol.erhsU[0] = 0\n",
"AU_Cart[0] - mwevol.AU_Cart[0] = 0\n",
"EU_Cart[0] - mwevol.EU_Cart[0] = 0\n",
"arhsU[1] - mwevol.arhsU[1] = 0\n",
"erhsU[1] - mwevol.erhsU[1] = 0\n",
"AU_Cart[1] - mwevol.AU_Cart[1] = 0\n",
"EU_Cart[1] - mwevol.EU_Cart[1] = 0\n",
"arhsU[2] - mwevol.arhsU[2] = 0\n",
"erhsU[2] - mwevol.erhsU[2] = 0\n",
"AU_Cart[2] - mwevol.AU_Cart[2] = 0\n",
"EU_Cart[2] - mwevol.EU_Cart[2] = 0\n"
]
}
],
"source": [
"# Reset the list of gridfunctions, as registering a gridfunction\n",
"# twice will spawn an error.\n",
"gri.glb_gridfcs_list = []\n",
"\n",
"# Call the VacuumMaxwellRHSs_rescaled() function from within the\n",
"# Maxwell/VacuumMaxwell_Flat_Evol_Curvilinear_rescaled.py module,\n",
"# which should do exactly the same as the above.\n",
"\n",
"# Set which system to use, which are defined in Maxwell/InitialData.py\n",
"par.initialize_param(par.glb_param(\"char\",\"Maxwell.InitialData\",\"System_to_use\",\"System_II\"))\n",
"\n",
"import Maxwell.VacuumMaxwell_Flat_Evol_Curvilinear_rescaled as mwevol\n",
"mwevol.VacuumMaxwellRHSs_rescaled()\n",
"\n",
"print(\"Consistency check between Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling tutorial and NRPy+ module: ALL SHOULD BE ZERO.\")\n",
"\n",
"print(\"C - mwevol.C = \" + str(C - mwevol.C))\n",
"print(\"G - mwevol.G = \" + str(G - mwevol.G))\n",
"print(\"psi_rhs - mwevol.psi_rhs = \" + str(psi_rhs - mwevol.psi_rhs))\n",
"print(\"Gamma_rhs - mwevol.Gamma_rhs = \" + str(Gamma_rhs - mwevol.Gamma_rhs))\n",
"\n",
"for i in range(DIM):\n",
" print(\"arhsU[\"+str(i)+\"] - mwevol.arhsU[\"+str(i)+\"] = \" + str(arhsU[i] - mwevol.arhsU[i]))\n",
" print(\"erhsU[\"+str(i)+\"] - mwevol.erhsU[\"+str(i)+\"] = \" + str(erhsU[i] - mwevol.erhsU[i]))\n",
" print(\"AU_Cart[\"+str(i)+\"] - mwevol.AU_Cart[\"+str(i)+\"] = \" + str(AU_Cart[i] - mwevol.AU_Cart[i]))\n",
" print(\"EU_Cart[\"+str(i)+\"] - mwevol.EU_Cart[\"+str(i)+\"] = \" + str(EU_Cart[i] - mwevol.EU_Cart[i]))\n"
]
},
{
"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](#top)\\]\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-VacuumMaxwell_Curvilinear_RHS-Rescaling.pdf](Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling.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": 15,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:18:54.875945Z",
"iopub.status.busy": "2021-03-07T17:18:54.874960Z",
"iopub.status.idle": "2021-03-07T17:18:56.391637Z",
"shell.execute_reply": "2021-03-07T17:18:56.392536Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[NbConvertApp] WARNING | pattern 'Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling.ipynb' matched no files\n",
"Created Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling.tex, and compiled\n",
" LaTeX file to PDF file Tutorial-VacuumMaxwell_Curvilinear_RHS-\n",
" Rescaling.pdf\n"
]
}
],
"source": [
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling\")"
]
}
],
"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.10.11"
}
},
"nbformat": 4,
"nbformat_minor": 4
}