{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Two Formulations of Maxwell's equations in Cartesian Coordinates\n", "\n", "## Authors: Patrick Nelson & Zach Etienne\n", "### Formatting improvements courtesy Brandon Clark\n", "\n", "[comment]: <> (Abstract: TODO)\n", "\n", "**Notebook Status:** Self-Validated \n", "\n", "**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented below, [here](#code_validation_sys1) and [here](#code_validation). **Additional validation tests may have been performed, but are as yet, undocumented. (TODO)**\n", "\n", "### NRPy+ Source Code for this module: \n", "* [Maxwell/MaxwellCartesian_Evol.py](../edit/Maxwell/MaxwellCartesian_Evol.py)\n", "* [Maxwell/MaxwellCartesian_ID.py](../edit/Maxwell/MaxwellCartesian_ID.py)\n", "\n", "## Introduction:\n", "\n", "This tutorial will draw on previous work done by Ian Ruchlin on [Maxwell's equations in Curvilinear Coordinates](Tutorial-MaxwellCurvilinear.ipynb), which itself drew on the two formulations described in [Illustrating Stability Properties of Numerical Relativity in Electrodynamics](https://arxiv.org/abs/gr-qc/0201051). This will be done to aid construction of an Einstein Toolkit thorn, which will itself be built in the [Tutorial-ETK_thorn-Maxwell](Tutorial-ETK_thorn-Maxwell.ipynb) tutorial notebook.\n", "\n", "The construction of our equations here will be nearly identical; however, by assuming Cartesian coordinates, we are able to make several simplifications and eliminate the need for [reference_metric.py](../edit/reference_metric.py) as a dependency.\n", "\n", "We will begin with their System I. While the Curvilinear version of this code assumed flat spacetime, we will be constructing the equations in a general spacetime. This allows us to simply replace the reference metric \"hatted\" quantities with their general counterparts." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Table of Contents:\n", "$$\\label{toc}$$\n", "\n", "This notebook is organized as follows\n", "\n", "1. [Step 1](#initializenrpy): Import core NRPy+ modules and set NRPy+ parameters\n", "1. [Step 2](#laplacian): Constructing the Covariant Laplacian \n", "1. [Step 3](#violation): Measuring the Constraint violation\n", "1. [Step 4](#code_validation_sys1): Code Validation against `Maxwell.MaxwellCartesian_Evol` NRPy+ Module (System I)\n", "1. [Step 5](#code_validation_sys2): Code Validation against `Maxwell.MaxwellCartesian_Evol` NRPy+ Module (System II)\n", "1. [Step 6](#id): Constructing the Initial Data\n", "1. [Step 7](#code_validation): Code Validation against `Maxwell.MaxwellCartesian_ID` NRPy+ Module\n", "1. [Step 8](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 1: Import core NRPy+ modules and set NRPy+ parameters \\[Back to [top](#toc)\\]\n", "$$\\label{initializenrpy}$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:42.379600Z", "iopub.status.busy": "2021-03-07T17:14:42.378579Z", "iopub.status.idle": "2021-03-07T17:14:42.857709Z", "shell.execute_reply": "2021-03-07T17:14:42.858202Z" } }, "outputs": [], "source": [ "#Step P1: Import needed Python modules\n", "import NRPy_param_funcs as par # NRPy+: Parameter interface\n", "import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends\n", "import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n", "import finite_difference as fin # NRPy+: Finite difference C code generation module\n", "import grid as gri # NRPy+: Functions having to do with numerical grids\n", "\n", "#Step P2: Name this Python module\n", "thismodule = \"MaxwellCartesian\"\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", "# Step 1: Set the finite differencing order to 4.\n", "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", 4)\n", "\n", "# Step 2: Register gridfunctions that are needed as input.\n", "psi = gri.register_gridfunctions(\"EVOL\", [\"psi\"])\n", "\n", "# Step 3a: Declare the rank-1 indexed expressions E_{i}, A_{i},\n", "# and \\partial_{i} \\psi. Derivative variables like these\n", "# must have an underscore in them, so the finite\n", "# difference module can parse the variable name properly.\n", "ED = ixp.register_gridfunctions_for_single_rank1(\"EVOL\", \"ED\")\n", "AD = ixp.register_gridfunctions_for_single_rank1(\"EVOL\", \"AD\")\n", "psi_dD = ixp.declarerank1(\"psi_dD\")\n", "x,y,z = gri.register_gridfunctions(\"AUX\",[\"x\",\"y\",\"z\"])\n", "\n", "## Step 3b: Declare the conformal metric tensor and its first\n", "# derivative. These are needed to find the Christoffel\n", "# symbols, which we need for covariant derivatives.\n", "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUX\",\"gammaDD\", \"sym01\") # The AUX or EVOL designation is *not*\n", " # used in diagnostic modules.\n", "gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\")\n", "gammaDD_dDD = ixp.declarerank4(\"gammaDD_dDD\",\"sym01_sym23\")\n", "\n", "gammaUU, detgamma = ixp.symm_matrix_inverter3x3(gammaDD)\n", "gammaUU_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\")\n", "\n", "# Define the Christoffel symbols\n", "GammaUDD = ixp.zerorank3(DIM)\n", "for i in range(DIM):\n", " for k in range(DIM):\n", " for l in range(DIM):\n", " for m in range(DIM):\n", " GammaUDD[i][k][l] += (sp.Rational(1,2))*gammaUU[i][m]*\\\n", " (gammaDD_dD[m][k][l] + gammaDD_dD[m][l][k] - gammaDD_dD[k][l][m])\n", "\n", "# Step 3b: Declare the rank-2 indexed expression \\partial_{j} A_{i},\n", "# which is not symmetric in its indices.\n", "# Derivative variables like these must have an underscore\n", "# in them, so the finite difference module can parse the\n", "# variable name properly.\n", "AD_dD = ixp.declarerank2(\"AD_dD\", \"nosym\")\n", "\n", "# Step 3c: Declare the rank-3 indexed expression \\partial_{jk} A_{i},\n", "# which is symmetric in the two {jk} indices.\n", "AD_dDD = ixp.declarerank3(\"AD_dDD\", \"sym12\")\n", "\n", "# Step 4: Calculate first and second covariant derivatives, and the\n", "# necessary contractions.\n", "# First covariant derivative\n", "# D_{j} A_{i} = A_{i,j} - \\Gamma^{k}_{ij} A_{k}\n", "AD_dcovD = ixp.zerorank2()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " AD_dcovD[i][j] = AD_dD[i][j]\n", " for k in range(DIM):\n", " AD_dcovD[i][j] -= GammaUDD[k][i][j] * AD[k]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: Constructing the Covariant Laplacian \\[Back to [top](#toc)\\]\n", "$$\\label{laplacian}$$\n", "\n", "One particular difficulty we will encounter here is taking the covariant Laplacian of a vector. This will here take the form of $D_j D^j A_i$. We will start with the outer derivative after lowering the index on the second operator. So, we see that \n", "\\begin{align}\n", "D_j D^j A_i &= D_j (\\gamma^{jk} D_k A_i) \\\\\n", "&= (D_k A_i) D_j \\gamma^{jk} + \\gamma^{jk} D_j D_k A_i \\\\\n", "&= \\gamma^{jk} [\\partial_j (D_k A_i) - \\Gamma^l_{ij} D_k A_l - \\Gamma^l_{jk} D_l A_i],\n", "\\end{align}\n", "dropping the first term from the second line because $D_j \\gamma^{jk} = 0$ Next, we will again apply the covariant derivative to $A_i$. First, however, we should consider that \n", "\\begin{align}\n", "D_k A_i &= \\partial_k A_i - \\Gamma^l_{ik} A_l \\\\\n", "& = \\partial_k A_i - \\gamma^{lm} \\Gamma_{mik} A_l \\\\\n", "& = \\partial_k A_i - \\Gamma_{mik} A^m, \\\\\n", "\\end{align}\n", "where $\\Gamma_{ljk} = \\frac{1}{2} (\\partial_k \\gamma_{lj} + \\partial_j \\gamma_{kl} - \\partial_l \\gamma_{jk})$ is the Christoffel symbol of the first kind. Note how we were able to use the raising operator to switch the height of two indices; we will use this in upcoming steps. Thus, our expression becomes\n", "\\begin{align}\n", "D_j D^j A_i &= \\gamma^{jk} [\\partial_j \\partial_k A_i - \\partial_j (\\Gamma_{lik} A^l) - \\Gamma^l_{ij} \\partial_k A_l + \\Gamma^m_{ij} \\Gamma_{lmk} A^l - \\Gamma^l_{jk} \\partial_l A_i + \\Gamma^m_{jk} \\Gamma_{lim} A^l] \\\\\n", "&= \\gamma^{jk} [\\partial_j \\partial_k A_i - \\underbrace{\\partial_j (\\Gamma_{lik} A^l)}_{\\text{sub-term}} - \\Gamma^l_{ij} \\partial_k A_l + \\Gamma^m_{ij} \\Gamma^l_{mk} A_l - \\Gamma^l_{jk} \\partial_l A_i + \\Gamma^m_{jk} \\Gamma^l_{im} A_l].\n", "\\end{align}\n", "Let's focus on the underbraced sub-term for a moment. Expanding this using the product rule and the definition of the Christoffel symbol, \n", "\\begin{align}\n", "\\partial_j (\\Gamma_{lik} A^l) &= A^l \\partial_j \\Gamma_{lik} + \\Gamma_{lik} \\partial_j A^l \\\\\n", "&= A^l \\partial_j (\\partial_k \\gamma_{li} + \\partial_i \\gamma_{kl} - \\partial_l \\gamma_{ik}) + \\Gamma_{lik} \\partial_j (\\gamma^{lm} A_m) \\\\\n", "&= A^l (\\gamma_{li,kj} + \\gamma_{kl,ij} - \\gamma_{ik,lj}) + \\Gamma_{lik} (\\gamma^{lm} A_{m,j} + A_m \\gamma^{lm}{}_{,j}), \\\\\n", "\\end{align}\n", "where commas in subscripts denote partial derivatives.\n", "\n", "So, the Laplacian becomes \n", "\\begin{align}\n", "D_j D^j A_i &= \\gamma^{jk} [A_{i,jk} - \n", "\\underbrace{A^l (\\gamma_{li,kj} + \\gamma_{kl,ij} - \\gamma_{ik,lj})}_{\\text{Term 1}} + \n", "\\underbrace{\\Gamma_{lik} (\\gamma^{lm} A_{m,j} + A_m \\gamma^{lm}{}_{,j})}_{\\text{Term 2}} - \n", "\\underbrace{(\\Gamma^l_{ij} A_{l,k} + \\Gamma^l_{jk} A_{i,l})}_{\\text{Term 3}} + \n", "\\underbrace{(\\Gamma^m_{ij} \\Gamma^l_{mk} A_l + \\Gamma ^m_{jk} \\Gamma^l_{im} A_l)}_{\\text{Term 4}}]; \\\\\n", "\\end{align}\n", "we will now begin to contruct these terms individually." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:42.922852Z", "iopub.status.busy": "2021-03-07T17:14:42.886964Z", "iopub.status.idle": "2021-03-07T17:14:43.251201Z", "shell.execute_reply": "2021-03-07T17:14:43.251734Z" } }, "outputs": [], "source": [ "# First, we must construct the lowered Christoffel symbols:\n", "# \\Gamma_{ijk} = \\gamma_{il} \\Gamma^l_{jk}\n", "# And raise the index on A:\n", "# A^j = \\gamma^{ij} A_i\n", "GammaDDD = ixp.zerorank3()\n", "AU = ixp.zerorank1()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " AU[j] += gammaUU[i][j] * AD[i]\n", " for k in range(DIM):\n", " for l in range(DIM):\n", " GammaDDD[i][j][k] += gammaDD[i][l] * GammaUDD[l][j][k]\n", "\n", "# Covariant second derivative (the bracketed terms):\n", "# D_j D^j A_i = \\gamma^{jk} [A_{i,jk} - A^l (\\gamma_{li,kj} + \\gamma_{kl,ij} - \\gamma_{ik,lj})\n", "# + \\Gamma_{lik} (\\gamma^{lm} A_{m,j} + A_m \\gamma^{lm}{}_{,j})\n", "# - (\\Gamma^l_{ij} A_{l,k} + \\Gamma^l_{jk} A_{i,l})\n", "# + (\\Gamma^m_{ij} \\Gamma^l_{mk} A_l + \\Gamma ^m_{jk} \\Gamma^l_{im} A_l)\n", "AD_dcovDD = ixp.zerorank3()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " for k in range(DIM):\n", " AD_dcovDD[i][j][k] = AD_dDD[i][j][k]\n", " for l in range(DIM):\n", " # Terms 1 and 3\n", " AD_dcovDD[i][j][k] -= AU[l] * (gammaDD_dDD[l][i][k][j] + gammaDD_dDD[k][l][i][j] - \\\n", " gammaDD_dDD[i][k][l][j]) \\\n", " + GammaUDD[l][i][j] * AD_dD[l][k] + GammaUDD[l][j][k] * AD_dD[i][l]\n", " for m in range(DIM):\n", " # Terms 2 and 4\n", " AD_dcovDD[i][j][k] += GammaDDD[l][i][k] * (gammaUU[l][m] * AD_dD[m][j] + AD[m] * gammaUU_dD[l][m][j]) \\\n", " + GammaUDD[m][i][j] * GammaUDD[l][m][k] * AD[l] \\\n", " + GammaUDD[m][j][k] * GammaUDD[l][i][m] * AD[l]\n", "\n", "# Covariant divergence\n", "# D_{i} A^{i} = \\gamma^{ij} D_{j} A_{i}\n", "DivA = 0\n", "# Gradient of covariant divergence\n", "# DivA_dD_{i} = \\gamma^{jk} A_{k;\\hat{j}\\hat{i}}\n", "DivA_dD = ixp.zerorank1()\n", "# Covariant Laplacian\n", "# LapAD_{i} = \\gamma^{jk} A_{i;\\hat{j}\\hat{k}}\n", "LapAD = ixp.zerorank1()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " DivA += gammaUU[i][j] * AD_dcovD[i][j]\n", " for k in range(DIM):\n", " DivA_dD[i] += gammaUU[j][k] * AD_dcovDD[k][j][i]\n", " LapAD[i] += gammaUU[j][k] * AD_dcovDD[i][j][k]\n", "\n", "# Step 5: Define right-hand sides for the evolution.\n", "ArhsD = ixp.zerorank1()\n", "ErhsD = ixp.zerorank1()\n", "for i in range(DIM):\n", " ArhsD[i] = -ED[i] - psi_dD[i]\n", " ErhsD[i] = -LapAD[i] + DivA_dD[i]\n", "psi_rhs = -DivA\n", "\n", "# Step 6: Generate C code for System I Maxwell's evolution equations,\n", "# print output to the screen (standard out, or stdout).\n", "#lhrh_list = []\n", "#for i in range(DIM):\n", "# lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"AD\" + str(i)), rhs=ArhsD[i]))\n", "# lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"ED\" + str(i)), rhs=ErhsD[i]))\n", "#lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"psi\"), rhs=psi_rhs))\n", "\n", "#fin.FD_outputC(\"stdout\", lhrh_list)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 3: Measuring the Constraint violation \\[Back to [top](#toc)\\]\n", "$$\\label{violation}$$\n", "\n", "To evaluate our results, we will need to measure how much our simulation differs from what should be physically possible. We will do this with the constraint equation $D_i E^i = 4 \\pi \\rho_e$; specifically, we will measure the constraint violation \n", "\\begin{align}\n", "\\mathcal{C} &= D_i E^i - 4 \\pi \\rho_e \\\\\n", "&= D_i (\\gamma^{ij} E_j) - 4 \\pi \\rho_e \\\\\n", "&= \\gamma^{ij} D_i E_j + E_j D_i \\gamma^{ij} - 4 \\pi \\rho_e \\\\\n", "&= \\gamma^{ij} D_i E_j,\n", "\\end{align}\n", "since the covariant derivative of the metric tensor is $0$ and $\\rho_e=0$ in free space. So, $\\mathcal{C} = \\gamma^{ij} (E_{j,i} - \\Gamma^b_{ij} E_b)$, which will be valid for both systems." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:43.284925Z", "iopub.status.busy": "2021-03-07T17:14:43.284266Z", "iopub.status.idle": "2021-03-07T17:14:43.287135Z", "shell.execute_reply": "2021-03-07T17:14:43.286601Z" } }, "outputs": [], "source": [ "ED_dD = ixp.declarerank2(\"ED_dD\",\"nosym\")\n", "Cviolation = gri.register_gridfunctions(\"AUX\", [\"Cviolation\"])\n", "Cviolation = sp.sympify(0)\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " Cviolation += gammaUU[i][j] * ED_dD[j][i]\n", " for b in range(DIM):\n", " Cviolation -= gammaUU[i][j] * GammaUDD[b][i][j] * ED[b]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 4: Code Validation against `Maxwell.MaxwellCartesian_Evol` NRPy+ Module (System I) \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation_sys1}$$\n", "\n", "Here, as a code validation check, we verify agreement in the SymPy expressions for the RHSs of Maxwell's equations (in System I) between\n", "\n", "1. this tutorial and \n", "2. the NRPy+ [Maxwell.MaxwellCartesian](../edit/Maxwell/MaxwellCartesian_Evol.py) module.\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:43.294584Z", "iopub.status.busy": "2021-03-07T17:14:43.293918Z", "iopub.status.idle": "2021-03-07T17:14:44.094691Z", "shell.execute_reply": "2021-03-07T17:14:44.095194Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Warning: System I is less stable!\n", "Consistency check between MaxwellCartesian tutorial and NRPy+ module: ALL SHOULD BE ZERO.\n", "psi_rhs - mwevol.psi_rhs = 0\n", "ArhsD[0] - mwevol.ArhsD[0] = 0\n", "ErhsD[0] - mwevol.ErhsD[0] = 0\n", "ArhsD[1] - mwevol.ArhsD[1] = 0\n", "ErhsD[1] - mwevol.ErhsD[1] = 0\n", "ArhsD[2] - mwevol.ArhsD[2] = 0\n", "ErhsD[2] - mwevol.ErhsD[2] = 0\n", "Cviolation - mwevol.Cviolation = 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", "# Step 18: Call the MaxwellCartesian_Evol() function from within the\n", "# Maxwell/MaxwellCartesian_Evol.py module,\n", "# which should do exactly the same as in Steps 1-16 above.\n", "import Maxwell.MaxwellCartesian_Evol as mwevol\n", "par.set_parval_from_str(\"System_to_use\",\"System_I\")\n", "mwevol.MaxwellCartesian_Evol()\n", "\n", "print(\"Consistency check between MaxwellCartesian tutorial and NRPy+ module: ALL SHOULD BE ZERO.\")\n", "\n", "print(\"psi_rhs - mwevol.psi_rhs = \" + str(psi_rhs - mwevol.psi_rhs))\n", "for i in range(DIM):\n", "\n", " print(\"ArhsD[\"+str(i)+\"] - mwevol.ArhsD[\"+str(i)+\"] = \" + str(ArhsD[i] - mwevol.ArhsD[i]))\n", " print(\"ErhsD[\"+str(i)+\"] - mwevol.ErhsD[\"+str(i)+\"] = \" + str(ErhsD[i] - mwevol.ErhsD[i]))\n", "\n", "print(\"Cviolation - mwevol.Cviolation = \" + str(Cviolation - mwevol.Cviolation))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now build the equations for System II." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:44.125326Z", "iopub.status.busy": "2021-03-07T17:14:44.120073Z", "iopub.status.idle": "2021-03-07T17:14:44.151640Z", "shell.execute_reply": "2021-03-07T17:14:44.151086Z" } }, "outputs": [], "source": [ "# We inherit here all of the definitions from System I, above\n", "\n", "# Step 7a: Register the scalar auxiliary variable \\Gamma\n", "Gamma = gri.register_gridfunctions(\"EVOL\", [\"Gamma\"])\n", "\n", "# Step 7b: Declare the ordinary gradient \\partial_{i} \\Gamma\n", "Gamma_dD = ixp.declarerank1(\"Gamma_dD\")\n", "\n", "# Step 8a: Construct the second covariant derivative of the scalar \\psi\n", "# \\psi_{;\\hat{i}\\hat{j}} = \\psi_{,i;\\hat{j}}\n", "# = \\psi_{,ij} - \\Gamma^{k}_{ij} \\psi_{,k}\n", "psi_dDD = ixp.declarerank2(\"psi_dDD\", \"sym01\")\n", "psi_dcovDD = ixp.zerorank2()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " psi_dcovDD[i][j] = psi_dDD[i][j]\n", " for k in range(DIM):\n", " psi_dcovDD[i][j] += - GammaUDD[k][i][j] * psi_dD[k]\n", "\n", "# Step 8b: Construct the covariant Laplacian of \\psi\n", "# Lappsi = ghat^{ij} D_{j} D_{i} \\psi\n", "Lappsi = 0\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " Lappsi += gammaUU[i][j] * psi_dcovDD[i][j]\n", "\n", "# Step 9: Define right-hand sides for the evolution.\n", "ArhsD = ixp.zerorank1()\n", "ErhsD = ixp.zerorank1()\n", "for i in range(DIM):\n", " ArhsD[i] = -ED[i] - psi_dD[i]\n", " ErhsD[i] = -LapAD[i] + Gamma_dD[i]\n", "psi_rhs = -Gamma\n", "Gamma_rhs = -Lappsi\n", "\n", "# Step 10: Generate C code for System II Maxwell's evolution equations,\n", "# print output to the screen (standard out, or stdout).\n", "#lhrh_list = []\n", "#for i in range(DIM):\n", "# lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"AD\" + str(i)), rhs=ArhsD[i]))\n", "# lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"ED\" + str(i)), rhs=ErhsD[i]))\n", "#lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"psi\"), rhs=psi_rhs))\n", "#lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"Gamma\"), rhs=Gamma_rhs))\n", "\n", "#fin.FD_outputC(\"stdout\", lhrh_list)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 5: Code Validation against `Maxwell.MaxwellCartesian_Evol` NRPy+ Module (System II) \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation_sys2}$$\n", "\n", "Here, as a code validation check, we verify agreement in the SymPy expressions for the RHSs of Maxwell's equations (in System II) between\n", "\n", "1. this tutorial and \n", "2. the NRPy+ [Maxwell.MaxwellCartesian](../edit/Maxwell/MaxwellCartesian_Evol.py) module.\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:44.191724Z", "iopub.status.busy": "2021-03-07T17:14:44.190673Z", "iopub.status.idle": "2021-03-07T17:14:44.730409Z", "shell.execute_reply": "2021-03-07T17:14:44.730939Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistency check between MaxwellCartesian tutorial and NRPy+ module: ALL SHOULD BE ZERO.\n", "psi_rhs - mwevol.psi_rhs = 0\n", "Gamma_rhs - mwevol.Gamma_rhs = 0\n", "ArhsD[0] - mwevol.ArhsD[0] = 0\n", "ErhsD[0] - mwevol.ErhsD[0] = 0\n", "ArhsD[1] - mwevol.ArhsD[1] = 0\n", "ErhsD[1] - mwevol.ErhsD[1] = 0\n", "ArhsD[2] - mwevol.ArhsD[2] = 0\n", "ErhsD[2] - mwevol.ErhsD[2] = 0\n", "Cviolation - mwevol.Cviolation = 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", "# Step 18: Call the MaxwellCartesian_Evol() function from within the\n", "# Maxwell/MaxwellCartesian_Evol.py module,\n", "# which should do exactly the same as in Steps 1-16 above.\n", "\n", "par.set_parval_from_str(\"System_to_use\",\"System_II\")\n", "mwevol.MaxwellCartesian_Evol()\n", "\n", "print(\"Consistency check between MaxwellCartesian tutorial and NRPy+ module: ALL SHOULD BE ZERO.\")\n", "\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", "for i in range(DIM):\n", "\n", " print(\"ArhsD[\"+str(i)+\"] - mwevol.ArhsD[\"+str(i)+\"] = \" + str(ArhsD[i] - mwevol.ArhsD[i]))\n", " print(\"ErhsD[\"+str(i)+\"] - mwevol.ErhsD[\"+str(i)+\"] = \" + str(ErhsD[i] - mwevol.ErhsD[i]))\n", "\n", "print(\"Cviolation - mwevol.Cviolation = \" + str(Cviolation - mwevol.Cviolation))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 6: Constructing the Initial Data \\[Back to [top](#toc)\\]\n", "$$\\label{id}$$\n", "\n", "Now that we have evolution equations in place, we must construct the initial data that will be evolved by the solver of our choice. We will start from the analytic solution to this system of equations, given in [Illustrating Stability Properties of Numerical Relativity in Electrodynamics](https://arxiv.org/abs/gr-qc/0201051) as\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", "for vanishing scalar potential $\\psi$, where $\\mathcal{A}$ gives the amplitude, $\\lambda$ describes the size of the wavepacket, $u = t+r$, and $v = t-r$. Other components of this field are $0$.To get initial data, then, we simply set $t=0$; since $\\psi=0$, $E_i = \\partial_t A_i$. Thus, our initial data becomes the equations\n", "\\begin{align}\n", "A^{\\hat{\\phi}} &= 0 \\\\\n", "E^{\\hat{\\phi}} &= 8 \\mathcal{A} r \\sin \\theta \\lambda^2 e^{-\\lambda r^2} \\\\\n", "\\psi &= 0\n", "\\end{align}\n", "where the non-$\\hat{\\phi}$ components are set to 0. We still will need to convert $E^i$ in spherical-like coordinates and then lower its index. Using the standard transformations for coordinates and unit vectors, \n", "\\begin{align}\n", "E^{\\hat{x}} &= -\\frac{y E^{\\hat{\\phi}}(x,y,z)}{\\sqrt{x^2+y^2}} \\\\\n", "E^{\\hat{y}} &= \\frac{x E^{\\hat{\\phi}}(x,y,z)}{\\sqrt{x^2+y^2}} \\\\\n", "E^{\\hat{z}} &= 0. \\\\\n", "\\end{align}\n", "We can lower the index in the usual way.\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$. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:44.759961Z", "iopub.status.busy": "2021-03-07T17:14:44.746005Z", "iopub.status.idle": "2021-03-07T17:14:44.762818Z", "shell.execute_reply": "2021-03-07T17:14:44.762120Z" } }, "outputs": [], "source": [ "# Step 1: Declare free parameters intrinsic to these initial data\n", "amp,lam = par.Cparameters(\"REAL\",thismodule,\n", " [\"amp\",\"lam\"],\n", " [1.0,1.0]) # __name__ = \"MaxwellCartesian_ID\", this module's name\n", "\n", "# Step 2: Set the initial data\n", "AidD = ixp.zerorank1()\n", "\n", "EidD = ixp.zerorank1()\n", "EidU = ixp.zerorank1()\n", "# Set the coordinate transformations:\n", "radial = sp.sqrt(x*x + y*y + z*z)\n", "polar = sp.atan2(sp.sqrt(x*x + y*y),z)\n", "EU_phi = 8*amp*radial*sp.sin(polar)*lam*lam*sp.exp(-lam*radial*radial)\n", "EidU[0] = -(y * EU_phi)/sp.sqrt(x*x + y*y)\n", "EidU[1] = (x * EU_phi)/sp.sqrt(x*x + y*y)\n", "# The z component (2)is zero.\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " EidD[i] += gammaDD[i][j] * EidU[j]\n", "\n", "psi_ID = sp.sympify(0)\n", "Gamma_ID = sp.sympify(0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 7: Code Validation against `Maxwell.MaxwellCartesian_ID` NRPy+ Module \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation}$$\n", "\n", "Here, as a code validation check, we verify agreement in the SymPy expressions for the initial data we intend to use between\n", "\n", "1. this tutorial and \n", "2. the NRPy+ [Maxwell.MaxwellCartesian_ID](../edit/Maxwell/MaxwellCartesian_ID.py) module.\n", "\n", "Since the initial data is identical between the two systems for $E_i$, $A_i$, and $\\psi$, so checking system I should be redundant; we will do it anyways, to be sure." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:44.770226Z", "iopub.status.busy": "2021-03-07T17:14:44.769265Z", "iopub.status.idle": "2021-03-07T17:14:44.777126Z", "shell.execute_reply": "2021-03-07T17:14:44.777693Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "System I consistency check between MaxwellCartesian tutorial and NRPy+ module;\n", " ALL SHOULD BE ZERO:\n", "psi_ID - mwid.psi_ID = 0\n", "AidD[0] - mwid.AidD[0] = 0\n", "EidD[0] - mwid.EidD[0] = 0\n", "AidD[1] - mwid.AidD[1] = 0\n", "EidD[1] - mwid.EidD[1] = 0\n", "AidD[2] - mwid.AidD[2] = 0\n", "EidD[2] - mwid.EidD[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", "par.set_parval_from_str(\"System_to_use\",\"System_I\")\n", "\n", "import Maxwell.MaxwellCartesian_ID as mwid\n", "mwid.MaxwellCartesian_ID()\n", "\n", "print(\"System I consistency check between MaxwellCartesian tutorial and NRPy+ module;\\n ALL SHOULD BE ZERO:\")\n", "\n", "print(\"psi_ID - mwid.psi_ID = \" + str(psi_ID - mwid.psi_ID))\n", "for i in range(DIM):\n", " print(\"AidD[\"+str(i)+\"] - mwid.AidD[\"+str(i)+\"] = \" + str(AidD[i] - mwid.AidD[i]))\n", " print(\"EidD[\"+str(i)+\"] - mwid.EidD[\"+str(i)+\"] = \" + str(EidD[i] - mwid.EidD[i]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we will repeat the check with system II initial data." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:44.786625Z", "iopub.status.busy": "2021-03-07T17:14:44.785839Z", "iopub.status.idle": "2021-03-07T17:14:44.789043Z", "shell.execute_reply": "2021-03-07T17:14:44.789549Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "System II consistency check between MaxwellCartesian tutorial and NRPy+ module;\n", " ALL SHOULD BE ZERO:\n", "psi_ID - mwid.psi_ID = 0\n", "Gamma_ID - mwid.Gamma_ID = 0\n", "AidD[0] - mwid.AidD[0] = 0\n", "EidD[0] - mwid.EidD[0] = 0\n", "AidD[1] - mwid.AidD[1] = 0\n", "EidD[1] - mwid.EidD[1] = 0\n", "AidD[2] - mwid.AidD[2] = 0\n", "EidD[2] - mwid.EidD[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", "par.set_parval_from_str(\"System_to_use\",\"System_II\")\n", "\n", "mwid.MaxwellCartesian_ID()\n", "\n", "print(\"System II consistency check between MaxwellCartesian tutorial and NRPy+ module;\\n ALL SHOULD BE ZERO:\")\n", "\n", "print(\"psi_ID - mwid.psi_ID = \" + str(psi_ID - mwid.psi_ID))\n", "print(\"Gamma_ID - mwid.Gamma_ID = \" + str(Gamma_ID - mwid.Gamma_ID))\n", "for i in range(DIM):\n", " print(\"AidD[\"+str(i)+\"] - mwid.AidD[\"+str(i)+\"] = \" + str(AidD[i] - mwid.AidD[i]))\n", " print(\"EidD[\"+str(i)+\"] - mwid.EidD[\"+str(i)+\"] = \" + str(EidD[i] - mwid.EidD[i]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 8: 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-MaxwellCartesian.pdf](Tutorial-MaxwellCartesian.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": 10, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:44.794251Z", "iopub.status.busy": "2021-03-07T17:14:44.793531Z", "iopub.status.idle": "2021-03-07T17:14:48.289556Z", "shell.execute_reply": "2021-03-07T17:14:48.288814Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Created Tutorial-MaxwellCartesian.tex, and compiled LaTeX file to PDF file\n", " Tutorial-MaxwellCartesian.pdf\n" ] } ], "source": [ "import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n", "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-MaxwellCartesian\")" ] } ], "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 }