diff --git a/Project.toml b/Project.toml index a60a54f..a068e89 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.2.0" [deps] JSServe = "824d6782-a2ef-11e9-3a09-e5662e0c26f9" LaserTypes = "e07c0bfa-524c-4f35-a151-c3dd916fa2f0" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" PICAnalysisTools = "049d77de-c438-40ac-b6d1-f1d179569bd2" PICDataStructures = "cb18667b-f199-4b8b-b886-9916fef47a26" @@ -25,7 +26,7 @@ WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008" [compat] JSServe = "1.2" LaserTypes = "0.1" -PICDataStructures = "0.3" +PICDataStructures = "0.4" PhysicalConstants = "0.2" Plots = "1.10" ProgressLogging = "0.1" diff --git a/src/SDFResultViewer.jl b/src/SDFResultViewer.jl index e3634b9..fb74952 100644 --- a/src/SDFResultViewer.jl +++ b/src/SDFResultViewer.jl @@ -1,12 +1,14 @@ module SDFResultViewer -export add_λ_units, unit_l, unit_t, unit_E, unit_B, unit_L, pₑ, +export add_λ_units, unit_l, unit_t, unit_E, unit_B, unit_L, unit_S, pₑ, batch_plot, mean_Lx_plot, Lx_section_plot, phase_space_summary_plot, average_linear_density_x, npart_plot, - section_widget, - compare_E_slice_with_analytic, similar_E, similar_laser + section_widget, approx_plot, + compare_E_slice_with_analytic, compare_S_slice_with_analytic, + poynting_summary, + apply_analytic, analytic_laser using SDFResults using PICDataStructures @@ -18,8 +20,9 @@ using Rotations using RecursiveArrayTools: recursive_bottom_eltype using Unitful using Unitful: superscript -using PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e +using PhysicalConstants.CODATA2018: c_0, ε_0, μ_0, m_e, e using Statistics, StatsBase +using LinearAlgebra using ProgressLogging # plots import Plots @@ -30,6 +33,8 @@ using JSServe using JSServe: DOM, Slider using Markdown +Unitful.register(@__MODULE__) + include("units.jl") include("batch.jl") include("plots/angular_momentum.jl") diff --git a/src/plots/fields.jl b/src/plots/fields.jl index 95ccb7c..48af174 100644 --- a/src/plots/fields.jl +++ b/src/plots/fields.jl @@ -1,33 +1,219 @@ -function E_slice_plot(E::ScalarField, dir, slice_location, t; - issliced=false, - save=false, - slice_dir=:x, +function approx_plot(f; + fig = Figure(), + figurepos = fig[1,1], + xlabel = "", + ylabel = "", + title = "", + downsample_size = nothing, + cbar_label="", kwargs...) - if !issliced - E_slice = slice(E, slice_dir, slice_location) + ax = Axis(figurepos[1,1]; + xlabel, ylabel, + title, + aspect = DataAspect() + ) + + if isnothing(downsample_size) + f_approx = downsample_approx(f) else - E_slice = E + f_approx = downsample_approx(f, downsample_size) end + @debug "Downsampled to $(size(f_approx))" + + plt = fieldplot!(ax, f_approx; kwargs...) + + Colorbar(figurepos[1,2], plt; width=20, label=cbar_label) + + # TODO: Implement save + # if save + # save && Plots.savefig(plt, "E$(dir_slice)_$loc.png") + # end + + return fig +end + +get_unit_label(f) = " (" * unitname(f) * ")" + +function compare_slice_with_analytic(numeric_f, f, laser, profile, file; + fig = Figure(resolution = (800, 700)), + component=1, + slice_location=0.0unit_l, + slice_dir=:x, + fieldname="", + f_unit) - grid = getdomain(E_slice) - cl = ustrip(max(abs.(extrema(E_slice))...)) + fx = get_parameter(file, :constant, :f_x)*u"m" |> unit_l + t = get_time(file) |> unit_t + s = analytic_laser(laser, profile, file) + + grid = getdomain(numeric_f) + + analytic_f = apply_analytic(f, grid, t, s; + component, + propagation_dir=:x, + x₀=fx) |> f_unit + + analytic_slice, loc = dynamic_slice(analytic_f, slice_dir, slice_location) + numeric_slice, loc = dynamic_slice(numeric_f, slice_dir, slice_location) + + units = get_unit_label(grid) + xlabel, ylabel = filter(x->x≠string(slice_dir), + ["x", "y", "z"] .* units) ts = sprint(show, t, context=:compact => true) - loc = sprint(show, slice_location, context=:compact => true) - xlabel, ylabel = filter(x->x≠string(slice_dir), ["x", "y", "z"]) - - plt = Plots.heatmap(grid..., E_slice; - title = "E$dir at t=$ts, x=$loc", - xlabel, ylabel, colorbar_title = "E$dir", - seriescolor = :jet1, - aspect_ratio = 1, - framestyle = :box, - clims = (-cl, cl), - kwargs...) - - if save - save && Plots.savefig(plt, "E$(dir_slice)_$loc.png") - end + loc_str = @lift sprint(show, + $loc, context=:compact => true) * " " * unitname(slice_location) + title = @lift fieldname * " slice at t = " * ts * " and " * + string(slice_dir) * " = " * $loc_str + + ax1 = Axis(fig[1,1][1,1]; + xlabel, ylabel, + title = "Analytic", + aspect = DataAspect() + ) + ax2 = Axis(fig[1,2][1,1]; + xlabel, ylabel, + title = "Numeric", + aspect = DataAspect() + ) + + plt1 = fieldplot!(ax1, analytic_slice) + plt2 = fieldplot!(ax2, numeric_slice) + + Colorbar(fig[1,1][2,1], plt1; + width = Relative(3/4), + height = 20, + vertical = false, flipaxis = false, + tellheight = true, + label = fieldname * get_unit_label(numeric_f), + ) + Colorbar(fig[1,2][2,1], plt2; + width = Relative(3/4), + height = 20, + vertical = false, flipaxis = false, + tellheight = true, + label = fieldname * get_unit_label(numeric_f), + ) + + Label(fig[0, :], title, textsize = 30, color = :black) + + return fig +end + +function compare_E_slice_with_analytic(laser, profile, file, dir; + slice_location=0.0unit_l, + slice_dir=:x) + + E = uconvert(unit_E, unit_l, file["e$dir"]) + + compare_slice_with_analytic(E, LaserTypes.E, laser, profile, file; + component=dir_to_idx(dir), + slice_location, + slice_dir, + f_unit=unit_E, + fieldname="E$dir" + ) +end + +function compare_S_slice_with_analytic(laser, profile, file, dir; + slice_location=0.0unit_l, + slice_dir=:x) + + Ex, Ey, Ez = uconvert.(unit_E, unit_l, file[:ex, :ey, :ez]) + Bx, By, Bz = uconvert.(unit_B, unit_l, file[:bx, :by, :bz]) + + E = build_vector((Ex, Ey, Ez), (:x, :y, :z)) + B = build_vector((Bx, By, Bz), (:x, :y, :z)) + S = ustrip(E × B) / ustrip(μ_0) + + compare_slice_with_analytic(S, LaserTypes.S, laser, profile, file; + component=dir_to_idx(dir), + slice_location=ustrip(slice_location), + slice_dir, + f_unit=unit_E, + ) +end + +function poynting_summary(file, laser, profile, slice_location) + t = get_time(file) |> unit_t + + # get data + fx = get_parameter(file, :constant, :f_x)*u"m" |> unit_l + Ex, Ey, Ez, E, Bx, By, Bz, B = read_EM_fields(file) + + Ex_grid = getdomain(Ex) + Bx_grid = getdomain(Bx) + + s = analytic_laser(laser, profile, file) + + analytic_E = get_analytic_E(Ex_grid, t, s, x₀=fx) + analytic_B = get_analytic_B(Bx_grid, t, s, x₀=fx) + analytic_S = uconvert(unit_S, analytic_E × analytic_B / μ_0) + + S = uconvert(unit_S, E × B / μ_0) + + S_approx = downsample_approx(S, 15) + analytic_S_approx = downsample_approx(analytic_S, 15) + S_slice = downsample_approx( + slice(S, :x, slice_location), + 50 + ) + analytic_S_slice = downsample_approx( + slice(analytic_S, :x, slice_location), + 50 + ) + + fig = Figure(resolution=(800,800)) + # Layout + lscene1 = LScene(fig[1, 1], scenekw = (camera = cam3d!, raw = false)) + lscene2 = LScene(fig[1, 2], scenekw = (camera = cam3d!, raw = false)) + ax1 = Axis(fig[2,1][1,1], + aspect = DataAspect(), + title="Numeric", + xlabel="y ($unit_l)", + ylabel="z ($unit_l)" + ) + ax2 = Axis(fig[2,2][1,1], + aspect = DataAspect(), + title="Analytic", + xlabel="y ($unit_l)", + ylabel="z ($unit_l)" + ) + # Plots + # 3D + arrowsize_factor = 0.02 + lengthscale_factor = 10 + fieldplot!(lscene1, S_approx; + arrowsize_factor, + lengthscale_factor + ) + fieldplot!(lscene2, analytic_S_approx; + arrowsize_factor, + lengthscale_factor + ) + # 2D + arrowsize_factor = 0.15 + lengthscale_factor = 6 + plt1 = fieldplot!(ax1, S_slice; + arrowsize_factor, + lengthscale_factor + ) + plt2 = fieldplot!(ax2, analytic_S_slice; + arrowsize_factor, + lengthscale_factor + ) + # Colorbars + Colorbar(fig[2,1][2,1], plt1, + width = Relative(3/4), height = 15, + vertical = false, flipaxis = false, + label = "S ($unit_S)", + ) + Colorbar(fig[2,2][2,1], plt2, + width = Relative(3/4), height = 15, + vertical = false, flipaxis = false, + label = "S ($unit_S)", + ) - return plt + WGLMakie.trim!(fig.layout) + fig end diff --git a/src/plots/utils.jl b/src/plots/utils.jl new file mode 100644 index 0000000..2814f07 --- /dev/null +++ b/src/plots/utils.jl @@ -0,0 +1,8 @@ +struct MinMaxTicks end + +function MakieLayout.get_tickvalues(::MinMaxTicks, min, max) + ticks = MakieLayout.get_tickvalues(WilkinsonTicks(5, k_min = 3), min, max) + pushfirst!(ticks, min) + push!(ticks, max) + ticks +end diff --git a/src/units.jl b/src/units.jl index 0999897..732e904 100644 --- a/src/units.jl +++ b/src/units.jl @@ -2,6 +2,7 @@ function add_λ_units(λ) expr = quote using Unitful + using Unitful: MW using PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e local λ = $λ @@ -12,6 +13,7 @@ function add_λ_units(λ) @unit unit_B "mₑ ω/e" UnitMagneticField m_e * ω / e false @unit unit_E "mₑ c ω/e" UnitElectricField m_e * ω / e * c_0 false @unit unit_L "mₑ c²/ω" UnitAngularMomenta m_e * c_0^2 / ω false + @unit unit_S "MW ω²/c²" UnitPoyntingFlux MW * ω^2 / c_0^2 false @unit pₑ "mₑ c" ElectronMomenta m_e * c_0 false Unitful.register(@__MODULE__) diff --git a/src/utils.jl b/src/utils.jl index ddcd144..e0e9855 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -7,7 +7,36 @@ function sdfcount(root_folder, pattern) Plots.plot(basename.(folders[idx]), n[idx], rotation=45, legend=false) end -function similar_laser(laser::Type{LaguerreGaussLaser}, t_profile, file) +function read_EM_fields(file) + Ex, Ey, Ez = uconvert.(unit_E, unit_l, file[:ex, :ey, :ez]) + Bx, By, Bz = uconvert.(unit_B, unit_l, file[:bx, :by, :bz]) + + E = build_vector((Ex, Ey, Ez), (:x, :y, :z)) + B = build_vector((Bx, By, Bz), (:x, :y, :z)) + + E, B +end + +function dynamic_slice(f, slice_dir, slice_location) + if f isa Observable + obs_f = ustrip(f) + else + @debug "Converted field to Observable" + obs_f = Observable(ustrip(f)) + end + if slice_location isa Observable + loc = @lift ustrip($slice_location) + else + @debug "Converted slice_location to Observable" + loc = Observable(ustrip(slice_location)) + end + + f_slice = @lift slice($obs_f, slice_dir, $loc) + + return f_slice, loc +end + +function analytic_laser(laser::Type{LaguerreGaussLaser}, t_profile, file) fx = get_parameter(file, :constant, :f_x)*u"m" |> unit_l λ = get_parameter(file, :laser, :lambda) m = get_parameter(file, :constant, :m) @@ -42,62 +71,94 @@ function field_rotation(dir) end end -function similar_E(::Val{true}, laser, t, propagation_dir, component::Int, grid; x₀, y₀, z₀) - map(Iterators.product(grid...)) do (x,y,z) +function apply_analytic(::Val{true}, f, laser, t, propagation_dir, component, grid; x₀, y₀, z₀) + mapgrid(grid) do (x,y,z) r = SVector{3}(x-x₀, y-y₀, z-z₀) R = field_rotation(propagation_dir) - analytic_E = R * (LaserTypes.E(R \ r, laser) * LaserTypes.g((R \ r)[3], t, laser)) - analytic_E[component] + analytic_E = R * (f(R \ r, laser) * LaserTypes.g((R \ r)[3], t, laser)) + + if component isa Number + analytic_E[component] + else + analytic_E + end end end -function similar_E(::Val{false}, laser, t, propagation_dir, component::Int, grid; x₀, y₀, z₀) - map(Iterators.product(grid...)) do (x,y,z) +function apply_analytic(::Val{false}, f, laser, t, propagation_dir, component, grid; x₀, y₀, z₀) + mapgrid(grid) do (x,y,z) r = SVector{3}(x-x₀, y-y₀, z-z₀) R = field_rotation(propagation_dir) - analytic_E = R * LaserTypes.E(R \ r, t, laser) - analytic_E[component] + analytic_E = R * f(R \ r, t, laser) + + if component isa Number + analytic_E[component] + else + analytic_E + end end end -function similar_E(E, t, laser; +function apply_analytic(f, grid, t, laser; component=1, propagation_dir=:x, complex=false, - x₀=zero(recursive_bottom_eltype(getdomain(E))), - y₀=zero(recursive_bottom_eltype(getdomain(E))), - z₀=zero(recursive_bottom_eltype(getdomain(E)))) + downsample_grid=true, + x₀=zero(recursive_bottom_eltype(grid)), + y₀=zero(recursive_bottom_eltype(grid)), + z₀=zero(recursive_bottom_eltype(grid))) + + grid = downsample_grid ? downsample_approx(grid) : grid + data = apply_analytic(Val(complex), + f, + laser, + t, + propagation_dir, + component, + grid; + x₀,y₀,z₀) + if eltype(data) <: Number + ScalarField(data, grid) + else + # broken + @error "Use analytic_E or analytic_B" + VectorField(data, grid) + end +end - E_approx = approximate_field(E) +# Stopgap solution untill we can greate vector fields correctly with the above +function get_analytic_E(grid, t, s; + downsample_grid=true, + x₀=zero(recursive_bottom_eltype(grid)), + y₀=zero(recursive_bottom_eltype(grid)), + z₀=zero(recursive_bottom_eltype(grid))) - grid = AxisGrid((getdomain(E_approx).*unit_l...,)) - data = similar_E(Val(complex), laser, t, propagation_dir, component, grid; x₀,y₀,z₀) + grid = downsample_grid ? downsample_approx(grid) : grid - ScalarField(data, grid) + analytic_Ex = apply_analytic(LaserTypes.E, grid, t, s; + x₀,y₀,z₀, component=1) |> unit_E + analytic_Ey = apply_analytic(LaserTypes.E, grid, t, s; + x₀,y₀,z₀, component=2) |> unit_E + analytic_Ez = apply_analytic(LaserTypes.E, grid, t, s; + x₀,y₀,z₀, component=3) |> unit_E + + build_vector((analytic_Ex, analytic_Ey, analytic_Ez), (:x,:y,:z)) end -function compare_E_slice_with_analytic(laser, profile, file, dir; - slice_location=0unit_l, - slice_dir=:x) +function get_analytic_B(grid, t, s; + downsample_grid=true, + x₀=zero(recursive_bottom_eltype(grid)), + y₀=zero(recursive_bottom_eltype(grid)), + z₀=zero(recursive_bottom_eltype(grid))) - fx = get_parameter(file, :constant, :f_x)*u"m" |> unit_l - t = get_time(file) |> unit_t - s = similar_laser(laser, profile, file) - - E = uconvert(unit_E, unit_l, file["e$dir"]) - - analytic_E = similar_E(E, t, s, - component=dir_to_idx(dir), - propagation_dir=:x, - x₀=fx) |> unit_E - - plt1 = E_slice_plot(analytic_E, dir, slice_location, t; - legend_title = "Analytic", - slice_dir, - issliced=false) - plt2 = E_slice_plot(E, dir, slice_location, t; - legend_title = "Numeric", - slice_dir, - issliced=false) - Plots.plot(plt1, plt2, size=(400,600), layout=(2,1)) + grid = downsample_grid ? downsample_approx(grid) : grid + + analytic_Bx = apply_analytic(LaserTypes.B, grid, t, s; + x₀,y₀,z₀, component=1) |> unit_B + analytic_By = apply_analytic(LaserTypes.B, grid, t, s; + x₀,y₀,z₀, component=2) |> unit_B + analytic_Bz = apply_analytic(LaserTypes.B, grid, t, s; + x₀,y₀,z₀, component=3) |> unit_B + + build_vector((analytic_Bx, analytic_By, analytic_Bz), (:x,:y,:z)) end