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.
# 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()
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:
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.
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)
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)
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}# 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])
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).
# 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:
# 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)
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}$.
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])
# \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
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}# 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]
# 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]
# 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]
# 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
Here, as a code validation check, we verify agreement in the SymPy expressions for the initial data we intend to use between
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$.
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
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.)
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