diff --git a/scri/asymptotic_bondi_data/constraints.py b/scri/asymptotic_bondi_data/constraints.py index 03777c0d..cbbc5ecc 100644 --- a/scri/asymptotic_bondi_data/constraints.py +++ b/scri/asymptotic_bondi_data/constraints.py @@ -11,10 +11,10 @@ def bondi_constraints(self, lhs=True, rhs=True): Bondi gauge establishes some relations that the data must satisfy: - ψ̇₀ = -ðψ₁ + 3 σ ψ₂ - ψ̇₁ = -ðψ₂ + 2 σ ψ₃ - ψ̇₂ = -ðψ₃ + 1 σ ψ₄ - ψ₃ = ∂ðσ̄/∂u + ψ̇₀ = ðψ₁ + 3 σ ψ₂ + ψ̇₁ = ðψ₂ + 2 σ ψ₃ + ψ̇₂ = ðψ₃ + 1 σ ψ₄ + ψ₃ = -∂ðσ̄/∂u ψ₄ = -∂²σ̄/∂u² Im[ψ₂] = -Im[ð²σ̄ + σ ∂σ̄/∂u] @@ -38,10 +38,10 @@ def bondi_violations(self): Bondi gauge establishes some relations that the data must satisfy: - ψ̇₀ = -ðψ₁ + 3 σ ψ₂ - ψ̇₁ = -ðψ₂ + 2 σ ψ₃ - ψ̇₂ = -ðψ₃ + 1 σ ψ₄ - ψ₃ = ∂ðσ̄/∂u + ψ̇₀ = ðψ₁ + 3 σ ψ₂ + ψ̇₁ = ðψ₂ + 2 σ ψ₃ + ψ̇₂ = ðψ₃ + 1 σ ψ₄ + ψ₃ = -∂ðσ̄/∂u ψ₄ = -∂²σ̄/∂u² Im[ψ₂] = -Im[ð²σ̄ + σ ∂σ̄/∂u] @@ -60,10 +60,10 @@ def bondi_violation_norms(self): Bondi gauge establishes some relations that the data must satisfy: - ψ̇₀ = -ðψ₁ + 3 σ ψ₂ - ψ̇₁ = -ðψ₂ + 2 σ ψ₃ - ψ̇₂ = -ðψ₃ + 1 σ ψ₄ - ψ₃ = ∂ðσ̄/∂u + ψ̇₀ = ðψ₁ + 3 σ ψ₂ + ψ̇₁ = ðψ₂ + 2 σ ψ₃ + ψ̇₂ = ðψ₃ + 1 σ ψ₄ + ψ₃ = -∂ðσ̄/∂u ψ₄ = -∂²σ̄/∂u² Im[ψ₂] = -Im[ð²σ̄ + σ ∂σ̄/∂u] @@ -83,7 +83,7 @@ def bianchi_0(self, lhs=True, rhs=True): In Bondi coordinates, the Bianchi identities simplify, resulting in this expression (among others) for the time derivative of ψ₀: - ψ̇₀ = -ðψ₁ + 3 σ ψ₂ + ψ̇₀ = ðψ₁ + 3 σ ψ₂ Parameters ========== @@ -99,7 +99,7 @@ def bianchi_0(self, lhs=True, rhs=True): if lhs: lhs_value = self.psi0.dot if rhs: - rhs_value = -self.psi1.eth + 3 * self.sigma * self.psi2 + rhs_value = self.psi1.eth_GHP + 3 * self.sigma * self.psi2 if lhs and rhs: return (lhs_value, rhs_value) elif lhs: @@ -114,7 +114,7 @@ def bianchi_1(self, lhs=True, rhs=True): In Bondi coordinates, the Bianchi identities simplify, resulting in this expression (among others) for the time derivative of ψ₁: - ψ̇₁ = -ðψ₂ + 2 σ ψ₃ + ψ̇₁ = ðψ₂ + 2 σ ψ₃ Parameters ========== @@ -130,7 +130,7 @@ def bianchi_1(self, lhs=True, rhs=True): if lhs: lhs_value = self.psi1.dot if rhs: - rhs_value = -self.psi2.eth + 2 * self.sigma * self.psi3 + rhs_value = self.psi2.eth_GHP + 2 * self.sigma * self.psi3 if lhs and rhs: return (lhs_value, rhs_value) elif lhs: @@ -145,7 +145,7 @@ def bianchi_2(self, lhs=True, rhs=True): In Bondi coordinates, the Bianchi identities simplify, resulting in this expression (among others) for the time derivative of ψ₂: - ψ̇₂ = -ðψ₃ + σ ψ₄ + ψ̇₂ = ðψ₃ + σ ψ₄ Parameters ========== @@ -161,7 +161,7 @@ def bianchi_2(self, lhs=True, rhs=True): if lhs: lhs_value = self.psi2.dot if rhs: - rhs_value = -self.psi3.eth + self.sigma * self.psi4 + rhs_value = self.psi3.eth_GHP + self.sigma * self.psi4 if lhs and rhs: return (lhs_value, rhs_value) elif lhs: @@ -176,7 +176,7 @@ def constraint_3(self, lhs=True, rhs=True): In Bondi coordinates, the value of ψ₃ is given by a time derivative and an angular derivative of the (conjugate) shear: - ψ₃ = ∂ðσ̄/∂u + ψ₃ = -∂ðσ̄/∂u Parameters ========== @@ -192,7 +192,7 @@ def constraint_3(self, lhs=True, rhs=True): if lhs: lhs_value = self.psi3 if rhs: - rhs_value = self.sigma.bar.dot.eth + rhs_value = -self.sigma.bar.dot.eth_GHP if lhs and rhs: return (lhs_value, rhs_value) elif lhs: @@ -256,7 +256,7 @@ def constraint_mass_aspect(self, lhs=True, rhs=True): if lhs: lhs_value = np.imag(self.psi2) if rhs: - rhs_value = -np.imag(self.sigma.bar.eth.eth + self.sigma * self.sigma.bar.dot) + rhs_value = -np.imag(self.sigma.bar.eth_GHP.eth_GHP + self.sigma * self.sigma.bar.dot) if lhs and rhs: return (lhs_value, rhs_value) elif lhs: diff --git a/scri/asymptotic_bondi_data/from_initial_values.py b/scri/asymptotic_bondi_data/from_initial_values.py index 3189a11c..2188e37c 100644 --- a/scri/asymptotic_bondi_data/from_initial_values.py +++ b/scri/asymptotic_bondi_data/from_initial_values.py @@ -4,10 +4,10 @@ def from_initial_values(cls, time, ell_max=8, sigma0=0.0, sigmadot0=0.0, sigmadd The initial-value formulation for Bondi gauge is determined by these relations: - ψ̇₀ = -ðψ₁ + 3 σ ψ₂ - ψ̇₁ = -ðψ₂ + 2 σ ψ₃ - ψ̇₂ = -ðψ₃ + 1 σ ψ₄ - ψ₃ = ∂ðσ̄/∂u + ψ̇₀ = ðψ₁ + 3 σ ψ₂ + ψ̇₁ = ðψ₂ + 2 σ ψ₃ + ψ̇₂ = ðψ₃ + 1 σ ψ₄ + ψ₃ = -∂ðσ̄/∂u ψ₄ = -∂²σ̄/∂u² We also have a constraint on the initial value of ψ₂: @@ -87,7 +87,7 @@ def asany_atleast2d_complex(a): # constant in time. If this is true, assumes sigmadot0 and sigmaddot0 # are constants in time, and just integrates them. u = abd.time - ð = lambda x: x.eth + ð = lambda x: x.eth_GHP conjugate = lambda x: x.bar σ_0 = ModesTimeSeries(asany_atleast2d_complex(sigma0), u, spin_weight=2, multiplication_truncator=max) σ_1 = ModesTimeSeries(asany_atleast2d_complex(sigmadot0), u, spin_weight=2, multiplication_truncator=max) @@ -96,33 +96,33 @@ def asany_atleast2d_complex(a): # ψ₄ = -∂²σ̄/∂u² ψ4_0 = -2 * conjugate(σ_2) abd.psi4 = ψ4_0 - # ψ₃ = ð ∂σ̄/∂u - ψ3_0 = ð(conjugate(σ_1)) - ψ3_1 = 2 * ð(conjugate(σ_2)) + # ψ₃ = -ð ∂σ̄/∂u + ψ3_0 = -ð(conjugate(σ_1)) + ψ3_1 = -2 * ð(conjugate(σ_2)) abd.psi3 = u * ψ3_1 + ψ3_0 # ψ₂ = ∫ (-ðψ₃ + σψ₄) du ψ2_0 = ( ModesTimeSeries(psi2, u, spin_weight=0, multiplication_truncator=max).real - 1j * (σ_0.bar.eth.eth + σ_0 * σ_1.bar).imag ) - ψ2_1 = σ_0 * ψ4_0 - ð(ψ3_0) - ψ2_2 = (σ_1 * ψ4_0 - ð(ψ3_1)) / 2 + ψ2_1 = σ_0 * ψ4_0 + ð(ψ3_0) + ψ2_2 = (σ_1 * ψ4_0 + ð(ψ3_1)) / 2 ψ2_3 = (1 / 3) * σ_2 * ψ4_0 abd.psi2 = u * (u * (u * ψ2_3 + ψ2_2) + ψ2_1) + ψ2_0 - # ψ₁ = ∫ (-ðψ₂ + 2σψ₃) du + # ψ₁ = ∫ (ðψ₂ + 2σψ₃) du ψ1_0 = ModesTimeSeries(psi1, u, spin_weight=1, multiplication_truncator=max) - ψ1_1 = 2 * σ_0 * ψ3_0 - ð(ψ2_0) - ψ1_2 = σ_0 * ψ3_1 + σ_1 * ψ3_0 - ð(ψ2_1) / 2 - ψ1_3 = (2 * σ_1 * ψ3_1 + 2 * σ_2 * ψ3_0 - ð(ψ2_2)) / 3 - ψ1_4 = (2 * σ_2 * ψ3_1 - ð(ψ2_3)) / 4 + ψ1_1 = 2 * σ_0 * ψ3_0 + ð(ψ2_0) + ψ1_2 = σ_0 * ψ3_1 + σ_1 * ψ3_0 + ð(ψ2_1) / 2 + ψ1_3 = (2 * σ_1 * ψ3_1 + 2 * σ_2 * ψ3_0 + ð(ψ2_2)) / 3 + ψ1_4 = (2 * σ_2 * ψ3_1 + ð(ψ2_3)) / 4 abd.psi1 = u * (u * (u * (u * ψ1_4 + ψ1_3) + ψ1_2) + ψ1_1) + ψ1_0 - # ψ₀ = ∫ (-ðψ₁ + 3σψ₂) du + # ψ₀ = ∫ (ðψ₁ + 3σψ₂) du ψ0_0 = ModesTimeSeries(psi0, u, spin_weight=2, multiplication_truncator=max) - ψ0_1 = 3 * σ_0 * ψ2_0 - ð(ψ1_0) - ψ0_2 = (3 * σ_0 * ψ2_1 + 3 * σ_1 * ψ2_0 - ð(ψ1_1)) / 2 - ψ0_3 = σ_0 * ψ2_2 + σ_1 * ψ2_1 + σ_2 * ψ2_0 - ð(ψ1_2) / 3 - ψ0_4 = (3 * σ_0 * ψ2_3 + 3 * σ_1 * ψ2_2 + 3 * σ_2 * ψ2_1 - ð(ψ1_3)) / 4 - ψ0_5 = (3 * σ_1 * ψ2_3 + 3 * σ_2 * ψ2_2 - ð(ψ1_4)) / 5 + ψ0_1 = 3 * σ_0 * ψ2_0 + ð(ψ1_0) + ψ0_2 = (3 * σ_0 * ψ2_1 + 3 * σ_1 * ψ2_0 + ð(ψ1_1)) / 2 + ψ0_3 = σ_0 * ψ2_2 + σ_1 * ψ2_1 + σ_2 * ψ2_0 + ð(ψ1_2) / 3 + ψ0_4 = (3 * σ_0 * ψ2_3 + 3 * σ_1 * ψ2_2 + 3 * σ_2 * ψ2_1 + ð(ψ1_3)) / 4 + ψ0_5 = (3 * σ_1 * ψ2_3 + 3 * σ_2 * ψ2_2 + ð(ψ1_4)) / 5 ψ0_6 = σ_2 * ψ2_3 / 2 abd.psi0 = u * (u * (u * (u * (u * (u * ψ0_6 + ψ0_5) + ψ0_4) + ψ0_3) + ψ0_2) + ψ0_1) + ψ0_0 elif np.ndim(sigma0) == 2: @@ -140,14 +140,14 @@ def adjust_psi2_imaginary_part(psi2, abd): abd.sigma = sigma0 # ψ₄ = -∂²σ̄/∂u² abd.psi4 = -abd.sigma.ddot.bar - # ψ₃ = ð ∂σ̄/∂u - abd.psi3 = abd.sigma.dot.bar.eth - # ψ₂ = ∫ (-ðψ₃ + σψ₄) du - abd.psi2 = (-abd.psi3.eth + abd.sigma * abd.psi4).int + adjust_psi2_imaginary_part(psi2, abd) - # ψ₁ = ∫ -ðψ₂ + σψ₃ du - abd.psi1 = (-abd.psi2.eth + 2 * abd.sigma * abd.psi3).int + psi1 - # ψ₀ = ∫ -ðψ₁ + σψ₂ dt - abd.psi0 = (-abd.psi1.eth + 3 * abd.sigma * abd.psi2).int + psi0 + # ψ₃ = -ð ∂σ̄/∂u + abd.psi3 = -abd.sigma.dot.bar.eth_GHP + # ψ₂ = ∫ (ðψ₃ + σψ₄) du + abd.psi2 = (abd.psi3.eth_GHP + abd.sigma * abd.psi4).int + adjust_psi2_imaginary_part(psi2, abd) + # ψ₁ = ∫ ðψ₂ + σψ₃ du + abd.psi1 = (abd.psi2.eth_GHP + 2 * abd.sigma * abd.psi3).int + psi1 + # ψ₀ = ∫ ðψ₁ + σψ₂ dt + abd.psi0 = (abd.psi1.eth_GHP + 3 * abd.sigma * abd.psi2).int + psi0 else: raise ValueError(f"Input `sigma0` must have 1 or 2 dimensions; it has {np.ndim(sigma0)}") diff --git a/scri/asymptotic_bondi_data/transformations.py b/scri/asymptotic_bondi_data/transformations.py index 8e76b045..40d0304b 100644 --- a/scri/asymptotic_bondi_data/transformations.py +++ b/scri/asymptotic_bondi_data/transformations.py @@ -310,8 +310,9 @@ def transform(self, **kwargs): # variation in u, the second axis variation in θ', and the third axis variation in ϕ'. u = self.u α = sf.Grid(supertranslation.evaluate(distorted_grid_rotors), spin_weight=0).real[np.newaxis, :, :] - ðα = sf.Grid(supertranslation.eth.evaluate(distorted_grid_rotors), spin_weight=α.s + 1)[np.newaxis, :, :] - ððα = sf.Grid(supertranslation.eth.eth.evaluate(distorted_grid_rotors), spin_weight=α.s + 2)[np.newaxis, :, :] + # The factors of 1/sqrt(2) and 1/2 come from using the GHP eth instead of the NP eth. + ðα = sf.Grid(supertranslation.eth.evaluate(distorted_grid_rotors)/np.sqrt(2), spin_weight=α.s + 1)[np.newaxis, :, :] + ððα = sf.Grid(0.5*supertranslation.eth.eth.evaluate(distorted_grid_rotors), spin_weight=α.s + 2)[np.newaxis, :, :] k, ðk_over_k, one_over_k, one_over_k_cubed = conformal_factors(boost_velocity, distorted_grid_rotors) # ðu'(u, θ', ϕ') exp(iλ) / k(θ', ϕ') diff --git a/scri/modes_time_series.py b/scri/modes_time_series.py index c93a0481..7f6677f8 100644 --- a/scri/modes_time_series.py +++ b/scri/modes_time_series.py @@ -119,3 +119,13 @@ def int(self): def iint(self): """Integrate modes twice with respect to time""" return self.antiderivative(2) + + @property + def eth_GHP(self): + """Raise spin-weight with GHP convention""" + return self.eth / np.sqrt(2) + + @property + def ethbar_GHP(self): + """Lower spin-weight with GHP convention""" + return self.ethbar / np.sqrt(2)