{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Generating C code for the right-hand sides of Maxwell's equations, in ***curvilinear*** coordinates, using a reference metric formalism\n", "\n", "## Author: Ian Ruchlin\n", "### Formatting improvements courtesy Brandon Clark\n", "\n", "[comment]: <> (Abstract: TODO)\n", "\n", "### The following formulations of Maxwell's equations, called System I and System II, are described in [Illustrating Stability Properties of Numerical Relativity in Electrodynamics](https://arxiv.org/abs/gr-qc/0201051) by Knapp et al.\n", "\n", "**Notebook Status:** In progress \n", "\n", "**Validation Notes:** This module has not yet undergone validation testing. Do ***not*** use it until after appropriate validation testing has been performed.\n", "\n", "## Introduction:\n", "[Maxwell's equations](https://en.wikipedia.org/wiki/Maxwell%27s_equations) are subject to the Gauss' law constraint\n", "$$\\mathcal{C} \\equiv \\hat{D}_{i} E^{i} - 4 \\pi \\rho = 0 \\; ,$$\n", "where $E^{i}$ is the electric vector field, $\\hat{D}_{i}$ is the [covariant derivative](https://en.wikipedia.org/wiki/Covariant_derivative) associated with the reference metric $\\hat{\\gamma}_{i j}$ (which is taken to represent flat space), and $\\rho$ is the electric charge density. We use $\\mathcal{C}$ as a measure of numerical error. Maxwell's equations are also required to satisfy $\\hat{D}_{i} B^{i} = 0$, where $B^{i}$ is the magnetic vector field. The magnetic constraint implies that the magnetic field can be expressed as\n", "$$B_{i} = \\epsilon_{i j k} \\hat{D}^{j} A^{k} \\; ,$$\n", "where $\\epsilon_{i j k}$ is the totally antisymmetric [Levi-Civita tensor](https://en.wikipedia.org/wiki/Levi-Civita_symbol) and $A^{i}$ is the vector potential field. Together with the scalar potential $\\psi$, the electric field can be expressed in terms of the potential fields as\n", "$$E_{i} = -\\hat{D}_{i} \\psi - \\partial_{t} A_{i} \\; .$$\n", "For now, we work in vacuum, where the electric charge density and the electric current density vector both vanish ($\\rho = 0$ and $j_{i} = 0$).\n", "\n", "In addition to the Gauss constraints, the electric and magnetic fields obey two independent [electromagnetic invariants](https://en.wikipedia.org/wiki/Classification_of_electromagnetic_fields#Invariants)\n", "\\begin{align}\n", "\\mathcal{P} &\\equiv B_{i} B^{i} - E_{i} E^{i} \\; , \\\\\n", "\\mathcal{Q} &\\equiv E_{i} B^{i} \\; .\n", "\\end{align}\n", "In vacuum, these satisfy $\\mathcal{P} = \\mathcal{Q} = 0$." ] }, { "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](#sys1): System I\n", "1. [Step 2](#sys2): System II\n", "1. [Step 3](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 1: System I \\[Back to [top](#toc)\\]\n", "$$\\label{sys1}$$\n", "\n", "In terms of the above definitions, the evolution Maxwell's equations take the form\n", "\\begin{align}\n", "\\partial_{t} A_{i} &= -E_{i} - \\hat{D}_{i} \\psi \\; , \\\\\n", "\\partial_{t} E_{i} &= -\\hat{D}_{j} \\hat{D}^{j} A_{i} + \\hat{D}_{i} \\hat{D}_{j} A^{j}\\; , \\\\\n", "\\partial_{t} \\psi &= -\\hat{D}_{i} A^{i} \\; .\n", "\\end{align}\n", "Note that this coupled system contains mixed second derivatives in the second term on the right hand side of the $E^{i}$ evolution equation. We will revisit this fact when building System II.\n", "\n", "It can be shown that the Gauss constraint satisfies the evolution equation\n", "$$\\partial_{t} \\mathcal{C} = 0 \\; .$$\n", "This implies that any constraint violating numerical error remains fixed in place during the evolution. This becomes problematic when the violations grow large and spoil the physics of the simulation." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:52.471716Z", "iopub.status.busy": "2021-03-07T17:14:52.466315Z", "iopub.status.idle": "2021-03-07T17:14:54.475005Z", "shell.execute_reply": "2021-03-07T17:14:54.475489Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{\n", " /*\n", " * NRPy+ Finite Difference Code Generation, Step 1 of 2: Read from main memory and compute finite difference stencils:\n", " */\n", " /*\n", " * Original SymPy expressions:\n", " * \"[const double AD_dD00 = invdx0*(-2*AD0_i0m1_i1_i2/3 + AD0_i0m2_i1_i2/12 + 2*AD0_i0p1_i1_i2/3 - AD0_i0p2_i1_i2/12),\n", " * const double AD_dD01 = invdx1*(-2*AD0_i0_i1m1_i2/3 + AD0_i0_i1m2_i2/12 + 2*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),\n", " * const double AD_dD02 = invdx2*(-2*AD0_i0_i1_i2m1/3 + AD0_i0_i1_i2m2/12 + 2*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),\n", " * const double AD_dD10 = invdx0*(-2*AD1_i0m1_i1_i2/3 + AD1_i0m2_i1_i2/12 + 2*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),\n", " * const double AD_dD11 = invdx1*(-2*AD1_i0_i1m1_i2/3 + AD1_i0_i1m2_i2/12 + 2*AD1_i0_i1p1_i2/3 - AD1_i0_i1p2_i2/12),\n", " * const double AD_dD12 = invdx2*(-2*AD1_i0_i1_i2m1/3 + AD1_i0_i1_i2m2/12 + 2*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),\n", " * const double AD_dD20 = invdx0*(-2*AD2_i0m1_i1_i2/3 + AD2_i0m2_i1_i2/12 + 2*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),\n", " * const double AD_dD21 = invdx1*(-2*AD2_i0_i1m1_i2/3 + AD2_i0_i1m2_i2/12 + 2*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),\n", " * const double AD_dD22 = invdx2*(-2*AD2_i0_i1_i2m1/3 + AD2_i0_i1_i2m2/12 + 2*AD2_i0_i1_i2p1/3 - AD2_i0_i1_i2p2/12),\n", " * const double AD_dDD001 = invdx0*invdx1*(4*AD0_i0m1_i1m1_i2/9 - AD0_i0m1_i1m2_i2/18 - 4*AD0_i0m1_i1p1_i2/9 + AD0_i0m1_i1p2_i2/18 - AD0_i0m2_i1m1_i2/18 + AD0_i0m2_i1m2_i2/144 + AD0_i0m2_i1p1_i2/18 - AD0_i0m2_i1p2_i2/144 - 4*AD0_i0p1_i1m1_i2/9 + AD0_i0p1_i1m2_i2/18 + 4*AD0_i0p1_i1p1_i2/9 - AD0_i0p1_i1p2_i2/18 + AD0_i0p2_i1m1_i2/18 - AD0_i0p2_i1m2_i2/144 - AD0_i0p2_i1p1_i2/18 + AD0_i0p2_i1p2_i2/144),\n", " * const double AD_dDD002 = invdx0*invdx2*(4*AD0_i0m1_i1_i2m1/9 - AD0_i0m1_i1_i2m2/18 - 4*AD0_i0m1_i1_i2p1/9 + AD0_i0m1_i1_i2p2/18 - AD0_i0m2_i1_i2m1/18 + AD0_i0m2_i1_i2m2/144 + AD0_i0m2_i1_i2p1/18 - AD0_i0m2_i1_i2p2/144 - 4*AD0_i0p1_i1_i2m1/9 + AD0_i0p1_i1_i2m2/18 + 4*AD0_i0p1_i1_i2p1/9 - AD0_i0p1_i1_i2p2/18 + AD0_i0p2_i1_i2m1/18 - AD0_i0p2_i1_i2m2/144 - AD0_i0p2_i1_i2p1/18 + AD0_i0p2_i1_i2p2/144),\n", " * const double AD_dDD011 = invdx1**2*(-5*AD0/2 + 4*AD0_i0_i1m1_i2/3 - AD0_i0_i1m2_i2/12 + 4*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),\n", " * const double AD_dDD022 = invdx2**2*(-5*AD0/2 + 4*AD0_i0_i1_i2m1/3 - AD0_i0_i1_i2m2/12 + 4*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),\n", " * const double AD_dDD100 = invdx0**2*(-5*AD1/2 + 4*AD1_i0m1_i1_i2/3 - AD1_i0m2_i1_i2/12 + 4*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),\n", " * const double AD_dDD101 = invdx0*invdx1*(4*AD1_i0m1_i1m1_i2/9 - AD1_i0m1_i1m2_i2/18 - 4*AD1_i0m1_i1p1_i2/9 + AD1_i0m1_i1p2_i2/18 - AD1_i0m2_i1m1_i2/18 + AD1_i0m2_i1m2_i2/144 + AD1_i0m2_i1p1_i2/18 - AD1_i0m2_i1p2_i2/144 - 4*AD1_i0p1_i1m1_i2/9 + AD1_i0p1_i1m2_i2/18 + 4*AD1_i0p1_i1p1_i2/9 - AD1_i0p1_i1p2_i2/18 + AD1_i0p2_i1m1_i2/18 - AD1_i0p2_i1m2_i2/144 - AD1_i0p2_i1p1_i2/18 + AD1_i0p2_i1p2_i2/144),\n", " * const double AD_dDD112 = invdx1*invdx2*(4*AD1_i0_i1m1_i2m1/9 - AD1_i0_i1m1_i2m2/18 - 4*AD1_i0_i1m1_i2p1/9 + AD1_i0_i1m1_i2p2/18 - AD1_i0_i1m2_i2m1/18 + AD1_i0_i1m2_i2m2/144 + AD1_i0_i1m2_i2p1/18 - AD1_i0_i1m2_i2p2/144 - 4*AD1_i0_i1p1_i2m1/9 + AD1_i0_i1p1_i2m2/18 + 4*AD1_i0_i1p1_i2p1/9 - AD1_i0_i1p1_i2p2/18 + AD1_i0_i1p2_i2m1/18 - AD1_i0_i1p2_i2m2/144 - AD1_i0_i1p2_i2p1/18 + AD1_i0_i1p2_i2p2/144),\n", " * const double AD_dDD122 = invdx2**2*(-5*AD1/2 + 4*AD1_i0_i1_i2m1/3 - AD1_i0_i1_i2m2/12 + 4*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),\n", " * const double AD_dDD200 = invdx0**2*(-5*AD2/2 + 4*AD2_i0m1_i1_i2/3 - AD2_i0m2_i1_i2/12 + 4*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),\n", " * const double AD_dDD202 = invdx0*invdx2*(4*AD2_i0m1_i1_i2m1/9 - AD2_i0m1_i1_i2m2/18 - 4*AD2_i0m1_i1_i2p1/9 + AD2_i0m1_i1_i2p2/18 - AD2_i0m2_i1_i2m1/18 + AD2_i0m2_i1_i2m2/144 + AD2_i0m2_i1_i2p1/18 - AD2_i0m2_i1_i2p2/144 - 4*AD2_i0p1_i1_i2m1/9 + AD2_i0p1_i1_i2m2/18 + 4*AD2_i0p1_i1_i2p1/9 - AD2_i0p1_i1_i2p2/18 + AD2_i0p2_i1_i2m1/18 - AD2_i0p2_i1_i2m2/144 - AD2_i0p2_i1_i2p1/18 + AD2_i0p2_i1_i2p2/144),\n", " * const double AD_dDD211 = invdx1**2*(-5*AD2/2 + 4*AD2_i0_i1m1_i2/3 - AD2_i0_i1m2_i2/12 + 4*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),\n", " * const double AD_dDD212 = invdx1*invdx2*(4*AD2_i0_i1m1_i2m1/9 - AD2_i0_i1m1_i2m2/18 - 4*AD2_i0_i1m1_i2p1/9 + AD2_i0_i1m1_i2p2/18 - AD2_i0_i1m2_i2m1/18 + AD2_i0_i1m2_i2m2/144 + AD2_i0_i1m2_i2p1/18 - AD2_i0_i1m2_i2p2/144 - 4*AD2_i0_i1p1_i2m1/9 + AD2_i0_i1p1_i2m2/18 + 4*AD2_i0_i1p1_i2p1/9 - AD2_i0_i1p1_i2p2/18 + AD2_i0_i1p2_i2m1/18 - AD2_i0_i1p2_i2m2/144 - AD2_i0_i1p2_i2p1/18 + AD2_i0_i1p2_i2p2/144),\n", " * const double psi_dD0 = invdx0*(-2*psi_i0m1_i1_i2/3 + psi_i0m2_i1_i2/12 + 2*psi_i0p1_i1_i2/3 - psi_i0p2_i1_i2/12),\n", " * const double psi_dD1 = invdx1*(-2*psi_i0_i1m1_i2/3 + psi_i0_i1m2_i2/12 + 2*psi_i0_i1p1_i2/3 - psi_i0_i1p2_i2/12),\n", " * const double psi_dD2 = invdx2*(-2*psi_i0_i1_i2m1/3 + psi_i0_i1_i2m2/12 + 2*psi_i0_i1_i2p1/3 - psi_i0_i1_i2p2/12)]\"\n", " */\n", " const double psi_i0_i1_i2m2 = in_gfs[IDX4(PSIGF, i0,i1,i2-2)];\n", " const double psi_i0_i1_i2m1 = in_gfs[IDX4(PSIGF, i0,i1,i2-1)];\n", " const double psi_i0_i1m2_i2 = in_gfs[IDX4(PSIGF, i0,i1-2,i2)];\n", " const double psi_i0_i1m1_i2 = in_gfs[IDX4(PSIGF, i0,i1-1,i2)];\n", " const double psi_i0m2_i1_i2 = in_gfs[IDX4(PSIGF, i0-2,i1,i2)];\n", " const double psi_i0m1_i1_i2 = in_gfs[IDX4(PSIGF, i0-1,i1,i2)];\n", " const double psi_i0p1_i1_i2 = in_gfs[IDX4(PSIGF, i0+1,i1,i2)];\n", " const double psi_i0p2_i1_i2 = in_gfs[IDX4(PSIGF, i0+2,i1,i2)];\n", " const double psi_i0_i1p1_i2 = in_gfs[IDX4(PSIGF, i0,i1+1,i2)];\n", " const double psi_i0_i1p2_i2 = in_gfs[IDX4(PSIGF, i0,i1+2,i2)];\n", " const double psi_i0_i1_i2p1 = in_gfs[IDX4(PSIGF, i0,i1,i2+1)];\n", " const double psi_i0_i1_i2p2 = in_gfs[IDX4(PSIGF, i0,i1,i2+2)];\n", " const double ED0 = in_gfs[IDX4(ED0GF, i0,i1,i2)];\n", " const double ED1 = in_gfs[IDX4(ED1GF, i0,i1,i2)];\n", " const double ED2 = in_gfs[IDX4(ED2GF, i0,i1,i2)];\n", " const double AD0_i0m2_i1_i2m2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2-2)];\n", " const double AD0_i0m1_i1_i2m2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2-2)];\n", " const double AD0_i0_i1_i2m2 = in_gfs[IDX4(AD0GF, i0,i1,i2-2)];\n", " const double AD0_i0p1_i1_i2m2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2-2)];\n", " const double AD0_i0p2_i1_i2m2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2-2)];\n", " const double AD0_i0m2_i1_i2m1 = in_gfs[IDX4(AD0GF, i0-2,i1,i2-1)];\n", " const double AD0_i0m1_i1_i2m1 = in_gfs[IDX4(AD0GF, i0-1,i1,i2-1)];\n", " const double AD0_i0_i1_i2m1 = in_gfs[IDX4(AD0GF, i0,i1,i2-1)];\n", " const double AD0_i0p1_i1_i2m1 = in_gfs[IDX4(AD0GF, i0+1,i1,i2-1)];\n", " const double AD0_i0p2_i1_i2m1 = in_gfs[IDX4(AD0GF, i0+2,i1,i2-1)];\n", " const double AD0_i0m2_i1m2_i2 = in_gfs[IDX4(AD0GF, i0-2,i1-2,i2)];\n", " const double AD0_i0m1_i1m2_i2 = in_gfs[IDX4(AD0GF, i0-1,i1-2,i2)];\n", " const double AD0_i0_i1m2_i2 = in_gfs[IDX4(AD0GF, i0,i1-2,i2)];\n", " const double AD0_i0p1_i1m2_i2 = in_gfs[IDX4(AD0GF, i0+1,i1-2,i2)];\n", " const double AD0_i0p2_i1m2_i2 = in_gfs[IDX4(AD0GF, i0+2,i1-2,i2)];\n", " const double AD0_i0m2_i1m1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1-1,i2)];\n", " const double AD0_i0m1_i1m1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1-1,i2)];\n", " const double AD0_i0_i1m1_i2 = in_gfs[IDX4(AD0GF, i0,i1-1,i2)];\n", " const double AD0_i0p1_i1m1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1-1,i2)];\n", " const double AD0_i0p2_i1m1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1-1,i2)];\n", " const double AD0_i0m2_i1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2)];\n", " const double AD0_i0m1_i1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2)];\n", " const double AD0 = in_gfs[IDX4(AD0GF, i0,i1,i2)];\n", " const double AD0_i0p1_i1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2)];\n", " const double AD0_i0p2_i1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2)];\n", " const double AD0_i0m2_i1p1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1+1,i2)];\n", " const double AD0_i0m1_i1p1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1+1,i2)];\n", " const double AD0_i0_i1p1_i2 = in_gfs[IDX4(AD0GF, i0,i1+1,i2)];\n", " const double AD0_i0p1_i1p1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1+1,i2)];\n", " const double AD0_i0p2_i1p1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1+1,i2)];\n", " const double AD0_i0m2_i1p2_i2 = in_gfs[IDX4(AD0GF, i0-2,i1+2,i2)];\n", " const double AD0_i0m1_i1p2_i2 = in_gfs[IDX4(AD0GF, i0-1,i1+2,i2)];\n", " const double AD0_i0_i1p2_i2 = in_gfs[IDX4(AD0GF, i0,i1+2,i2)];\n", " const double AD0_i0p1_i1p2_i2 = in_gfs[IDX4(AD0GF, i0+1,i1+2,i2)];\n", " const double AD0_i0p2_i1p2_i2 = in_gfs[IDX4(AD0GF, i0+2,i1+2,i2)];\n", " const double AD0_i0m2_i1_i2p1 = in_gfs[IDX4(AD0GF, i0-2,i1,i2+1)];\n", " const double AD0_i0m1_i1_i2p1 = in_gfs[IDX4(AD0GF, i0-1,i1,i2+1)];\n", " const double AD0_i0_i1_i2p1 = in_gfs[IDX4(AD0GF, i0,i1,i2+1)];\n", " const double AD0_i0p1_i1_i2p1 = in_gfs[IDX4(AD0GF, i0+1,i1,i2+1)];\n", " const double AD0_i0p2_i1_i2p1 = in_gfs[IDX4(AD0GF, i0+2,i1,i2+1)];\n", " const double AD0_i0m2_i1_i2p2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2+2)];\n", " const double AD0_i0m1_i1_i2p2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2+2)];\n", " const double AD0_i0_i1_i2p2 = in_gfs[IDX4(AD0GF, i0,i1,i2+2)];\n", " const double AD0_i0p1_i1_i2p2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2+2)];\n", " const double AD0_i0p2_i1_i2p2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2+2)];\n", " const double AD1_i0_i1m2_i2m2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2-2)];\n", " const double AD1_i0_i1m1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2-2)];\n", " const double AD1_i0_i1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1,i2-2)];\n", " const double AD1_i0_i1p1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2-2)];\n", " const double AD1_i0_i1p2_i2m2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2-2)];\n", " const double AD1_i0_i1m2_i2m1 = in_gfs[IDX4(AD1GF, i0,i1-2,i2-1)];\n", " const double AD1_i0_i1m1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1-1,i2-1)];\n", " const double AD1_i0_i1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1,i2-1)];\n", " const double AD1_i0_i1p1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1+1,i2-1)];\n", " const double AD1_i0_i1p2_i2m1 = in_gfs[IDX4(AD1GF, i0,i1+2,i2-1)];\n", " const double AD1_i0m2_i1m2_i2 = in_gfs[IDX4(AD1GF, i0-2,i1-2,i2)];\n", " const double AD1_i0m1_i1m2_i2 = in_gfs[IDX4(AD1GF, i0-1,i1-2,i2)];\n", " const double AD1_i0_i1m2_i2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2)];\n", " const double AD1_i0p1_i1m2_i2 = in_gfs[IDX4(AD1GF, i0+1,i1-2,i2)];\n", " const double AD1_i0p2_i1m2_i2 = in_gfs[IDX4(AD1GF, i0+2,i1-2,i2)];\n", " const double AD1_i0m2_i1m1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1-1,i2)];\n", " const double AD1_i0m1_i1m1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1-1,i2)];\n", " const double AD1_i0_i1m1_i2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2)];\n", " const double AD1_i0p1_i1m1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1-1,i2)];\n", " const double AD1_i0p2_i1m1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1-1,i2)];\n", " const double AD1_i0m2_i1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1,i2)];\n", " const double AD1_i0m1_i1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1,i2)];\n", " const double AD1 = in_gfs[IDX4(AD1GF, i0,i1,i2)];\n", " const double AD1_i0p1_i1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1,i2)];\n", " const double AD1_i0p2_i1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1,i2)];\n", " const double AD1_i0m2_i1p1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1+1,i2)];\n", " const double AD1_i0m1_i1p1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1+1,i2)];\n", " const double AD1_i0_i1p1_i2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2)];\n", " const double AD1_i0p1_i1p1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1+1,i2)];\n", " const double AD1_i0p2_i1p1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1+1,i2)];\n", " const double AD1_i0m2_i1p2_i2 = in_gfs[IDX4(AD1GF, i0-2,i1+2,i2)];\n", " const double AD1_i0m1_i1p2_i2 = in_gfs[IDX4(AD1GF, i0-1,i1+2,i2)];\n", " const double AD1_i0_i1p2_i2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2)];\n", " const double AD1_i0p1_i1p2_i2 = in_gfs[IDX4(AD1GF, i0+1,i1+2,i2)];\n", " const double AD1_i0p2_i1p2_i2 = in_gfs[IDX4(AD1GF, i0+2,i1+2,i2)];\n", " const double AD1_i0_i1m2_i2p1 = in_gfs[IDX4(AD1GF, i0,i1-2,i2+1)];\n", " const double AD1_i0_i1m1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1-1,i2+1)];\n", " const double AD1_i0_i1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1,i2+1)];\n", " const double AD1_i0_i1p1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1+1,i2+1)];\n", " const double AD1_i0_i1p2_i2p1 = in_gfs[IDX4(AD1GF, i0,i1+2,i2+1)];\n", " const double AD1_i0_i1m2_i2p2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2+2)];\n", " const double AD1_i0_i1m1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2+2)];\n", " const double AD1_i0_i1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1,i2+2)];\n", " const double AD1_i0_i1p1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2+2)];\n", " const double AD1_i0_i1p2_i2p2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2+2)];\n", " const double AD2_i0_i1m2_i2m2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2-2)];\n", " const double AD2_i0_i1m1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2-2)];\n", " const double AD2_i0m2_i1_i2m2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2-2)];\n", " const double AD2_i0m1_i1_i2m2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2-2)];\n", " const double AD2_i0_i1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1,i2-2)];\n", " const double AD2_i0p1_i1_i2m2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2-2)];\n", " const double AD2_i0p2_i1_i2m2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2-2)];\n", " const double AD2_i0_i1p1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2-2)];\n", " const double AD2_i0_i1p2_i2m2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2-2)];\n", " const double AD2_i0_i1m2_i2m1 = in_gfs[IDX4(AD2GF, i0,i1-2,i2-1)];\n", " const double AD2_i0_i1m1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1-1,i2-1)];\n", " const double AD2_i0m2_i1_i2m1 = in_gfs[IDX4(AD2GF, i0-2,i1,i2-1)];\n", " const double AD2_i0m1_i1_i2m1 = in_gfs[IDX4(AD2GF, i0-1,i1,i2-1)];\n", " const double AD2_i0_i1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1,i2-1)];\n", " const double AD2_i0p1_i1_i2m1 = in_gfs[IDX4(AD2GF, i0+1,i1,i2-1)];\n", " const double AD2_i0p2_i1_i2m1 = in_gfs[IDX4(AD2GF, i0+2,i1,i2-1)];\n", " const double AD2_i0_i1p1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1+1,i2-1)];\n", " const double AD2_i0_i1p2_i2m1 = in_gfs[IDX4(AD2GF, i0,i1+2,i2-1)];\n", " const double AD2_i0_i1m2_i2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2)];\n", " const double AD2_i0_i1m1_i2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2)];\n", " const double AD2_i0m2_i1_i2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2)];\n", " const double AD2_i0m1_i1_i2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2)];\n", " const double AD2 = in_gfs[IDX4(AD2GF, i0,i1,i2)];\n", " const double AD2_i0p1_i1_i2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2)];\n", " const double AD2_i0p2_i1_i2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2)];\n", " const double AD2_i0_i1p1_i2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2)];\n", " const double AD2_i0_i1p2_i2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2)];\n", " const double AD2_i0_i1m2_i2p1 = in_gfs[IDX4(AD2GF, i0,i1-2,i2+1)];\n", " const double AD2_i0_i1m1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1-1,i2+1)];\n", " const double AD2_i0m2_i1_i2p1 = in_gfs[IDX4(AD2GF, i0-2,i1,i2+1)];\n", " const double AD2_i0m1_i1_i2p1 = in_gfs[IDX4(AD2GF, i0-1,i1,i2+1)];\n", " const double AD2_i0_i1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1,i2+1)];\n", " const double AD2_i0p1_i1_i2p1 = in_gfs[IDX4(AD2GF, i0+1,i1,i2+1)];\n", " const double AD2_i0p2_i1_i2p1 = in_gfs[IDX4(AD2GF, i0+2,i1,i2+1)];\n", " const double AD2_i0_i1p1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1+1,i2+1)];\n", " const double AD2_i0_i1p2_i2p1 = in_gfs[IDX4(AD2GF, i0,i1+2,i2+1)];\n", " const double AD2_i0_i1m2_i2p2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2+2)];\n", " const double AD2_i0_i1m1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2+2)];\n", " const double AD2_i0m2_i1_i2p2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2+2)];\n", " const double AD2_i0m1_i1_i2p2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2+2)];\n", " const double AD2_i0_i1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1,i2+2)];\n", " const double AD2_i0p1_i1_i2p2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2+2)];\n", " const double AD2_i0p2_i1_i2p2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2+2)];\n", " const double AD2_i0_i1p1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2+2)];\n", " const double AD2_i0_i1p2_i2p2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2+2)];\n", " const double FDPart1_Rational_2_3 = 2.0/3.0;\n", " const double FDPart1_Rational_1_12 = 1.0/12.0;\n", " const double FDPart1_Rational_4_9 = 4.0/9.0;\n", " const double FDPart1_Rational_1_18 = 1.0/18.0;\n", " const double FDPart1_Rational_1_144 = 1.0/144.0;\n", " const double FDPart1_Rational_5_2 = 5.0/2.0;\n", " const double FDPart1_Rational_4_3 = 4.0/3.0;\n", " const double FDPart1_1 = -AD0_i0_i1_i2p2;\n", " const double FDPart1_9 = -AD0*FDPart1_Rational_5_2;\n", " const double FDPart1_12 = -AD1*FDPart1_Rational_5_2;\n", " const double FDPart1_14 = -AD2*FDPart1_Rational_5_2;\n", " const double AD_dD00 = invdx0*(FDPart1_Rational_1_12*(AD0_i0m2_i1_i2 - AD0_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD0_i0m1_i1_i2 + AD0_i0p1_i1_i2));\n", " const double AD_dD01 = invdx1*(FDPart1_Rational_1_12*(AD0_i0_i1m2_i2 - AD0_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));\n", " const double AD_dD02 = invdx2*(FDPart1_Rational_1_12*(AD0_i0_i1_i2m2 + FDPart1_1) + FDPart1_Rational_2_3*(-AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));\n", " const double AD_dD10 = invdx0*(FDPart1_Rational_1_12*(AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));\n", " const double AD_dD11 = invdx1*(FDPart1_Rational_1_12*(AD1_i0_i1m2_i2 - AD1_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD1_i0_i1m1_i2 + AD1_i0_i1p1_i2));\n", " const double AD_dD12 = invdx2*(FDPart1_Rational_1_12*(AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));\n", " const double AD_dD20 = invdx0*(FDPart1_Rational_1_12*(AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));\n", " const double AD_dD21 = invdx1*(FDPart1_Rational_1_12*(AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));\n", " const double AD_dD22 = invdx2*(FDPart1_Rational_1_12*(AD2_i0_i1_i2m2 - AD2_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD2_i0_i1_i2m1 + AD2_i0_i1_i2p1));\n", " const double AD_dDD001 = invdx0*invdx1*(FDPart1_Rational_1_144*(AD0_i0m2_i1m2_i2 - AD0_i0m2_i1p2_i2 - AD0_i0p2_i1m2_i2 + AD0_i0p2_i1p2_i2) + FDPart1_Rational_1_18*(-AD0_i0m1_i1m2_i2 + AD0_i0m1_i1p2_i2 - AD0_i0m2_i1m1_i2 + AD0_i0m2_i1p1_i2 + AD0_i0p1_i1m2_i2 - AD0_i0p1_i1p2_i2 + AD0_i0p2_i1m1_i2 - AD0_i0p2_i1p1_i2) + FDPart1_Rational_4_9*(AD0_i0m1_i1m1_i2 - AD0_i0m1_i1p1_i2 - AD0_i0p1_i1m1_i2 + AD0_i0p1_i1p1_i2));\n", " const double AD_dDD002 = invdx0*invdx2*(FDPart1_Rational_1_144*(AD0_i0m2_i1_i2m2 - AD0_i0m2_i1_i2p2 - AD0_i0p2_i1_i2m2 + AD0_i0p2_i1_i2p2) + FDPart1_Rational_1_18*(-AD0_i0m1_i1_i2m2 + AD0_i0m1_i1_i2p2 - AD0_i0m2_i1_i2m1 + AD0_i0m2_i1_i2p1 + AD0_i0p1_i1_i2m2 - AD0_i0p1_i1_i2p2 + AD0_i0p2_i1_i2m1 - AD0_i0p2_i1_i2p1) + FDPart1_Rational_4_9*(AD0_i0m1_i1_i2m1 - AD0_i0m1_i1_i2p1 - AD0_i0p1_i1_i2m1 + AD0_i0p1_i1_i2p1));\n", " const double AD_dDD011 = ((invdx1)*(invdx1))*(FDPart1_9 + FDPart1_Rational_1_12*(-AD0_i0_i1m2_i2 - AD0_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));\n", " const double AD_dDD022 = ((invdx2)*(invdx2))*(FDPart1_9 + FDPart1_Rational_1_12*(-AD0_i0_i1_i2m2 + FDPart1_1) + FDPart1_Rational_4_3*(AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));\n", " const double AD_dDD100 = ((invdx0)*(invdx0))*(FDPart1_12 + FDPart1_Rational_1_12*(-AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));\n", " const double AD_dDD101 = invdx0*invdx1*(FDPart1_Rational_1_144*(AD1_i0m2_i1m2_i2 - AD1_i0m2_i1p2_i2 - AD1_i0p2_i1m2_i2 + AD1_i0p2_i1p2_i2) + FDPart1_Rational_1_18*(-AD1_i0m1_i1m2_i2 + AD1_i0m1_i1p2_i2 - AD1_i0m2_i1m1_i2 + AD1_i0m2_i1p1_i2 + AD1_i0p1_i1m2_i2 - AD1_i0p1_i1p2_i2 + AD1_i0p2_i1m1_i2 - AD1_i0p2_i1p1_i2) + FDPart1_Rational_4_9*(AD1_i0m1_i1m1_i2 - AD1_i0m1_i1p1_i2 - AD1_i0p1_i1m1_i2 + AD1_i0p1_i1p1_i2));\n", " const double AD_dDD112 = invdx1*invdx2*(FDPart1_Rational_1_144*(AD1_i0_i1m2_i2m2 - AD1_i0_i1m2_i2p2 - AD1_i0_i1p2_i2m2 + AD1_i0_i1p2_i2p2) + FDPart1_Rational_1_18*(-AD1_i0_i1m1_i2m2 + AD1_i0_i1m1_i2p2 - AD1_i0_i1m2_i2m1 + AD1_i0_i1m2_i2p1 + AD1_i0_i1p1_i2m2 - AD1_i0_i1p1_i2p2 + AD1_i0_i1p2_i2m1 - AD1_i0_i1p2_i2p1) + FDPart1_Rational_4_9*(AD1_i0_i1m1_i2m1 - AD1_i0_i1m1_i2p1 - AD1_i0_i1p1_i2m1 + AD1_i0_i1p1_i2p1));\n", " const double AD_dDD122 = ((invdx2)*(invdx2))*(FDPart1_12 + FDPart1_Rational_1_12*(-AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));\n", " const double AD_dDD200 = ((invdx0)*(invdx0))*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));\n", " const double AD_dDD202 = invdx0*invdx2*(FDPart1_Rational_1_144*(AD2_i0m2_i1_i2m2 - AD2_i0m2_i1_i2p2 - AD2_i0p2_i1_i2m2 + AD2_i0p2_i1_i2p2) + FDPart1_Rational_1_18*(-AD2_i0m1_i1_i2m2 + AD2_i0m1_i1_i2p2 - AD2_i0m2_i1_i2m1 + AD2_i0m2_i1_i2p1 + AD2_i0p1_i1_i2m2 - AD2_i0p1_i1_i2p2 + AD2_i0p2_i1_i2m1 - AD2_i0p2_i1_i2p1) + FDPart1_Rational_4_9*(AD2_i0m1_i1_i2m1 - AD2_i0m1_i1_i2p1 - AD2_i0p1_i1_i2m1 + AD2_i0p1_i1_i2p1));\n", " const double AD_dDD211 = ((invdx1)*(invdx1))*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));\n", " const double AD_dDD212 = invdx1*invdx2*(FDPart1_Rational_1_144*(AD2_i0_i1m2_i2m2 - AD2_i0_i1m2_i2p2 - AD2_i0_i1p2_i2m2 + AD2_i0_i1p2_i2p2) + FDPart1_Rational_1_18*(-AD2_i0_i1m1_i2m2 + AD2_i0_i1m1_i2p2 - AD2_i0_i1m2_i2m1 + AD2_i0_i1m2_i2p1 + AD2_i0_i1p1_i2m2 - AD2_i0_i1p1_i2p2 + AD2_i0_i1p2_i2m1 - AD2_i0_i1p2_i2p1) + FDPart1_Rational_4_9*(AD2_i0_i1m1_i2m1 - AD2_i0_i1m1_i2p1 - AD2_i0_i1p1_i2m1 + AD2_i0_i1p1_i2p1));\n", " const double psi_dD0 = invdx0*(FDPart1_Rational_1_12*(psi_i0m2_i1_i2 - psi_i0p2_i1_i2) + FDPart1_Rational_2_3*(-psi_i0m1_i1_i2 + psi_i0p1_i1_i2));\n", " const double psi_dD1 = invdx1*(FDPart1_Rational_1_12*(psi_i0_i1m2_i2 - psi_i0_i1p2_i2) + FDPart1_Rational_2_3*(-psi_i0_i1m1_i2 + psi_i0_i1p1_i2));\n", " const double psi_dD2 = invdx2*(FDPart1_Rational_1_12*(psi_i0_i1_i2m2 - psi_i0_i1_i2p2) + FDPart1_Rational_2_3*(-psi_i0_i1_i2m1 + psi_i0_i1_i2p1));\n", " /*\n", " * NRPy+ Finite Difference Code Generation, Step 2 of 2: Evaluate SymPy expressions and write to main memory:\n", " */\n", " /*\n", " * Original SymPy expressions:\n", " * \"[rhs_gfs[IDX4(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0,\n", " * rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = (AD0 + AD_dD00*xx0 + AD_dDD101 - 2*(AD0*xx0 + AD_dD11)/xx0)/xx0**2 - (AD_dD00*xx0 - AD_dD11/xx0 + AD_dDD011 - (AD0*xx0 + AD_dD11)/xx0)/xx0**2 + (AD0*sin(xx1)**2 + AD_dD00*xx0*sin(xx1)**2 + AD_dD10*sin(2*xx1)/2 + AD_dDD202 - 2*(AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/xx0)/(xx0**2*sin(xx1)**2) - (AD_dD00*xx0*sin(xx1)**2 - AD_dD22/xx0 + AD_dDD022 + (-AD1/xx0 + AD_dD01)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/xx0)/(xx0**2*sin(xx1)**2),\n", " * rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1,\n", " * rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1/xx0**2 + AD_dD10/xx0 + AD_dDD001 - AD_dDD100 - (-AD1/xx0 + AD_dD01)/xx0 - (-AD_dD22*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD122 + xx0*(-AD1/xx0 + AD_dD10)*sin(xx1)**2 + (AD0*xx0 + AD_dD11)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)*sin(2*xx1)/(2*sin(xx1)**2))/(xx0**2*sin(xx1)**2) + (2*AD0*xx0*sin(xx1)*cos(xx1) + AD1*cos(2*xx1) + AD_dD01*xx0*sin(xx1)**2 + AD_dD11*sin(2*xx1)/2 + AD_dDD212 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)*sin(2*xx1)/sin(xx1)**2)/(xx0**2*sin(xx1)**2),\n", " * rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2,\n", " * rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2/xx0**2 + AD_dD20/xx0 + AD_dDD002 - AD_dDD200 - (-AD2/xx0 + AD_dD02)/xx0 + (AD_dD02*xx0 + AD_dDD112 - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD12)*sin(2*xx1)/(2*sin(xx1)**2) - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/(2*sin(xx1)**2))/xx0**2 - (AD2*(-cos(2*xx1)/sin(xx1)**2 + sin(2*xx1)*cos(xx1)/sin(xx1)**3) - AD_dD21*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD211 + xx0*(-AD2/xx0 + AD_dD20) - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/(2*sin(xx1)**2))/xx0**2,\n", " * rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -AD_dD00 - (AD0*xx0 + AD_dD11)/xx0**2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/(xx0**2*sin(xx1)**2)]\"\n", " */\n", " const double FDPart3_0 = (1.0/((xx0)*(xx0)));\n", " const double FDPart3_1 = AD_dD00*xx0;\n", " const double FDPart3_2 = (1.0/(xx0));\n", " const double FDPart3_3 = AD0*xx0;\n", " const double FDPart3_4 = AD_dD11 + FDPart3_3;\n", " const double FDPart3_6 = sin(xx1);\n", " const double FDPart3_7 = ((FDPart3_6)*(FDPart3_6));\n", " const double FDPart3_9 = sin(2*xx1);\n", " const double FDPart3_10 = (1.0/2.0)*FDPart3_9;\n", " const double FDPart3_12 = AD1*FDPart3_10 + AD_dD22 + FDPart3_3*FDPart3_7;\n", " const double FDPart3_14 = (1.0/(FDPart3_7));\n", " const double FDPart3_15 = FDPart3_0*FDPart3_14;\n", " const double FDPart3_16 = -AD1*FDPart3_2;\n", " const double FDPart3_18 = cos(2*xx1);\n", " const double FDPart3_20 = cos(xx1);\n", " const double FDPart3_22 = FDPart3_10*FDPart3_14;\n", " const double FDPart3_23 = -AD2*FDPart3_2;\n", " const double FDPart3_24 = -AD2*FDPart3_22;\n", " const double FDPart3_25 = -FDPart3_22*(AD_dD21 + FDPart3_24);\n", " rhs_gfs[IDX4(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0;\n", " rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = FDPart3_0*(AD0 + AD_dDD101 + FDPart3_1 - 2*FDPart3_2*FDPart3_4) - FDPart3_0*(-AD_dD11*FDPart3_2 + AD_dDD011 + FDPart3_1 - FDPart3_2*FDPart3_4) + FDPart3_15*(AD0*FDPart3_7 + AD_dD10*FDPart3_10 + AD_dDD202 + FDPart3_1*FDPart3_7 - 2*FDPart3_12*FDPart3_2) - FDPart3_15*(-AD_dD22*FDPart3_2 + AD_dDD022 + FDPart3_1*FDPart3_7 + FDPart3_10*(AD_dD01 + FDPart3_16) - FDPart3_12*FDPart3_2);\n", " rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1;\n", " rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1*FDPart3_0 + AD_dD10*FDPart3_2 + AD_dDD001 - AD_dDD100 - FDPart3_15*(-AD_dD22*FDPart3_22 + AD_dDD122 - FDPart3_10*FDPart3_12*FDPart3_14 + FDPart3_10*FDPart3_4 + FDPart3_7*xx0*(AD_dD10 + FDPart3_16)) + FDPart3_15*(AD1*FDPart3_18 + AD_dD01*FDPart3_7*xx0 + AD_dD11*FDPart3_10 + AD_dDD212 - FDPart3_12*FDPart3_14*FDPart3_9 + 2*FDPart3_20*FDPart3_3*FDPart3_6) - FDPart3_2*(AD_dD01 + FDPart3_16);\n", " rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2;\n", " rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2*FDPart3_0 + AD_dD20*FDPart3_2 + AD_dDD002 - AD_dDD200 + FDPart3_0*(AD_dD02*xx0 + AD_dDD112 - FDPart3_22*(AD_dD12 + FDPart3_24) + FDPart3_25) - FDPart3_0*(AD2*(-FDPart3_14*FDPart3_18 + FDPart3_20*FDPart3_9/((FDPart3_6)*(FDPart3_6)*(FDPart3_6))) - AD_dD21*FDPart3_22 + AD_dDD211 + FDPart3_25 + xx0*(AD_dD20 + FDPart3_23)) - FDPart3_2*(AD_dD02 + FDPart3_23);\n", " rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -AD_dD00 - FDPart3_0*FDPart3_4 - FDPart3_12*FDPart3_15;\n", "}\n" ] } ], "source": [ "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 grid as gri # NRPy+: Functions having to do with numerical grids\n", "import finite_difference as fin # NRPy+: Finite difference C code generation module\n", "import reference_metric as rfm # NRPy+: Reference metric support\n", "from outputC import lhrh # NRPy+: Core C code output module\n", "\n", "par.set_parval_from_str(\"reference_metric::CoordSystem\", \"Spherical\")\n", "par.set_parval_from_str(\"grid::DIM\", 3)\n", "\n", "rfm.reference_metric()\n", "\n", "# The name of this module (\"maxwell\") is given by __name__:\n", "thismodule = __name__\n", "\n", "# Step 0: Read the spatial dimension parameter as DIM.\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", "\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_dHatD = ixp.zerorank2()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " AD_dHatD[i][j] = AD_dD[i][j]\n", " for k in range(DIM):\n", " AD_dHatD[i][j] -= rfm.GammahatUDD[k][i][j] * AD[k]\n", "\n", "# Second covariant derivative\n", "# D_{k} D_{j} A_{i} = \\partial_{k} D_{j} A_{i} - \\Gamma^{l}_{jk} D_{l} A_{i}\n", "# - \\Gamma^{l}_{ik} D_{j} A_{l}\n", "# = A_{i,jk}\n", "# - \\Gamma^{l}_{ij,k} A_{l}\n", "# - \\Gamma^{l}_{ij} A_{l,k}\n", "# - \\Gamma^{l}_{jk} A_{i;\\hat{l}}\n", "# - \\Gamma^{l}_{ik} A_{l;\\hat{j}}\n", "AD_dHatDD = ixp.zerorank3()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " for k in range(DIM):\n", " AD_dHatDD[i][j][k] = AD_dDD[i][j][k]\n", " for l in range(DIM):\n", " AD_dHatDD[i][j][k] += - rfm.GammahatUDDdD[l][i][j][k] * AD[l] \\\n", " - rfm.GammahatUDD[l][i][j] * AD_dD[l][k] \\\n", " - rfm.GammahatUDD[l][j][k] * AD_dHatD[i][l] \\\n", " - rfm.GammahatUDD[l][i][k] * AD_dHatD[l][j]\n", "\n", "# Covariant divergence\n", "# D_{i} A^{i} = ghat^{ij} D_{j} A_{i}\n", "DivA = 0\n", "# Gradient of covariant divergence\n", "# DivA_dD_{i} = ghat^{jk} A_{k;\\hat{j}\\hat{i}}\n", "DivA_dD = ixp.zerorank1()\n", "# Covariant Laplacian\n", "# LapAD_{i} = ghat^{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 += rfm.ghatUU[i][j] * AD_dHatD[i][j]\n", " for k in range(DIM):\n", " DivA_dD[i] += rfm.ghatUU[j][k] * AD_dHatDD[k][j][i]\n", " LapAD[i] += rfm.ghatUU[j][k] * AD_dHatDD[i][j][k]\n", "\n", "# Step 5: Define right-hand sides for the evolution.\n", "AD_rhs = ixp.zerorank1()\n", "ED_rhs = ixp.zerorank1()\n", "for i in range(DIM):\n", " AD_rhs[i] = -ED[i] - psi_dD[i]\n", " ED_rhs[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=AD_rhs[i]))\n", " lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"ED\" + str(i)), rhs=ED_rhs[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 2: System II \\[Back to [top](#toc)\\]\n", "$$\\label{sys2}$$\n", "\n", "Define the auxiliary variable\n", "$$\\Gamma \\equiv \\hat{D}_{i} A^{i} \\; .$$\n", "Substituting this into Maxwell's equations yields the system\n", "\\begin{align}\n", "\\partial_{t} A_{i} &= -E_{i} - \\hat{D}_{i} \\psi \\; , \\\\\n", "\\partial_{t} E_{i} &= -\\hat{D}_{j} \\hat{D}^{j} A_{i} + \\hat{D}_{i} \\Gamma \\; , \\\\\n", "\\partial_{t} \\psi &= -\\Gamma \\; , \\\\\n", "\\partial_{t} \\Gamma &= -\\hat{D}_{i} \\hat{D}^{i} \\psi \\; .\n", "\\end{align}\n", "\n", "\n", "\n", "It can be shown that the Gauss constraint now satisfies the wave equation\n", "$$\\partial_{t}^{2} \\mathcal{C} = \\hat{D}_{i} \\hat{D}^{i} \\mathcal{C} \\; .$$\n", "Thus, any constraint violation introduced by numerical error propagates away at the speed of light. This property increases the stability of of the simulation, compared to System I above. A similar trick is used in the [BSSN formulation](Tutorial-BSSNCurvilinear.ipynb) of Einstein's equations." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:54.549764Z", "iopub.status.busy": "2021-03-07T17:14:54.513988Z", "iopub.status.idle": "2021-03-07T17:14:55.234470Z", "shell.execute_reply": "2021-03-07T17:14:55.233887Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{\n", " /*\n", " * NRPy+ Finite Difference Code Generation, Step 1 of 2: Read from main memory and compute finite difference stencils:\n", " */\n", " /*\n", " * Original SymPy expressions:\n", " * \"[const double AD_dD00 = invdx0*(-2*AD0_i0m1_i1_i2/3 + AD0_i0m2_i1_i2/12 + 2*AD0_i0p1_i1_i2/3 - AD0_i0p2_i1_i2/12),\n", " * const double AD_dD01 = invdx1*(-2*AD0_i0_i1m1_i2/3 + AD0_i0_i1m2_i2/12 + 2*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),\n", " * const double AD_dD02 = invdx2*(-2*AD0_i0_i1_i2m1/3 + AD0_i0_i1_i2m2/12 + 2*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),\n", " * const double AD_dD10 = invdx0*(-2*AD1_i0m1_i1_i2/3 + AD1_i0m2_i1_i2/12 + 2*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),\n", " * const double AD_dD11 = invdx1*(-2*AD1_i0_i1m1_i2/3 + AD1_i0_i1m2_i2/12 + 2*AD1_i0_i1p1_i2/3 - AD1_i0_i1p2_i2/12),\n", " * const double AD_dD12 = invdx2*(-2*AD1_i0_i1_i2m1/3 + AD1_i0_i1_i2m2/12 + 2*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),\n", " * const double AD_dD20 = invdx0*(-2*AD2_i0m1_i1_i2/3 + AD2_i0m2_i1_i2/12 + 2*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),\n", " * const double AD_dD21 = invdx1*(-2*AD2_i0_i1m1_i2/3 + AD2_i0_i1m2_i2/12 + 2*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),\n", " * const double AD_dD22 = invdx2*(-2*AD2_i0_i1_i2m1/3 + AD2_i0_i1_i2m2/12 + 2*AD2_i0_i1_i2p1/3 - AD2_i0_i1_i2p2/12),\n", " * const double AD_dDD000 = invdx0**2*(-5*AD0/2 + 4*AD0_i0m1_i1_i2/3 - AD0_i0m2_i1_i2/12 + 4*AD0_i0p1_i1_i2/3 - AD0_i0p2_i1_i2/12),\n", " * const double AD_dDD011 = invdx1**2*(-5*AD0/2 + 4*AD0_i0_i1m1_i2/3 - AD0_i0_i1m2_i2/12 + 4*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),\n", " * const double AD_dDD022 = invdx2**2*(-5*AD0/2 + 4*AD0_i0_i1_i2m1/3 - AD0_i0_i1_i2m2/12 + 4*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),\n", " * const double AD_dDD100 = invdx0**2*(-5*AD1/2 + 4*AD1_i0m1_i1_i2/3 - AD1_i0m2_i1_i2/12 + 4*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),\n", " * const double AD_dDD111 = invdx1**2*(-5*AD1/2 + 4*AD1_i0_i1m1_i2/3 - AD1_i0_i1m2_i2/12 + 4*AD1_i0_i1p1_i2/3 - AD1_i0_i1p2_i2/12),\n", " * const double AD_dDD122 = invdx2**2*(-5*AD1/2 + 4*AD1_i0_i1_i2m1/3 - AD1_i0_i1_i2m2/12 + 4*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),\n", " * const double AD_dDD200 = invdx0**2*(-5*AD2/2 + 4*AD2_i0m1_i1_i2/3 - AD2_i0m2_i1_i2/12 + 4*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),\n", " * const double AD_dDD211 = invdx1**2*(-5*AD2/2 + 4*AD2_i0_i1m1_i2/3 - AD2_i0_i1m2_i2/12 + 4*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),\n", " * const double AD_dDD222 = invdx2**2*(-5*AD2/2 + 4*AD2_i0_i1_i2m1/3 - AD2_i0_i1_i2m2/12 + 4*AD2_i0_i1_i2p1/3 - AD2_i0_i1_i2p2/12),\n", " * const double Gamma_dD0 = invdx0*(-2*Gamma_i0m1_i1_i2/3 + Gamma_i0m2_i1_i2/12 + 2*Gamma_i0p1_i1_i2/3 - Gamma_i0p2_i1_i2/12),\n", " * const double Gamma_dD1 = invdx1*(-2*Gamma_i0_i1m1_i2/3 + Gamma_i0_i1m2_i2/12 + 2*Gamma_i0_i1p1_i2/3 - Gamma_i0_i1p2_i2/12),\n", " * const double Gamma_dD2 = invdx2*(-2*Gamma_i0_i1_i2m1/3 + Gamma_i0_i1_i2m2/12 + 2*Gamma_i0_i1_i2p1/3 - Gamma_i0_i1_i2p2/12),\n", " * const double psi_dD0 = invdx0*(-2*psi_i0m1_i1_i2/3 + psi_i0m2_i1_i2/12 + 2*psi_i0p1_i1_i2/3 - psi_i0p2_i1_i2/12),\n", " * const double psi_dD1 = invdx1*(-2*psi_i0_i1m1_i2/3 + psi_i0_i1m2_i2/12 + 2*psi_i0_i1p1_i2/3 - psi_i0_i1p2_i2/12),\n", " * const double psi_dD2 = invdx2*(-2*psi_i0_i1_i2m1/3 + psi_i0_i1_i2m2/12 + 2*psi_i0_i1_i2p1/3 - psi_i0_i1_i2p2/12),\n", " * const double psi_dDD00 = invdx0**2*(-5*psi/2 + 4*psi_i0m1_i1_i2/3 - psi_i0m2_i1_i2/12 + 4*psi_i0p1_i1_i2/3 - psi_i0p2_i1_i2/12),\n", " * const double psi_dDD11 = invdx1**2*(-5*psi/2 + 4*psi_i0_i1m1_i2/3 - psi_i0_i1m2_i2/12 + 4*psi_i0_i1p1_i2/3 - psi_i0_i1p2_i2/12),\n", " * const double psi_dDD22 = invdx2**2*(-5*psi/2 + 4*psi_i0_i1_i2m1/3 - psi_i0_i1_i2m2/12 + 4*psi_i0_i1_i2p1/3 - psi_i0_i1_i2p2/12)]\"\n", " */\n", " const double psi_i0_i1_i2m2 = in_gfs[IDX4(PSIGF, i0,i1,i2-2)];\n", " const double psi_i0_i1_i2m1 = in_gfs[IDX4(PSIGF, i0,i1,i2-1)];\n", " const double psi_i0_i1m2_i2 = in_gfs[IDX4(PSIGF, i0,i1-2,i2)];\n", " const double psi_i0_i1m1_i2 = in_gfs[IDX4(PSIGF, i0,i1-1,i2)];\n", " const double psi_i0m2_i1_i2 = in_gfs[IDX4(PSIGF, i0-2,i1,i2)];\n", " const double psi_i0m1_i1_i2 = in_gfs[IDX4(PSIGF, i0-1,i1,i2)];\n", " const double psi = in_gfs[IDX4(PSIGF, i0,i1,i2)];\n", " const double psi_i0p1_i1_i2 = in_gfs[IDX4(PSIGF, i0+1,i1,i2)];\n", " const double psi_i0p2_i1_i2 = in_gfs[IDX4(PSIGF, i0+2,i1,i2)];\n", " const double psi_i0_i1p1_i2 = in_gfs[IDX4(PSIGF, i0,i1+1,i2)];\n", " const double psi_i0_i1p2_i2 = in_gfs[IDX4(PSIGF, i0,i1+2,i2)];\n", " const double psi_i0_i1_i2p1 = in_gfs[IDX4(PSIGF, i0,i1,i2+1)];\n", " const double psi_i0_i1_i2p2 = in_gfs[IDX4(PSIGF, i0,i1,i2+2)];\n", " const double ED0 = in_gfs[IDX4(ED0GF, i0,i1,i2)];\n", " const double ED1 = in_gfs[IDX4(ED1GF, i0,i1,i2)];\n", " const double ED2 = in_gfs[IDX4(ED2GF, i0,i1,i2)];\n", " const double AD0_i0_i1_i2m2 = in_gfs[IDX4(AD0GF, i0,i1,i2-2)];\n", " const double AD0_i0_i1_i2m1 = in_gfs[IDX4(AD0GF, i0,i1,i2-1)];\n", " const double AD0_i0_i1m2_i2 = in_gfs[IDX4(AD0GF, i0,i1-2,i2)];\n", " const double AD0_i0_i1m1_i2 = in_gfs[IDX4(AD0GF, i0,i1-1,i2)];\n", " const double AD0_i0m2_i1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2)];\n", " const double AD0_i0m1_i1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2)];\n", " const double AD0 = in_gfs[IDX4(AD0GF, i0,i1,i2)];\n", " const double AD0_i0p1_i1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2)];\n", " const double AD0_i0p2_i1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2)];\n", " const double AD0_i0_i1p1_i2 = in_gfs[IDX4(AD0GF, i0,i1+1,i2)];\n", " const double AD0_i0_i1p2_i2 = in_gfs[IDX4(AD0GF, i0,i1+2,i2)];\n", " const double AD0_i0_i1_i2p1 = in_gfs[IDX4(AD0GF, i0,i1,i2+1)];\n", " const double AD0_i0_i1_i2p2 = in_gfs[IDX4(AD0GF, i0,i1,i2+2)];\n", " const double AD1_i0_i1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1,i2-2)];\n", " const double AD1_i0_i1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1,i2-1)];\n", " const double AD1_i0_i1m2_i2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2)];\n", " const double AD1_i0_i1m1_i2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2)];\n", " const double AD1_i0m2_i1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1,i2)];\n", " const double AD1_i0m1_i1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1,i2)];\n", " const double AD1 = in_gfs[IDX4(AD1GF, i0,i1,i2)];\n", " const double AD1_i0p1_i1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1,i2)];\n", " const double AD1_i0p2_i1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1,i2)];\n", " const double AD1_i0_i1p1_i2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2)];\n", " const double AD1_i0_i1p2_i2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2)];\n", " const double AD1_i0_i1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1,i2+1)];\n", " const double AD1_i0_i1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1,i2+2)];\n", " const double AD2_i0_i1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1,i2-2)];\n", " const double AD2_i0_i1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1,i2-1)];\n", " const double AD2_i0_i1m2_i2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2)];\n", " const double AD2_i0_i1m1_i2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2)];\n", " const double AD2_i0m2_i1_i2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2)];\n", " const double AD2_i0m1_i1_i2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2)];\n", " const double AD2 = in_gfs[IDX4(AD2GF, i0,i1,i2)];\n", " const double AD2_i0p1_i1_i2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2)];\n", " const double AD2_i0p2_i1_i2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2)];\n", " const double AD2_i0_i1p1_i2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2)];\n", " const double AD2_i0_i1p2_i2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2)];\n", " const double AD2_i0_i1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1,i2+1)];\n", " const double AD2_i0_i1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1,i2+2)];\n", " const double Gamma_i0_i1_i2m2 = in_gfs[IDX4(GAMMAGF, i0,i1,i2-2)];\n", " const double Gamma_i0_i1_i2m1 = in_gfs[IDX4(GAMMAGF, i0,i1,i2-1)];\n", " const double Gamma_i0_i1m2_i2 = in_gfs[IDX4(GAMMAGF, i0,i1-2,i2)];\n", " const double Gamma_i0_i1m1_i2 = in_gfs[IDX4(GAMMAGF, i0,i1-1,i2)];\n", " const double Gamma_i0m2_i1_i2 = in_gfs[IDX4(GAMMAGF, i0-2,i1,i2)];\n", " const double Gamma_i0m1_i1_i2 = in_gfs[IDX4(GAMMAGF, i0-1,i1,i2)];\n", " const double Gamma = in_gfs[IDX4(GAMMAGF, i0,i1,i2)];\n", " const double Gamma_i0p1_i1_i2 = in_gfs[IDX4(GAMMAGF, i0+1,i1,i2)];\n", " const double Gamma_i0p2_i1_i2 = in_gfs[IDX4(GAMMAGF, i0+2,i1,i2)];\n", " const double Gamma_i0_i1p1_i2 = in_gfs[IDX4(GAMMAGF, i0,i1+1,i2)];\n", " const double Gamma_i0_i1p2_i2 = in_gfs[IDX4(GAMMAGF, i0,i1+2,i2)];\n", " const double Gamma_i0_i1_i2p1 = in_gfs[IDX4(GAMMAGF, i0,i1,i2+1)];\n", " const double Gamma_i0_i1_i2p2 = in_gfs[IDX4(GAMMAGF, i0,i1,i2+2)];\n", " const double FDPart1_Rational_2_3 = 2.0/3.0;\n", " const double FDPart1_Rational_1_12 = 1.0/12.0;\n", " const double FDPart1_Rational_5_2 = 5.0/2.0;\n", " const double FDPart1_Rational_4_3 = 4.0/3.0;\n", " const double FDPart1_1 = -AD0_i0_i1p2_i2;\n", " const double FDPart1_9 = ((invdx0)*(invdx0));\n", " const double FDPart1_10 = -AD0*FDPart1_Rational_5_2;\n", " const double FDPart1_11 = ((invdx1)*(invdx1));\n", " const double FDPart1_12 = ((invdx2)*(invdx2));\n", " const double FDPart1_13 = -AD1*FDPart1_Rational_5_2;\n", " const double FDPart1_14 = -AD2*FDPart1_Rational_5_2;\n", " const double FDPart1_18 = -FDPart1_Rational_5_2*psi;\n", " const double AD_dD00 = invdx0*(FDPart1_Rational_1_12*(AD0_i0m2_i1_i2 - AD0_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD0_i0m1_i1_i2 + AD0_i0p1_i1_i2));\n", " const double AD_dD01 = invdx1*(FDPart1_Rational_1_12*(AD0_i0_i1m2_i2 + FDPart1_1) + FDPart1_Rational_2_3*(-AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));\n", " const double AD_dD02 = invdx2*(FDPart1_Rational_1_12*(AD0_i0_i1_i2m2 - AD0_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));\n", " const double AD_dD10 = invdx0*(FDPart1_Rational_1_12*(AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));\n", " const double AD_dD11 = invdx1*(FDPart1_Rational_1_12*(AD1_i0_i1m2_i2 - AD1_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD1_i0_i1m1_i2 + AD1_i0_i1p1_i2));\n", " const double AD_dD12 = invdx2*(FDPart1_Rational_1_12*(AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));\n", " const double AD_dD20 = invdx0*(FDPart1_Rational_1_12*(AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));\n", " const double AD_dD21 = invdx1*(FDPart1_Rational_1_12*(AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));\n", " const double AD_dD22 = invdx2*(FDPart1_Rational_1_12*(AD2_i0_i1_i2m2 - AD2_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD2_i0_i1_i2m1 + AD2_i0_i1_i2p1));\n", " const double AD_dDD000 = FDPart1_9*(FDPart1_10 + FDPart1_Rational_1_12*(-AD0_i0m2_i1_i2 - AD0_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD0_i0m1_i1_i2 + AD0_i0p1_i1_i2));\n", " const double AD_dDD011 = FDPart1_11*(FDPart1_10 + FDPart1_Rational_1_12*(-AD0_i0_i1m2_i2 + FDPart1_1) + FDPart1_Rational_4_3*(AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));\n", " const double AD_dDD022 = FDPart1_12*(FDPart1_10 + FDPart1_Rational_1_12*(-AD0_i0_i1_i2m2 - AD0_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));\n", " const double AD_dDD100 = FDPart1_9*(FDPart1_13 + FDPart1_Rational_1_12*(-AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));\n", " const double AD_dDD111 = FDPart1_11*(FDPart1_13 + FDPart1_Rational_1_12*(-AD1_i0_i1m2_i2 - AD1_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD1_i0_i1m1_i2 + AD1_i0_i1p1_i2));\n", " const double AD_dDD122 = FDPart1_12*(FDPart1_13 + FDPart1_Rational_1_12*(-AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));\n", " const double AD_dDD200 = FDPart1_9*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));\n", " const double AD_dDD211 = FDPart1_11*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));\n", " const double AD_dDD222 = FDPart1_12*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0_i1_i2m2 - AD2_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD2_i0_i1_i2m1 + AD2_i0_i1_i2p1));\n", " const double Gamma_dD0 = invdx0*(FDPart1_Rational_1_12*(Gamma_i0m2_i1_i2 - Gamma_i0p2_i1_i2) + FDPart1_Rational_2_3*(-Gamma_i0m1_i1_i2 + Gamma_i0p1_i1_i2));\n", " const double Gamma_dD1 = invdx1*(FDPart1_Rational_1_12*(Gamma_i0_i1m2_i2 - Gamma_i0_i1p2_i2) + FDPart1_Rational_2_3*(-Gamma_i0_i1m1_i2 + Gamma_i0_i1p1_i2));\n", " const double Gamma_dD2 = invdx2*(FDPart1_Rational_1_12*(Gamma_i0_i1_i2m2 - Gamma_i0_i1_i2p2) + FDPart1_Rational_2_3*(-Gamma_i0_i1_i2m1 + Gamma_i0_i1_i2p1));\n", " const double psi_dD0 = invdx0*(FDPart1_Rational_1_12*(psi_i0m2_i1_i2 - psi_i0p2_i1_i2) + FDPart1_Rational_2_3*(-psi_i0m1_i1_i2 + psi_i0p1_i1_i2));\n", " const double psi_dD1 = invdx1*(FDPart1_Rational_1_12*(psi_i0_i1m2_i2 - psi_i0_i1p2_i2) + FDPart1_Rational_2_3*(-psi_i0_i1m1_i2 + psi_i0_i1p1_i2));\n", " const double psi_dD2 = invdx2*(FDPart1_Rational_1_12*(psi_i0_i1_i2m2 - psi_i0_i1_i2p2) + FDPart1_Rational_2_3*(-psi_i0_i1_i2m1 + psi_i0_i1_i2p1));\n", " const double psi_dDD00 = FDPart1_9*(FDPart1_18 + FDPart1_Rational_1_12*(-psi_i0m2_i1_i2 - psi_i0p2_i1_i2) + FDPart1_Rational_4_3*(psi_i0m1_i1_i2 + psi_i0p1_i1_i2));\n", " const double psi_dDD11 = FDPart1_11*(FDPart1_18 + FDPart1_Rational_1_12*(-psi_i0_i1m2_i2 - psi_i0_i1p2_i2) + FDPart1_Rational_4_3*(psi_i0_i1m1_i2 + psi_i0_i1p1_i2));\n", " const double psi_dDD22 = FDPart1_12*(FDPart1_18 + FDPart1_Rational_1_12*(-psi_i0_i1_i2m2 - psi_i0_i1_i2p2) + FDPart1_Rational_4_3*(psi_i0_i1_i2m1 + psi_i0_i1_i2p1));\n", " /*\n", " * NRPy+ Finite Difference Code Generation, Step 2 of 2: Evaluate SymPy expressions and write to main memory:\n", " */\n", " /*\n", " * Original SymPy expressions:\n", " * \"[rhs_gfs[IDX4(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0,\n", " * rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = -AD_dDD000 + Gamma_dD0 - (AD_dD00*xx0 - AD_dD11/xx0 + AD_dDD011 - (AD0*xx0 + AD_dD11)/xx0)/xx0**2 - (AD_dD00*xx0*sin(xx1)**2 - AD_dD22/xx0 + AD_dDD022 + (-AD1/xx0 + AD_dD01)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/xx0)/(xx0**2*sin(xx1)**2),\n", " * rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1,\n", " * rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1/xx0**2 + AD_dD10/xx0 - AD_dDD100 + Gamma_dD1 + (-AD1/xx0 + AD_dD10)/xx0 - (AD_dD01*xx0 + AD_dDD111 + xx0*(-AD1/xx0 + AD_dD01) + xx0*(-AD1/xx0 + AD_dD10))/xx0**2 - (-AD_dD22*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD122 + xx0*(-AD1/xx0 + AD_dD10)*sin(xx1)**2 + (AD0*xx0 + AD_dD11)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)*sin(2*xx1)/(2*sin(xx1)**2))/(xx0**2*sin(xx1)**2),\n", " * rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2,\n", " * rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2/xx0**2 + AD_dD20/xx0 - AD_dDD200 + Gamma_dD2 + (-AD2/xx0 + AD_dD20)/xx0 - (AD2*(-cos(2*xx1)/sin(xx1)**2 + sin(2*xx1)*cos(xx1)/sin(xx1)**3) - AD_dD21*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD211 + xx0*(-AD2/xx0 + AD_dD20) - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/(2*sin(xx1)**2))/xx0**2 - (AD_dD02*xx0*sin(xx1)**2 + AD_dD12*sin(2*xx1)/2 + AD_dDD222 + xx0*(-AD2/xx0 + AD_dD02)*sin(xx1)**2 + xx0*(-AD2/xx0 + AD_dD20)*sin(xx1)**2 + (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD12)*sin(2*xx1)/2 + (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/2)/(xx0**2*sin(xx1)**2),\n", " * rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -Gamma,\n", " * rhs_gfs[IDX4(GAMMAGF, i0, i1, i2)] = -psi_dDD00 - (psi_dD0*xx0 + psi_dDD11)/xx0**2 - (psi_dD0*xx0*sin(xx1)**2 + psi_dD1*sin(2*xx1)/2 + psi_dDD22)/(xx0**2*sin(xx1)**2)]\"\n", " */\n", " const double FDPart3_0 = (1.0/((xx0)*(xx0)));\n", " const double FDPart3_2 = (1.0/(xx0));\n", " const double FDPart3_4 = AD0*xx0 + AD_dD11;\n", " const double FDPart3_5 = sin(xx1);\n", " const double FDPart3_6 = ((FDPart3_5)*(FDPart3_5));\n", " const double FDPart3_7 = -AD1*FDPart3_2;\n", " const double FDPart3_10 = sin(2*xx1);\n", " const double FDPart3_11 = (1.0/2.0)*FDPart3_10;\n", " const double FDPart3_12 = AD0*FDPart3_6*xx0 + AD1*FDPart3_11 + AD_dD22;\n", " const double FDPart3_13 = (1.0/(FDPart3_6));\n", " const double FDPart3_14 = FDPart3_0*FDPart3_13;\n", " const double FDPart3_16 = xx0*(AD_dD10 + FDPart3_7);\n", " const double FDPart3_17 = FDPart3_11*FDPart3_13;\n", " const double FDPart3_18 = -AD2*FDPart3_2;\n", " const double FDPart3_20 = xx0*(AD_dD20 + FDPart3_18);\n", " const double FDPart3_21 = -AD2*FDPart3_17;\n", " const double FDPart3_22 = FDPart3_11*(AD_dD21 + FDPart3_21);\n", " rhs_gfs[IDX4(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0;\n", " rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = -AD_dDD000 - FDPart3_0*(AD_dD00*xx0 - AD_dD11*FDPart3_2 + AD_dDD011 - FDPart3_2*FDPart3_4) - FDPart3_14*(AD_dD00*FDPart3_6*xx0 - AD_dD22*FDPart3_2 + AD_dDD022 + FDPart3_11*(AD_dD01 + FDPart3_7) - FDPart3_12*FDPart3_2) + Gamma_dD0;\n", " rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1;\n", " rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1*FDPart3_0 + AD_dD10*FDPart3_2 - AD_dDD100 - FDPart3_0*(AD_dD01*xx0 + AD_dDD111 + FDPart3_16 + xx0*(AD_dD01 + FDPart3_7)) - FDPart3_14*(-AD_dD22*FDPart3_17 + AD_dDD122 + FDPart3_11*FDPart3_4 - FDPart3_12*FDPart3_17 + FDPart3_16*FDPart3_6) + FDPart3_2*(AD_dD10 + FDPart3_7) + Gamma_dD1;\n", " rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2;\n", " rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2*FDPart3_0 + AD_dD20*FDPart3_2 - AD_dDD200 - FDPart3_0*(AD2*(FDPart3_10*cos(xx1)/((FDPart3_5)*(FDPart3_5)*(FDPart3_5)) - FDPart3_13*cos(2*xx1)) - AD_dD21*FDPart3_17 + AD_dDD211 - FDPart3_13*FDPart3_22 + FDPart3_20) - FDPart3_14*(AD_dD02*FDPart3_6*xx0 + AD_dD12*FDPart3_11 + AD_dDD222 + FDPart3_11*(AD_dD12 + FDPart3_21) + FDPart3_20*FDPart3_6 + FDPart3_22 + FDPart3_6*xx0*(AD_dD02 + FDPart3_18)) + FDPart3_2*(AD_dD20 + FDPart3_18) + Gamma_dD2;\n", " rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -Gamma;\n", " rhs_gfs[IDX4(GAMMAGF, i0, i1, i2)] = -FDPart3_0*(psi_dD0*xx0 + psi_dDD11) - FDPart3_14*(FDPart3_11*psi_dD1 + FDPart3_6*psi_dD0*xx0 + psi_dDD22) - psi_dDD00;\n", "}\n" ] } ], "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_dHatDD = ixp.zerorank2()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " psi_dHatDD[i][j] = psi_dDD[i][j]\n", " for k in range(DIM):\n", " psi_dHatDD[i][j] += - rfm.GammahatUDD[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 += rfm.ghatUU[i][j] * psi_dHatDD[i][j]\n", "\n", "# Step 9: Define right-hand sides for the evolution.\n", "AD_rhs = ixp.zerorank1()\n", "ED_rhs = ixp.zerorank1()\n", "for i in range(DIM):\n", " AD_rhs[i] = -ED[i] - psi_dD[i]\n", " ED_rhs[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=AD_rhs[i]))\n", " lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"ED\" + str(i)), rhs=ED_rhs[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 3: 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-MaxwellCurvilinear.pdf](Tutorial-MaxwellCurvilinear.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": 3, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:14:55.238869Z", "iopub.status.busy": "2021-03-07T17:14:55.238152Z", "iopub.status.idle": "2021-03-07T17:14:58.527854Z", "shell.execute_reply": "2021-03-07T17:14:58.527060Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Created Tutorial-MaxwellCurvilinear.tex, and compiled LaTeX file to PDF\n", " file Tutorial-MaxwellCurvilinear.pdf\n" ] } ], "source": [ "import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n", "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-MaxwellCurvilinear\")" ] } ], "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 }