WIP: remodeling PNSystem infrastructure

docs/src/internals/pn_systems.md

@@ -2,9 +2,9 @@
 
 ```@docs
 PNSystem
+Quasispherical
 BBH
 BHNS
 NSNS
 FDPNSystem
-fd_pnsystem
 ```

src/PostNewtonian.jl

@@ -1,12 +1,15 @@
 module PostNewtonian
 
-# Always explicitly address functions similar to functions defined in this package,
-# which come from these packages:
+# We must always explicitly qualify functions similar to functions defined in this package
+# by the name of the package.  We will use such functions from these packages:
 using MacroTools: MacroTools
 using FastDifferentiation: FastDifferentiation
-using RuntimeGeneratedFunctions: RuntimeGeneratedFunctions
 
-# Otherwise, we just explicitly import specific functions:
+# Otherwise, we just explicitly import specific functions / types.  Note that the difference
+# between `using` and `import` in the following lines is that `using` will only allow us to
+# call the functions, while `import` would also allow us to specialize them (define new
+# methods of the imported functions).
+using FastDifferentiation: Node as FDNode
 using DataInterpolations: CubicSpline
 using InteractiveUtils: methodswith
 using LinearAlgebra: mul!

src/dynamics/right_hand_sides.jl

@@ -21,6 +21,31 @@ function TaylorT5_v̇(p)
     return inv(truncated_series_ratio(v̇_denominator_coeffs(p), v̇_numerator_coeffs(p)))
 end
 
+"""
+    causes_domain_error!(u̇, p)
+
+Ensure that these parameters correspond to a physically valid set of PN parameters.
+
+If the parameters are not valid, this function should modify `u̇` to indicate that the
+current step is invalid.  This is done by filling `u̇` with `NaN`s, which will be detected
+by the ODE solver and cause it to try a different (smaller) step size.
+
+Currently, the only check that is done is to test that these parameters result in a PN
+parameter v>0.  In the future, this function may be expanded to include other checks, or it
+may be specialized for specific `PNSystem` subtypes.
+"""
+function causes_domain_error!(u̇::ST, p::PNSystem{NT}) where {ST,NT}
+    if !ismutabletype(ST)
+        error("`causes_domain_error!` cannot modify input `u̇` because it is immutable")
+    end
+    if v(p) ≤ 0  # If this is expanded, document the change in the docstring.
+        u̇ .= convert(eltype(NT), NaN)
+        true
+    else
+        false
+    end
+end
+
 @pn_expression function TaylorTn!(pnsystem, u̇, TaylorTn_v̇::V̇) where {V̇}
     # If these parameters result in v≤0, fill u̇ with NaNs so that `solve` will
     # know that this was a bad step and try again.

src/pn_systems/BBH.jl

@@ -0,0 +1,117 @@
+"""
+    BBH{NT, ST, PNOrder} <: PNSystem{NT, ST, PNOrder}
+
+The [`PNSystem`](@ref) subtype describing a binary black hole system.
+
+The `state` vector here holds the fundamental state variables characterizing the masses,
+spins, orientation, velocity, and orbital phase of the system.  The spins unpacked into
+three components each.  The orientation is described by the four components of the `Rotor`
+`R`.  This gives us a total of 14 elements:
+
+    M₁, M₂, χ⃗₁ˣ, χ⃗₁ʸ, χ⃗₁ᶻ, χ⃗₂ˣ, χ⃗₂ʸ, χ⃗₂ᶻ, Rʷ, Rˣ, Rʸ, Rᶻ, v, Φ
+
+The "orbital phase" `Φ` is tracked as the 14th element of the `state` vector.  This is just
+the integral of the (scalar) orbital angular frequency `Ω`, and holds little interest for
+general systems beyond a convenient description of how "far" the system has evolved.  For
+nonprecessing systems, `Φ` would be sufficient to describe the system's position, which is
+more completely described by the `Rotor` `R`.  However, for precessing systems, it is
+difficult to extract this quantity from `R`.
+"""
+struct BBH{NT,ST,PNOrder} <: PNSystem{NT,ST,PNOrder}
+    state::ST
+
+    function BBH{NT,ST,PNOrder}(state) where {NT,ST,PNOrder}
+        if eachindex(state) != Base.OneTo(14)
+            error(
+                "The `state` vector for `BBH` must be indexed from 1 to 14; " *
+                "input is indexed `$(eachindex(state))`.",
+            )
+        end
+        new{NT,ST,PNOrder}(state)
+    end
+    function BBH(; M₁, M₂, χ⃗₁, χ⃗₂, v, R=Rotor(1), Φ=0, PNOrder=typemax(Int), kwargs...)
+        (NT, ST, PNOrder, state) = prepare_system(; M₁, M₂, χ⃗₁, χ⃗₂, R, v, Φ, PNOrder)
+        return new{NT,ST,PNOrder}(state)
+    end
+    function BBH(state; Λ₁=0, Λ₂=0, PNOrder=typemax(Int))
+        if eachindex(state) != Base.OneTo(14)
+            error(
+                "The `state` vector for `BBH` must be indexed from 1 to 14; " *
+                "input is indexed `$(eachindex(state))`.",
+            )
+        end
+        @assert Λ₁ == 0
+        @assert Λ₂ == 0
+        return new{eltype(state),typeof(state),prepare_pn_order(PNOrder)}(state)
+    end
+end
+const BHBH = BBH
+
+# The following are methods of functions defined in `state_variables.jl`, specialized for
+# `BBH` systems.
+state(pnsystem::BBH) = pnsystem.state
+function symbols(::Type{<:BBH})
+    (:M₁, :M₂, :χ⃗₁ˣ, :χ⃗₁ʸ, :χ⃗₁ᶻ, :χ⃗₂ˣ, :χ⃗₂ʸ, :χ⃗₂ᶻ, :Rʷ, :Rˣ, :Rʸ, :Rᶻ, :v, :Φ)
+end
+function ascii_symbols(::Type{<:BBH})
+    (:M1, :M2, :chi1x, :chi1y, :chi1z, :chi2x, :chi2y, :chi2z, :Rw, :Rx, :Ry, :Rz, :v, :Phi)
+end
+for (i, symbol) ∈ enumerate(symbols(BBH))
+    # This will define, e.g., `M₁(pnsystem::BBH) = pnsystem.state[1]`.  We
+    # could do this manually, but this is more concise and less error-prone.
+    @eval begin
+        $(symbol)(pnsystem::BBH) = @inbounds pnsystem.state[$i]
+        function symbol_index(::Type{T}, ::Val{Symbol($symbol)}) where {T<:BBH}
+            $i
+        end
+    end
+end
+
+Λ₁(pnsystem::BBH) = zero(pnsystem)
+Λ₂(pnsystem::BBH) = zero(pnsystem)
+
+@testitem "BBH constructors" begin
+    using Quaternionic
+
+    pnA = BBH(;
+        M₁=1.0f0, M₂=2.0f0, χ⃗₁=Float32[3.0, 4.0, 5.0], χ⃗₂=Float32[6.0, 7.0, 8.0], v=0.23f0
+    )
+    @test pnA.state ==
+        Float32[1.0; 2.0; 3.0; 4.0; 5.0; 6.0; 7.0; 8.0; 1.0; 0.0; 0.0; 0.0; 0.23; 0.0]
+
+    pnB = BBH(;
+        M₁=1.0f0,
+        M₂=2.0f0,
+        χ⃗₁=Float32[3.0, 4.0, 5.0],
+        χ⃗₂=Float32[6.0, 7.0, 8.0],
+        v=0.23f0,
+        Φ=9.0f0,
+    )
+    @test pnB.state ==
+        Float32[1.0; 2.0; 3.0; 4.0; 5.0; 6.0; 7.0; 8.0; 1.0; 0.0; 0.0; 0.0; 0.23; 9.0]
+
+    R = randn(RotorF32)
+    pn1 = BBH(;
+        M₁=1.0f0,
+        M₂=2.0f0,
+        χ⃗₁=Float32[3.0, 4.0, 5.0],
+        χ⃗₂=Float32[6.0, 7.0, 8.0],
+        R=R,
+        v=0.23f0,
+    )
+    @test pn1.state ≈ [1.0; 2.0; 3.0; 4.0; 5.0; 6.0; 7.0; 8.0; components(R)...; 0.23; 0.0]
+
+    pn2 = BBH(;
+        M₁=1.0f0,
+        M₂=2.0f0,
+        χ⃗₁=Float32[3.0, 4.0, 5.0],
+        χ⃗₂=Float32[6.0, 7.0, 8.0],
+        R=R,
+        v=0.23f0,
+        Φ=9.0f0,
+    )
+    @test pn2.state ≈ [1.0; 2.0; 3.0; 4.0; 5.0; 6.0; 7.0; 8.0; components(R)...; 0.23; 9.0]
+
+    pn1.state[end] = 9.0f0
+    @test pn1.state == pn2.state
+end

src/pn_systems/BHNS.jl

@@ -0,0 +1,80 @@
+"""
+    BHNS{NT, ST, PNOrder} <: PNSystem{NT, ST, PNOrder}
+
+The [`PNSystem`](@ref) subtype describing a black-hole—neutron-star binary system.
+
+The `state` vector is the same as for a [`BBH`](@ref), with an additional field `Λ₂` holding
+the (constant) tidal-coupling parameter of the neutron star.
+
+Note that the neutron star is *always* object 2 — meaning that `M₂`, `χ⃗₂`, and `Λ₂` always
+refer to it; `M₁` and `χ⃗₁` always refer to the black hole.  (It's "BHNS", not "NSBH".)  See
+also [`NSNS`](@ref).
+"""
+struct BHNS{NT,ST,PNOrder} <: PNSystem{NT,ST,PNOrder}
+    state::ST
+
+    function BHNS{NT,ST,PNOrder}(state) where {NT,ST,PNOrder}
+        if eachindex(state) != Base.OneTo(15)
+            error(
+                "The `state` vector for `BHNS` must be indexed from 1 to 15; " *
+                "input is indexed `$(eachindex(state))`.",
+            )
+        end
+        new{NT,ST,PNOrder}(state)
+    end
+    function BHNS(;
+        M₁, M₂, χ⃗₁, χ⃗₂, v, R=Rotor(1), Φ=0, Λ₂, PNOrder=typemax(Int), kwargs...
+    )
+        NT, ST, PNOrder, state = prepare_system(; M₁, M₂, χ⃗₁, χ⃗₂, R, v, Φ, Λ₂, PNOrder)
+        return new{NT,ST,PNOrder}(state)
+    end
+    function BHNS(state; PNOrder=typemax(Int))
+        if eachindex(state) != Base.OneTo(15)
+            error(
+                "The `state` vector for `BHNS` must be indexed from 1 to 15; " *
+                "input is indexed `$(eachindex(state))`.",
+            )
+        end
+        NT, ST, PNOrder = eltype(state), typeof(state), prepare_pn_order(PNOrder)
+        return new{NT,ST,PNOrder}(state)
+    end
+end
+
+# The following are methods of functions defined in `state_variables.jl`, specialized for
+# `BHNS` systems.
+state(pnsystem::BHNS) = pnsystem.state
+function symbols(::Type{<:BHNS})
+    (:M₁, :M₂, :χ⃗₁ˣ, :χ⃗₁ʸ, :χ⃗₁ᶻ, :χ⃗₂ˣ, :χ⃗₂ʸ, :χ⃗₂ᶻ, :Rʷ, :Rˣ, :Rʸ, :Rᶻ, :v, :Φ, :Λ₂)
+end
+function ascii_symbols(::Type{<:BHNS})
+    (
+        :M1,
+        :M2,
+        :chi1x,
+        :chi1y,
+        :chi1z,
+        :chi2x,
+        :chi2y,
+        :chi2z,
+        :Rw,
+        :Rx,
+        :Ry,
+        :Rz,
+        :v,
+        :Phi,
+        :Lambda2,
+    )
+end
+for (i, symbol) ∈ enumerate(symbols(BHNS))
+    # This will define, e.g., `M₁(pnsystem::BHNS) = pnsystem.state[1]`.  We
+    # could do this manually, but this is more concise and less error-prone.
+    @eval begin
+        $(symbol)(pnsystem::BHNS) = @inbounds pnsystem.state[$i]
+        function symbol_index(::Type{T}, ::Val{Symbol($symbol)}) where {T<:BHNS}
+            $i
+        end
+    end
+end
+
+Λ₁(pnsystem::BHNS) = zero(pnsystem)
+Λ₂(pnsystem::BHNS) = @inbounds pnsystem.state[15]

src/pn_systems/FDPNsystem.jl

@@ -0,0 +1,34 @@
+"""
+    FDPNSystem{NT, PN, PNOrder} <: PNSystem{FDNode, Vector{FDNode}, PNOrder}
+
+A `PNSystem` that contains information as variables from
+[`FastDifferentiation.jl`](https://docs.juliahub.com/General/FastDifferentiation/stable/).
+
+Note that this type also involves the type parameter `PN`, which is actually the type of a
+`PNSystem`, and its type parameter `NT`, which will be the number type of actual numbers
+that eventually get fed into (and will be passed out from) functions that use this system.
+
+One important example of what this type is used for is computing the derivative of the
+orbital binding energy, `𝓔′` — and in particular, for generating the corresponding function
+method to apply to a given `PNSystem`.
+"""
+struct FDPNSystem{NT,PN<:PNSystem{NT},PNOrder} <: PNSystem{FDNode,Vector{FDNode},PNOrder}
+    state::Vector{FDNode}
+
+    function FDPNSystem(::Type{PN}, PNOrder=typemax(Int)) where {NT,PN<:PNSystem{NT}}
+        return new{NT,prepare_pn_order(PNOrder)}([FDNode(s) for s ∈ symbols(PN)])
+    end
+end
+
+symbols(pnsystem::FDPNSystem{NT,PN}) where {NT,PN} = symbols(PN)
+
+function symbol_index(pnsystem::FDPNSystem{NT,PN}, s::Symbol) where {NT,PN}
+    symbol_index(PN, Val(s))
+end
+
+## TODO: See if this method is needed
+
+## The old code had this, but I think it would probably just cause errors.  It might be
+## relied upon in the functions where we take derivatives — 𝓔′code and γₚₙ₀′ — but even if
+## so, maybe we could work around it with another function.
+#Base.eltype(::FDPNSystem{FT}) where {FT} = FT

src/pn_systems/NSNS.jl

@@ -0,0 +1,98 @@
+"""
+    NSNS{NT, ST, PNOrder} <: PNSystem{NT, ST, PNOrder}
+
+The [`PNSystem`](@ref) subtype describing a neutron-star—neutron-star binary system.
+
+The `state` vector is the same as for a [`BBH`](@ref), with two additional fields `Λ₁`
+and `Λ₂` holding the (constant) tidal-coupling parameters of the neutron stars.  See also
+[`BHNS`](@ref).
+"""
+struct NSNS{NT,ST,PNOrder} <: PNSystem{NT,ST,PNOrder}
+    state::ST
+
+    function NSNS{NT,ST,PNOrder}(state) where {NT,ST,PNOrder}
+        if eachindex(state) != Base.OneTo(16)
+            error(
+                "The `state` vector for `NSNS` must be indexed from 1 to 16; " *
+                "input is indexed `$(eachindex(state))`.",
+            )
+        end
+        new{NT,ST,PNOrder}(state)
+    end
+    function NSNS(;
+        M₁, M₂, χ⃗₁, χ⃗₂, v, R=Rotor(1), Φ=0, Λ₁, Λ₂, PNOrder=typemax(Int), kwargs...
+    )
+        NT, ST, PNOrder, state = prepare_system(;
+            M₁, M₂, χ⃗₁, χ⃗₂, R, v, Φ, Λ₁, Λ₂, PNOrder
+        )
+        return new{NT,ST,PNOrder}(state)
+    end
+    function NSNS(state; PNOrder=typemax(Int))
+        if eachindex(state) != Base.OneTo(16)
+            error(
+                "The `state` vector for `NSNS` must be indexed from 1 to 16; " *
+                "input is indexed `$(eachindex(state))`.",
+            )
+        end
+        NT, ST, PNOrder = eltype(state), typeof(state), prepare_pn_order(PNOrder)
+        return new{NT,ST,PNOrder}(state)
+    end
+end
+const BNS = NSNS
+
+# The following are methods of functions defined in `state_variables.jl`, specialized for
+# `NSNS` systems.
+state(pnsystem::NSNS) = pnsystem.state
+function symbols(::Type{<:NSNS})
+    (
+        :M₁,
+        :M₂,
+        :χ⃗₁ˣ,
+        :χ⃗₁ʸ,
+        :χ⃗₁ᶻ,
+        :χ⃗₂ˣ,
+        :χ⃗₂ʸ,
+        :χ⃗₂ᶻ,
+        :Rʷ,
+        :Rˣ,
+        :Rʸ,
+        :Rᶻ,
+        :v,
+        :Φ,
+        :Λ₁,
+        :Λ₂,
+    )
+end
+function ascii_symbols(::Type{<:NSNS})
+    (
+        :M1,
+        :M2,
+        :chi1x,
+        :chi1y,
+        :chi1z,
+        :chi2x,
+        :chi2y,
+        :chi2z,
+        :Rw,
+        :Rx,
+        :Ry,
+        :Rz,
+        :v,
+        :Phi,
+        :Lambda1,
+        :Lambda2,
+    )
+end
+for (i, symbol) ∈ enumerate(symbols(NSNS))
+    # This will define, e.g., `M₁(pnsystem::NSNS) = pnsystem.state[1]`.  We
+    # could do this manually, but this is more concise and less error-prone.
+    @eval begin
+        $(symbol)(pnsystem::NSNS) = @inbounds pnsystem.state[$i]
+        function symbol_index(::Type{T}, ::Val{Symbol($symbol)}) where {T<:NSNS}
+            $i
+        end
+    end
+end
+
+Λ₁(pnsystem::NSNS) = @inbounds pnsystem.state[15]
+Λ₂(pnsystem::NSNS) = @inbounds pnsystem.state[16]