{
"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",
"# 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:** <font color='green'><b> Validated </b></font>\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": [
"<a id='toc'></a>\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": [
"<a id='initializenrpy'></a>\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": [
"<a id='step2'></a>\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": [
"<a id='cart_basis'></a>\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": [
"<a id='dst_basis'></a>\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": [
"<a id='step3'></a>\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": [
"<a id='lorentz'></a>\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": [
"<a id='wave_eq'></a>\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": [
"<a id='step4'></a>\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": [
"<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-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
}