Notebook

Initial Data for Solving Maxwell's Equations in Flat Spacetime¶

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

This module constructs the initial data for Maxwell's equations as symbolic (SymPy) expressions, for a purely toriodal dipole field, as defined in Knapp, Walker & Baumgarte (2002).¶

Notebook Status: Validated

Validation Notes: All expressions generated in this module have been validated, against the Dendro code Maxwell initial data, and have satisfied the contraints given in Tutorial-VacuumMaxwell_Curvilinear_RHS-Rescaling, as well as the wave equation for the electric field and the vector potential.

NRPy+ Source Code for this module: Maxwell/InitialData.py, reference_metric.py ¶

Table of Contents¶

$\label{toc}$

Step 1: Initialize needed Python/NRPy+ modules and set destination basis
Step 2: A Purely Toriodal Dipole Field
1. Step 2.a: Converting to Cartesian Basis
2. Step 2.b: Converting to Destination Basis
Step 3: Checks
1. Step 3.a: Lorentz Gauge Condition & Divergence Constraint
2. Step 3.b: Check that $A^i$ satisfies the wave equation
Step 4: Code Validation
Step 5: Output this notebook to $\LaTeX$ -formatted PDF file

Step 1: Initialize needed Python/NRPy+ modules and set destination basis [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 reference_metric as rfm   # NRPy+: Reference metric support
import sympy as sp                # SymPy: The Python computer algebra package upon which NRPy+ depends

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

# Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical,
#              SymTP, SinhSymTP
dst_basis = "Cylindrical"

# To help with simplifications, we tell Sympy that
#  the coordinate xx0 is radial like (positive real)
radial_like_dst_xx0 = True

# Set coordinate system to Cartesian
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()

Step 2: A Purely Toriodal Dipole Field [Back to top]¶

$\label{step2}$

Having the evolution equations from Knapp, Walker & Baumgarte (2002) written in Tutorial-VacuumMaxwell_Cartesian_RHSs, we must construct the initial data that will then be time evolved. Beginning from the analytic solution to this system of equation given by equation 16 of Knapp, Walker & Baumgarte (2002),

$\begin{align} 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), \\ \end{align}$

where $A^{\hat{\phi}} = A^{\phi} r\sin\theta$ , $\mathcal{A}$ gives the amplitude, $\lambda$ describes the size of the wavepacket, $u = t+r$ , and $v = t-r$ . Other components of the vector potential are $0$ . Note that these expressions repesent the exact solution to both systems of equations at any time $t \geq 0$ , at all points on our numerical grid. Thus, to get initial data we set $t=0$ .

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$ .

We may calculate $E^i$ using

$\begin{align} E^i = -\partial_t A^i. \end{align}$

Inputs for initial data:

amp - $A$
lam - $\lambda$
time - $t$

Below we define the Cartesian coordinates $x, y$ and $z$ . We then define the vector potential $A^i$ in spherical coordinates, but each component is written in terms of Cartesian coordinates. This makes the subsequent basis changes easier.

In [2]:

x = rfm.xx_to_Cart[0]
y = rfm.xx_to_Cart[1]
z = rfm.xx_to_Cart[2]

# Step 1: Declare free parameters intrinsic to these initial data
# Amplitude
amp = par.Cparameters("REAL",__name__,"amp", default_vals=1.0)
# lambda
lam = par.Cparameters("REAL",__name__,"lam", default_vals=1.0)
time = par.Cparameters("REAL",__name__,"time", default_vals=0.0)
wavespeed = par.Cparameters("REAL",__name__,"wavespeed", default_vals=1.0)

psi_ID = sp.sympify(0)

Gamma_ID = sp.sympify(0)

$\begin{align} 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), \\ A^{\phi} &= \frac{A^{\hat{\phi}}} {r\sin\theta} \end{align}$

In [3]:

AidU_Sph = ixp.zerorank1()
# Set coordinate transformations:
r = sp.sqrt(x*x + y*y + z*z)
sin_theta = z / r

u = time + r
v = time - r
e_lam_u = sp.exp(-lam*u**2)
e_lam_v = sp.exp(-lam*v**2)

# Equation 16 from https://arxiv.org/abs/gr-qc/0201051
AU_phi_hat = (amp*sin_theta)*( ((e_lam_v - e_lam_u)/r**2) - \
                        2*lam*(v*e_lam_v + u*e_lam_u)/r )

AidU_Sph[2] = AU_phi_hat/(r*sin_theta)

Step 2.a: Converting to Cartesian Basis [Back to top]¶

$\label{cart_basis}$

Note that $A^i$ is defined in sperical coordinates, so we must therefore transform to Cartesian coordinates using the Jacobian. Here 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 Sph}^j}, \end{align}$

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

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

In [4]:

# Coordinate transformation from spherical to Cartesian
AidU_Cart = ixp.zerorank1()
Jac_dxSphU_dxCartD = ixp.zerorank2()
for i in range(DIM):
    for j in range(DIM):
        Jac_dxSphU_dxCartD[i][j] = sp.diff(rfm.xxSph[i],rfm.xx_to_Cart[j])

#         Jac_dxCartU_dxSphD[i][j] = sp.diff(rfm.xx_to_Cart[i],rfm.xx[j])
Jac_dxCartU_dxSphD,dummy = ixp.generic_matrix_inverter3x3(Jac_dxSphU_dxCartD)

for i in range(DIM):
    for j in range(DIM):
        AidU_Cart[i] += Jac_dxCartU_dxSphD[i][j]*AidU_Sph[j]
for i in range(DIM):
    AidU_Cart[i] = sp.simplify(AidU_Cart[i])

Step 2.b: Converting to Destination Basis [Back to top]¶

$\label{dst_basis}$

Here we prepare to convert $A^i$ from the Cartesian basis to the destination basis. To do so, we first rewrite each component of $A^i$ in terms of the destination coordinates. This is done by first re-labelling the NRPy+ coordinates $xx0, xx1, xx2$ as $cart_{xx0}, cart_{xx1}, cart_{xx2}$ . Then, each $cart_{xxi}$ is replaced by its counterpart expression in the destination basis using reference_metric.py.

Note that for algebraic simplification, we tell sympy that the coordinate $xx0$ is radial like and thus positive and real (if the destination coordinates are curvilinear).

In [5]:

# rfm is still defined in Cartesian coordinates
cart_xx = ixp.declarerank1("cart_xx")
for i in range(DIM):
    for k in range(DIM):
        AidU_Cart[i] = AidU_Cart[i].subs(rfm.xx[k], cart_xx[k])

# Set coordinate system to dst_basis
par.set_parval_from_str("reference_metric::CoordSystem",dst_basis)
rfm.reference_metric()

for i in range(DIM):
    for k in range(DIM):
        AidU_Cart[i] = AidU_Cart[i].subs(cart_xx[k], rfm.xx_to_Cart[k])

if radial_like_dst_xx0:
    for j in range(DIM):
        AidU_Cart[j] =  sp.refine(sp.simplify(AidU_Cart[j]), sp.Q.positive(rfm.xx[0]))

We define Jacobians relative to the center of the destination grid, at a point $x^j_{\rm dst}=$ (xx0,xx1,xx2) ${}_{\rm dst}$ on the destination grid:

${\rm Jac\_dUCart\_dDdstUD[i][j]} = \frac{\partial x^i_{\rm Cart}}{\partial x^j_{\rm dst}},$

via exact differentiation (courtesy SymPy), and the inverse Jacobian

${\rm Jac\_dUdst\_dDCartUD[i][j]} = \frac{\partial x^i_{\rm dst}}{\partial x^j_{\rm Cart}},$

using NRPy+'s generic_matrix_inverter3x3() function. In terms of these, the transformation of BSSN tensors from Cartesian to the destination grid's "reference_metric::CoordSystem" coordinates may be written:

$A^i_{\rm dst} = \frac{\partial x^i_{\rm dst}}{\partial x^\ell_{\rm Cart}} A^\ell_{\rm Cart}$

In [6]:

# Step 3: Transform BSSN tensors in Cartesian basis to destination grid basis, using center of dest. grid as origin

# Step 3.a: Next construct Jacobian and inverse Jacobian matrices:
Jac_dUCart_dDrfmUD,Jac_dUrfm_dDCartUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Cartesian()

# Step 3.b: Convert basis of all BSSN *vectors* from Cartesian to destination basis
AidU = rfm.basis_transform_vectorU_from_Cartesian_to_rfmbasis(Jac_dUrfm_dDCartUD, AidU_Cart)

# Define electric field --> E^i = -\partial_t A^i
EidU = ixp.zerorank1()
for j in range(DIM):
    EidU[j] = -sp.diff(AidU[j], time)

Step 3: Checks [Back to top]¶

$\label{step3}$

Here we validate the initial data. Specifically, we check that the constraints from Tutorial-VacuumMaxwell_formulation_Curvilinear are satisfied;

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

Note that the above simply to their usual forms in Cartesian coordinates.

Finally, we check that $A^i$ satisfies the covariant wave equation,

$\begin{align} \partial_t^2 A^i - \hat{\gamma}^{jk} \hat{\nabla}_j \hat{\nabla}_k A^i = 0, \end{align}$

where $\hat{\nabla}_j$ is the covariant derivative associated with the spatial metric $\hat{\gamma}_{jk}$ .

In [7]:

AidU_dD = ixp.zerorank2()
for i in range(DIM):
    for j in range(DIM):
        AidU_dD[i][j] += sp.diff(AidU[i], rfm.xx[j])

AidU_dDD = ixp.zerorank3()
for i in range(DIM):
    for j in range(DIM):
        for k in range(DIM):
            AidU_dDD[i][j][k] += sp.diff(AidU[i], rfm.xx[j], rfm.xx[k])

Step 3.a: Lorentz Gauge Condition & Divergence Constraint [Back to top]¶

$\label{lorentz}$

In [8]:

# \mathcal{G} \equiv \Gamma - \partial_i A^i - \hat{\Gamma}^i_{ji} A^j
G = Gamma_ID
for i in range(DIM):
    G -= AidU_dD[i][i]
    for j in range(DIM):
        G -= rfm.GammahatUDD[i][j][i]*AidU[j]

print('G should evaluate to zero:', sp.simplify(G), '\n')

# \mathcal{C} \equiv \partial_i E^i + \hat{\Gamma}^i_{ji} E^j
C = sp.sympify(0)
for i in range(DIM):
    C += sp.diff(EidU[i], rfm.xx[i], 1)
    for j in range(DIM):
        C += rfm.GammahatUDD[i][j][i]*EidU[j]
print('C should evaluate to zero:', sp.simplify(C))

G should evaluate to zero: 0 

C should evaluate to zero: 0

Step 3.b: Check that $A^i$ satisfies the wave equation [Back to top]¶

$\label{wave_eq}$

Based on the definition of covariant derivative, we have

$\hat{\nabla}_{k} A^{i} = A^i_{,k} + \hat{\Gamma}^i_{mk} A^m$

Since $\hat{\nabla}_{k} A^{i}$ is a tensor, the covariant derivative of this will have the same indexing as a tensor $T_k^i$ :

$\hat{\nabla}_{j} T^i_k = T^i_{k,j} + \hat{\Gamma}^i_{dj} T^d_k - \hat{\Gamma}^d_{kj} T^i_d.$

Therefore,

$\begin{align} \hat{\nabla}_{j} \left(\hat{\nabla}_{k} A^{i}\right) &= \left(A^i_{,k} + \hat{\Gamma}^i_{mk} A^m\right)_{,j} + \hat{\Gamma}^i_{dj} \left(A^d_{,k} + \hat{\Gamma}^d_{mk} A^m\right) - \hat{\Gamma}^d_{kj} \left(A^i_{,d} + \hat{\Gamma}^i_{md} A^m\right) \\ &= A^i_{,kj} + \hat{\Gamma}^i_{mk,j} A^m + \hat{\Gamma}^i_{mk} A^m_{,j} + \hat{\Gamma}^i_{dj}A^d_{,k} + \hat{\Gamma}^i_{dj}\hat{\Gamma}^d_{mk} A^m - \hat{\Gamma}^d_{kj} A^i_{,d} - \hat{\Gamma}^d_{kj} \hat{\Gamma}^i_{md} A^m \\ &= {\underbrace {\textstyle A^i_{,kj}}_{\text{Term 1}}}+ {\underbrace {\textstyle \hat{\Gamma}^i_{mk,j} A^m + \hat{\Gamma}^i_{mk} A^m_{,j} + \hat{\Gamma}^i_{dj} A^d_{,k} - \hat{\Gamma}^d_{kj} A^i_{,d}}_{\text{Term 2}}} + {\underbrace {\textstyle \hat{\Gamma}^i_{dj}\hat{\Gamma}^d_{mk} A^m - \hat{\Gamma}^d_{kj} \hat{\Gamma}^i_{md} A^m}_{\text{Term 3}}}. \end{align}$

Thus

$\hat{\gamma}^{jk} \hat{\nabla}_{j} \left(\hat{\nabla}_{k} A^{i}\right) = \hat{\gamma}^{jk} \left(\text{Term 1} + \text{Term 2} + \text{Term 3}\right).$

We use the above to confirm

$\begin{align} \partial_t^2 A^i - \hat{\gamma}^{jk} \hat{\nabla}_j \hat{\nabla}_k A^i = 0, \end{align}$

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

In [9]:

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

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

In [10]:

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

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

In [11]:

# Term 3: \hat{\Gamma}^i_{dj}\hat{\Gamma}^d_{mk} A^m -
#         \hat{\Gamma}^d_{kj} \hat{\Gamma}^i_{md} A^m
Term3UDD = 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):
                    Term3UDD[i][j][k] += ( rfm.GammahatUDD[i][d][j]*rfm.GammahatUDD[d][m][k] \
                                           -rfm.GammahatUDD[d][k][j]*rfm.GammahatUDD[i][m][d])*AidU[m]

$\hat{\gamma}^{jk} \hat{\nabla}_{j} \left(\hat{\nabla}_{k} A^{i}\right) = \hat{\gamma}^{jk} \left(\text{Term 1} + \text{Term 2} + \text{Term 3}\right),$

$\partial_t^2 A^i - \hat{\gamma}^{jk} \hat{\nabla}_j \hat{\nabla}_k A^i = 0$

In [12]:

# A^i_{,kj} + \hat{\Gamma}^i_{mk,j} A^m + \hat{\Gamma}^i_{mk} A^m_{,j} +
# \hat{\Gamma}^i_{dj}A^d_{,k} + \hat{\Gamma}^i_{dj}\hat{\Gamma}^d_{mk} A^m -
# \hat{\Gamma}^d_{kj} A^i_{,d} - \hat{\Gamma}^d_{kj} \hat{\Gamma}^i_{md} A^m

Difference = ixp.zerorank1()
for i in range(DIM):
    Difference[i] = sp.diff(AidU[i], time, 2)
    for j in range(DIM):
        for k in range(DIM):
            Difference[i] += -rfm.ghatUU[k][j]*(Term1UDD[i][j][k] + Term2UDD[i][j][k] + Term3UDD[i][j][k])

for i in range(DIM):
    print(str(i)+"th component of A-field equation. Should be zero (takes a bit, please be patient): ")
    print("                "+str(sp.simplify(Difference[i])))

0th component of A-field equation. Should be zero (takes a bit, please be patient): 
                0
1th component of A-field equation. Should be zero (takes a bit, please be patient): 
                0
2th component of A-field equation. Should be zero (takes a bit, please be patient): 
                0

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 initial data we intend to use between

this tutorial and
the NRPy+ InitialData module.

Since the initial data is identical between the two systems for $E^i, A^i$ , and $\psi$ , we also set and validate initial data for $\Gamma$ .

In [13]:

import Maxwell.InitialData as mwid
par.set_parval_from_str("Maxwell.InitialData::System_to_use","System_II")

mwid.InitialData()

# Again, to help sympy with simplifications
if radial_like_dst_xx0:
    for j in range(DIM):
        mwid.AidU[j] =  sp.refine(sp.simplify(mwid.AidU[j]), sp.Q.positive(rfm.xx[0]))
        mwid.EidU[j] =  sp.refine(sp.simplify(mwid.EidU[j]), sp.Q.positive(rfm.xx[0]))

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

print("psi_ID - mwid.psi_ID = " + str(sp.simplify(psi_ID) - mwid.psi_ID))
print("Gamma_ID - mwid.Gamma_ID = " + str(Gamma_ID - mwid.Gamma_ID))

for i in range(DIM):
    print("AidU["+str(i)+"] - mwid.AidU["+str(i)+"] = " + str(sp.simplify(AidU[i] - mwid.AidU[i])))
    print("EidU["+str(i)+"] - mwid.EidU["+str(i)+"] = " + str(sp.simplify(EidU[i] - mwid.EidU[i])))

Currently using System_II initial data
Consistency check between this tutorial and NRPy+ module: ALL SHOULD BE ZERO.
psi_ID - mwid.psi_ID = 0
Gamma_ID - mwid.Gamma_ID = 0
AidU[0] - mwid.AidU[0] = 0
EidU[0] - mwid.EidU[0] = 0
AidU[1] - mwid.AidU[1] = 0
EidU[1] - mwid.EidU[1] = 0
AidU[2] - mwid.AidU[2] = 0
EidU[2] - mwid.EidU[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_InitialData.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

In [14]:

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

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