{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\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:** Validated \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": [ "\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": [ "\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": [ "\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": [ " \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": [ " \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": [ "\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", "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.8.2" } }, "nbformat": 4, "nbformat_minor": 2 }