WIP: Major update to free-space simulations

Project.toml

@@ -46,6 +46,7 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 Scratch = "6c6a2e73-6563-6170-7368-637461726353"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -96,6 +97,7 @@ QuadGK = "2.4"
 Random = "1.9"
 Reexport = "1.2"
 Roots = "2"
+Rotations = "1.7.1"
 Scratch = "1"
 SpecialFunctions = "0.10.3, 1, 2"
 StaticArrays = "1"

examples/low_level_interface/freespace/chi2_benchmarking.jl

@@ -0,0 +1,37 @@
+using Luna
+using BenchmarkTools
+import Rotations: RotZY, RotYZ, RotMatrix, RotMatrix3
+import LinearAlgebra: mul!
+using Metal
+
+θ = deg2rad(23.3717)
+ϕ = deg2rad(30)
+material = :BBO
+
+Et = rand(2^12, 2)
+out = zero(Et)
+
+c = Nonlinear.Chi2Field(θ, ϕ, PhysData.χ2(material))
+# r = Nonlinear.Kerr_field(PhysData.χ3(material))
+
+function toprofile(n)
+    for i in 1:n
+        c(out, Et, 1)
+    end
+end
+##
+@profview toprofile(1)
+@profview toprofile(100)
+@btime c(out, Et, 1)
+
+##
+Enl2 = rand(2^12, 6)'
+Et2 = rand(2^12, 3)'
+##
+Et2_m = MtlArray(convert(Matrix{Float32}, Et2))
+χ2_m = MtlArray(convert(Matrix{Float32}, c.χ2_toLab))
+Enl2_m = MtlArray(convert(Matrix{Float32}, Enl2))
+@btime mul!(Et2_m, χ2_m, Enl2_m)
+
+##
+@btime mul!(Et2, c.χ2_toLab, Enl2)

examples/low_level_interface/freespace/free2D.jl

@@ -0,0 +1,85 @@
+using Luna
+Luna.set_fftw_mode(:estimate)
+import FFTW
+import PyPlot:pygui, plt
+pygui(true)
+
+gas = :Ar
+pres = 4
+
+τ = 30e-15
+λ0 = 800e-9
+
+w0 = 2e-3
+energy = 1.5e-3
+L = 2
+
+R = 6e-3
+N = 128
+
+grid = Grid.RealGrid(L, 800e-9, (400e-9, 2000e-9), 0.2e-12)
+xgrid = Grid.Free2DGrid(R, N)
+
+x = xgrid.x
+energyfun, energyfunω = Fields.energyfuncs(grid, xgrid)
+
+densityfun = let dens0=PhysData.density(gas, pres)
+    z -> dens0
+end
+
+responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),)
+
+linop = LinearOps.make_const_linop(grid, xgrid, PhysData.ref_index_fun(gas, pres))
+normfun = NonlinearRHS.const_norm_free2D(grid, xgrid, PhysData.ref_index_fun(gas, pres))
+
+inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0)
+
+Eω, transform, FT = Luna.setup(grid, xgrid, densityfun, normfun, responses, inputs)
+
+# statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid))
+output = Output.MemoryOutput(0, grid.zmax, 21)
+
+Luna.run(Eω, grid, linop, transform, FT, output, max_dz=Inf, init_dz=1e-1)
+
+ω = grid.ω
+t = grid.t
+
+zout = output.data["z"]
+Eout = output.data["Eω"]
+
+println("Transforming...")
+Eωx = FFTW.ifft(Eout, 3)
+Etx = FFTW.irfft(Eout, length(grid.t), (1, 3))
+println("...done")
+
+
+Ilog = log10.(Maths.normbymax(abs2.(Eωx)))
+
+Iωx = abs2.(Eωx);
+
+Ix = zeros(Float64, (length(x), length(zout)));
+energy = zeros(length(zout));
+for ii in axes(Etx, 4)
+    energy[ii] = energyfun(Etx[:, 1, :, ii]);
+    Ix[:, ii] = (grid.ω[2]-grid.ω[1]) .* sum(Iωx[:, 1, :, ii], dims=1);
+end
+
+ω0idx = argmin(abs.(grid.ω .- 2π*PhysData.c/λ0))
+
+E0ωyx = FFTW.ifft(Eω[ω0idx, 1, :, :], (1, 2))
+
+Iωxlog = log10.(Maths.normbymax(Iωx))
+
+plt.figure()
+plt.pcolormesh(zout, ω.*1e-15/2π, Iωxlog[:, 1, N÷2+1, :])
+plt.xlabel("Z (m)")
+plt.ylabel("f (PHz)")
+plt.title("I(ω, x=0, z)")
+plt.clim(-6, 0)
+plt.colorbar()
+
+
+plt.figure()
+plt.plot(zout.*1e2, energy.*1e6)
+plt.xlabel("Distance [cm]")
+plt.ylabel("Energy [μJ]")

examples/low_level_interface/freespace/free2D_bbo.jl

@@ -0,0 +1,192 @@
+#=
+Here we simulate type I SHG in a BBO crystal driven at 800 nm. The focus is small and the crystal
+is reasonably thick, so we observe strong temporal and spatial walk-off effects.
+=#
+using Luna
+import FFTW
+import Luna: Hankel
+import PyPlot: plt
+import NumericalIntegration: integrate
+
+λ0 = 800e-9 # driving wavelength
+τfwhm = 30e-15 # pulse duration
+w0 = 20e-6 # 1/e² beam radius
+energy = 10e-9 # pulse energy
+
+material = :BBO
+thickness = 200e-6 # BBO thickness
+
+R = 4*w0 # radius of the spatial window
+N = 2^7 # number of spatial points
+
+grid = Grid.RealGrid(thickness, λ0, (250e-9, 2e-6), 120e-15)
+xgrid = Grid.Free2DGrid(R, N)
+
+θ = deg2rad(29.2) # type I phase-matching angle
+ϕ = deg2rad(30) # ϕ for type I phase-matching
+
+##
+responses = (Nonlinear.Chi2Field(θ, ϕ, PhysData.χ2(material)),
+             Nonlinear.Kerr_field(PhysData.χ3(material)))
+
+#=
+Ref index functions. Note that here nfunx takes (λ, δθ=0) as arguments,
+which allows make_const_linop to calculate the actual internal angle depending
+on frequency and transverse k-vector component.
+=#
+nfunx, nfuny = PhysData.ref_index_fun_xy(material, θ)
+linop = LinearOps.make_const_linop(grid, xgrid, (nfunx, nfuny))
+
+normfun = NonlinearRHS.const_norm_free2D(grid, xgrid, (nfunx, nfuny))
+densityfun = z -> 1 # density is unity because we're considering a solid.
+##
+# scaling factor in front of the energy corresponds to the integral over y
+inputs = Fields.GaussGaussField(;λ0, τfwhm, energy=energy/(sqrt(π/2)*w0), w0)
+Eω, transform, FT = Luna.setup(grid, xgrid, densityfun, normfun, responses, inputs)
+
+##
+output = Output.MemoryOutput(0, grid.zmax, 101)
+Luna.run(Eω, grid, linop, transform, FT, output; init_dz=1e-6)
+
+##
+z = output["z"]
+Eωk = output["Eω"] # (ω, pol, k, z)
+x = xgrid.x
+
+# normalisation prefactor for spectral intensity
+ωprefac = PhysData.c*PhysData.ε_0/2 * 2π/(grid.ω[end]^2) * sqrt(π/2)*w0
+
+Eωr = FFTW.ifft(Eωk, 3) # (ω, pol, x, z)
+Etr = FFTW.irfft(Eωr, 2*(length(grid.ω)-1), 1) # (t, pol, x, z)
+EtrH = Maths.hilbert(Etr)
+Iωr = abs2.(Eωr) # (ω, pol, x, z)
+Itr = 0.5*PhysData.c*PhysData.ε_0*abs2.(EtrH) # (t, pol, x, z)
+
+Irxy = dropdims(sum(Iωr; dims=1); dims=1) # (pol, x, z)
+Iωxy = dropdims(Maths.integrateNd(x, Iωr; dim=3); dims=3)*ωprefac # (ω, pol, z)
+Ir = dropdims(sum(Iωr; dims=(1, 2)); dims=(1, 2)) # (x, z)
+
+Itxy = dropdims(Maths.integrateNd(x, Itr; dim=3); dims=3)*sqrt(π/2)*w0 # (t, pol, z)
+
+Eω0 = Eωr[:, :, length(x)÷2+1, :]
+Et0 = FFTW.irfft(Eω0, 2*(length(grid.ω)-1), 1) # (t, pol, z)
+Et0 = Maths.hilbert(Et0)
+It0 = 0.5*PhysData.c*PhysData.ε_0*abs2.(Et0)
+
+et, eω = Fields.energyfuncs(grid, xgrid)
+
+energy_out = dropdims(mapslices(eω, Eωk; dims=(1, 3)); dims=(1, 3))*sqrt(π/2)*w0
+
+ω = grid.ω
+
+
+##
+plt.figure()
+plt.pcolormesh(z*1e3, x*1e6, Ir)
+plt.xlabel("Distance (mm)")
+plt.ylabel("X (μm)")
+
+##
+plt.figure()
+plt.plot(z*1e6, 1e9energy_out[1, :]; label="X polarisation")
+plt.plot(z*1e6, 1e9energy_out[2, :]; label="Y polarisation")
+plt.xlabel("Distance (μm)")
+plt.ylabel("Energy (nJ)")
+plt.legend()
+
+##
+fig = plt.figure()
+fig.set_size_inches(12, 6)
+plt.subplot(1, 2, 1)
+plt.pcolormesh(z*1e3, x*1e6, Irxy[1, :, :])
+plt.xlabel("Distance (mm)")
+plt.ylabel("radius (μm)")
+plt.title("X Polarisation")
+plt.subplot(1, 2, 2)
+plt.pcolormesh(z*1e3, x*1e6, Irxy[2, :, :])
+plt.xlabel("Distance (mm)")
+plt.ylabel("radius (μm)")
+plt.title("Y Polarisation")
+fig.tight_layout()
+
+
+##
+plt.figure()
+plt.subplot(1, 2, 1)
+plt.pcolormesh(z*1e3, grid.t*1e15, Itxy[:, 1, :])
+plt.title("X Pol")
+plt.xlabel("Distance (mm)")
+plt.ylabel("Time (fs)")
+plt.subplot(1, 2, 2)
+plt.pcolormesh(z*1e3, grid.t*1e15, Itxy[:, 2, :])
+plt.title("Y Pol")
+plt.xlabel("Distance (mm)")
+plt.suptitle("Time domain")
+
+##
+plt.figure()
+plt.subplot(2, 1, 1)
+plt.plot(grid.t*1e15, Itxy[:, 1, 1]; label="input X")
+plt.plot(grid.t*1e15, Itxy[:, 2, 1]; label="input Y")
+plt.legend()
+plt.subplot(2, 1, 2)
+plt.plot(grid.t*1e15, Itxy[:, 1, end]; label="output X")
+plt.plot(grid.t*1e15, Itxy[:, 2, end]; label="output Y")
+plt.xlabel("Time (fs)")
+plt.legend()
+
+##
+mm = 2π*maximum(Iωxy)*1e12
+plt.figure()
+plt.subplot(2, 1, 1)
+plt.semilogy(ω*1e-12/2π, 2π*1e12*Iωxy[:, 1, 1]; label="input X")
+plt.semilogy(ω*1e-12/2π, 2π*1e12*Iωxy[:, 2, 1]; label="input Y")
+plt.ylim(1e-5mm, 5mm)
+plt.xlim(200, 900)
+plt.legend()
+plt.subplot(2, 1, 2)
+plt.semilogy(ω*1e-12/2π, 2π*1e12*Iωxy[:, 1, end]; label="output X")
+plt.semilogy(ω*1e-12/2π, 2π*1e12*Iωxy[:, 2, end]; label="output Y")
+# plt.semilogy(ω*1e-12/2π, 2π*1e12*Iωxy[:, 2, 1]; c="C1", linestyle="--", label="input Y")
+plt.ylim(1e-5mm, 5mm)
+plt.xlim(200, 900)
+plt.xlabel("Frequency (THz)")
+plt.ylabel("SED (J/Hz)")
+plt.legend()
+
+##
+fig = plt.figure()
+fig.set_size_inches(12, 7)
+plt.subplot(1, 2, 1)
+plt.pcolormesh(grid.t*1e15, xgrid.x*1e6, (Etr[:, 1, :, end])'; cmap="seismic")
+plt.clim([-1, 1].*maximum(abs.(Etr[:, 1, :, end])))
+plt.xlim(-50, 50)
+plt.ylim(-100, 100)
+plt.xlabel("Time (fs)")
+plt.ylabel("X (μm)")
+plt.title("X polarisation")
+
+plt.subplot(1, 2, 2)
+plt.pcolormesh(grid.t*1e15, xgrid.x*1e6, (Etr[:, 2, :, end])'; cmap="seismic")
+plt.clim([-1, 1].*maximum(abs.(Etr[:, 2, :, end])))
+plt.xlim(-50, 50)
+plt.ylim(-100, 100)
+plt.xlabel("Time (fs)")
+plt.title("Y polarisation")
+
+##
+fig = plt.figure()
+fig.set_size_inches(12, 6)
+plt.subplot(1, 2, 1)
+plt.plot(grid.t*1e15, 1e-9*(Etr[:, 1, length(xgrid.x)÷2+1, end]);
+        label="Max: $(maximum(1e-9*(Etr[:, 1, length(xgrid.x)÷2+1, end])))")
+plt.ylabel("E (GV/m)")
+plt.legend()
+plt.xlabel("Time (fs)")
+plt.title("X Polarisation")
+plt.subplot(1, 2, 2)
+plt.plot(grid.t*1e15, 1e-9*(Etr[:, 2, length(xgrid.x)÷2+1, end]);
+          label="Max: $(maximum(1e-9*(Etr[:, 2, length(xgrid.x)÷2+1, end])))")
+plt.legend()
+plt.xlabel("Time (fs)")
+plt.title("Y Polarisation")