Notebook Status: Validated
Validation Notes: The initial data generated by the NRPy+ module corresponding to this tutorial notebook are used shown to satisfy Einstein's equations as expected in this tutorial notebook.
# Step 0: Load all needed Python/NRPy+ modules
import os,sys,shutil # Standard Python modules for multiplatform OS-level functions
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import numpy as np # NumPy: A large collection of mathematical functions for Python
from scipy.sparse import spdiags # SciPy: Sparse, tri-diagonal matrix setup function
from scipy.sparse import csc_matrix # SciPy: Sparse matrix optimization function
from scipy.sparse.linalg import spsolve # SciPy: Solver of linear systems involving sparse matrices
import outputC as outC # NRPy+: Core C code output module
import reference_metric as rfm # NRPy+: Reference metric support
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
# Step 0.a: Create the output directory
Ccodesdir = "ScalarFieldID_validation"
shutil.rmtree(Ccodesdir,ignore_errors=True)
cmd.mkdir(Ccodesdir)
In this section we will set up time symmetric initial data for the gravitational collapse of a massless scalar field, in spherical coordinates. Our discussion will follow closely section III.A of Akbarian & Choptuik (2015) (henceforth A&C). We will be using a uniform radial sampling. All initial data quantities will be written with tildes over them, meaning that, for example, ˜α≡α(0,r).
We are here considering a spherically symmetric problem, so that f=f(t,r), for every function discussed in this tutorial. The demand for time-symmetric initial data then imples that
˜Kij=0 ,˜K=0 ,˜βi=0 ,˜Bi=0 .For the scalar field, φ, it also demands
∂tφ(0,r)=0 ,which we discuss below.
We will be implementing the following options for the initial profile of the scalar field
˜φI=φ0exp(−r2σ2) ,˜φII=φ0r3exp[−(r−r0σ)2] ,˜φIII=φ0{1−tanh[(r−r0σ)2]}.We introduce the two auxiliary fields
˜Φ≡∂r˜φandΠ≡−1α(∂tφ−βi∂iφ) ,of which ˜Φ will only be used as an auxiliary variable for setting the initial data, but Π is a dynamical variable which will be evolved in time. Because we are setting time-symmetric initial data, ∂t=0=βi, and thus ˜Π=0.
# Step 1: Setting up time-symmetric initial data
# Step 1.a: Define basic parameters
# Step 1.a.i: Domain size
RMAX = 50
# Step 1.a.ii: Number of gridpoints in the radial direction
NR = 30000
# Step 1.a.iii: Initial data family. Available options are:
# Gaussian_pulse, Gaussian_pulsev2, and Tanh_pulse
ID_Family = "Gaussian_pulsev2"
# Step 1.a.iv: Coordinate system. Available options are:
# Spherical and SinhSpherical
CoordSystem = "Spherical"
# Step 1.a.v: SinhSpherical parameters
sinhA = RMAX
sinhW = 0.1
# Step 1.b: Set the radial array
if CoordSystem == "Spherical":
r = np.linspace(0,RMAX,NR+1) # Set the r array
dr = np.zeros(NR)
for i in range(NR):
dr[i] = r[1]-r[0]
r = np.delete(r-dr[0]/2,0) # Shift the vector by -dr/2 and remove the negative entry
elif CoordSystem == "SinhSpherical":
if sinhA is None or sinhW is None:
print("Error: SinhSpherical coordinates require initialization of both sinhA and sinhW")
sys.exit(1)
else:
x = np.linspace(0,1.0,NR+1)
dx = 1.0/(NR+1)
x = np.delete(x-dx/2,0) # Shift the vector by -dx/2 and remove the negative entry
r = sinhA * np.sinh( x/sinhW ) / np.sinh( 1.0/sinhW )
dr = sinhA * np.cosh( x/sinhW ) / np.sinh( 1.0/sinhW ) * dx
else:
print("Error: Unknown coordinate system")
sys.exit(1)
# Step 1.c: Step size squared
dr2 = dr**2
# Step 1.d: Set SymPy variables for the initial condition
phi0,rr,rr0,sigma = sp.symbols("phi0 rr rr0 sigma",real=True)
# Step 1.e: Now set the initial profile of the scalar field
if ID_Family == "Gaussian_pulse":
phiID = phi0 * sp.exp( -r**2/sigma**2 )
elif ID_Family == "Gaussian_pulsev2":
phiID = phi0 * rr**3 * sp.exp( -(rr-rr0)**2/sigma**2 )
elif ID_Family == "Tanh_pulse":
phiID = phi0 * ( 1 - sp.tanh( (rr-rr0)**2/sigma**2 ) )
else:
print("Unkown initial data family: ",ID_Family)
print("Available options are: Gaussian_pulse, Gaussian_pulsev2, and Tanh_pulse")
sys.exit(1)
# Step 1.f: Compute Phi := \partial_{r}phi
PhiID = sp.diff(phiID,rr)
# Step 1.g: Generate NumPy functions for phi
# and Phi from the SymPy variables.
phi = sp.lambdify((phi0,rr,rr0,sigma),phiID)
Phi = sp.lambdify((phi0,rr,rr0,sigma),PhiID)
# Step 1.h: populating the varphi(0,r) array
phi0 = 0.1
r0 = 0
sigma = 1
ID_sf = phi(phi0,r,r0,sigma)
To set up the physical metric initial data, ˜γij, we will start by considering the conformal transformation
γij=e4ϕˉγij ,where ˉγij is the conformal metric and eϕ is the conformal factor. We then fix the initial value of ˉγij according to eqs. (32) and (43) of A&C
ˉγij=ˆγij ,where ˆγij is the reference metric, which is the flat metric in spherical symmetry
ˆγij=(1000r2000r2sin2θ) .To determine the physical metric, we must then determine the conformal factor eϕ. This is done by solving the Hamiltonian constraint (cf. eq. (12) of Baumgarte)
ˆγijˆDiˆDjψ=−2πψ5ρ ,where ψ≡e˜ϕ. For a massless scalar field, we know that
Tμν=∂μφ∂νφ−12gμν(∂λφ∂λφ) .where gμν is the inverse of the ADM 4-metric given by eq. (2.119) of Baumgarte & Shapiro's Numerical Relativity,
gμν=(−α−2α−2βiα−2βjγij−α−2βiβj) .We know that (see Step 2 in this tutorial module for the details)
∂tφ=α−1Π ,∂λφ∂λφ=−Π2+γij∂iφ∂jφ .The tt-component of the energy-momentum tensor at the initial time is then given by (we will ommit the "tildes" below to avoid cluttering the equation, but keep in mind that all quantities are considered at t=0)
Ttt=(∂tφ)2−12gtt(∂λφ∂λφ)=(Πα)2−12(−1α2)(−Π2+γij∂iφ∂jφ)=γij∂iφ∂jφ2α2=e−4ϕˉγij∂iφ∂jφ2α2=e−4ϕˆγij∂iφ∂jφ2α2=e−4ϕˆγrr∂rφ∂rφ2α2=e−4ϕΦ22α2By remembering the definition of the normal vector nμ=(−α,0,0,0) (eq. (2.117) of B&S), we can then evaluate the energy density ρ given by eq. (24) of A&C
˜ρ=˜nμ˜nν˜Tμν=e−4˜ϕ2˜Φ2 .Plugging this result in the Hamiltonian constraint, remembering that ψ≡e˜ϕ, we have
∂2rψ+2r∂rψ+πψΦ2=0 .This is a linear elliptic equation which will solve using the procedure described in detail in section 6.2.2 of B&S.
We will discretize the Hamiltonian constraint using second-order accurate finite differences. We get
ψi+1−2ψi+ψi−1Δr2+2ri(ψi+1−ψi−12Δr)+πψiΦ2i=0 ,or, by multiplying the entire equation by Δr2 and then grouping the coefficients of each ψj:
(1−Δrri)ψi−1+(πΔr2Φ2i−2)ψi+(1+Δrri)ψi+1=0 .We choose to set up a grid that is cell-centered, with:
ri=(i−12)Δr ,so that r0=−Δr2. This is a two-point boundary value problem, which we solve using the same strategy as A&C, described in eqs. (48)-(50):
∂rψ|r=0=0 ,limr→∞ψ=1 .In terms of our grid structure, the first boundary condition (regularity at the origin) is written to second-order in Δr as:
∂rψ|r=0=ψ1−ψ0Δr=0⇒ψ0=ψ1 .The second boundary condition (asymptotic flatness) can be interpreted as
ψN=1+CrN (rN≫1) ,which then implies
∂rψN=−Cr2N=−1rN(CrN)=−1rN(ψN−1)=1−ψNrN ,which can then be written as
ψN+1−ψN−12Δr=1−ψNrN⇒ψN+1=ψN−1−2ΔrrNψN+2ΔrrN .Substituting the boundary conditions at the boxed equations above, we end up with
(πΔr2Φ21−1−Δrr1)ψ1+(1+Δrr1)ψ2=0(i=1) ,(1−Δrri)ψi−1+(πΔr2Φ2i−2)ψi+(1+Δrri)ψi+1=0(1<i<N) ,2ψN−1+[πΔr2Φ2N−2−2ΔrrN(1+ΔrrN)]ψN=−2ΔrrN(1+ΔrrN)(i=N) .This results in the following tridiagonal system of linear equations
A⋅→ψ=→s⇒→ψ=A−1⋅→s ,where
A=((πΔr2Φ21−1−Δrr1)(1+Δrr1)00000(1−Δrr2)(πΔr2Φ22−2)(1+Δrr2)00000⋱⋱⋱00000(1−Δrri)(πΔr2Φ2i−2)(1+Δrri)00000⋱⋱⋱00000(1−ΔrrN−1)(πΔr2Φ2N−1−2)(1+ΔrrN−1)000002[πΔr2Φ2N−2−2ΔrrN(1+ΔrrN)]) ,and
→s=(00⋮0⋮0−2ΔrrN(1+ΔrrN))We now start solving the tridiagonal linear system. We start by implementing the tridiagonal matrix A defined above. We break down it down by implementing each diagonal into an array. We start by looking at the main diagonal:
diagmain=((πΔr2Φ21−1−Δrr1)(πΔr2Φ22−2)⋮(πΔr2Φ2i−2)⋮(πΔr2Φ2N−1−2)[πΔr2Φ2N−2−2ΔrrN(1+ΔrrN)])=((πΔr2Φ21−2)(πΔr2Φ22−2)⋮(πΔr2Φ2i−2)⋮(πΔr2Φ2N−1−2)(πΔr2Φ2N−2))+(1−Δrr10⋮0⋮0−2ΔrrN(1+ΔrrN))}N elements# Set the main diagonal
main_diag = np.pi * dr2 * Phi(phi0,r,r0,sigma)**2 - 2
# Update the first element of the main diagonal
main_diag[0] += 1 - dr[0]/r[0]
# Update the last element of the main diagonal
main_diag[NR-1] += - (2 * dr[NR-1] / r[NR-1])*(1 + dr[NR-1] / r[NR-1])
Then we look at the upper diagonal of the A matrix:
diagupper=(1+Δrr11+Δrr2⋮1+Δrri⋮1+ΔrrN−21+ΔrrN−1)}N-1 elements# Set the upper diagonal, ignoring the last point in the r array
upper_diag = np.zeros(NR)
upper_diag[1:] = 1 + dr[:-1]/r[:-1]
Finally, we look at the lower diagonal of the A matrix:
diaglower=(1−Δrr21−Δrr3⋮1−Δrri+1⋮1−ΔrrN−12)}N-1 elements# Set the lower diagonal, start counting the r array at the second element
lower_diag = np.zeros(NR)
lower_diag[:-1] = 1 - dr[1:]/r[1:]
# Change the last term in the lower diagonal to its correct value
lower_diag[NR-2] = 2
Finally, we construct the tridiagonal matrix by adding the three diagonals, while shifting the upper and lower diagonals to the right and left, respectively. Because A is a sparse matrix, we will also use scipy to solve the linear system faster.
!pip install scipy >/dev/null
# Set the sparse matrix A by adding up the three diagonals
A = spdiags([main_diag,upper_diag,lower_diag],[0,1,-1],NR,NR)
# Then compress the sparse matrix A column wise, so that SciPy can invert it later
A = csc_matrix(A)
# Set up the right-hand side of the linear system: s
s = np.zeros(NR)
# Update the last entry of the vector s
s[NR-1] = - (2 * dr[NR-1] / r[NR-1])*(1 + dr[NR-1] / r[NR-1])
# Compress the vector s column-wise
s = csc_matrix(s)
# Solve the sparse linear system using scipy
# https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.spsolve.html
psi = spsolve(A, s.T)
We then show useful plots of the conformal factor ψ and of the evolved conformal factors
ϕ=logψ ,W=ψ−2 ,χ=ψ−4 .import matplotlib.pyplot as plt
# Compute phi
phi = np.log(psi)
# Compute W
W = psi**(-2)
# Compute chi
chi = psi**(-4)
f = plt.figure(figsize=(12,8),dpi=100)
ax = f.add_subplot(221)
ax.set_title(r"Conformal factor $\psi(0,r)$")
ax.set_ylabel(r"$\psi(0,r)$")
ax.plot(r,psi,'k-')
ax.grid()
ax2 = f.add_subplot(222)
ax2.set_title(r"Evolved conformal factor $\phi(0,r)$")
ax2.set_ylabel(r"$\phi(0,r)$")
ax2.plot(r,phi,'r-')
ax2.grid()
ax3 = f.add_subplot(223)
ax3.set_title(r"Evolved conformal factor $W(0,r)$")
ax3.set_xlabel(r"$r$")
ax3.set_ylabel(r"$W(0,r)$")
ax3.plot(r,W,'b-')
ax3.grid()
ax4 = f.add_subplot(224)
ax4.set_title(r"Evolved conformal factor $\chi(0,r)$")
ax4.set_xlabel(r"$r$")
ax4.set_ylabel(r"$\chi(0,r)$")
ax4.plot(r,chi,'c-')
ax4.grid()
outfile = os.path.join(Ccodesdir,"cfs_scalarfield_id.png")
plt.savefig(outfile)
plt.close(f)
# Display the figure
from IPython.display import Image
Image(outfile)
# Set the unity lapse initial condition
alpha_unity = np.ones(NR)
The second one is discussed in the last paragraph of section II.B in Baumgarte, which is to set the "pre-collapsed lapse"
˜α=ψ−2 .# Set the "pre-collapsed lapse" initial condition
alpha_precollapsed = psi**(-2)
# Check to see which version of Python is being used
# For a machine running the final release of Python 3.7.1,
# sys.version_info should return the tuple [3,7,1,'final',0]
if sys.version_info[0] == 3:
np.savetxt(os.path.join(Ccodesdir,"outputSFID_unity_lapse.txt"), list(zip( r, ID_sf, psi**4, alpha_unity )),
fmt="%.15e")
np.savetxt(os.path.join(Ccodesdir,"outputSFID_precollapsed_lapse.txt"), list(zip( r, ID_sf, psi**4, alpha_precollapsed )),
fmt="%.15e")
elif sys.version_info[0] == 2:
np.savetxt(os.path.join(Ccodesdir,"outputSFID_unity_lapse.txt"), zip( r, ID_sf, psi**4, alpha_unity ),
fmt="%.15e")
np.savetxt(os.path.join(Ccodesdir,"outputSFID_precollapsed_lapse.txt"), zip( r, ID_sf, psi**4, alpha_precollapsed ),
fmt="%.15e")
In order to use the initial data file properly, we must tell the program how to interpolate the values we just computed to the values of r in our numerical grid. We do this by creating two C functions: one that interpolates the ADM quantities, {γij,Kij,α,βi,Bi}, and one that interpolates the scalar field quantities, {φ,Π}. The two files written below use the scalarfield_interpolate_1D( ) function, which is defined in the ScalarField/ScalarField_interp.h file. This function performs a Lagrange polynomial interpolation between the initial data file and the numerical grid used during the simulation.
def ID_scalarfield_ADM_quantities(Ccodesdir=".",new_way=False):
includes = ["NRPy_basic_defines.h", "NRPy_function_prototypes.h"]
desc = """(c) 2021 Leo Werneck
This function takes as input either (x,y,z) or (r,th,ph) and outputs
all ADM quantities in the Cartesian or Spherical basis, respectively.
"""
c_type = "void"
name = "ID_scalarfield_ADM_quantities"
params = """const REAL xyz_or_rthph[3],const ID_inputs other_inputs,
REAL *restrict gammaDD00,REAL *restrict gammaDD01,REAL *restrict gammaDD02,
REAL *restrict gammaDD11,REAL *restrict gammaDD12,REAL *restrict gammaDD22,
REAL *restrict KDD00,REAL *restrict KDD01,REAL *restrict KDD02,
REAL *restrict KDD11,REAL *restrict KDD12,REAL *restrict KDD22,
REAL *restrict alpha,
REAL *restrict betaU0,REAL *restrict betaU1,REAL *restrict betaU2,
REAL *restrict BU0,REAL *restrict BU1,REAL *restrict BU2"""
body = """
const REAL r = xyz_or_rthph[0];
const REAL th = xyz_or_rthph[1];
const REAL ph = xyz_or_rthph[2];
REAL sf_star,psi4_star,alpha_star;
scalarfield_interpolate_1D(r,
other_inputs.interp_stencil_size,
other_inputs.numlines_in_file,
other_inputs.r_arr,
other_inputs.sf_arr,
other_inputs.psi4_arr,
other_inputs.alpha_arr,
&sf_star,&psi4_star,&alpha_star);
// Update alpha
*alpha = alpha_star;
// gamma_{rr} = psi^4
*gammaDD00 = psi4_star;
// gamma_{thth} = psi^4 r^2
*gammaDD11 = psi4_star*r*r;
// gamma_{phph} = psi^4 r^2 sin^2(th)
*gammaDD22 = psi4_star*r*r*sin(th)*sin(th);
// All other quantities ARE ZERO:
*gammaDD01 = 0.0; *gammaDD02 = 0.0;
/**/ *gammaDD12 = 0.0;
*KDD00 = 0.0; *KDD01 = 0.0; *KDD02 = 0.0;
/**/ *KDD11 = 0.0; *KDD12 = 0.0;
/**/ *KDD22 = 0.0;
*betaU0 = 0.0; *betaU1 = 0.0; *betaU2 = 0.0;
*BU0 = 0.0; *BU1 = 0.0; *BU2 = 0.0;
"""
if new_way == True:
outC.add_to_Cfunction_dict(includes=includes,desc=desc,c_type=c_type,name=name,
params=params,body=body,enableCparameters=False)
else:
outfile = os.path.join(Ccodesdir,"ID_scalarfield_ADM_quantities-validation.h")
outC.outCfunction(outfile=outfile,
includes=None,desc=desc,c_type=c_type,name=name,
params=params,body=body,enableCparameters=False)
def ID_scalarfield_spherical(Ccodesdir=".",new_way=False):
includes = ["NRPy_basic_defines.h", "NRPy_function_prototypes.h"]
desc = """(c) 2021 Leo Werneck
This function takes as input either (x,y,z) or (r,th,ph) and outputs all
scalar field quantities in the Cartesian or Spherical basis, respectively.
"""
c_type = "void"
name = "ID_scalarfield_spherical"
params = "const REAL xyz_or_rthph[3],const ID_inputs other_inputs,REAL *restrict sf,REAL *restrict sfM"
body = """
const REAL r = xyz_or_rthph[0];
const REAL th = xyz_or_rthph[1];
const REAL ph = xyz_or_rthph[2];
REAL sf_star,psi4_star,alpha_star;
scalarfield_interpolate_1D(r,
other_inputs.interp_stencil_size,
other_inputs.numlines_in_file,
other_inputs.r_arr,
other_inputs.sf_arr,
other_inputs.psi4_arr,
other_inputs.alpha_arr,
&sf_star,&psi4_star,&alpha_star);
// Update varphi
*sf = sf_star;
// Update Pi
*sfM = 0;
"""
if new_way == True:
outC.add_to_Cfunction_dict(includes=includes,desc=desc,c_type=c_type,name=name,
params=params,body=body,enableCparameters=False)
else:
outfile = os.path.join(Ccodesdir,"ID_scalarfield_spherical-validation.h")
outC.outCfunction(outfile=outfile,
includes=None,desc=desc,c_type=c_type,name=name,
params=params,body=body,enableCparameters=False)
In this tutorial module we have explained how to obtain spherically symmetric, time-symmetric initial data for the collapse of a massless scalar field in Spherical coordinates (see Step 1). We have also explained how to interpolate the initial data file to the numerical grid we will use during the simulation (see Step 2).
NRPy+ is capable of generating the BSSN evolution equations in many different Curvilinear coordinates (for example SinhSpherical coordinates, which are of particular interest for this problem). Therefore, it is essential that we convert the Spherical initial data generated here to any Curvilinear system supported by NRPy+.
We start by calling the reference_metric() function within the reference_metric.py NRPy+ module. This will set up a variety of useful quantities for us.
Then the code below interpolate the values of the Spherical grid {r,θ,ϕ} to the Curvilinear grid {xx0,xx1,xx2}.
def ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2(Ccodesdir=".",pointer_to_ID_inputs=False,new_way=False):
rfm.reference_metric()
rthph = outC.outputC(rfm.xxSph[0:3],["rthph[0]", "rthph[1]", "rthph[2]"],
"returnstring", "includebraces=False,outCverbose=False,preindent=1")
includes = ["NRPy_basic_defines.h", "NRPy_function_prototypes.h"]
desc = """(c) 2021 Leo Werneck
This function takes as input either (x,y,z) or (r,th,ph) and outputs all
scalar field quantities in the Cartesian or Spherical basis, respectively.
"""
c_type = "void"
name = "ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2"
params = "const paramstruct *restrict params,const REAL xx0xx1xx2[3],\n"
if pointer_to_ID_inputs == True:
params += "ID_inputs *other_inputs,\n"
else:
params += "ID_inputs other_inputs,\n"
params += "REAL *restrict sf, REAL *restrict sfM"
body = """
const REAL xx0 = xx0xx1xx2[0];
const REAL xx1 = xx0xx1xx2[1];
const REAL xx2 = xx0xx1xx2[2];
REAL rthph[3];
"""+rthph+"""
ID_scalarfield_spherical(rthph,other_inputs,sf,sfM);
"""
if new_way == True:
outC.add_to_Cfunction_dict(includes=includes,desc=desc,c_type=c_type,name=name,
params=params,body=body)
else:
outfile = os.path.join(Ccodesdir,"ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2-validation.h")
outC.outCfunction(outfile=outfile,
includes=None,desc=desc,c_type=c_type,name=name,
params=params,body=body)
Finally, we create the driver function which puts everything together using OpenMP.
def ID_scalarfield(Ccodesdir=".",new_way=False):
includes = ["NRPy_basic_defines.h", "NRPy_function_prototypes.h"]
desc = """(c) 2021 Leo Werneck
This is the scalar field initial data driver functiono.
"""
c_type = "void"
name = "ID_scalarfield"
params = """const paramstruct *restrict params,REAL *restrict xx[3],
ID_inputs other_inputs,REAL *restrict in_gfs"""
body = """
const int idx = IDX3S(i0,i1,i2);
const REAL xx0xx1xx2[3] = {xx0,xx1,xx2};
ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2(params,xx0xx1xx2,other_inputs,
&in_gfs[IDX4ptS(SFGF,idx)],
&in_gfs[IDX4ptS(SFMGF,idx)]);
"""
loopopts = "AllPoints,Read_xxs"
if new_way == True:
outC.add_to_Cfunction_dict(includes=includes,desc=desc,c_type=c_type,name=name,
params=params,body=body,loopopts=loopopts)
else:
outfile = os.path.join(Ccodesdir,"ID_scalarfield-validation.h")
outC.outCfunction(outfile=outfile,
includes=None,desc=desc,c_type=c_type,name=name,
params=params,body=body,loopopts=loopopts)
def NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(Ccodesdir=".",pointer_to_ID_inputs=False,new_way=False):
ID_scalarfield_ADM_quantities(Ccodesdir=Ccodesdir,new_way=new_way)
ID_scalarfield_spherical(Ccodesdir=Ccodesdir,new_way=new_way)
ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2(Ccodesdir=Ccodesdir,pointer_to_ID_inputs=pointer_to_ID_inputs,new_way=new_way)
ID_scalarfield(Ccodesdir=Ccodesdir,new_way=new_way)
NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(Ccodesdir=Ccodesdir)
Output C function ID_scalarfield_ADM_quantities() to file ScalarFieldID_validation/ID_scalarfield_ADM_quantities-validation.h Output C function ID_scalarfield_spherical() to file ScalarFieldID_validation/ID_scalarfield_spherical-validation.h Output C function ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2() to file ScalarFieldID_validation/ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2-validation.h Output C function ID_scalarfield() to file ScalarFieldID_validation/ID_scalarfield-validation.h
First we load the ScalarField/ScalarField_InitialData.py module and compute everything by using the scalarfield_initial_data( ) function, which should do exactly the same as we have done in this tutorial.
# Import the ScalarField.ScalarField_InitialData NRPy module
import ScalarField.ScalarField_InitialData as sfid
# Output the unity lapse initial data file
outputname = os.path.join(Ccodesdir,"outputSFID_unity_lapse-validation.txt")
sfid.ScalarField_InitialData(outputname,ID_Family,
phi0,r0,sigma,NR,RMAX,CoordSystem=CoordSystem,
sinhA=sinhA,sinhW=sinhW,lapse_condition="Unity")
# Output the "pre-collapsed" lapse initial data file
outputname = os.path.join(Ccodesdir,"outputSFID_precollapsed_lapse-validation.txt")
sfid.ScalarField_InitialData(outputname,ID_Family,
phi0,r0,sigma,NR,RMAX,CoordSystem=CoordSystem,
sinhA=sinhA,sinhW=sinhW,lapse_condition="Pre-collapsed")
# Output C codes
sfid.NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(Ccodesdir=Ccodesdir)
Generated the ADM initial data for the gravitational collapse of a massless scalar field in Spherical coordinates. Type of initial condition: Scalar field: "Gaussian" Shell ADM quantities: Time-symmetric Lapse condition: Unity Parameters: amplitude = 0.1, center = 0, width = 1, domain size = 50, number of points = 30000, Initial data file = ScalarFieldID_validation/outputSFID_unity_lapse-validation.txt. Generated the ADM initial data for the gravitational collapse of a massless scalar field in Spherical coordinates. Type of initial condition: Scalar field: "Gaussian" Shell ADM quantities: Time-symmetric Lapse condition: Pre-collapsed Parameters: amplitude = 0.1, center = 0, width = 1, domain size = 50, number of points = 30000, Initial data file = ScalarFieldID_validation/outputSFID_precollapsed_lapse-validation.txt. Output C function ID_scalarfield_ADM_quantities() to file ScalarFieldID_validation/ID_scalarfield_ADM_quantities.h Output C function ID_scalarfield_spherical() to file ScalarFieldID_validation/ID_scalarfield_spherical.h Output C function ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2() to file ScalarFieldID_validation/ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2.h Output C function ID_scalarfield() to file ScalarFieldID_validation/ID_scalarfield.h
import filecmp
if filecmp.cmp(os.path.join(Ccodesdir,'outputSFID_unity_lapse.txt'),
os.path.join(Ccodesdir,'outputSFID_unity_lapse-validation.txt')) == False:
print("ERROR: Unity lapse initial data test FAILED!")
sys.exit(1)
else:
print(" Unity lapse initial data test: PASSED!")
if filecmp.cmp(os.path.join(Ccodesdir,'outputSFID_precollapsed_lapse.txt'),
os.path.join(Ccodesdir,'outputSFID_precollapsed_lapse-validation.txt')) == False:
print("ERROR: \"Pre-collapsed\" lapse initial data test FAILED!")
sys.exit(1)
else:
print(" \"Pre-collapsed\" lapse initial data test: PASSED!")
if filecmp.cmp(os.path.join(Ccodesdir,'ID_scalarfield_ADM_quantities.h'),
os.path.join(Ccodesdir,'ID_scalarfield_ADM_quantities-validation.h')) == False:
print("ERROR: ADM quantities interpolation file test FAILED!")
sys.exit(1)
else:
print(" ADM quantities interpolation file test: PASSED!")
if filecmp.cmp(os.path.join(Ccodesdir,'ID_scalarfield_spherical.h'),
os.path.join(Ccodesdir,'ID_scalarfield_spherical-validation.h')) == False:
print("ERROR: Scalar field interpolation file test FAILED!")
sys.exit(1)
else:
print(" Scalar field interpolation file test: PASSED!")
if filecmp.cmp(os.path.join(Ccodesdir,'ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2.h'),
os.path.join(Ccodesdir,'ID_scalarfield_xx0xx1xx2_to_BSSN_xx0xx1xx2-validation.h')) == False:
print("ERROR: Scalar field Spherical to Curvilinear test FAILED!")
sys.exit(1)
else:
print("Scalar field Spherical to Curvilinear test: PASSED!")
if filecmp.cmp(os.path.join(Ccodesdir,'ID_scalarfield.h'),
os.path.join(Ccodesdir,'ID_scalarfield-validation.h')) == False:
print("ERROR: Scalar field driver test: FAILED!")
sys.exit(1)
else:
print(" Scalar field driver test: PASSED!")
Unity lapse initial data test: PASSED! "Pre-collapsed" lapse initial data test: PASSED! ADM quantities interpolation file test: PASSED! Scalar field interpolation file test: PASSED! Scalar field Spherical to Curvilinear test: PASSED! Scalar field driver test: PASSED!
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-ADM_Initial_Data-ScalarField.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-ADM_Initial_Data-ScalarField")
Created Tutorial-ADM_Initial_Data-ScalarField.tex, and compiled LaTeX file to PDF file Tutorial-ADM_Initial_Data-ScalarField.pdf