update to use SVectors inside of property structs

Project.toml

@@ -1,15 +1,15 @@
 name = "ShockwaveProperties"
 uuid = "77d2bf28-a3e9-4b9c-9fcf-b85f74cc8a50"
 authors = ["Alex Fleming <[email protected]> and contributors"]
-version = "0.1.11"
+version = "0.2.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulChainRules = "f31437dd-25a7-4345-875f-756556e6935d"
 
-
 [compat]
 Unitful = "1"
 UnitfulChainRules = "0.1"

README.md

@@ -18,23 +18,23 @@ Property computation and conversion can be done via a myriad of functions. Any a
 Interrogating `ConservedProps` and `PrimitiveProps` is easily done via:
 
 - `density(u)`: Density of state `u`.
-- `momentum(u; gas)`: Momentum of state `u`
-- `velocity(u; gas)`: Velocity of state `u`
-- `mach_number(u; gas)`: Mach number of state `u`
-- `temperature(u; gas)`: Temperature of state `u`
-- `total_internal_energy_density(u; gas)`: Total (internal + kinetic) energy density of state `u`
+- `momentum(u, gas)`: Momentum of state `u`
+- `velocity(u, gas)`: Velocity of state `u`
+- `mach_number(u, gas)`: Mach number of state `u`
+- `temperature(u, gas)`: Temperature of state `u`
+- `total_internal_energy_density(u, gas)`: Total (internal + kinetic) energy density of state `u`
 
 Other properties can be directly computed:
 
-- `pressure(u; gas)`: also offers overloads for certain special cases.
-- `static_enthalpy_density(u; gas)`: Computes the static enthalpy density of `u`. Does **NOT** include kinetic energy!
-- `total_enthalpy_density(u; gas)`: Computes total enthalpy density of `u`. **DOES** include kinetic energy!
+- `pressure(u, gas)`: also offers overloads for certain special cases.
+- `static_enthalpy_density(u, gas)`: Computes the static enthalpy density of `u`. Does **NOT** include kinetic energy!
+- `total_enthalpy_density(u, gas)`: Computes total enthalpy density of `u`. **DOES** include kinetic energy!
 
 "Specific" properties are related to the mass of a state, rather than its volume.
 
-- `specific_internal_energy(u; gas)`
-- `specific_static_enthalpy(u; gas)`
-- `specific_total_enthalpy(u; gas)`
+- `specific_internal_energy(u, gas)`
+- `specific_static_enthalpy(u, gas)`
+- `specific_total_enthalpy(u, gas)`
 
 We can compute the change in properties across a shock wave from the shock normal $\hat n$ and shock tangent $\hat t$ via
 

src/ShockwaveProperties.jl

@@ -1,6 +1,7 @@
 module ShockwaveProperties
 
 using LinearAlgebra
+using StaticArrays
 using Unitful
 using UnitfulChainRules
 

src/cpg.jl

@@ -65,35 +65,39 @@ These completely determine the state of a calorically perfect gas.
 """
 struct PrimitiveProps{N,DTYPE,Q1<:Density{DTYPE},Q2<:Temperature{DTYPE}}
     ρ::Q1
-    M::NTuple{N,Float64}
+    M::SVector{N,DTYPE}
     T::Q2
 end
 
 """
-    PrimitiveProps(ρ::Density, M::{Real, DimensionlessQuantity}, T::Temperature)
+    PrimitiveProps(ρ::Real, M, T::Real)
+Construct a PrimitiveProps and assign the default units.
+"""
+function PrimitiveProps(ρ, M, T)
+    return PrimitiveProps(
+        Quantity(ρ, _units_ρ),
+        SVector{length(M)}(M),
+        Quantity(T, _units_T),
+    )
+end
+
+"""
+    PrimitiveProps(ρ::Density, M, T::Temperature)
 
 Construct a `PrimitiveProps` and strip any vestigal units from the mach number.
 """
 function PrimitiveProps(ρ::Density, M, T::Temperature)
-    return PrimitiveProps(ρ, tuple(uconvert.(NoUnits, M)...), T)
+    return PrimitiveProps(ρ, SVector{length(M)}(uconvert.(NoUnits, M)), T)
 end
 
 """
-    PrimitiveProps(ρ::Density, M::Vector{NoDims}, P::Pressure, gas::CaloricallyPerfectGas)
+    PrimitiveProps(ρ::Density, M, P::Pressure, gas::CaloricallyPerfectGas)
 
 Construct a PrimitiveProps from `[ρ, M, P]`. 
 """
 function PrimitiveProps(ρ::Density, M, P::Pressure, gas::CaloricallyPerfectGas)
     T = P / (ρ * gas.R)
-    return PrimitiveProps(ρ, tuple(uconvert.(NoUnits, M)...), T)
-end
-
-"""
-    PrimitiveProps(ρ::Real, M, T::Real)
-Construct a PrimitiveProps and assign the default units.
-"""
-function PrimitiveProps(ρ::Real, M, T::Real)
-    return PrimitiveProps(Quantity(ρ, _units_ρ), tuple(M...), Quantity(T, _units_T))
+    return PrimitiveProps(ρ, SVector{length(M)}(uconvert.(NoUnits, M)), T)
 end
 
 """
@@ -103,7 +107,7 @@ Construct a PrimitiveProps from a vector and assign default units.
 function PrimitiveProps(s::AbstractVector)
     return PrimitiveProps(
         Quantity(s[1], _units_ρ),
-        tuple(s[2:end-1]...),
+        s[2:end-1],
         Quantity(s[end], _units_T),
     )
 end
@@ -124,30 +128,38 @@ struct ConservedProps{
     U3<:EnergyDensity{DTYPE},
 }
     ρ::U1
-    ρv::NTuple{N,U2}
+    ρv::SVector{N,U2}
     ρE::U3
 end
 
 """
-    ConservedProps(ρ::Float64, ρv, ρE::Float64)
-Construct a ConservedState and assign the default units.
+    ConservedProps(ρ, ρv, ρE)
+Construct a ConservedState and assign the default units, if none are provided.
 """
-function ConservedProps(ρ::Real, ρv, ρE::Real)
+function ConservedProps(ρ, ρv, ρE)
     return ConservedProps(
         Quantity(Float64(ρ), _units_ρ),
-        tuple(Quantity.(ρv, _units_ρv)...),
+        SVector{length(ρv)}(Quantity.(ρv, _units_ρv)),
         Quantity(Float64(ρE), _units_ρE),
     )
 end
 
+function ConservedProps(
+    ρ::Density,
+    ρv::Union{Tuple{Vararg{U}},AbstractVector{U}},
+    ρE::EnergyDensity,
+) where {U<:MomentumDensity}
+    return ConservedProps(ρ, SVector{length(ρv)}(ρv), ρE)
+end
+
 """
     ConservedProps(p::AbstractVector)
 Construct a ConservedProps from a vector and assign default units.
 """
 function ConservedProps(u::AbstractVector)
     return ConservedProps(
         Quantity(u[1], _units_ρ),
-        tuple(Quantity.(u[2:end-1], _units_ρv)...),
+        Quantity.(u[2:end-1], _units_ρv),
         Quantity(u[end], _units_ρE),
     )
 end
@@ -404,21 +416,11 @@ speed_of_sound(ρ::Density, P::Pressure, gas::CaloricallyPerfectGas) = sqrt(gas.
 
 Compute the speed of sound from conserved state properties.
 """
-function speed_of_sound(
-    ρ::Real,
-    ρv,
-    ρE::Real,
-    gas::CaloricallyPerfectGas,
-)
+function speed_of_sound(ρ::Real, ρv, ρE::Real, gas::CaloricallyPerfectGas)
     return speed_of_sound(ρ, ustrip(pressure(ρ, ρv, ρE, gas)), gas)
 end
 
-function speed_of_sound(
-    ρ::Density,
-    ρv,
-    ρE::EnergyDensity;
-    gas::CaloricallyPerfectGas,
-)
+function speed_of_sound(ρ::Density, ρv, ρE::EnergyDensity, gas::CaloricallyPerfectGas)
     return speed_of_sound(ρ, pressure(ρ, ρv, ρE, gas), gas)
 end
 
@@ -462,8 +464,9 @@ end
 
 ### DISRESPECT UNITS AND WORK WITH STATES AS COLLECTIONS ###
 
-function state_to_vector(state::T) where {T}
-    return vcat(map(sym -> ustrip.(getfield(state, sym)), fieldnames(T))...)
+function state_to_vector(state)
+    v = map(sym -> ustrip.(getfield(state, sym)), fieldnames(typeof(state)))
+    return vcat(v[1], v[2]..., v[3])
 end
 
 """

test/runtests.jl

@@ -4,7 +4,139 @@ using Test
 using Unitful
 
 @testset verbose = true "ShockwaveProperties.jl" begin
-    include("test_property_correctness.jl")
+    @testset "Dimensional Analysis" begin
+        # giving different units shouldn't mess with the actual results
+        # these are the same state
+        s1 = PrimitiveProps(1.225, [2.0, 0.0], 300.0)
+        s2 = PrimitiveProps(0.002376892407u"slug/ft^3", [2.0, 0.0], 540.0u"Ra")
+        u1 = ConservedProps(s1, DRY_AIR)
+        u2 = ConservedProps(s2, DRY_AIR)
+        #check dimensions
+        for v ∈ (s1, s2, u1, u2)
+            @test pressure(v, DRY_AIR) isa Unitful.Pressure
+            @test temperature(v, DRY_AIR) isa Unitful.Temperature
+            @test (
+                specific_static_internal_energy(v, DRY_AIR) isa
+                ShockwaveProperties.SpecificEnergy
+            )
+            @test (
+                total_internal_energy_density(v, DRY_AIR) isa
+                ShockwaveProperties.EnergyDensity
+            )
+            @test speed_of_sound(v, DRY_AIR) isa Unitful.Velocity
+        end
+        # check equivalency between primitive/conserved
+        # cross test with unit change as well
+        test_items = ((s1, u1), (s2, u2), (s1, u2), (s2, u1))
+        @testset "Pairwise Equivalency $i" for i ∈ eachindex(test_items)
+            (v1, v2) = test_items[i]
+            @test pressure(v1, DRY_AIR) ≈ pressure(v2, DRY_AIR)
+            @test temperature(v1, DRY_AIR) ≈ temperature(v2, DRY_AIR)
+            @test all(
+                momentum_density(v1, DRY_AIR) ≈ momentum_density(v2, DRY_AIR),
+            )
+            @test all(velocity(v1, DRY_AIR) ≈ velocity(v2, DRY_AIR))
+            @test (
+                specific_static_internal_energy(v1, DRY_AIR) ≈
+                specific_static_internal_energy(v2, DRY_AIR)
+            )
+            @test speed_of_sound(v1, DRY_AIR) ≈ speed_of_sound(v2, DRY_AIR)
+        end
+    end
+
+    @testset "Convert Primitve ↔ Conserved" begin
+        s1 = PrimitiveProps(1.225, [2.0, 0.0], 300.0)
+        u = ConservedProps(s1, DRY_AIR)
+        s2 = PrimitiveProps(u, DRY_AIR)
+        @test s1.ρ ≈ s2.ρ
+        @test all(s1.M ≈ s2.M)
+        @test s1.T ≈ s2.T
+    end
+
+    @testset "Equivalency Between Unitful and Unitless" begin
+        n = [-1.0, 0]
+        t = [0.0, 1.0]
+        sL = PrimitiveProps(1.225, [2.0, 0.0], 300.0)
+        sR = state_behind(sL, n, t, DRY_AIR)
+        sL_nounits = state_to_vector(sL)
+        sR_nounits = primitive_state_behind(sL_nounits, n, t, DRY_AIR)
+        @testset "Primitive State Quantities" begin
+            @test sR_nounits[1] ≈ ustrip(sR.ρ)
+            @test all(sR_nounits[2:end-1] .≈ sR.M)
+            @test sR_nounits[end] ≈ ustrip(sR.T)
+            @test pressure(sL_nounits[1], sL_nounits[end], DRY_AIR) ≈
+                  pressure(sL, DRY_AIR)
+            @test pressure(sR_nounits[1], sR_nounits[end], DRY_AIR) ≈
+                  pressure(sR, DRY_AIR)
+        end
+
+        uL = ConservedProps(sL, DRY_AIR)
+        uL_nounits = state_to_vector(uL)
+        uR = state_behind(uL, n, t, DRY_AIR)
+        uR_nounits = conserved_state_behind(uL_nounits, n, t, DRY_AIR)
+
+        @testset "Conserved State Quantities" begin
+            @test uR_nounits[1] ≈ ustrip(uR.ρ)
+            @test all(uR_nounits[2:end-1] .≈ ustrip.(uR.ρv))
+            @test uR_nounits[end] ≈ ustrip(uR.ρE)
+            ρe_nounits = static_internal_energy_density(
+                uR_nounits[1],
+                uR_nounits[2:end-1],
+                uR_nounits[end],
+            )
+            @test ρe_nounits ≈ ustrip(static_internal_energy_density(uR))
+            @test pressure(ρe_nounits, DRY_AIR) ≈ pressure(uR, DRY_AIR)
+        end
+
+        @testset "Conversion" begin
+            uR_nounits2 = conserved_state_vector(sR_nounits, DRY_AIR)
+            sR_nounits2 = primitive_state_vector(uR_nounits, DRY_AIR)
+            @test uR_nounits[1] ≈ uR_nounits2[1]
+            @test all(uR_nounits .≈ uR_nounits2)
+            @test all(sR_nounits .≈ sR_nounits2)
+        end
+    end
+
+    # flux for the euler equations
+    function F(u::ConservedProps{N, T, U1, U2, U3}) where {N, T, U1, U2, U3}
+        v = velocity(u)
+        P = pressure(u, DRY_AIR)
+        return vcat(u.ρv', u.ρv * v' + I * P, (v * (u.ρE + P))')
+    end
+
+    @testset "Rankine-Hugoniot Condition" begin
+        free_stream = PrimitiveProps(1.225, [2.0, 0.0], 300.0)
+        u_L = ConservedProps(free_stream, DRY_AIR)
+        # restricted domain here because formula from Anderson&Anderson
+        # breaks down near β = 0
+        # TODO investigate this.
+        for θ ∈ range(0, π / 3; length = 40)
+            n = [-cos(θ), sin(θ)]
+            t = [0 1; -1 0] * n
+            u_R = state_behind(u_L, n, t, DRY_AIR)
+            # F(u_l)⋅n̂ - F(u_r)⋅n̂ = 0 ⟹ F(u_l)⋅n̂ = F(u_r)⋅n̂ 
+            @test all(F(u_L) * n .≈ F(u_R) * n)
+        end
+    end
+
+    @testset "Speed of Sound" begin
+        gas = DRY_AIR
+        s = PrimitiveProps(1.225, [2.0, 0.0], 300.0)
+        u = ConservedProps(s, gas)
+
+        a_s = speed_of_sound(s, gas)
+        a_u = speed_of_sound(u, gas)
+
+        a_s_pressure = speed_of_sound(density(s), pressure(s, gas), gas)
+        a_u_pressure = speed_of_sound(density(u), pressure(u, gas), gas)
+
+        a_s_ie = speed_of_sound(density(s), momentum_density(s, gas), total_internal_energy_density(s, gas), gas)
+        a_u_ie = speed_of_sound(density(u), momentum_density(u), total_internal_energy_density(u), gas)
+
+        @test a_s ≈ a_u
+        @test a_s_pressure ≈ a_u_pressure
+        @test a_s_ie ≈ a_u_ie
+    end
 
     @testset verbose = true "Billig Shockwave Parametrization" begin
         using .BilligShockParametrization