evaluate_bondi_norms.py

#%%
"""
This script shows how to load CCE data into a `scri.AsymptoticBondiData` object
and evaluate a set of constraints that should be satisfied involving the strain
and the Newman-Penrose quantities (assuming the waveform is in Bondi gauge).
These constraints are listed in the plot legend (see below) as equations to be
satisfied at each instant.  The violation is just one side minus the other, and
the norm is that difference squared integrated over the sphere (and
square-rooted).
"""

import time
import numpy as np
import matplotlib.pyplot as plt
import scri


# The data output by SpEC make certain choices for tetrad
# normalizations, which I assume are shared by SpECTRE
convention = "SpEC"


# This is for the current SpECTRE CCE output format
file_format = "spectrecce1"


# Absolute or relative path to the file containing the CCE data
file_name = "/home/fs01/nd357/Runs/GhDgFdPaper/8Points/CceOutputR200.h5"


# Now we load up the files.  This may take several minutes (probably due to weird chunking)
t1 = time.perf_counter()
abd1 = scri.SpEC.file_io.create_abd_from_h5(
    file_format,
    convention=convention,
    file_name=file_name
)
t2 = time.perf_counter()
print(f"Time to load file: {t2-t1} s")


# The next step will take forever if we don't downsample the data
abd2 = abd1.interpolate(abd1.time[::50])


# Now, just compute the norms of the Bondi constraints
constraint_norms = abd2.bondi_violation_norms

#%% Plot the results
plt.semilogy(abd2.t, np.array(constraint_norms).T)
plt.xlabel("Time (code units)")
plt.ylabel(r"Bondi-constraint violation $L^2$ norms")
plt.legend([
    "ψ̇₀ = ðψ₁ + 3 σ ψ₂",
    "ψ̇₁ = ðψ₂ + 2 σ ψ₃",
    "ψ̇₂ = ðψ₃ + 1 σ ψ₄",
    "ψ₃ = -∂ðσ̄/∂u",
    "ψ₄ = -∂²σ̄/∂u²",
    "Im[ψ₂] = -Im[ð²σ̄ + σ ∂σ̄/∂u]"
])
# plt.title("")
plt.savefig("BondiViolationNorms.png");
plt.show()