From 99251b1f44ce91be9c74bbf9764cbca914cb91ba Mon Sep 17 00:00:00 2001 From: Severin Denk Date: Tue, 13 May 2025 13:32:15 -0700 Subject: [PATCH 01/10] This adds a routine for ray intersection with a torus --- src/geometry.jl | 104 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) diff --git a/src/geometry.jl b/src/geometry.jl index e2028da3..cbcef41d 100644 --- a/src/geometry.jl +++ b/src/geometry.jl @@ -3,6 +3,7 @@ document[Symbol("Geometry")] = Symbol[] import StaticArrays import MillerExtendedHarmonic: MXH import LinearAlgebra +import Roots """ centroid(x::AbstractVector{<:T}, y::AbstractVector{<:T}) where {T<:Real} @@ -82,6 +83,109 @@ end @compat public revolution_volume push!(document[Symbol("Geometry")], :revolution_volume) +""" + ray_rect_intersection(O, D, rect_min, rect_max) + +Compute the closest intersection point between the ray + + P(t) = O + t*D, t ≥ 0 + +and the axis-aligned rectangle with corners `rect_min = (x_min,y_min)` and +`rect_max = (x_max,y_max)`. Returns `(x,y)` of the hit point, or `nothing` +if there is no intersection. +""" +function ray_rect_intersection( + O::NTuple{2,Float64}, + D::NTuple{2,Float64}, + rect_min::NTuple{2,Float64}, + rect_max::NTuple{2,Float64}, +) + x0, y0 = O + dx, dy = D + x_min, y_min = rect_min + x_max, y_max = rect_max + + # Compute slab intersections for x + if dx != 0.0 + tx1, tx2 = (x_min - x0)/dx, (x_max - x0)/dx + t_x_min, t_x_max = min(tx1,tx2), max(tx1,tx2) + else + # Ray parallel to x-planes: must lie within the slab + if x0 < x_min || x0 > x_max + return nothing + end + t_x_min, t_x_max = -Inf, Inf + end + + # Compute slab intersections for y + if dy != 0.0 + ty1, ty2 = (y_min - y0)/dy, (y_max - y0)/dy + t_y_min, t_y_max = min(ty1,ty2), max(ty1,ty2) + else + if y0 < y_min || y0 > y_max + return nothing + end + t_y_min, t_y_max = -Inf, Inf + end + + # Find overlap of the two t-intervals + t_enter = max(t_x_min, t_y_min) + t_exit = min(t_x_max, t_y_max) + + # No intersection if empty interval or entirely behind the ray + if t_enter > t_exit || t_exit < 0.0 + return nothing + end + + # Use the entry point if in front of origin, else the exit point + t_hit = t_enter ≥ 0.0 ? t_enter : t_exit + + return (x0 + t_hit*dx, y0 + t_hit*dy) +end + +""" + ray_rect_torus_intersect(origin, direction, ρ_bounds, z_bounds) + +Finds intersection of a ray with a rectangular torus defined in cylindrical coordinates. + +- `origin`: (x, y, z) vector +- `direction`: normalized (dx, dy, dz) +- `ρ_bounds`: (ρ_min, ρ_max) +- `z_bounds`: (z_min, z_max) + +Returns (t, point) or `nothing` if no intersection. +""" +function ray_rect_torus_intersect(origin, direction, ρ_bounds, z_bounds) + x0, y0, z0 = origin + dx, dy, dz = direction + ρ_min, ρ_max = ρ_bounds + z_min, z_max = z_bounds + + # Parametrize ρ(t) and z(t) + ρ(t) = hypot(x0 + dx * t, y0 + dy * t) + z(t) = z0 + dz * t + + # Solve for t where z(t) in [z_min, z_max] + t_z1 = dz != 0 ? (z_min - z0) / dz : -Inf + t_z2 = dz != 0 ? (z_max - z0) / dz : Inf + tmin, tmax = sort([t_z1, t_z2]) + + # Search for a t in [tmin, tmax] such that ρ(t) ∈ [ρ_min, ρ_max] + ts = range(tmin, tmax; length=1000) + + for t in ts + ρt = ρ(t) + zt = z(t) + if ρt ≥ ρ_min && ρt ≤ ρ_max && zt ≥ z_min && zt ≤ z_max + point = origin .+ t .* direction + return t, point + end + end + + return nothing +end + + """ intersection_angles( path1_r::AbstractVector{T}, From 211599093b0b546c3401e6b6c846d031d55f4562 Mon Sep 17 00:00:00 2001 From: Severin Denk Date: Wed, 14 May 2025 11:32:48 -0700 Subject: [PATCH 02/10] Sneak in little addition for new source --- src/physics/sources.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/physics/sources.jl b/src/physics/sources.jl index e101bed0..7ceed980 100644 --- a/src/physics/sources.jl +++ b/src/physics/sources.jl @@ -625,6 +625,7 @@ function new_source( source.identifier.name = name source.identifier.index = index cs1d = resize!(source.profiles_1d) + csglbl = resize!(source.global_quantities) cs1d.grid.rho_tor_norm = rho cs1d.grid.volume = volume cs1d.grid.area = area @@ -634,11 +635,13 @@ function new_source( cs1d.electrons.power_inside = cumtrapz(volume, value) elseif electrons_power_inside !== missing cs1d.electrons.power_inside = value = electrons_power_inside + cs1d.electrons.energy = gradient(volume, value) else cs1d.electrons.energy = zero(volume) cs1d.electrons.power_inside = zero(volume) end + csglbl.electrons.power = cs1d.electrons.power_inside[end] if total_ion_energy !== missing cs1d.total_ion_energy = value = total_ion_energy @@ -650,6 +653,8 @@ function new_source( cs1d.total_ion_energy = zero(volume) cs1d.total_ion_power_inside = zero(volume) end + csglbl.total_ion_power = cs1d.total_ion_power_inside[end] + csglbl.power = csglbl.total_ion_power + csglbl.electrons.power if electrons_particles !== missing cs1d.electrons.particles = value = electrons_particles @@ -661,6 +666,7 @@ function new_source( cs1d.electrons.particles = zero(volume) cs1d.electrons.particles_inside = zero(volume) end + csglbl.electrons.particles = cs1d.electrons.particles_inside[end] if j_parallel !== missing cs1d.j_parallel = value = j_parallel @@ -672,6 +678,7 @@ function new_source( cs1d.j_parallel = zero(area) cs1d.current_parallel_inside = zero(area) end + csglbl.current_parallel = cs1d.current_parallel_inside[end] if momentum_tor !== missing cs1d.momentum_tor = value = momentum_tor @@ -683,6 +690,7 @@ function new_source( cs1d.momentum_tor = zero(volume) cs1d.torque_tor_inside = zero(volume) end + csglbl.torque_tor = cs1d.torque_tor_inside[end] return source end From 4fd08d158ca072c82bc0ab67a333070b40efcea8 Mon Sep 17 00:00:00 2001 From: Severin Denk Date: Mon, 2 Jun 2025 16:45:54 -0700 Subject: [PATCH 03/10] Add tests and fix bugs --- src/geometry.jl | 129 ++++++++++++++++++---- test/geometry/test_ray_torus_intersect.jl | 16 +++ test/geometry/test_solve_r_intersect.jl | 24 ++++ 3 files changed, 145 insertions(+), 24 deletions(-) create mode 100644 test/geometry/test_ray_torus_intersect.jl create mode 100644 test/geometry/test_solve_r_intersect.jl diff --git a/src/geometry.jl b/src/geometry.jl index aaaf8fd0..2f0ad622 100644 --- a/src/geometry.jl +++ b/src/geometry.jl @@ -142,46 +142,127 @@ function ray_rect_intersection( return (x0 + t_hit*dx, y0 + t_hit*dy) end +""" + r_intersect_interval(x0::T, y0::T, dx::T, dy::T, r_min::T, r_max::T) where {T<:Real} + +Finds intersection of ray starting at x0 and y0 with direction (dx,dy) + with reference major radius value r_ref + +- `x0``, `y0``: Origin +- `dx``, `dy``: Normalized direction in x-y plane +- `r_min`, `r_max`: Edges of r_range + +Returns interval [t1, t2] for which r_min <= r(t) <= r_max with r(t)^2 = (x+dx*t)^2 + (y+dy*t)^2 """ - ray_rect_torus_intersect(origin, direction, ρ_bounds, z_bounds) + +function solve_r_intersect(x0::T, y0::T, dx::T, dy::T, r_min::T, r_max::T) where {T<:Real} + r0 = sqrt(x0^2 + y0^2) + t_crossings = zeros(Float64, 4) + crossing = zeros(Float64, 4) # +1 -> into domain + # -1 -> out of + # 0 -> no crossing + into_domain = [1.0, -1.0] + if dx == 0 && dy == 0 + # Handle vertical ray + if (r0 > r_min && r0 < r_max) + return [[0.0 Inf];] + else + return zeros(Float64, 0) + end + end + for (i_ref, r_ref) in enumerate([r_min, r_max]) + Δ = r_ref^2 * (dx^2 + dy^2) - x0^2 * dy^2 - y0^2 * dx^2 + 2 * x0 * y0* dx*dy + if Δ < 0 + # println("No intersection for ", r_ref) + t_crossings[2*(i_ref-1) + 1] = Inf + t_crossings[2*(i_ref-1) + 2] = Inf + crossing[2*(i_ref-1) + 1] = 0.0 + crossing[2*(i_ref-1) + 2] = 0.0 + continue + end + t1 = -(x0*dx + y0 * dy) - sqrt(Δ) + t2 = -(x0*dx + y0 * dy) + sqrt(Δ) + t_crossings[2*(i_ref-1) + 1] = t1 >= 0 ? t1/ (dx^2 + dy^2) : Inf + t_crossings[2*(i_ref-1) + 2] = t2 >= 0 ? t2/ (dx^2 + dy^2) : Inf + crossing[2*(i_ref-1) + 1] = t1 >= 0 ? into_domain[i_ref]*sign(x0*dx + dx^2*t1 + y0 * dy + dy^2*t1) : 0.0 + crossing[2*(i_ref-1) + 2] = t2 >= 0 ? into_domain[i_ref]*sign(x0*dx + dx^2*t2 + y0 * dy + dy^2*t2) : 0.0 + end + last = 0 # 1 in, 0 out + # If the first point is in we add it to the crossings + if (r0 > r_min && r0 < r_max) + append!(t_crossings, 0.0) + append!(crossing, 1.0) + end + i_sort = sortperm(t_crossings) + intervals = zeros(Float64, 0) + for i in i_sort + if t_crossings[i] == Inf + continue + end + if last == 0 && crossing[i] > 0.0 # looking for a crossing into the rectangle + append!(intervals, t_crossings[i]) + last = 1 + elseif last == 1 && crossing[i] < 0.0 # looking for a crossing leaving the rectangle + append!(intervals, t_crossings[i]) + last = 0 + end + end + if length(intervals) % 2 != 0 + throw(ErrorException("Somehow only found odd number of intersections which is impossible.")) + end + if length(intervals) > 2 # Return as shape (2,2). Need to transpose because column-major... + return permutedims(reshape(intervals, :, 2)) + elseif length(intervals) > 0 + return reshape(intervals, :, 2) # Return as shape (1,2) + else + return intervals # Return empty + end +end + + +""" + ray_torus_intersect(origin, direction, ρ_bounds, z_bounds) Finds intersection of a ray with a rectangular torus defined in cylindrical coordinates. - `origin`: (x, y, z) vector - `direction`: normalized (dx, dy, dz) -- `ρ_bounds`: (ρ_min, ρ_max) +- `r_bounds`: (r_min, r_max) - `z_bounds`: (z_min, z_max) -Returns (t, point) or `nothing` if no intersection. +Returns intersection point between ray closest to the origin or `nothing` if no intersection. """ -function ray_rect_torus_intersect(origin, direction, ρ_bounds, z_bounds) +function ray_torus_intersect(origin, direction, r_bounds, z_bounds) x0, y0, z0 = origin - dx, dy, dz = direction - ρ_min, ρ_max = ρ_bounds + if norm(direction) == 0 + throw(ArgumentError("Direction most not have norm zero")) + end + dx, dy, dz = direction/norm(direction) + r_min, r_max = r_bounds z_min, z_max = z_bounds + # Intervals in between the ray passes throught the r limits of the box + t_intervals_R = solve_r_intersect(x0, y0, dx, dy, r_min, r_max) + if dz == 0 + if z0 < z_min || z0 > z_max + return nothing, nothing + else + t_interval_z = [0.0 Inf] + end + else + t_interval_z = sort![(z_min - z0) / dz, (z_max - z0) / dz] + end + t_interval_z[1] = t_interval_z[1] > 0.0 ? t_interval_z[1] : 0.0 # Remove intersection at negative values + # Parametrize ρ(t) and z(t) - ρ(t) = hypot(x0 + dx * t, y0 + dy * t) - z(t) = z0 + dz * t - - # Solve for t where z(t) in [z_min, z_max] - t_z1 = dz != 0 ? (z_min - z0) / dz : -Inf - t_z2 = dz != 0 ? (z_max - z0) / dz : Inf - tmin, tmax = sort([t_z1, t_z2]) - - # Search for a t in [tmin, tmax] such that ρ(t) ∈ [ρ_min, ρ_max] - ts = range(tmin, tmax; length=1000) - - for t in ts - ρt = ρ(t) - zt = z(t) - if ρt ≥ ρ_min && ρt ≤ ρ_max && zt ≥ z_min && zt ≤ z_max - point = origin .+ t .* direction - return t, point + for t_interval_R in eachrow(t_intervals_R) + t_min = max(t_interval_R[1], t_interval_z[1]) + t_max = min(t_interval_R[2], t_interval_z[2]) + if t_min < t_max + return [x0 y0 z0] + [dx dy dz] * t_min end end - return nothing end diff --git a/test/geometry/test_ray_torus_intersect.jl b/test/geometry/test_ray_torus_intersect.jl new file mode 100644 index 00000000..be575cdf --- /dev/null +++ b/test/geometry/test_ray_torus_intersect.jl @@ -0,0 +1,16 @@ +using IMAS +using Test + +@testset "Test ray_torus_intersect" begin + test_cases = [ + ([6.0,0.0,0.0], [-1.0,0.0,0.0], [1.0, 4.0], [-1.0, 1.0], [4.0 0.0 0.0], ) + ] + for (origin, direction, r_bounds, z_bounds, expected) in test_cases + @testset "origin=$origin, direction=$direction, r_bounds=$r_bounds, z_bounds=$z_bounds" begin + # println("--- x0=$x0, y0=$y0, dx=$dx, dy=$dy, r_min=$r_min, r_max=$r_max ---") + test_result = IMAS.ray_torus_intersect(origin, direction, r_bounds, z_bounds) + println("--- Result | expected: $test_result | $expected ---") + @test test_result ≈ expected + end + end +end \ No newline at end of file diff --git a/test/geometry/test_solve_r_intersect.jl b/test/geometry/test_solve_r_intersect.jl new file mode 100644 index 00000000..34b39fc6 --- /dev/null +++ b/test/geometry/test_solve_r_intersect.jl @@ -0,0 +1,24 @@ +using IMAS +using Test + +@testset "Test solve_r_intersect" begin + test_cases = [ + (3.0, 0.0, 0.0, 0.0, 1.0, 4.0, [[0.0 Inf];]), # ray is not moving, but starts on the grid + (6.0, 0.0, 0.0, 0.0, 1.0, 4.0, zeros(Float64, 0)), # ray is not moving, and starts off the grid + (6.0, 0.0, -1.0, 0.0, 1.0, 4.0, [2. 5.; 7. 10.]), # ray in x-direction + (0.0, 6.0, 0.0, -1.0, 2.0, 4.0, [2. 4.; 8. 10.]), # ray in y-direction + (3.0, 0.0, -1.0, 0.0, 1.0, 4.0, [0.0 2.; 4.0 7.0]), # ray in x-direction with start in torus + (6.0, 6.0, 0.0, -1.0, 2.0, 4.0, zeros(Float64, 0)), # ray in x-direction with offset causing no intersection + (3.0, 6.0, -3.0/5.0, -4.0/5.0, 1.0, 2.0, [[5.0 41.0/5.0];]), # Single crossing with x-y ray + (3.0, 8.0, -3.0/5.0, -4.0/5.0, 3.0, 4.0, [5.0 32.0/5.0; 10.0 57.0/5.0]) # double crossing with x-y ray + ] + + for (x0, y0, dx, dy, r_min, r_max, expected) in test_cases + @testset "x0=$x0, y0=$y0, dx=$dx, dy=$dy, r_min=$r_min, r_max=$r_max" begin + # println("--- x0=$x0, y0=$y0, dx=$dx, dy=$dy, r_min=$r_min, r_max=$r_max ---") + test_result = IMAS.solve_r_intersect(x0, y0, dx, dy, r_min, r_max) + println("--- Result | expected: $test_result | $expected ---") + @test test_result ≈ expected + end + end +end \ No newline at end of file From da464d5591f26aa00cb0ed5a1283f4eb5b922637 Mon Sep 17 00:00:00 2001 From: Severin Denk Date: Tue, 17 Jun 2025 11:38:29 -0700 Subject: [PATCH 04/10] Remove uneeded, untested routine --- src/geometry.jl | 60 ------------------------------------------------- 1 file changed, 60 deletions(-) diff --git a/src/geometry.jl b/src/geometry.jl index 81d2aec0..df55cc6b 100644 --- a/src/geometry.jl +++ b/src/geometry.jl @@ -82,66 +82,6 @@ end @compat public revolution_volume push!(document[Symbol("Geometry")], :revolution_volume) - -""" - ray_rect_intersection(O, D, rect_min, rect_max) - -Compute the closest intersection point between the ray - - P(t) = O + t*D, t ≥ 0 - -and the axis-aligned rectangle with corners `rect_min = (x_min,y_min)` and -`rect_max = (x_max,y_max)`. Returns `(x,y)` of the hit point, or `nothing` -if there is no intersection. -""" -function ray_rect_intersection( - O::NTuple{2,Float64}, - D::NTuple{2,Float64}, - rect_min::NTuple{2,Float64}, - rect_max::NTuple{2,Float64}, -) - x0, y0 = O - dx, dy = D - x_min, y_min = rect_min - x_max, y_max = rect_max - - # Compute slab intersections for x - if dx != 0.0 - tx1, tx2 = (x_min - x0)/dx, (x_max - x0)/dx - t_x_min, t_x_max = min(tx1,tx2), max(tx1,tx2) - else - # Ray parallel to x-planes: must lie within the slab - if x0 < x_min || x0 > x_max - return nothing - end - t_x_min, t_x_max = -Inf, Inf - end - - # Compute slab intersections for y - if dy != 0.0 - ty1, ty2 = (y_min - y0)/dy, (y_max - y0)/dy - t_y_min, t_y_max = min(ty1,ty2), max(ty1,ty2) - else - if y0 < y_min || y0 > y_max - return nothing - end - t_y_min, t_y_max = -Inf, Inf - end - - # Find overlap of the two t-intervals - t_enter = max(t_x_min, t_y_min) - t_exit = min(t_x_max, t_y_max) - - # No intersection if empty interval or entirely behind the ray - if t_enter > t_exit || t_exit < 0.0 - return nothing - end - - # Use the entry point if in front of origin, else the exit point - t_hit = t_enter ≥ 0.0 ? t_enter : t_exit - - return (x0 + t_hit*dx, y0 + t_hit*dy) -end """ r_intersect_interval(x0::T, y0::T, dx::T, dy::T, r_min::T, r_max::T) where {T<:Real} From eb2edce8ce8827befa9c2445e8db592c54fb2cc6 Mon Sep 17 00:00:00 2001 From: Severin Denk Date: Tue, 17 Jun 2025 11:49:58 -0700 Subject: [PATCH 05/10] Make sure geometry tests are run --- test/runtests.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index e5db99ed..67c260d6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,4 +11,8 @@ else include("runtests_extract.jl") include("runtests_plot_recipes.jl") + + include("geometry/test_ray_torus_intersect.jl") + + include("geometry/test_solve_r_intersect.jl") end From 544e4ac8e9535c1c0f27786d5bada36fdada7086 Mon Sep 17 00:00:00 2001 From: Severin Denk Date: Tue, 17 Jun 2025 16:54:51 -0700 Subject: [PATCH 06/10] Welp theere goes my contribution --- src/geometry.jl | 124 ------------------ src/physics/particles.jl | 14 ++ test/geometry/test_ray_torus_intersect.jl | 16 --- test/geometry/test_solve_r_intersect.jl | 24 ---- .../particles/test_toroidal_intersection.jl | 16 +++ test/runtests.jl | 6 +- 6 files changed, 32 insertions(+), 168 deletions(-) delete mode 100644 test/geometry/test_ray_torus_intersect.jl delete mode 100644 test/geometry/test_solve_r_intersect.jl create mode 100644 test/physics/particles/test_toroidal_intersection.jl diff --git a/src/geometry.jl b/src/geometry.jl index df55cc6b..862f88cf 100644 --- a/src/geometry.jl +++ b/src/geometry.jl @@ -82,130 +82,6 @@ end @compat public revolution_volume push!(document[Symbol("Geometry")], :revolution_volume) -""" - r_intersect_interval(x0::T, y0::T, dx::T, dy::T, r_min::T, r_max::T) where {T<:Real} - -Finds intersection of ray starting at x0 and y0 with direction (dx,dy) - with reference major radius value r_ref - -- `x0``, `y0``: Origin -- `dx``, `dy``: Normalized direction in x-y plane -- `r_min`, `r_max`: Edges of r_range - -Returns interval [t1, t2] for which r_min <= r(t) <= r_max with r(t)^2 = (x+dx*t)^2 + (y+dy*t)^2 - -""" - -function solve_r_intersect(x0::T, y0::T, dx::T, dy::T, r_min::T, r_max::T) where {T<:Real} - r0 = sqrt(x0^2 + y0^2) - t_crossings = zeros(Float64, 4) - crossing = zeros(Float64, 4) # +1 -> into domain - # -1 -> out of - # 0 -> no crossing - into_domain = [1.0, -1.0] - if dx == 0 && dy == 0 - # Handle vertical ray - if (r0 > r_min && r0 < r_max) - return [[0.0 Inf];] - else - return zeros(Float64, 0) - end - end - for (i_ref, r_ref) in enumerate([r_min, r_max]) - Δ = r_ref^2 * (dx^2 + dy^2) - x0^2 * dy^2 - y0^2 * dx^2 + 2 * x0 * y0* dx*dy - if Δ < 0 - # println("No intersection for ", r_ref) - t_crossings[2*(i_ref-1) + 1] = Inf - t_crossings[2*(i_ref-1) + 2] = Inf - crossing[2*(i_ref-1) + 1] = 0.0 - crossing[2*(i_ref-1) + 2] = 0.0 - continue - end - t1 = -(x0*dx + y0 * dy) - sqrt(Δ) - t2 = -(x0*dx + y0 * dy) + sqrt(Δ) - t_crossings[2*(i_ref-1) + 1] = t1 >= 0 ? t1/ (dx^2 + dy^2) : Inf - t_crossings[2*(i_ref-1) + 2] = t2 >= 0 ? t2/ (dx^2 + dy^2) : Inf - crossing[2*(i_ref-1) + 1] = t1 >= 0 ? into_domain[i_ref]*sign(x0*dx + dx^2*t1 + y0 * dy + dy^2*t1) : 0.0 - crossing[2*(i_ref-1) + 2] = t2 >= 0 ? into_domain[i_ref]*sign(x0*dx + dx^2*t2 + y0 * dy + dy^2*t2) : 0.0 - end - last = 0 # 1 in, 0 out - # If the first point is in we add it to the crossings - if (r0 > r_min && r0 < r_max) - append!(t_crossings, 0.0) - append!(crossing, 1.0) - end - i_sort = sortperm(t_crossings) - intervals = zeros(Float64, 0) - for i in i_sort - if t_crossings[i] == Inf - continue - end - if last == 0 && crossing[i] > 0.0 # looking for a crossing into the rectangle - append!(intervals, t_crossings[i]) - last = 1 - elseif last == 1 && crossing[i] < 0.0 # looking for a crossing leaving the rectangle - append!(intervals, t_crossings[i]) - last = 0 - end - end - if length(intervals) % 2 != 0 - throw(ErrorException("Somehow only found odd number of intersections which is impossible.")) - end - if length(intervals) > 2 # Return as shape (2,2). Need to transpose because column-major... - return permutedims(reshape(intervals, :, 2)) - elseif length(intervals) > 0 - return reshape(intervals, :, 2) # Return as shape (1,2) - else - return intervals # Return empty - end -end - - -""" - ray_torus_intersect(origin, direction, ρ_bounds, z_bounds) - -Finds intersection of a ray with a rectangular torus defined in cylindrical coordinates. - -- `origin`: (x, y, z) vector -- `direction`: normalized (dx, dy, dz) -- `r_bounds`: (r_min, r_max) -- `z_bounds`: (z_min, z_max) - -Returns intersection point between ray closest to the origin or `nothing` if no intersection. -""" -function ray_torus_intersect(origin, direction, r_bounds, z_bounds) - x0, y0, z0 = origin - if norm(direction) == 0 - throw(ArgumentError("Direction most not have norm zero")) - end - dx, dy, dz = direction/norm(direction) - r_min, r_max = r_bounds - z_min, z_max = z_bounds - # Intervals in between the ray passes throught the r limits of the box - t_intervals_R = solve_r_intersect(x0, y0, dx, dy, r_min, r_max) - if dz == 0 - if z0 < z_min || z0 > z_max - return nothing, nothing - else - t_interval_z = [0.0 Inf] - end - else - t_interval_z = sort![(z_min - z0) / dz, (z_max - z0) / dz] - end - t_interval_z[1] = t_interval_z[1] > 0.0 ? t_interval_z[1] : 0.0 # Remove intersection at negative values - - - # Parametrize ρ(t) and z(t) - for t_interval_R in eachrow(t_intervals_R) - t_min = max(t_interval_R[1], t_interval_z[1]) - t_max = min(t_interval_R[2], t_interval_z[2]) - if t_min < t_max - return [x0 y0 z0] + [dx dy dz] * t_min - end - end - return nothing -end - """ intersection_angles( diff --git a/src/physics/particles.jl b/src/physics/particles.jl index 1091d31e..2a917287 100644 --- a/src/physics/particles.jl +++ b/src/physics/particles.jl @@ -378,6 +378,20 @@ function toroidal_intersection(r1::Real, z1::Real, r2::Real, z2::Real, px::Real, return t end +""" + toroidal_intersections(wallr::Vector{T}, wallz::vector{T}, p::vector{T}, v:vector{T}) where {T<:Real} + +Returns the first time of intersection between a moving particle and a wall in toroidal geometry + + - `wallr`: Vector with r coordinates of the wall (must be closed) + - `wallz`: Vector with z coordinates of the wall (must be closed) + - `p`: Vector with current position (3D) of the particle in Cartesian coordinates. + - `v`: Vector with current velocity components of the particle +""" +function toroidal_intersection(wallr::Vector{T}, wallz::Vector{T}, p::Vector{T}, v::Vector{T}) where {T<:Real} + return toroidal_intersection(wallr, wallz, p[1], p[2], p[3], v[1], v[2], v[3]) +end + """ toroidal_intersections(wallr::Vector{T}, wallz::vector{T}, px::Real, py::Real, pz::Real, vx::Real, vy::Real, vz::Real) where {T<:Real} diff --git a/test/geometry/test_ray_torus_intersect.jl b/test/geometry/test_ray_torus_intersect.jl deleted file mode 100644 index be575cdf..00000000 --- a/test/geometry/test_ray_torus_intersect.jl +++ /dev/null @@ -1,16 +0,0 @@ -using IMAS -using Test - -@testset "Test ray_torus_intersect" begin - test_cases = [ - ([6.0,0.0,0.0], [-1.0,0.0,0.0], [1.0, 4.0], [-1.0, 1.0], [4.0 0.0 0.0], ) - ] - for (origin, direction, r_bounds, z_bounds, expected) in test_cases - @testset "origin=$origin, direction=$direction, r_bounds=$r_bounds, z_bounds=$z_bounds" begin - # println("--- x0=$x0, y0=$y0, dx=$dx, dy=$dy, r_min=$r_min, r_max=$r_max ---") - test_result = IMAS.ray_torus_intersect(origin, direction, r_bounds, z_bounds) - println("--- Result | expected: $test_result | $expected ---") - @test test_result ≈ expected - end - end -end \ No newline at end of file diff --git a/test/geometry/test_solve_r_intersect.jl b/test/geometry/test_solve_r_intersect.jl deleted file mode 100644 index 34b39fc6..00000000 --- a/test/geometry/test_solve_r_intersect.jl +++ /dev/null @@ -1,24 +0,0 @@ -using IMAS -using Test - -@testset "Test solve_r_intersect" begin - test_cases = [ - (3.0, 0.0, 0.0, 0.0, 1.0, 4.0, [[0.0 Inf];]), # ray is not moving, but starts on the grid - (6.0, 0.0, 0.0, 0.0, 1.0, 4.0, zeros(Float64, 0)), # ray is not moving, and starts off the grid - (6.0, 0.0, -1.0, 0.0, 1.0, 4.0, [2. 5.; 7. 10.]), # ray in x-direction - (0.0, 6.0, 0.0, -1.0, 2.0, 4.0, [2. 4.; 8. 10.]), # ray in y-direction - (3.0, 0.0, -1.0, 0.0, 1.0, 4.0, [0.0 2.; 4.0 7.0]), # ray in x-direction with start in torus - (6.0, 6.0, 0.0, -1.0, 2.0, 4.0, zeros(Float64, 0)), # ray in x-direction with offset causing no intersection - (3.0, 6.0, -3.0/5.0, -4.0/5.0, 1.0, 2.0, [[5.0 41.0/5.0];]), # Single crossing with x-y ray - (3.0, 8.0, -3.0/5.0, -4.0/5.0, 3.0, 4.0, [5.0 32.0/5.0; 10.0 57.0/5.0]) # double crossing with x-y ray - ] - - for (x0, y0, dx, dy, r_min, r_max, expected) in test_cases - @testset "x0=$x0, y0=$y0, dx=$dx, dy=$dy, r_min=$r_min, r_max=$r_max" begin - # println("--- x0=$x0, y0=$y0, dx=$dx, dy=$dy, r_min=$r_min, r_max=$r_max ---") - test_result = IMAS.solve_r_intersect(x0, y0, dx, dy, r_min, r_max) - println("--- Result | expected: $test_result | $expected ---") - @test test_result ≈ expected - end - end -end \ No newline at end of file diff --git a/test/physics/particles/test_toroidal_intersection.jl b/test/physics/particles/test_toroidal_intersection.jl new file mode 100644 index 00000000..80a1ddc9 --- /dev/null +++ b/test/physics/particles/test_toroidal_intersection.jl @@ -0,0 +1,16 @@ +using IMAS +using Test +@testset "Test toroical_intersection" begin + test_cases = [ + ([1.0, 4.0, 4.0, 1.0, 1.0], [-1.0, -1.0, 1.0, 1.0, -1.0], [6.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [4.0,0.0,0.0]) + ] + for (r_box, z_box, origin, direction, expected) in test_cases + @testset "origin=$origin, direction=$direction, r_box=$r_box, z_box=$z_box" begin + # println("--- x0=$x0, y0=$y0, dx=$dx, dy=$dy, r_min=$r_min, r_max=$r_max ---") + t = IMAS.toroidal_intersection(r_box, z_box, origin, direction) + test_result = origin .+ direction .* t + println("--- Result | expected: $test_result | $expected ---") + @test test_result ≈ expected + end + end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 67c260d6..617167f0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,8 +11,6 @@ else include("runtests_extract.jl") include("runtests_plot_recipes.jl") - - include("geometry/test_ray_torus_intersect.jl") - - include("geometry/test_solve_r_intersect.jl") + + include("physics/particles/test_toroidal_intersection.jl") end From ef97206c8c5865596d24563666debe7247c9ba2e Mon Sep 17 00:00:00 2001 From: Severin Denk Date: Tue, 17 Jun 2025 17:00:09 -0700 Subject: [PATCH 07/10] Fix accidental deletion when merging --- src/physics/sources.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/physics/sources.jl b/src/physics/sources.jl index e1bc1656..fd4d7579 100644 --- a/src/physics/sources.jl +++ b/src/physics/sources.jl @@ -673,6 +673,7 @@ function new_source( setproperty!(cs1d, :total_ion_power_inside, value; error_on_missing_coordinates = false) setproperty!(cs1d, :total_ion_energy, gradient(volume, value); error_on_missing_coordinates = false) else + value = electrons_power_inside setproperty!(cs1d, :total_ion_energy, zero(volume); error_on_missing_coordinates = false) setproperty!(cs1d, :total_ion_power_inside, zero(volume); error_on_missing_coordinates = false) end From 91a5b95b87c0410a98c859a613f7453daaedbf02 Mon Sep 17 00:00:00 2001 From: Severin Denk Date: Tue, 17 Jun 2025 17:02:02 -0700 Subject: [PATCH 08/10] Wrong spot --- src/physics/sources.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics/sources.jl b/src/physics/sources.jl index fd4d7579..b2f889c7 100644 --- a/src/physics/sources.jl +++ b/src/physics/sources.jl @@ -656,6 +656,7 @@ function new_source( setproperty!(electrons, :energy, value; error_on_missing_coordinates = false) setproperty!(electrons, :power_inside, cumtrapz(volume, value); error_on_missing_coordinates = false) elseif electrons_power_inside !== missing + value = electrons_power_inside setproperty!(electrons, :power_inside, value; error_on_missing_coordinates = false) setproperty!(electrons, :energy, gradient(volume, value); error_on_missing_coordinates = false) else @@ -673,7 +674,6 @@ function new_source( setproperty!(cs1d, :total_ion_power_inside, value; error_on_missing_coordinates = false) setproperty!(cs1d, :total_ion_energy, gradient(volume, value); error_on_missing_coordinates = false) else - value = electrons_power_inside setproperty!(cs1d, :total_ion_energy, zero(volume); error_on_missing_coordinates = false) setproperty!(cs1d, :total_ion_power_inside, zero(volume); error_on_missing_coordinates = false) end From 14f7d46cc0f893e079be27443a332e4d08474749 Mon Sep 17 00:00:00 2001 From: Severin Denk Date: Tue, 17 Jun 2025 17:03:13 -0700 Subject: [PATCH 09/10] Remove cruft --- src/geometry.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/geometry.jl b/src/geometry.jl index 862f88cf..2f7ac233 100644 --- a/src/geometry.jl +++ b/src/geometry.jl @@ -3,7 +3,6 @@ document[Symbol("Geometry")] = Symbol[] import StaticArrays import MillerExtendedHarmonic: MXH import LinearAlgebra -import Roots """ centroid(x::AbstractVector{<:T}, y::AbstractVector{<:T}) where {T<:Real} From 381cc894a077bd9ccafc459bb50d1e63debe0fa5 Mon Sep 17 00:00:00 2001 From: Severin Denk <60154343+AreWeDreaming@users.noreply.github.com> Date: Tue, 12 Aug 2025 10:59:07 -0700 Subject: [PATCH 10/10] Update src/physics/particles.jl Add assertion --- src/physics/particles.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/physics/particles.jl b/src/physics/particles.jl index 2a917287..93ab61f2 100644 --- a/src/physics/particles.jl +++ b/src/physics/particles.jl @@ -389,6 +389,7 @@ Returns the first time of intersection between a moving particle and a wall in t - `v`: Vector with current velocity components of the particle """ function toroidal_intersection(wallr::Vector{T}, wallz::Vector{T}, p::Vector{T}, v::Vector{T}) where {T<:Real} + @assert length(p) == length(v) == 3 return toroidal_intersection(wallr, wallz, p[1], p[2], p[3], v[1], v[2], v[3]) end