Notebook

Time Evolution of Maxwell's Equations in Flat Spacetime and Curvilinear 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_Curvilinear.¶

Notebook Status: Validated

Validation Notes: All expressions generated in this here have been validated against the VacuumMaxwell_Flat_Evol_Curvilinear_rescaled module, as well as the Maxwell/VacuumMaxwell_Flat_Evol_Cartesian module when setting the coordinate system to Cartesian.

NRPy+ Source Code for this module: VacuumMaxwell_Flat_Evol_Curvilinear_rescaled ¶

$\label{top}$

Table of Contents:¶

Step 1: Set core NRPy+ parameters for numerical grids and reference metric
Step 2: System II in curvilinear coordinates, using the rescaled quantities $a^i$ and $e^i$
Step 3: Convert $A^i$ and $E^i$ to the Cartesian basis
Step 4: Code Validation
Step 5: Output this notebook to $\LaTeX$ -formatted PDF file

Step 1: Preliminaries [Back to top]¶

$\label{step1}$

Set up the needed NRPy+ infrastructure, such the number of dimensions and finite differencing order.

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 reference_metric as rfm   # NRPy+: Reference metric support
import grid as gri               # NRPy+: Functions having to do with numerical grids
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends

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

# Set coordinate system
# Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,
#              SymTP, SinhSymTP

CoordSystem = "Spherical"
par.set_parval_from_str("reference_metric::CoordSystem",CoordSystem)
# Set reference metric related quantities
rfm.reference_metric()

Step 2: System II in Curvilinear Coordinates, using the rescaled quantities $a^i$ and $e^i$ [Back to top]¶

$\label{step2}$

(Following discussion reproduced from Tutorial-VacuumMaxwell_formulation_Curvilinear)

Consider an arbitrary vector $\Lambda^i$ , with smooth (continous) Cartesian components $\Lambda^x$ , $\Lambda^y$ , and $\Lambda^z$ . Transforming $\Lambda^i$ to, e.g. spherical coordinates, introduces terms that spoil the smoothness of $\Lambda^i$ ;

$\Lambda^\phi = \frac{1}{r \sin \theta} \times \left[ \text{smooth part} \right].$

Evolving $\Lambda^\phi$ will introduce instabilities along the $z$ -axis. To avoid this, we instead evolve the rescaled quantity $\lambda^i$ , defined by

$\bar{\Lambda}^i = \frac{\lambda^i}{\text{scalefactor}[i]}.$

where we use the Hadamard product of matrices, and no sums are implied by the repeated indices.

Thus, we evolve the smoothed variable $\lambda^i$ , via

$\lambda^i = \bar{\Lambda}^i \text{scalefactor}[i].$

Within Nrpy+, ReU[i] = 1/scalefactor[i], giving

$\lambda^i = \frac{\bar{\Lambda}^i}{\text{ReU}[i]}.$

We now define the rescaled quantities $a^i$ and $e^i$ and rewrite our formulation of Maxwell's equations in curvilinear coordinates;

$\begin{align} a^i &= \frac{A^i}{\text{ReU}[i]},\\ \\ e^i &= \frac{E^i}{\text{ReU}[i]}, \end{align}$

Taking a time derivative on both sides,

$\begin{align} \partial_t a^i &= \frac{\partial_t A^i}{\text{ReU}[i]} = \frac{ -E^i - \hat{g}^{ij}\partial_j \varphi}{\text{ReU}[i]} = -e^i - \frac{\hat{g}^{ij}\partial_j \varphi}{\text{ReU}[i]},\\ \\ \partial_t e^i &= \frac{\partial_t E^i}{\text{ReU}[i]} = \frac{\hat{g}^{ij}\partial_j \Gamma - \hat{g}^{jk} \hat{\nabla}_{j} \left(\hat{\nabla}_{k} A^{i}\right)}{\text{ReU}[i]} = \frac{\hat{g}^{ij}\partial_j \Gamma}{\text{ReU}[i]} - \frac{\hat{g}^{jk} \hat{\nabla}_{j} \left(\hat{\nabla}_{k} a^{i} \text{ReU}[i] \right)}{\text{ReU}[i]}. \end{align}$

Given that

$\partial_t E^i = {\underbrace {\textstyle \hat{g}^{ij}\partial_j \Gamma}_{\text{Term 1}}} - \hat{\gamma}^{jk} \left({\underbrace {\textstyle A^i_{,kj}}_{\text{Term 2}}} + {\underbrace {\textstyle \hat{\Gamma}^i_{mk,j} A^m + \hat{\Gamma}^i_{mk} A^m_{,j} + \hat{\Gamma}^i_{dj} A^d_{,k} - \hat{\Gamma}^d_{kj} A^i_{,d}}_{\text{Term 3}}} + {\underbrace {\textstyle \hat{\Gamma}^i_{dj}\hat{\Gamma}^d_{mk} A^m - \hat{\Gamma}^d_{kj} \hat{\Gamma}^i_{md} A^m}_{\text{Term 4}}}\right),$

we can make the following replacements within the above, in terms of NRPy+ code;

$\begin{align} A^i = \text{AU[i]} &\to \text{aU[i] * rfm.ReU[i]} \\ \partial_j A^i = \text{AUdD[i][j]} &\to \text{aU_dD[i][j] * rfm.ReU[i]} +\text{aU[i] * rfm.ReUdD[i][j]} \\ \partial_k \partial_j A^i = \text{AUdDD[i][j][k]} &\to \text{aU_dDD[i][j][k] * rfm.ReU[i]} + \text{aU_dDD[i][j] * rfm.ReUdD[i][k]} \\ &+ \text{aU_dD[i][k] * rfm.ReUdD[i][j]} + \text{aU[i] * rfm.ReUdDD[i][j][k]} \end{align}$

The remainder of Maxwell's equations are unchanged;

$\partial_t \Gamma = -\hat{g}^{ij} \left( \partial_i \partial_j \varphi - \hat{\Gamma}^k_{ji} \partial_k \varphi \right),$

$\partial_t \varphi = -\Gamma,$

subject to constraints

$\begin{align} \mathcal{G} &\equiv \Gamma - \partial_i A^i + \hat{\Gamma}^i_{ji} A^j &= 0,\\ \mathcal{C} &\equiv \partial_i E^i + \hat{\Gamma}^i_{ji} E^j &= 0. \end{align}$

In [2]:

# Register gridfunctions that are needed as input.

#  Declare the rank-1 indexed expressions e^{i}, e^{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")

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

# a^i
aU = ixp.register_gridfunctions_for_single_rank1("EVOL", "aU")

# \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")

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

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

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

# \partial_i \Gamma
Gamma_dD = ixp.declarerank1("Gamma_dD")

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

ghatUU = rfm.ghatUU
GammahatUDD = rfm.GammahatUDD
GammahatUDDdD = rfm.GammahatUDDdD
ReU = rfm.ReU
ReUdD = rfm.ReUdD
ReUdDD = rfm.ReUdDD

$\partial_t a^i = -e^i - \frac{\hat{g}^{ij}\partial_j \varphi}{\text{ReU}[i]},$

In [3]:

# \partial_t a^i = -e^i - \frac{\hat{g}^{ij}\partial_j \varphi}{\text{ReU}[i]}
arhsU = ixp.zerorank1()
for i in range(DIM):
    arhsU[i] -= eU[i]
    for j in range(DIM):
        arhsU[i] -= (ghatUU[i][j]*psi_dD[j])/ReU[i]

$\partial_t e^i = \frac{\hat{g}^{ij}\partial_j \Gamma - \hat{g}^{jk} \hat{\nabla}_{j} \left(\hat{\nabla}_{k} A^{i}\right)}{\text{ReU}[i]} = \frac{\hat{g}^{ij}\partial_j \Gamma}{\text{ReU}[i]} - \frac{\hat{g}^{jk} \hat{\nabla}_{j} \left(\hat{\nabla}_{k} a^{i} \text{ReU}[i] \right)}{\text{ReU}[i]}.$

Given that

we can make the following replacements within the above, in terms of NRPy+ code;

$\begin{align} A^i = \text{AU[i]} &\to \text{aU[i] * rfm.ReU[i]} \\ \partial_j A^i = \text{AUdD[i][j]} &\to \text{aU_dD[i][j] * rfm.ReU[i]} +\text{aU[i] * rfm.ReUdD[i][j]} \\ \partial_k \partial_j A^i = \text{AUdDD[i][j][k]} &\to \text{aU_dDD[i][j][k] * rfm.ReU[i]} + \text{aU_dD[i][j] * rfm.ReUdD[i][k]} \\ &+ \text{aU_dD[i][k] * rfm.ReUdD[i][j]} + \text{aU[i] * rfm.ReUdDD[i][j][k]} \end{align}$

In [4]:

# A^i
AU = ixp.zerorank1()

# \partial_k ( A^i ) --> rank two tensor
AU_dD = ixp.zerorank2()

# \partial_k partial_m ( A^i ) --> rank three tensor
AU_dDD = ixp.zerorank3()

for i in range(DIM):
    AU[i] = aU[i]*ReU[i]
    for j in range(DIM):
        AU_dD[i][j] = aU_dD[i][j]*ReU[i] + aU[i]*ReUdD[i][j]
        for k in range(DIM):
            AU_dDD[i][j][k] = aU_dDD[i][j][k]*ReU[i] + aU_dD[i][j]*ReUdD[i][k] +\
                              aU_dD[i][k]*ReUdD[i][j] + aU[i]*ReUdDD[i][j][k]

$\text{Term 1} = \hat{g}^{ij}\partial_j \Gamma$

In [5]:

# Term 1 = \hat{g}^{ij}\partial_j \Gamma
Term1U = ixp.zerorank1()
for i in range(DIM):
    for j in range(DIM):
        Term1U[i] += ghatUU[i][j]*Gamma_dD[j]

$\text{Term 2} = A^i_{,kj}$

In [6]:

# Term 2: A^i_{,kj}
Term2UDD = ixp.zerorank3()
for i in range(DIM):
    for j in range(DIM):
        for k in range(DIM):
            Term2UDD[i][j][k] += AU_dDD[i][k][j]

$\text{Term 3} = \hat{\Gamma}^i_{mk,j} A^m + \hat{\Gamma}^i_{mk} A^m_{,j} + \hat{\Gamma}^i_{dj}A^d_{,k} - \hat{\Gamma}^d_{kj} A^i_{,d}$

In [7]:

# Term 3: \hat{\Gamma}^i_{mk,j} A^m + \hat{\Gamma}^i_{mk} A^m_{,j}
#        + \hat{\Gamma}^i_{dj}A^d_{,k} - \hat{\Gamma}^d_{kj} A^i_{,d}
Term3UDD = ixp.zerorank3()
for i in range(DIM):
    for j in range(DIM):
        for k in range(DIM):
            for m in range(DIM):
                Term3UDD[i][j][k] += GammahatUDDdD[i][m][k][j]*AU[m]    \
                                      + GammahatUDD[i][m][k]*AU_dD[m][j] \
                                      + GammahatUDD[i][m][j]*AU_dD[m][k] \
                                      - GammahatUDD[m][k][j]*AU_dD[i][m]

$\text{Term 4} = \hat{\Gamma}^i_{dj}\hat{\Gamma}^d_{mk} A^m - \hat{\Gamma}^d_{kj} \hat{\Gamma}^i_{md} A^m$

In [8]:

# Term 4: \hat{\Gamma}^i_{dj}\hat{\Gamma}^d_{mk} A^m -
#         \hat{\Gamma}^d_{kj} \hat{\Gamma}^i_{md} A^m
Term4UDD = ixp.zerorank3()
for i in range(DIM):
    for j in range(DIM):
        for k in range(DIM):
            for m in range(DIM):
                for d in range(DIM):
                    Term4UDD[i][j][k] += ( GammahatUDD[i][d][j]*GammahatUDD[d][m][k] \
                                           -GammahatUDD[d][k][j]*GammahatUDD[i][m][d])*AU[m]

Finally, we build up the RHS of $E^i$ ,

and divide through by ReU[i] to get $e^i$ .

In [9]:

# \partial_t E^i = \hat{g}^{ij}\partial_j \Gamma - \hat{\gamma}^{jk}*
#    (A^i_{,kj}
#   + \hat{\Gamma}^i_{mk,j} A^m + \hat{\Gamma}^i_{mk} A^m_{,j}
#   + \hat{\Gamma}^i_{dj} A^d_{,k} - \hat{\Gamma}^d_{kj} A^i_{,d}
#   + \hat{\Gamma}^i_{dj}\hat{\Gamma}^d_{mk} A^m
#   - \hat{\Gamma}^d_{kj} \hat{\Gamma}^i_{md} A^m)

ErhsU = ixp.zerorank1()
for i in range(DIM):
    ErhsU[i] += Term1U[i]
    for j in range(DIM):
        for k in range(DIM):
            ErhsU[i] -= ghatUU[j][k]*(Term2UDD[i][j][k] + Term3UDD[i][j][k] + Term4UDD[i][j][k])

erhsU = ixp.zerorank1()
for i in range(DIM):
    erhsU[i] = ErhsU[i]/ReU[i]

$\partial_t \Gamma = -\hat{g}^{ij} \left( \partial_i \partial_j \varphi - \hat{\Gamma}^k_{ji} \partial_k \varphi \right)$

In [10]:

# \partial_t \Gamma = -\hat{g}^{ij} (\partial_i \partial_j \varphi -
# \hat{\Gamma}^k_{ji} \partial_k \varphi)
Gamma_rhs = sp.sympify(0)
for i in range(DIM):
    for j in range(DIM):
        Gamma_rhs -= ghatUU[i][j]*psi_dDD[i][j]
        for k in range(DIM):
            Gamma_rhs += ghatUU[i][j]*GammahatUDD[k][j][i]*psi_dD[k]

$\partial_t \varphi = -\Gamma$

In [11]:

# \partial_t \varphi = -\Gamma
psi_rhs = -Gamma

Constraints:

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

In [12]:

# \mathcal{G} \equiv \Gamma - \partial_i A^i + \hat{\Gamma}^i_{ji} A^j
G = Gamma
for i in range(DIM):
    G -= AU_dD[i][i]
    for j in range(DIM):
        G += GammahatUDD[i][j][i]*AU[j]

# E^i
EU = ixp.zerorank1()

# \partial_k ( A^i ) --> rank two tensor
EU_dD = ixp.zerorank2()
for i in range(DIM):
    EU[i] = eU[i]*ReU[i]
    for j in range(DIM):
        EU_dD[i][j] = eU_dD[i][j]*ReU[i] + eU[i]*ReUdD[i][j]

C = sp.sympify(0)
for i in range(DIM):
    C += EU_dD[i][i]
    for j in range(DIM):
        C += GammahatUDD[i][j][i]*EU[j]

Step 3: Convert $A^i$ and $E^i$ to the Cartesian basis [Back to top]¶

$\label{cart_transform}$

Here we convert $A^i$ and $E^i$ to the Cartesian basis, to make convergence tests within Tutorial-Start_to_Finish-Solving_Maxwells_Equations_in_Vacuum-Curvilinear easier. Specifically, we will use the coordinate transformation definitions provided by reference_metric.py to build the Jacobian:

$\begin{align} \frac{\partial x_{\rm Cart}^i}{\partial x_{\rm Orig}^j}, \end{align}$

where $x_{\rm Cart}^i \in \{x,y,z\}$ . We then apply it to $A^i$ and $E^i$ to transform into Cartesian coordinates, via

$\begin{align} A^i_{\rm Cart} = \frac{\partial x_{\rm Cart}^i}{\partial x_{\rm Orig}^j} A^j_{\rm Orig}. \end{align}$

In [13]:

def Convert_to_Cartesian_basis(VU):
    # Coordinate transformation from original basis to Cartesian
    rfm.reference_metric()

    VU_Cart = ixp.zerorank1()
    Jac_dxCartU_dxOrigD = ixp.zerorank2()
    for i in range(DIM):
        for j in range(DIM):
            Jac_dxCartU_dxOrigD[i][j] = sp.diff(rfm.xx_to_Cart[i], rfm.xx[j])

    for i in range(DIM):
        for j in range(DIM):
            VU_Cart[i] += Jac_dxCartU_dxOrigD[i][j]*VU[j]
    return VU_Cart

AU_Cart = Convert_to_Cartesian_basis(AU)
EU_Cart = Convert_to_Cartesian_basis(EU)

Step 4: NRPy+ Module Code Validation [Back to top]¶

$\label{step4}$

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

this tutorial and
the NRPy+ VacuumMaxwell_Flat_Evol_Curvilinear_rescaled module.

In [14]:

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

#          Call the VacuumMaxwellRHSs_rescaled() function from within the
#          Maxwell/VacuumMaxwell_Flat_Evol_Curvilinear_rescaled.py module,
#          which should do exactly the same as the 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_II"))

import Maxwell.VacuumMaxwell_Flat_Evol_Curvilinear_rescaled as mwevol
mwevol.VacuumMaxwellRHSs_rescaled()

print("Consistency check between Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling 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]))
    print("AU_Cart["+str(i)+"] - mwevol.AU_Cart["+str(i)+"] = " + str(AU_Cart[i] - mwevol.AU_Cart[i]))
    print("EU_Cart["+str(i)+"] - mwevol.EU_Cart["+str(i)+"] = " + str(EU_Cart[i] - mwevol.EU_Cart[i]))

Consistency check between Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling 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
AU_Cart[0] - mwevol.AU_Cart[0] = 0
EU_Cart[0] - mwevol.EU_Cart[0] = 0
arhsU[1] - mwevol.arhsU[1] = 0
erhsU[1] - mwevol.erhsU[1] = 0
AU_Cart[1] - mwevol.AU_Cart[1] = 0
EU_Cart[1] - mwevol.EU_Cart[1] = 0
arhsU[2] - mwevol.arhsU[2] = 0
erhsU[2] - mwevol.erhsU[2] = 0
AU_Cart[2] - mwevol.AU_Cart[2] = 0
EU_Cart[2] - mwevol.EU_Cart[2] = 0

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

In [15]:

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

[NbConvertApp] WARNING | pattern 'Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling.ipynb' matched no files
Created Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling.tex, and compiled
    LaTeX file to PDF file Tutorial-VacuumMaxwell_Curvilinear_RHS-
    Rescaling.pdf