Merge pull request #37 from 10220/set_convention

Set ABD internal conventions to Moreschi-Boyle

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:

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)}")
 

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(θ', ϕ')

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)