From 6f82873baeb6085064bbfafdba78fd19decf8646 Mon Sep 17 00:00:00 2001 From: Aaron Kaw Date: Tue, 6 Aug 2024 22:23:52 +1000 Subject: [PATCH] Ensemble ray tracing --- src/01_oceanography/oceanography.jl | 2 + src/02_acoustics/03_propagation/tracing.jl | 56 ------------ .../03_propagation/tracing/03_beam.jl | 91 +++++++++++++++++++ .../03_propagation/tracing/04_fan.jl | 45 +++++++++ .../03_propagation/tracing/05_grid.jl | 0 src/OceanSonar.jl | 13 +++ 6 files changed, 151 insertions(+), 56 deletions(-) create mode 100644 src/01_oceanography/oceanography.jl delete mode 100644 src/02_acoustics/03_propagation/tracing.jl create mode 100644 src/02_acoustics/03_propagation/tracing/03_beam.jl create mode 100644 src/02_acoustics/03_propagation/tracing/04_fan.jl create mode 100644 src/02_acoustics/03_propagation/tracing/05_grid.jl diff --git a/src/01_oceanography/oceanography.jl b/src/01_oceanography/oceanography.jl new file mode 100644 index 0000000..9e1f5bd --- /dev/null +++ b/src/01_oceanography/oceanography.jl @@ -0,0 +1,2 @@ +"Circumference [m] of earth's equator." +const EQUATORIAL_EARTH_CIRCUMFERENCE = 40_075e3 \ No newline at end of file diff --git a/src/02_acoustics/03_propagation/tracing.jl b/src/02_acoustics/03_propagation/tracing.jl deleted file mode 100644 index 71ba7f6..0000000 --- a/src/02_acoustics/03_propagation/tracing.jl +++ /dev/null @@ -1,56 +0,0 @@ -function TracingODESystem( - c_ocn::Function, - r_max::Real; - name -) - @parameters begin - s, [description = "Ray arc length"] - end - - pars = @parameters begin - r₀ = 0.0 - z₀ - θ₀ - end - - @variables begin - θ(s), [description = "Ray angle, positive downwards"] - c(s), [description = "Ray celerity"] - end - - deps = @variables begin - r(s) = r₀ - z(s) = z₀ - ξ(s) = cos(θ₀) / c_ocn(r₀, z₀) - ζ(s) = sin(θ₀) / c_ocn(r₀, z₀) - end - - Ds = Differential(s) - - c²(r′, z′) = c_ocn(r′, z′)^2 - ∂c_∂r(r′, z′) = ModelingToolkit.derivative(c_ocn(r, z′), r′) - ∂c_∂z(r′, z′) = ModelingToolkit.derivative(c_ocn(r′, z), z′) - - eqns = [ - Ds(r) ~ c_ocn(r, z) * ξ - Ds(z) ~ c_ocn(r, z) * ζ - Ds(ξ) ~ -∂c_∂r(r, z) / c_ocn(r, z)^2 - Ds(ζ) ~ -∂c_∂z(r, z) / c_ocn(r, z)^2 - θ ~ atan(ζ, ξ) - c ~ c_ocn(r, z) - ] - - maximum_range = [r ~ r_max] => ( - (ntg, _, _, _) -> terminate!(ntg), - [], [], [] - ) - - events = [ - maximum_range - ] - - return ODESystem(eqns, s, deps, pars; - name = name, - continuous_events = events - ) -end diff --git a/src/02_acoustics/03_propagation/tracing/03_beam.jl b/src/02_acoustics/03_propagation/tracing/03_beam.jl new file mode 100644 index 0000000..8c1cd66 --- /dev/null +++ b/src/02_acoustics/03_propagation/tracing/03_beam.jl @@ -0,0 +1,91 @@ +export Beam + +const DEFAULT_RAY_ARC_LENGTH_SPAN = (0, EQUATORIAL_EARTH_CIRCUMFERENCE) + +function AcousticTracingODESystem( + c_ocn::Function, + r_max::Real; + name +) + @parameters begin + s, [description = "Ray arc length"] + end + + pars = @parameters begin + r₀ = 0.0 + z₀ + θ₀ + end + + deps = @variables begin + r(s) = r₀ + z(s) = z₀ + ξ(s) = cos(θ₀) / c_ocn(r₀, z₀) + ζ(s) = sin(θ₀) / c_ocn(r₀, z₀) + θ(s), [description = "Ray angle, positive downwards"] + c(s), [description = "Ray celerity"] + end + + Ds = Differential(s) + + c²(r′, z′) = c_ocn(r′, z′)^2 + ∂c_∂r(r′, z′) = ModelingToolkit.derivative(c_ocn(r, z′), r′) + ∂c_∂z(r′, z′) = ModelingToolkit.derivative(c_ocn(r′, z), z′) + + eqns = [ + Ds(r) ~ c_ocn(r, z) * ξ + Ds(z) ~ c_ocn(r, z) * ζ + Ds(ξ) ~ -∂c_∂r(r, z) / c_ocn(r, z)^2 + Ds(ζ) ~ -∂c_∂z(r, z) / c_ocn(r, z)^2 + θ ~ atan(ζ, ξ) + c ~ c_ocn(r, z) + ] + + maximum_range = [r ~ r_max] => ( + (ntg, _, _, _) -> terminate!(ntg), + [], [], [] + ) + + events = [ + maximum_range + ] + + return ODESystem(eqns, s, deps, pars; + name = name, + continuous_events = events + ) +end + +struct Beam{MN <: ModelName} + model::MN + + s_max::Float64 + + c::Function + θ::Function + r::Function + z::Function + ξ::Function + ζ::Function + + sol::ODESolution + + function Beam(model::ModelName, sys::ODESystem, sol::ODESolution) + r(s::Real) = sol(s, idxs = sys.r) + r(s::AbstractVector{<:Real}) = sol(s, idxs = sys.r) |> collect + z(s::Real) = sol(s, idxs = sys.z) + z(s::AbstractVector{<:Real}) = sol(s, idxs = sys.z) |> collect + ξ(s::Real) = sol(s, idxs = sys.ξ) + ξ(s::AbstractVector{<:Real}) = sol(s, idxs = sys.ξ) |> collect + ζ(s::Real) = sol(s, idxs = sys.ζ) + ζ(s::AbstractVector{<:Real}) = sol(s, idxs = sys.ζ) |> collect + c(s::Real) = sol(s, idxs = sys.c) + c(s::AbstractVector{<:Real}) = sol(s, idxs = sys.c) |> collect + θ(s::Real) = sol(s, idxs = sys.θ) + θ(s::AbstractVector{<:Real}) = sol(s, idxs = sys.θ) |> collect + + new{model |> typeof}(model, sol.t[end], c, θ, r, z, ξ, ζ, sol) + end +end + +show(io::IO, beam::Beam) = print(io, "Beam($(beam.θ(0) |> rad2deg)°)") diff --git a/src/02_acoustics/03_propagation/tracing/04_fan.jl b/src/02_acoustics/03_propagation/tracing/04_fan.jl new file mode 100644 index 0000000..7c70600 --- /dev/null +++ b/src/02_acoustics/03_propagation/tracing/04_fan.jl @@ -0,0 +1,45 @@ +export Fan + +struct Fan{MN <: ModelName} + beams::Vector{Beam{MN}} + sys::ODESystem + prob::EnsembleProblem + + function Fan( + model::ModelName, + θ₀s::AbstractVector{<:Real}, + z₀::Real, + r_max::Real, + c_ocn::Function + ) + @mtkbuild sys = AcousticTracingODESystem(c_ocn, r_max) + prob = ODEProblem(sys, [], DEFAULT_RAY_ARC_LENGTH_SPAN, [sys.θ₀ => 0.0, sys.z₀ => z₀]) + + # Attempt 1: Errors on `init`. + # tiltray! = setp(sys, sys.θ₀) + # prob_func(internal_prob, n, _) = tiltray!(internal_prob, θ₀s[n]) + + # Attempt 2: Doesn't change launch angle parameter + # function prob_func(internal_prob, n, _) + # new_prob = remake(internal_prob, p = [sys.θ₀ => θ₀s[n], sys.z₀ => z₀]) + # @show prob.p[1] + # return new_prob + # end + + # Attempt 3: Make an entirely new ODEProblem + prob_func(_, n, _) = ODEProblem(sys, [], DEFAULT_RAY_ARC_LENGTH_SPAN, [sys.θ₀ => θ₀s[n], sys.z₀ => z₀]) + + ens_prob = EnsembleProblem(prob, prob_func = prob_func) + ens_sol = solve(ens_prob, Tsit5(), EnsembleThreads(), + reltol = 1e-50, + trajectories = length(θ₀s) + ) + beams = [Beam(model, sys, sol) for sol in ens_sol] + new{model |> typeof}(beams, sys, ens_prob) + end +end + +function show(io::IO, fan::Fan) + xtr = [beam.θ(0) for beam in fan.beams] |> extrema .|> rad2deg + print(io, "Fan($(fan.beams |> length): $(xtr |> minimum)° .. $(xtr |> maximum)°)") +end diff --git a/src/02_acoustics/03_propagation/tracing/05_grid.jl b/src/02_acoustics/03_propagation/tracing/05_grid.jl new file mode 100644 index 0000000..e69de29 diff --git a/src/OceanSonar.jl b/src/OceanSonar.jl index d9dc947..84faf54 100644 --- a/src/OceanSonar.jl +++ b/src/OceanSonar.jl @@ -2,17 +2,30 @@ module OceanSonar export subtypes +import Base: getproperty, show + using InteractiveUtils: subtypes using ModelingToolkit: ModelingToolkit, Differential, + EnsembleProblem, + EnsembleSolution, + EnsembleThreads, + ODEProblem, + ODESolution, ODESystem, + remake, + # setp, terminate!, @mtkbuild, @parameters, @variables +using OrdinaryDiffEq: + solve, + Tsit5 + """ ``` include_subroots(path::AbstractString)