Notebook

Time Evolution of Maxwell's Equations & Constraints in Flat Spacetime and Cartesian Coordinates¶

Authors: Terrence Pierre Jacques, Zachariah Etienne and Ian Ruchlin¶

This module constructs the evolution equations for Maxwell's equations as symbolic (SymPy) expressions, for an electromagnetic field in vacuum, as defined in Tutorial-VacuumMaxwell_formulation_Cartesian.¶

Notebook Status: Validated

Validation Notes: All expressions generated in this module have been validated.

NRPy+ Source Code for this module: Maxwell/VacuumMaxwell_Flat_Evol_Cartesian.py ¶

This tutorial takes the expressions defined in Tutorial-VacuumMaxwell_formulation_Cartesian, and constructs them using NRPy+ in SymPy expressions.

Table of Contents¶

$\label{toc}$

Step 1: Initialize needed Python/NRPy+ modules
Step 2: System 1
1. Step 2.a: Construct evolution equations
2. Step 2.b: Code Validation
Step 3: System 2
1. Step 3.a: Construct evolution equations
2. Step 3.b: Code Validation
Step 4: Output this notebook to $\LaTeX$ -formatted PDF file

Step 1: Initialize needed Python/NRPy+ modules [Back to top]¶

$\label{initializenrpy}$

In [1]:

# Import needed Python modules
import NRPy_param_funcs as par   # NRPy+: Parameter interface
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import grid as gri               # NRPy+: Functions having to do with numerical grids

#Step 0: Set the spatial dimension parameter to 3.
par.set_parval_from_str("grid::DIM", 3)
DIM = par.parval_from_str("grid::DIM")

Step 2: System I [Back to top]¶

$\label{system1}$

Step 2.a: Construct evolution equations [Back to top]¶

$\label{system1_rhs}$

Following our derivations in Tutorial-VacuumMaxwell_formulation_Cartesian, we construct the system of equations for System I, in flat space and cartesian coordinates:

$\begin{align} \partial_t A^i &= -E^i - \partial^i \varphi, \\ \partial_t E^i &= \partial^i \partial_j A^j - \partial_j \partial^j A^i, \\ \partial_t \varphi &= -\partial_i A^i, \end{align}$

In [2]:

# Register gridfunctions that are needed as input.

#  Declare the rank-1 indexed expressions E_{i}, A_{i},
#  and \partial_{i} \psi, that are to be evolved in time.
#  Derivative variables like these must have an underscore
#  in them, so the finite difference module can parse
#  the variable name properly.

# E_i
EU = ixp.register_gridfunctions_for_single_rank1("EVOL", "EU")

# A_i, _AD is unused
_AU = ixp.register_gridfunctions_for_single_rank1("EVOL", "AU")

# \psi is a scalar function that is time evolved
# _psi is unused
_psi = gri.register_gridfunctions("EVOL", ["psi"])

# \partial_i \psi
psi_dD = ixp.declarerank1("psi_dD")

# \partial_k ( A_i ) --> rank two tensor
AU_dD = ixp.declarerank2("AU_dD", "nosym")

# \partial_k partial_m ( A_i ) --> rank three tensor
AU_dDD = ixp.declarerank3("AU_dDD", "sym12")

$\begin{align} \partial_t A^i &= -E^i - \partial^i \varphi, \\ \partial_t E^i &= \partial^i \partial_j A^j - \partial_j \partial^j A^i, \\ \partial_t \varphi &= -\partial_i A^i, \end{align}$

In [3]:

# \partial_t \psi = -\partial_i A_i
psi_rhs = 0

# \partial_t A_i = E_i - \partial_i \psi
ArhsU = ixp.zerorank1()

# \partial_t E_i = -\partial_j^2 A_i + \partial_j \partial_i A_j
ErhsU = ixp.zerorank1()

# RHSs - equations 10 and 11 from https://arxiv.org/abs/gr-qc/0201051
for i in range(DIM):
    ArhsU[i] = -EU[i] - psi_dD[i]
    psi_rhs += -AU_dD[i][i]
    for j in range(DIM):
        ErhsU[i] += -AU_dDD[i][j][j] + AU_dDD[j][j][i]

Constraint:

$\begin{align} \mathcal{C} \equiv \partial_i E^i \end{align}$

In [4]:

# \partial_k ( E^i ) --> rank two tensor
EU_dD = ixp.declarerank2("EU_dD", "nosym")

# C = \partial_i E^i
C = EU_dD[0][0] + EU_dD[1][1] + EU_dD[2][2]

In [5]:

print('C = ', C)
print('ArhsU[0] = ', ArhsU[0])
print('ArhsU[1] = ', ArhsU[1])
print('ArhsU[2] = ', ArhsU[2])
print('ErhsU[0] = ', ErhsU[0])
print('ErhsU[1] = ', ErhsU[1])
print('ErhsU[2] = ', ErhsU[2])
print('psi_rhs  = ', psi_rhs)

C =  EU_dD00 + EU_dD11 + EU_dD22
ArhsU[0] =  -EU0 - psi_dD0
ArhsU[1] =  -EU1 - psi_dD1
ArhsU[2] =  -EU2 - psi_dD2
ErhsU[0] =  -AU_dDD011 - AU_dDD022 + AU_dDD101 + AU_dDD202
ErhsU[1] =  AU_dDD001 - AU_dDD100 - AU_dDD122 + AU_dDD212
ErhsU[2] =  AU_dDD002 + AU_dDD112 - AU_dDD200 - AU_dDD211
psi_rhs  =  -AU_dD00 - AU_dD11 - AU_dD22

Step 2.b: NRPy+ Module Code Validation [Back to top]¶

$\label{system1_val}$

Here, as a code validation check, we verify agreement in the SymPy expressions for the RHSs of Maxwell's equations (in System I) between

this tutorial and
the NRPy+ Maxwell/VacuumMaxwell_Flat_Evol_Cartesian.py module.

In [6]:

# Reset the list of gridfunctions, as registering a gridfunction
#   twice will spawn an error.
gri.glb_gridfcs_list = []

#          Call the VacuumMaxwellRHSs() function from within the
#          Maxwell/VacuumMaxwell_Flat_Evol_Cartesian.py module,
#          which should do exactly the same as the steps above.

# Set which system to use, which are defined in Maxwell/InitialData.py
par.initialize_param(par.glb_param("char", "Maxwell.InitialData","System_to_use","System_I"))

import Maxwell.VacuumMaxwell_Flat_Evol_Cartesian as mwevol
mwevol.VacuumMaxwellRHSs()

print("Consistency check between Tutorial-VacuumMaxwell_Cartesian_RHSs tutorial and NRPy+ module: ALL SHOULD BE ZERO.")

print("C - mwevol.C = " + str(C - mwevol.C))
print("psi_rhs - mwevol.psi_rhs = " + str(psi_rhs - mwevol.psi_rhs))
for i in range(DIM):

    print("ArhsU["+str(i)+"] - mwevol.ArhsU["+str(i)+"] = " + str(ArhsU[i] - mwevol.ArhsU[i]))
    print("ErhsU["+str(i)+"] - mwevol.ErhsU["+str(i)+"] = " + str(ErhsU[i] - mwevol.ErhsU[i]))

Currently using System_I RHSs 

Consistency check between Tutorial-VacuumMaxwell_Cartesian_RHSs tutorial and NRPy+ module: ALL SHOULD BE ZERO.
C - mwevol.C = 0
psi_rhs - mwevol.psi_rhs = 0
ArhsU[0] - mwevol.ArhsU[0] = 0
ErhsU[0] - mwevol.ErhsU[0] = 0
ArhsU[1] - mwevol.ArhsU[1] = 0
ErhsU[1] - mwevol.ErhsU[1] = 0
ArhsU[2] - mwevol.ArhsU[2] = 0
ErhsU[2] - mwevol.ErhsU[2] = 0

Step 3: System II [Back to top]¶

$\label{system2}$

Step 3.a: Construct evolution equations [Back to top]¶

$\label{system2_rhs}$

Next, again following Tutorial-VacuumMaxwell_formulation_Cartesian, we construct the system of equations for System II, in flat space and cartesian coordinates:

$\begin{align} \partial_t A^i &= -E^i - \partial^i \varphi, \\ \partial_t E^i &= \partial^i \Gamma - \partial_j \partial^j A^i, \\ \partial_t \Gamma &= - \partial_i \partial^i \varphi, \\ \partial_t \varphi &= -\Gamma, \end{align}$

In [7]:

# We inherit here all of the definitions from System I, above

# Register the scalar auxiliary variable \Gamma
Gamma = gri.register_gridfunctions("EVOL", ["Gamma"])

# Declare the ordinary gradient \partial_{i} \Gamma
Gamma_dD = ixp.declarerank1("Gamma_dD")

# partial_i \partial_j \psi
psi_dDD = ixp.declarerank2("psi_dDD", "sym01")

psi_rhs = -Gamma

In [8]:

# Equation 15 https://arxiv.org/abs/gr-qc/0201051
# \partial_t \Gamma = -\partial_i^2 \psi
Gamma_rhs = 0

# \partial_t A_i = -E_i - \partial_i \psi
ArhsU = ixp.zerorank1()

# \partial_t E_i = -\partial_j^2 A_i + \partial_i \Gamma
ErhsU = ixp.zerorank1()

# RHSs - equations 10 and 14 https://arxiv.org/abs/gr-qc/0201051
for i in range(DIM):
    ArhsU[i] = -EU[i] - psi_dD[i]
    Gamma_rhs -= psi_dDD[i][i]
    ErhsU[i] += Gamma_dD[i]
    for j in range(DIM):
        ErhsU[i] -= AU_dDD[i][j][j]

Constraints:

$\begin{align} \mathcal{G} &\equiv \Gamma - \partial_i A^i, \\ \mathcal{C} &\equiv \partial_i E^i. \end{align}$

In [9]:

# \partial_k ( E^i ) --> rank two tensor
EU_dD = ixp.declarerank2("EU_dD", "nosym")

# C = \partial_i E^i
C = EU_dD[0][0] + EU_dD[1][1] + EU_dD[2][2]

# G = \Gamma - \partial_i A^i
G = Gamma - AU_dD[0][0] - AU_dD[1][1] - AU_dD[2][2]

In [10]:

print('C = ', C)
print('G = ', G)
print('ArhsU[0]  = ', ArhsU[0])
print('ArhsU[1]  = ', ArhsU[1])
print('ArhsU[2]  = ', ArhsU[2])
print('ErhsU[0]  = ', ErhsU[0])
print('ErhsU[1]  = ', ErhsU[1])
print('ErhsU[2]  = ', ErhsU[2])
print('psi_rhs   = ', psi_rhs)
print('Gamma_rhs =', Gamma_rhs)

C =  EU_dD00 + EU_dD11 + EU_dD22
G =  -AU_dD00 - AU_dD11 - AU_dD22 + Gamma
ArhsU[0]  =  -EU0 - psi_dD0
ArhsU[1]  =  -EU1 - psi_dD1
ArhsU[2]  =  -EU2 - psi_dD2
ErhsU[0]  =  -AU_dDD000 - AU_dDD011 - AU_dDD022 + Gamma_dD0
ErhsU[1]  =  -AU_dDD100 - AU_dDD111 - AU_dDD122 + Gamma_dD1
ErhsU[2]  =  -AU_dDD200 - AU_dDD211 - AU_dDD222 + Gamma_dD2
psi_rhs   =  -Gamma
Gamma_rhs = -psi_dDD00 - psi_dDD11 - psi_dDD22

Step 3.b: NRPy+ Module Code Validation [Back to top]¶

$\label{system2_val}$

Here, as a code validation check, we verify agreement in the SymPy expressions for the RHSs of Maxwell's equations (in System II) between

this tutorial and
the NRPy+ Maxwell/VacuumMaxwell_Flat_Evol_Cartesian.py module.

In [11]:

# Reset the list of gridfunctions, as registering a gridfunction
#   twice will spawn an error.
gri.glb_gridfcs_list = []

#          Call the VacuumMaxwellRHSs() function from within the
#          Maxwell/VacuumMaxwell_Flat_Evol_Cartesian.py module,
#          which should do exactly the same as the steps above.

par.set_paramsvals_value("Maxwell.InitialData::System_to_use=System_II")
mwevol.VacuumMaxwellRHSs()

print("Consistency check between Tutorial-VacuumMaxwell_Cartesian_RHSs tutorial and NRPy+ module: ALL SHOULD BE ZERO.")

print("C - mwevol.C = " + str(C - mwevol.C))
print("G - mwevol.G = " + str(G - mwevol.G))
print("psi_rhs - mwevol.psi_rhs = " + str(psi_rhs - mwevol.psi_rhs))
print("Gamma_rhs - mwevol.Gamma_rhs = " + str(Gamma_rhs - mwevol.Gamma_rhs))
for i in range(DIM):

    print("ArhsU["+str(i)+"] - mwevol.ArhsU["+str(i)+"] = " + str(ArhsU[i] - mwevol.ArhsU[i]))
    print("ErhsU["+str(i)+"] - mwevol.ErhsU["+str(i)+"] = " + str(ErhsU[i] - mwevol.ErhsU[i]))

Currently using System_II RHSs 

Consistency check between Tutorial-VacuumMaxwell_Cartesian_RHSs tutorial and NRPy+ module: ALL SHOULD BE ZERO.
C - mwevol.C = 0
G - mwevol.G = 0
psi_rhs - mwevol.psi_rhs = 0
Gamma_rhs - mwevol.Gamma_rhs = 0
ArhsU[0] - mwevol.ArhsU[0] = 0
ErhsU[0] - mwevol.ErhsU[0] = 0
ArhsU[1] - mwevol.ArhsU[1] = 0
ErhsU[1] - mwevol.ErhsU[1] = 0
ArhsU[2] - mwevol.ArhsU[2] = 0
ErhsU[2] - mwevol.ErhsU[2] = 0

Step 4: Output this notebook to $\LaTeX$ -formatted PDF file [Back to top]¶

$\label{latex_pdf_output}$

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 Tutorial-VacuumMaxwell_Cartesian_RHSs.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

In [12]:

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-VacuumMaxwell_Cartesian_RHSs")

Created Tutorial-VacuumMaxwell_Cartesian_RHSs.tex, and compiled LaTeX file
    to PDF file Tutorial-VacuumMaxwell_Cartesian_RHSs.pdf