Add Liénard-Wiechert potential evaluation

.github/workflows/CI.yml

@@ -15,7 +15,7 @@ jobs:
   test:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
-    timeout-minutes: 60
+    timeout-minutes: 120
     permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created
       actions: write
       contents: read
@@ -37,6 +37,9 @@ jobs:
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
+      - name: Monitor memory
+        if: always()
+        run: free -h
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v5

Project.toml

@@ -4,6 +4,7 @@ version = "0.1.0-DEV"
 authors = ["Sebastian Micluța-Câmpeanu <sebastian.mc95@proton.me> and contributors"]
 
 [deps]
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
@@ -19,6 +20,8 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 [compat]
 Aqua = "0.8"
 CSV = "0.10"
+DataInterpolations = "8.9.0"
+FFTW = "1"
 HypergeometricFunctions = "0.3.28"
 JET = "0.9, 0.10, 0.11"
 LaserTypes = "0.2"
@@ -27,6 +30,7 @@ ModelingToolkit = "11"
 ModelingToolkitBase = "1.20.0"
 OrdinaryDiffEqNonlinearSolve = "1.10.0"
 OrdinaryDiffEqRosenbrock = "1.11.0"
+OrdinaryDiffEqTsit5 = "1.9"
 OrdinaryDiffEqVerner = "1.2.0"
 PhysicalConstants = "0.2.3"
 Plots = "1"
@@ -44,11 +48,13 @@ julia = "1.11"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 LaserTypes = "e07c0bfa-524c-4f35-a151-c3dd916fa2f0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -59,4 +65,4 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "CSV", "JET", "LaserTypes", "LinearAlgebra", "Test", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqVerner", "Plots", "Random", "SciMLBase", "StaticArrays", "Statistics", "SymbolicIndexingInterface"]
+test = ["Aqua", "CSV", "FFTW", "JET", "LaserTypes", "LinearAlgebra", "Test", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner", "Plots", "Random", "SciMLBase", "StaticArrays", "Statistics", "SymbolicIndexingInterface"]

scripts/Project.toml

@@ -1,10 +1,18 @@
 [deps]
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+DiffEqGPU = "071ae1c0-96b5-11e9-1965-c90190d839ea"
 ElectronDynamicsModels = "acecdaf2-97b2-47e1-90eb-3efa7bb274e5"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
+ProfileCanvas = "efd6af41-a80b-495e-886c-e51b0c7d77a3"
+ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"

scripts/thomson_scattering.jl

@@ -0,0 +1,149 @@
+using ElectronDynamicsModels
+using ModelingToolkit
+using OrdinaryDiffEqVerner, OrdinaryDiffEqTsit5
+using OrdinaryDiffEqNonlinearSolve
+using SciMLBase
+using StaticArrays
+using SymbolicIndexingInterface
+using LinearAlgebra
+using FFTW
+using CairoMakie
+
+# Atomic units
+const c = 137.03599908330932
+
+# Laser parameters
+ω = 0.057
+τ = 150 / ω
+λ = 2π * c / ω
+w₀ = 75λ
+Rmax = 3.25w₀
+
+a₀ = 10.0
+
+@named ref_frame = ProperFrame(:atomic)
+
+@named laser = LaguerreGaussLaser(;
+    wavelength = λ,
+    a0 = a₀,
+    beam_waist = w₀,
+    radial_index = 2,
+    azimuthal_index = -2,
+    ref_frame,
+    temporal_profile = :gaussian,
+    temporal_width = τ,
+    focus_position = 0.0,
+    polarization = :circular
+)
+@named elec = ClassicalElectron(; laser)
+sys = mtkcompile(elec)
+
+# Time span
+τi = -8τ
+τf = 8τ
+tspan = (τi, τf)
+
+# Single electron solve (for parameter access)
+x⁰ = [τi * c, 0.0, 0.0, 0.0]
+u⁰ = [c, 0.0, 0.0, 0.0]
+
+u0 = [
+    sys.x => x⁰,
+    sys.u => u⁰,
+]
+
+prob = ODEProblem{false, SciMLBase.FullSpecialize}(
+    sys, u0, tspan, u0_constructor = SVector{8}, fully_determined = true
+)
+sol0 = solve(prob, Vern9(), reltol = 1.0e-15, abstol = 1.0e-12)
+
+# Sunflower distribution for electron positions
+const ϕ = (1 + √5) / 2
+
+function radius(k, n, b)
+    if k > n - b
+        return 1.0
+    else
+        return sqrt(k - 0.5) / sqrt(n - (b + 1) / 2)
+    end
+end
+
+function sunflower(n, α)
+    points = Vector{Vector{Float64}}()
+    angle_stride = 2π / ϕ^2
+    b = round(Int, α * sqrt(n))
+    for k in 1:n
+        r = radius(k, n, b)
+        θ = k * angle_stride
+        push!(points, [r * cos(θ), r * sin(θ)])
+    end
+    return points
+end
+
+# Ensemble solve
+N = 300
+R₀ = Rmax * sunflower(N, 2)
+xμ = [[τi * c, r..., 0.0] for r in R₀]
+
+set_x = setsym_oop(prob, [Initial(sys.x); Initial(sys.u)])
+
+function prob_func(prob, i, repeat)
+    x_new = SVector{4}(xμ[i]...)
+    u_new = SVector{4}(c, 0.0, 0.0, 0.0)
+    u0, p = set_x(prob, SVector{8}(x_new..., u_new...))
+    return remake(prob; u0, p)
+end
+
+function abserr(a₀)
+    amp = log10(a₀)
+    expo = -amp^2 / 27 + 32amp / 27 - 220 / 27
+    return 10^expo
+end
+
+ensemble = EnsembleProblem(prob; prob_func, safetycopy = false)
+solution = solve(
+    ensemble, Vern9(), EnsembleThreads();
+    reltol = 1.0e-12, abstol = abserr(a₀), trajectories = N
+)
+
+# Radiation computation
+
+trajs = trajectory_interpolants(solution)
+
+# Screen parameters
+const Z = 2.0e5λ
+const δt = 2π / ω / 4
+const N_samples = floor(Int, (τf - τi) / δt)
+const x⁰_start = c * τi + hypot(Z, 25w₀ + Rmax)
+
+Nx = 50
+Ny = 50
+
+x⁰_samples = range(start = x⁰_start, step = c * δt, length = N_samples)
+
+screen = ObserverScreen(
+    LinRange(-25w₀, 25w₀, Nx),
+    LinRange(-25w₀, 25w₀, Ny),
+    Z,
+    x⁰_samples
+)
+
+A_s = accumulate_potential(trajs, screen, Tsit5())
+
+# GPU
+# accumulate_potential(trajs, screen, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()))
+A_ω = rfft(A_s, 1)
+
+# Find fundamental frequency bin
+freqs = rfftfreq(N_samples, 1 / δt)
+idx_f1 = findmin(x -> abs(x - ω / 2π), freqs)[2]
+
+# Plot the y-component of the potential at the fundamental frequency
+fig = Figure()
+ax = Axis(fig[1, 1], aspect = 1, xlabel = "x", ylabel = "y", title = "Thomson scattering (ω₁)")
+field = real.(A_ω[idx_f1, 3, :, :])
+heatmap!(
+    ax, collect(screen.x_grid), collect(screen.y_grid), field,
+    colorrange = maximum(abs, field) .* (-1, 1), colormap = :seismic
+)
+fig

src/ElectronDynamicsModels.jl

@@ -2,32 +2,32 @@ module ElectronDynamicsModels
 
 using ModelingToolkit
 using ModelingToolkitBase: AbstractSystem, SymbolicT, build_explicit_observed_function, get_systems
-using SymbolicIndexingInterface: getname, setsym_oop
+using SymbolicIndexingInterface: getname, setsym_oop, variable_index
 using PhysicalConstants, Unitful, UnitfulAtomic
 using PhysicalConstants.CODATA2018: c_0, e, m_e, ε_0
 using LinearAlgebra
 using Symbolics
 using HypergeometricFunctions: HypergeometricFunctions, _₁F₁, pochhammer
 using StaticArrays
+using SciMLBase
+using DataInterpolations
 
 m_dot(x, y) = x[1] * y[1] - x[2] * y[2] - x[3] * y[3] - x[4] * y[4]
 
 export GaussLaser, LaguerreGaussLaser
 export ReferenceFrame, ProperFrame, LabFrame,
-    UniformField,
-    PlaneWave,
+    UniformField, PlaneWave,
     ParticleDynamics,
-    LandauLifshitzRadiation,
-    AbrahamLorentzRadiation,
+    LandauLifshitzRadiation, AbrahamLorentzRadiation,
     ChargedParticle,
-    ClassicalElectron,
-    RadiatingElectron,
-    LandauLifshitzElectron,
-    FieldEvaluator
+    ClassicalElectron, RadiatingElectron, LandauLifshitzElectron,
+    FieldEvaluator,
+    ObserverScreen, trajectory_interpolants, TrajectoryInterpolant, accumulate_potential
 
 include("base.jl")
 include("dynamics.jl")
 include("fields.jl")
+include("radiation.jl")
 include("radiation_reaction.jl")
 include("external_fields.jl")
 include("systems.jl")

src/external_fields.jl

@@ -13,9 +13,11 @@ Parameters:
 - t₀: pulse center time (= n_cycles × 2π/ω)
 - z₀: focus position along z-axis
 """
-@component function GaussLaser(; name, wavelength = nothing, frequency = nothing,
+@component function GaussLaser(;
+        name, wavelength = nothing, frequency = nothing,
         a0 = 10.0, beam_waist = nothing, polarization = :linear,
-        n_cycles = 5, focus_position = nothing, ref_frame)
+        n_cycles = 5, focus_position = nothing, ref_frame
+    )
     if wavelength === nothing && frequency === nothing
         wavelength = 1.0
     end
@@ -82,8 +84,7 @@ Parameters:
         E[1] ~ real(ξx * E_g)
         E[2] ~ real(ξy * E_g)
         E[3] ~ real(
-            2im / (k * wz^2) *
-                (1 + im * (z / z_R)) *
+            2im / (k * wz^2) * (1 + im * (z / z_R)) *
                 (x[2] * (ξx * E_g) + x[3] * (ξy * E_g)),
         )
 
@@ -228,7 +229,7 @@ Reference: Sarachik & Schappert, Phys. Rev. D 1, 2738 (1970)
 
     vars = nameof(iv) == :τ ? [x, t] : [x]
 
-    sys = System(eqs, iv, vars, [A, ω, k_dir, pol, λ]; name, systems = [ref_frame], initialization_eqs, bindings, initial_conditions)
+    sys = System(eqs, iv; name, systems = [ref_frame], initialization_eqs, bindings, initial_conditions)
 
     extend(sys, field_dynamics)
 end