{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Initial Data for Solving Maxwell's Equations in Flat Spacetime\n", "\n", "## Authors: Terrence Pierre Jacques, Zachariah Etienne and Ian Ruchlin\n", "\n", "\n", "## This module constructs the initial data for Maxwell's equations as symbolic (SymPy) expressions, for a purely toriodal dipole field, as defined in [Knapp, Walker & Baumgarte (2002)](https://arxiv.org/abs/gr-qc/0201051).\n", "\n", "**Notebook Status:** Validated \n", "\n", "**Validation Notes:** All expressions generated in this module have been validated, against the [Dendro code Maxwell initial data](https://github.com/paralab/Dendro-GR), and have satisfied the contraints given in [Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling](Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling.ipynb), as well as the wave equation for the electric field and the vector potential.\n", "\n", "### NRPy+ Source Code for this module: [Maxwell/InitialData.py](../edit/Maxwell/InitialData.py), [reference_metric.py](../edit/reference_metric.py)\n", "\n", "\n", "[comment]: <> (Introduction: TODO)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Table of Contents\n", "$$\\label{toc}$$ \n", "\n", "1. [Step 1](#initializenrpy): Initialize needed Python/NRPy+ modules and set destination basis\n", "1. [Step 2](#step2): A Purely Toriodal Dipole Field\n", " 1. [Step 2.a](#cart_basis): Converting to Cartesian Basis\n", " 1. [Step 2.b](#dst_basis): Converting to Destination Basis\n", "1. [Step 3](#step3): Checks\n", " 1. [Step 3.a](#lorentz): Lorentz Gauge Condition & Divergence Constraint\n", " 1. [Step 3.b](#wave_eq): Check that $A^i$ satisfies the wave equation\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", "\n", "# Step 1: Initialize needed Python/NRPy+ modules and set destination basis \\[Back to [top](#toc)\\]\n", "\n", "$$\\label{initializenrpy}$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:15.368159Z", "iopub.status.busy": "2021-03-07T17:19:15.367393Z", "iopub.status.idle": "2021-03-07T17:19:15.855737Z", "shell.execute_reply": "2021-03-07T17:19:15.855156Z" } }, "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 sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends\n", "\n", "#Step 0: 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", "# Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,\n", "# SymTP, SinhSymTP\n", "dst_basis = \"Cylindrical\"\n", "\n", "# To help with simplifications, we tell Sympy that\n", "# the coordinate xx0 is radial like (positive real)\n", "radial_like_dst_xx0 = True\n", "\n", "# Set coordinate system to Cartesian\n", "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n", "rfm.reference_metric()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: A Purely Toriodal Dipole Field \\[Back to [top](#toc)\\]\n", "$$\\label{step2}$$\n", "\n", "Having the evolution equations from [Knapp, Walker & Baumgarte (2002)](https://arxiv.org/abs/gr-qc/0201051) written in [Tutorial-VacuumMaxwell_Cartesian_RHSs](Tutorial-VacuumMaxwell_Cartesian_RHSs.ipynb), we must construct the initial data that will then be time evolved. Beginning from the analytic solution to this system of equation given by equation 16 of [Knapp, Walker & Baumgarte (2002)](https://arxiv.org/abs/gr-qc/0201051),\n", "\n", "\\begin{align}\n", "A^{\\hat{\\phi}} &= \\mathcal{A} \\sin \\theta \\left( \\frac{e^{-\\lambda v^2}-e^{-\\lambda u^2}}{r^2} - 2 \\lambda \\frac{ve^{-\\lambda v^2}-ue^{-\\lambda u^2}}{r} \\right), \\\\\n", "\\end{align}\n", "\n", "where $A^{\\hat{\\phi}} = A^{\\phi} r\\sin\\theta$, $\\mathcal{A}$ gives the amplitude, $\\lambda$ describes the size of the wavepacket, $u = t+r$, and $v = t-r$. Other components of the vector potential are $0$. Note that these expressions repesent the exact solution to both systems of equations at any time $t \\geq 0$, at all points on our numerical grid. Thus, to get initial data we set $t=0$.\n", "\n", "For system II, we will also need to set initial data for $\\Gamma$. Since $\\Gamma = -\\partial_t \\psi$ and we have chosen $\\psi(t=0) = 0$, $\\Gamma(t=0) = 0$. \n", "\n", "We may calculate $E^i$ using \n", "\n", "\\begin{align}\n", "E^i = -\\partial_t A^i.\n", "\\end{align}\n", "\n", "\n", "**Inputs for initial data**:\n", "\n", "* amp - $A$\n", "* lam - $\\lambda$\n", "* time - $t$\n", "\n", "Below we define the Cartesian coordinates $x, y$ and $z$. We then define the vector potential $A^i$ in spherical coordinates, but each component is written in terms of Cartesian coordinates. This makes the subsequent basis changes easier." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:15.893846Z", "iopub.status.busy": "2021-03-07T17:19:15.893051Z", "iopub.status.idle": "2021-03-07T17:19:15.895285Z", "shell.execute_reply": "2021-03-07T17:19:15.895751Z" } }, "outputs": [], "source": [ "x = rfm.xx_to_Cart[0]\n", "y = rfm.xx_to_Cart[1]\n", "z = rfm.xx_to_Cart[2]\n", "\n", "# Step 1: Declare free parameters intrinsic to these initial data\n", "# Amplitude\n", "amp = par.Cparameters(\"REAL\",__name__,\"amp\", default_vals=1.0)\n", "# lambda\n", "lam = par.Cparameters(\"REAL\",__name__,\"lam\", default_vals=1.0)\n", "time = par.Cparameters(\"REAL\",__name__,\"time\", default_vals=0.0)\n", "wavespeed = par.Cparameters(\"REAL\",__name__,\"wavespeed\", default_vals=1.0)\n", "\n", "psi_ID = sp.sympify(0)\n", "\n", "Gamma_ID = sp.sympify(0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{align}\n", "A^{\\hat{\\phi}} &= \\mathcal{A} \\sin \\theta \\left( \\frac{e^{-\\lambda v^2}-e^{-\\lambda u^2}}{r^2} - 2 \\lambda \\frac{ve^{-\\lambda v^2}-ue^{-\\lambda u^2}}{r} \\right), \\\\\n", "A^{\\phi} &= \\frac{A^{\\hat{\\phi}}} {r\\sin\\theta}\n", "\\end{align}" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:15.903914Z", "iopub.status.busy": "2021-03-07T17:19:15.903210Z", "iopub.status.idle": "2021-03-07T17:19:15.924590Z", "shell.execute_reply": "2021-03-07T17:19:15.923952Z" } }, "outputs": [], "source": [ "AidU_Sph = ixp.zerorank1()\n", "# Set coordinate transformations:\n", "r = sp.sqrt(x*x + y*y + z*z)\n", "sin_theta = z / r\n", "\n", "u = time + r\n", "v = time - r\n", "e_lam_u = sp.exp(-lam*u**2)\n", "e_lam_v = sp.exp(-lam*v**2)\n", "\n", "# Equation 16 from https://arxiv.org/abs/gr-qc/0201051\n", "AU_phi_hat = (amp*sin_theta)*( ((e_lam_v - e_lam_u)/r**2) - \\\n", " 2*lam*(v*e_lam_v + u*e_lam_u)/r )\n", "\n", "AidU_Sph[2] = AU_phi_hat/(r*sin_theta)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.a: Converting to Cartesian Basis \\[Back to [top](#toc)\\]\n", "$$\\label{cart_basis}$$\n", "\n", "Note that $A^i$ is defined in sperical coordinates, so we must therefore transform to Cartesian coordinates using the [Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant#Example_3:_spherical-Cartesian_transformation). Here 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 Sph}^j},\n", "\\end{align}\n", "\n", "where $x_{\\rm Sph}^j \\in \\{r,\\theta,\\phi\\}$ and $x_{\\rm Cart}^i \\in \\{x,y,z\\}$. We then apply it to $A^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 Sph}^j} A^j_{\\rm Sph}.\n", "\\end{align}" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:15.998886Z", "iopub.status.busy": "2021-03-07T17:19:15.963298Z", "iopub.status.idle": "2021-03-07T17:19:19.451832Z", "shell.execute_reply": "2021-03-07T17:19:19.452719Z" } }, "outputs": [], "source": [ "# Coordinate transformation from spherical to Cartesian\n", "AidU_Cart = ixp.zerorank1()\n", "Jac_dxSphU_dxCartD = ixp.zerorank2()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " Jac_dxSphU_dxCartD[i][j] = sp.diff(rfm.xxSph[i],rfm.xx_to_Cart[j])\n", "\n", "# Jac_dxCartU_dxSphD[i][j] = sp.diff(rfm.xx_to_Cart[i],rfm.xx[j])\n", "Jac_dxCartU_dxSphD,dummy = ixp.generic_matrix_inverter3x3(Jac_dxSphU_dxCartD)\n", "\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " AidU_Cart[i] += Jac_dxCartU_dxSphD[i][j]*AidU_Sph[j]\n", "for i in range(DIM):\n", " AidU_Cart[i] = sp.simplify(AidU_Cart[i])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.b: Converting to Destination Basis \\[Back to [top](#toc)\\]\n", "$$\\label{dst_basis}$$\n", "\n", "Here we prepare to convert $A^i$ from the Cartesian basis to the destination basis. To do so, we first rewrite each component of $A^i$ in terms of the destination coordinates. This is done by first re-labelling the NRPy+ coordinates $xx0, xx1, xx2$ as $cart_{xx0}, cart_{xx1}, cart_{xx2}$. Then, each $cart_{xxi}$ is replaced by its counterpart expression in the destination basis using [reference_metric.py](../edit/reference_metric.py).\n", "\n", "Note that for algebraic simplification, we tell sympy that the coordinate $xx0$ is radial like and thus positive and real (if the destination coordinates are curvilinear)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:19.518257Z", "iopub.status.busy": "2021-03-07T17:19:19.482280Z", "iopub.status.idle": "2021-03-07T17:19:24.784045Z", "shell.execute_reply": "2021-03-07T17:19:24.783406Z" } }, "outputs": [], "source": [ "# rfm is still defined in Cartesian coordinates\n", "cart_xx = ixp.declarerank1(\"cart_xx\")\n", "for i in range(DIM):\n", " for k in range(DIM):\n", " AidU_Cart[i] = AidU_Cart[i].subs(rfm.xx[k], cart_xx[k])\n", "\n", "# Set coordinate system to dst_basis\n", "par.set_parval_from_str(\"reference_metric::CoordSystem\",dst_basis)\n", "rfm.reference_metric()\n", "\n", "for i in range(DIM):\n", " for k in range(DIM):\n", " AidU_Cart[i] = AidU_Cart[i].subs(cart_xx[k], rfm.xx_to_Cart[k])\n", "\n", "if radial_like_dst_xx0:\n", " for j in range(DIM):\n", " AidU_Cart[j] = sp.refine(sp.simplify(AidU_Cart[j]), sp.Q.positive(rfm.xx[0]))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define Jacobians relative to the center of the destination grid, at a point $x^j_{\\rm dst}=$(`xx0,xx1,xx2`)${}_{\\rm dst}$ on the destination grid:\n", "$$\n", "{\\rm Jac\\_dUCart\\_dDdstUD[i][j]} = \\frac{\\partial x^i_{\\rm Cart}}{\\partial x^j_{\\rm dst}},\n", "$$\n", "\n", "via exact differentiation (courtesy SymPy), and the inverse Jacobian\n", "$$\n", "{\\rm Jac\\_dUdst\\_dDCartUD[i][j]} = \\frac{\\partial x^i_{\\rm dst}}{\\partial x^j_{\\rm Cart}},\n", "$$\n", "\n", "using NRPy+'s `generic_matrix_inverter3x3()` function. In terms of these, the transformation of BSSN tensors from Cartesian to the destination grid's `\"reference_metric::CoordSystem\"` coordinates may be written:\n", "\n", "$$\n", "A^i_{\\rm dst} = \\frac{\\partial x^i_{\\rm dst}}{\\partial x^\\ell_{\\rm Cart}} A^\\ell_{\\rm Cart}\n", "$$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:24.843655Z", "iopub.status.busy": "2021-03-07T17:19:24.807219Z", "iopub.status.idle": "2021-03-07T17:19:24.861864Z", "shell.execute_reply": "2021-03-07T17:19:24.861189Z" } }, "outputs": [], "source": [ "# Step 3: Transform BSSN tensors in Cartesian basis to destination grid basis, using center of dest. grid as origin\n", "\n", "# Step 3.a: Next construct Jacobian and inverse Jacobian matrices:\n", "Jac_dUCart_dDrfmUD,Jac_dUrfm_dDCartUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Cartesian()\n", "\n", "# Step 3.b: Convert basis of all BSSN *vectors* from Cartesian to destination basis\n", "AidU = rfm.basis_transform_vectorU_from_Cartesian_to_rfmbasis(Jac_dUrfm_dDCartUD, AidU_Cart)\n", "\n", "# Define electric field --> E^i = -\\partial_t A^i\n", "EidU = ixp.zerorank1()\n", "for j in range(DIM):\n", " EidU[j] = -sp.diff(AidU[j], time)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 3: Checks \\[Back to [top](#toc)\\]\n", "$$\\label{step3}$$\n", "\n", "Here we validate the initial data. Specifically, we check that the constraints from [Tutorial-VacuumMaxwell_formulation_Curvilinear](Tutorial-VacuumMaxwell_formulation_Curvilinear.ipynb) are satisfied;\n", "\n", "\\begin{align}\n", "\\mathcal{G} &\\equiv \\Gamma - \\partial_i A^i - \\hat{\\Gamma}^i_{ji} A^j &= 0, \\quad &\\text{Lorenz gauge condition} \\\\\n", "\\mathcal{C} &\\equiv \\partial_i E^i + \\hat{\\Gamma}^i_{ji} E^j &= 0, \\quad &\\text{Divergence Constraint}.\n", "\\end{align}\n", "\n", "Note that the above simply to their usual forms in Cartesian coordinates.\n", "\n", "Finally, we check that $A^i$ satisfies the covariant wave equation,\n", "\n", "\\begin{align}\n", "\\partial_t^2 A^i - \\hat{\\gamma}^{jk} \\hat{\\nabla}_j \\hat{\\nabla}_k A^i = 0,\n", "\\end{align}\n", "\n", "where $\\hat{\\nabla}_j$ is the covariant derivative associated with the spatial metric $\\hat{\\gamma}_{jk}$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:24.922180Z", "iopub.status.busy": "2021-03-07T17:19:24.886444Z", "iopub.status.idle": "2021-03-07T17:19:30.717132Z", "shell.execute_reply": "2021-03-07T17:19:30.717629Z" } }, "outputs": [], "source": [ "AidU_dD = ixp.zerorank2()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " AidU_dD[i][j] += sp.diff(AidU[i], rfm.xx[j])\n", "\n", "AidU_dDD = ixp.zerorank3()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " for k in range(DIM):\n", " AidU_dDD[i][j][k] += sp.diff(AidU[i], rfm.xx[j], rfm.xx[k])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.a: Lorentz Gauge Condition & Divergence Constraint \\[Back to [top](#toc)\\]\n", "$$\\label{lorentz}$$\n", "\n", "\\begin{align}\n", "\\mathcal{G} &\\equiv \\Gamma - \\partial_i A^i - \\hat{\\Gamma}^i_{ji} A^j &= 0, \\quad &\\text{Lorenz gauge condition} \\\\\n", "\\mathcal{C} &\\equiv \\partial_i E^i + \\hat{\\Gamma}^i_{ji} E^j &= 0, \\quad &\\text{Divergence Constraint}\n", "\\end{align}" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:30.797840Z", "iopub.status.busy": "2021-03-07T17:19:30.761781Z", "iopub.status.idle": "2021-03-07T17:19:30.827150Z", "shell.execute_reply": "2021-03-07T17:19:30.827909Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "G should evaluate to zero: 0 \n", "\n", "C should evaluate to zero: 0\n" ] } ], "source": [ "# \\mathcal{G} \\equiv \\Gamma - \\partial_i A^i - \\hat{\\Gamma}^i_{ji} A^j\n", "G = Gamma_ID\n", "for i in range(DIM):\n", " G -= AidU_dD[i][i]\n", " for j in range(DIM):\n", " G -= rfm.GammahatUDD[i][j][i]*AidU[j]\n", "\n", "print('G should evaluate to zero:', sp.simplify(G), '\\n')\n", "\n", "# \\mathcal{C} \\equiv \\partial_i E^i + \\hat{\\Gamma}^i_{ji} E^j\n", "C = sp.sympify(0)\n", "for i in range(DIM):\n", " C += sp.diff(EidU[i], rfm.xx[i], 1)\n", " for j in range(DIM):\n", " C += rfm.GammahatUDD[i][j][i]*EidU[j]\n", "print('C should evaluate to zero:', sp.simplify(C))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.b: Check that $A^i$ satisfies the wave equation \\[Back to [top](#toc)\\]\n", "$$\\label{wave_eq}$$\n", "\n", "Based on the definition of covariant derivative, we have\n", "$$\n", "\\hat{\\nabla}_{k} A^{i} = A^i_{,k} + \\hat{\\Gamma}^i_{mk} A^m\n", "$$\n", "\n", "Since $\\hat{\\nabla}_{k} A^{i}$ is a tensor, the covariant derivative of this will have the same indexing as a tensor $T_k^i$:\n", "\n", "$$\n", "\\hat{\\nabla}_{j} T^i_k = T^i_{k,j} + \\hat{\\Gamma}^i_{dj} T^d_k - \\hat{\\Gamma}^d_{kj} T^i_d.\n", "$$\n", "\n", "Therefore,\n", "\\begin{align}\n", "\\hat{\\nabla}_{j} \\left(\\hat{\\nabla}_{k} A^{i}\\right) &= \\left(A^i_{,k} + \\hat{\\Gamma}^i_{mk} A^m\\right)_{,j} + \\hat{\\Gamma}^i_{dj} \\left(A^d_{,k} + \\hat{\\Gamma}^d_{mk} A^m\\right) - \\hat{\\Gamma}^d_{kj} \\left(A^i_{,d} + \\hat{\\Gamma}^i_{md} A^m\\right) \\\\\n", "&= A^i_{,kj} + \\hat{\\Gamma}^i_{mk,j} A^m + \\hat{\\Gamma}^i_{mk} A^m_{,j} + \\hat{\\Gamma}^i_{dj}A^d_{,k} + \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} A^m - \\hat{\\Gamma}^d_{kj} A^i_{,d} - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} A^m \\\\\n", "&= {\\underbrace {\\textstyle A^i_{,kj}}_{\\text{Term 1}}}+\n", "{\\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 2}}} +\n", "{\\underbrace {\\textstyle \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} A^m - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} A^m}_{\\text{Term 3}}}.\n", "\\end{align}\n", "\n", "Thus \n", "$$\n", "\\hat{\\gamma}^{jk} \\hat{\\nabla}_{j} \\left(\\hat{\\nabla}_{k} A^{i}\\right) = \\hat{\\gamma}^{jk} \\left(\\text{Term 1} + \\text{Term 2} + \\text{Term 3}\\right).\n", "$$\n", "\n", "We use the above to confirm\n", "\n", "\\begin{align}\n", "\\partial_t^2 A^i - \\hat{\\gamma}^{jk} \\hat{\\nabla}_j \\hat{\\nabla}_k A^i = 0,\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\text{Term 1} = A^i_{,kj}\n", "$$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:30.835228Z", "iopub.status.busy": "2021-03-07T17:19:30.834164Z", "iopub.status.idle": "2021-03-07T17:19:30.837527Z", "shell.execute_reply": "2021-03-07T17:19:30.836645Z" } }, "outputs": [], "source": [ "# Term 1: A^i_{,kj}\n", "Term1UDD = ixp.zerorank3()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " for k in range(DIM):\n", " Term1UDD[i][j][k] += AidU_dDD[i][k][j]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\text{Term 2} = \\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": 10, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:30.848503Z", "iopub.status.busy": "2021-03-07T17:19:30.847450Z", "iopub.status.idle": "2021-03-07T17:19:30.850228Z", "shell.execute_reply": "2021-03-07T17:19:30.850771Z" } }, "outputs": [], "source": [ "# Term 2: \\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", "Term2UDD = 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", " Term2UDD[i][j][k] += rfm.GammahatUDDdD[i][m][k][j]*AidU[m] \\\n", " + rfm.GammahatUDD[i][m][k]*AidU_dD[m][j] \\\n", " + rfm.GammahatUDD[i][m][j]*AidU_dD[m][k] \\\n", " - rfm.GammahatUDD[m][k][j]*AidU_dD[i][m]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\text{Term 3} = \\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": 11, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:30.880036Z", "iopub.status.busy": "2021-03-07T17:19:30.872938Z", "iopub.status.idle": "2021-03-07T17:19:30.887747Z", "shell.execute_reply": "2021-03-07T17:19:30.888526Z" } }, "outputs": [], "source": [ "# Term 3: \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} A^m -\n", "# \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} A^m\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", " for d in range(DIM):\n", " Term3UDD[i][j][k] += ( rfm.GammahatUDD[i][d][j]*rfm.GammahatUDD[d][m][k] \\\n", " -rfm.GammahatUDD[d][k][j]*rfm.GammahatUDD[i][m][d])*AidU[m]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\hat{\\gamma}^{jk} \\hat{\\nabla}_{j} \\left(\\hat{\\nabla}_{k} A^{i}\\right) = \\hat{\\gamma}^{jk} \\left(\\text{Term 1} + \\text{Term 2} + \\text{Term 3}\\right),\n", "$$\n", "\n", "$$\n", "\\partial_t^2 A^i - \\hat{\\gamma}^{jk} \\hat{\\nabla}_j \\hat{\\nabla}_k A^i = 0\n", "$$" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:19:30.902385Z", "iopub.status.busy": "2021-03-07T17:19:30.901614Z", "iopub.status.idle": "2021-03-07T17:20:05.258047Z", "shell.execute_reply": "2021-03-07T17:20:05.257528Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0th component of A-field equation. Should be zero (takes a bit, please be patient): \n", " 0\n", "1th component of A-field equation. Should be zero (takes a bit, please be patient): \n", " 0\n", "2th component of A-field equation. Should be zero (takes a bit, please be patient): \n", " 0\n" ] } ], "source": [ "# A^i_{,kj} + \\hat{\\Gamma}^i_{mk,j} A^m + \\hat{\\Gamma}^i_{mk} A^m_{,j} +\n", "# \\hat{\\Gamma}^i_{dj}A^d_{,k} + \\hat{\\Gamma}^i_{dj}\\hat{\\Gamma}^d_{mk} A^m -\n", "# \\hat{\\Gamma}^d_{kj} A^i_{,d} - \\hat{\\Gamma}^d_{kj} \\hat{\\Gamma}^i_{md} A^m\n", "\n", "Difference = ixp.zerorank1()\n", "for i in range(DIM):\n", " Difference[i] = sp.diff(AidU[i], time, 2)\n", " for j in range(DIM):\n", " for k in range(DIM):\n", " Difference[i] += -rfm.ghatUU[k][j]*(Term1UDD[i][j][k] + Term2UDD[i][j][k] + Term3UDD[i][j][k])\n", "\n", "for i in range(DIM):\n", " print(str(i)+\"th component of A-field equation. Should be zero (takes a bit, please be patient): \")\n", " print(\" \"+str(sp.simplify(Difference[i])))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 4: NRPy+ Module Code Validation \\[Back to [top](#toc)\\]\n", "$$\\label{step4}$$\n", "\n", "Here, as a code validation check, we verify agreement in the SymPy expressions for the initial data we intend to use between\n", "1. this tutorial and \n", "2. the NRPy+ [InitialData](../edit/Maxwell/InitialData.py) module.\n", "Since the initial data is identical between the two systems for $E^i, A^i$, and $\\psi$, we also set and validate initial data for $\\Gamma$." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:20:05.266526Z", "iopub.status.busy": "2021-03-07T17:20:05.265906Z", "iopub.status.idle": "2021-03-07T17:20:16.008366Z", "shell.execute_reply": "2021-03-07T17:20:16.007753Z" }, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Currently using System_II initial data\n", "Consistency check between this tutorial and NRPy+ module: ALL SHOULD BE ZERO.\n", "psi_ID - mwid.psi_ID = 0\n", "Gamma_ID - mwid.Gamma_ID = 0\n", "AidU[0] - mwid.AidU[0] = 0\n", "EidU[0] - mwid.EidU[0] = 0\n", "AidU[1] - mwid.AidU[1] = 0\n", "EidU[1] - mwid.EidU[1] = 0\n", "AidU[2] - mwid.AidU[2] = 0\n", "EidU[2] - mwid.EidU[2] = 0\n" ] } ], "source": [ "import Maxwell.InitialData as mwid\n", "par.set_parval_from_str(\"Maxwell.InitialData::System_to_use\",\"System_II\")\n", "\n", "mwid.InitialData()\n", "\n", "# Again, to help sympy with simplifications\n", "if radial_like_dst_xx0:\n", " for j in range(DIM):\n", " mwid.AidU[j] = sp.refine(sp.simplify(mwid.AidU[j]), sp.Q.positive(rfm.xx[0]))\n", " mwid.EidU[j] = sp.refine(sp.simplify(mwid.EidU[j]), sp.Q.positive(rfm.xx[0]))\n", "\n", "print(\"Consistency check between this tutorial and NRPy+ module: ALL SHOULD BE ZERO.\")\n", "\n", "print(\"psi_ID - mwid.psi_ID = \" + str(sp.simplify(psi_ID) - mwid.psi_ID))\n", "print(\"Gamma_ID - mwid.Gamma_ID = \" + str(Gamma_ID - mwid.Gamma_ID))\n", "\n", "for i in range(DIM):\n", " print(\"AidU[\"+str(i)+\"] - mwid.AidU[\"+str(i)+\"] = \" + str(sp.simplify(AidU[i] - mwid.AidU[i])))\n", " print(\"EidU[\"+str(i)+\"] - mwid.EidU[\"+str(i)+\"] = \" + str(sp.simplify(EidU[i] - mwid.EidU[i])))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\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-VacuumMaxwell_InitialData.pdf](Tutorial-VacuumMaxwell_InitialData.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": 14, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:20:16.013027Z", "iopub.status.busy": "2021-03-07T17:20:16.012303Z", "iopub.status.idle": "2021-03-07T17:20:19.711800Z", "shell.execute_reply": "2021-03-07T17:20:19.711134Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Created Tutorial-VacuumMaxwell_InitialData.tex, and compiled LaTeX file to\n", " PDF file Tutorial-VacuumMaxwell_InitialData.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_InitialData\")" ] } ], "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 }