Fix full modal nonlinear polarisation

Project.toml

@@ -1,7 +1,7 @@
 name = "Luna"
 uuid = "30eb0fb0-5147-11e9-3356-d75b018717ce"
 authors = ["chrisbrahms <[email protected]>"]
-version = "0.1.0"
+version = "0.1.2"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"

src/Capillary.jl

@@ -261,15 +261,16 @@ dimlimits(m::MarcatiliMode; z=0) = (:polar, (0.0, 0.0), (radius(m, z), 2π))
 
 # we use polar coords, so xs = (r, θ)
 function field(m::MarcatiliMode, xs; z=0)
+    r, θ = xs
     if m.kind == :HE
-        return (besselj(m.n-1, xs[1]*m.unm/radius(m, z)) .* SVector(
-            cos(xs[2])*sin(m.n*(xs[2] + m.ϕ)) - sin(xs[2])*cos(m.n*(xs[2] + m.ϕ)),
-            sin(xs[2])*sin(m.n*(xs[2] + m.ϕ)) + cos(xs[2])*cos(m.n*(xs[2] + m.ϕ))
+        return (besselj(m.n-1, r*m.unm/radius(m, z)) .* SVector(
+            cos(θ)*sin(m.n*(θ + m.ϕ)) - sin(θ)*cos(m.n*(θ + m.ϕ)),
+            sin(θ)*sin(m.n*(θ + m.ϕ)) + cos(θ)*cos(m.n*(θ + m.ϕ))
             ))
     elseif m.kind == :TE
-        return besselj(1, xs[1]*m.unm/radius(m, z)) .* SVector(-sin(xs[2]), cos(xs[2]))
+        return besselj(1, r*m.unm/radius(m, z)) .* SVector(-sin(θ), cos(θ))
     elseif m.kind == :TM
-        return besselj(1, xs[1]*m.unm/radius(m, z)) .* SVector(cos(xs[2]), sin(xs[2]))
+        return besselj(1, r*m.unm/radius(m, z)) .* SVector(cos(θ), sin(θ))
     end
 end
 

src/NonlinearRHS.jl

@@ -239,7 +239,7 @@ function pointcalc!(fval, xs, t::TransModal)
             fval[:, i] .= 0.0
             continue
         end
-        if size(xs, 1) > 1
+        if size(xs, 1) > 1 # full 2-D mode integral
             x2 = xs[2, i]
             if t.dimlimits[1] == :polar
                 pre = x1
@@ -264,6 +264,7 @@ function pointcalc!(fval, xs, t::TransModal)
         to_time!(t.Er, t.Erω, t.Erωo, inv(t.FT))
         # get nonlinear pol at r,θ
         Et_to_Pt!(t.Pr, t.Er, t.resp, t.density)
+        t.ncalls += 1
         @. t.Pr *= t.grid.towin
         to_freq!(t.Prω, t.Prωo, t.Pr, t.FT)
         @. t.Prω *= t.grid.ωwin
@@ -277,17 +278,18 @@ end
 
 function (t::TransModal)(nl, Eω, z)
     reset!(t, Eω, z)
+    _, ll, ul = t.dimlimits
     if t.full
-        val, err = Cubature.pcubature_v(
+        val, err = Cubature.hcubature_v(
             length(Eω)*2,
             (x, fval) -> pointcalc!(fval, x, t),
-            t.dimlimits[2], t.dimlimits[3], 
+            ll, ul, 
             reltol=t.rtol, abstol=t.atol, maxevals=t.mfcn, error_norm=Cubature.L2)
     else
         val, err = Cubature.pcubature_v(
             length(Eω)*2,
             (x, fval) -> pointcalc!(fval, x, t),
-            (t.dimlimits[2][1],), (t.dimlimits[3][1],), 
+            (ll[1],), (ul[1],), 
             reltol=t.rtol, abstol=t.atol, maxevals=t.mfcn, error_norm=Cubature.L2)
     end
     nl .= reshape(reinterpret(ComplexF64, val), size(nl))

test/test_multimode.jl

@@ -58,7 +58,9 @@ end
     energyfun, energyfunω = Fields.energyfuncs(grid)
 
     dens0 = PhysData.density(gas, pres)
-    densityfun(z) = dens0
+    densityfun = let dens0=PhysData.density(gas, pres)
+        z -> dens0
+    end
     responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
     linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)
     inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=1e-6)
@@ -86,6 +88,58 @@ end
     @test norm(Iω - Iωf)/norm(Iω) < 0.003
 end
 
+@testset "Full, LP11 vs TM01" begin
+    # TM01 and LP11 (which is TM01+HE21) have the same dispersion and a ratio of effective
+    # areas of 2/3 -- mode averaged propagation in TM01 with 2/3x Aeff should be identical
+    # to full modal propagation in LP11, as long as only Kerr is considered.
+    using Luna
+    import LinearAlgebra: norm
+    a = 13e-6
+    gas = :Ar
+    pres = 5
+    τ = 30e-15
+    λ0 = 800e-9
+    energy = 1e-6
+    grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 0.4e-12)
+    m = Capillary.MarcatiliMode(a, gas, pres; kind=:TM, n=0, m=1, loss=false)
+    aeff(z) = 2/3*Modes.Aeff(m, z=z) # Aeff of LP11 is 2/3x smaller than that of TM01
+    energyfun, energyfunω = Fields.energyfuncs(grid)
+
+    dens0 = PhysData.density(gas, pres)
+    densityfun = let dens0=PhysData.density(gas, pres)
+        z -> dens0
+    end
+    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
+    linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0)
+    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy)
+    Eω, transform, FT = Luna.setup(
+        grid, densityfun, responses, inputs, βfun!, aeff)
+    statsfun = Stats.collect_stats(grid, Eω,
+                               Stats.ω0(grid),
+                               Stats.energy(grid, energyfunω))
+    output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
+    Luna.run(Eω, grid, linop, transform, FT, output, status_period=5)
+
+    # vertical LP11 double-lobe
+    modes = (
+        Capillary.MarcatiliMode(a, gas, pres, n=0, m=1, kind=:TM, ϕ=0.0, loss=false),
+        Capillary.MarcatiliMode(a, gas, pres, n=2, m=1, kind=:HE, ϕ=-π/4, loss=false),
+    )
+
+    field = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy/2)
+    inputs = ((mode=1, fields=(field,)), (mode=2, fields=(field,)))
+    # inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy)
+    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
+                                   modes, :xy; full=true)
+    outputf = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
+    linop = LinearOps.make_const_linop(grid, modes, λ0)
+    Luna.run(Eω, grid, linop, transform, FT, outputf, status_period=10)
+
+    Iω = abs2.(output["Eω"])
+    Iωf = dropdims(sum(abs2.(outputf["Eω"]); dims=2), dims=2)
+    @test norm(Iω - Iωf)/norm(Iω) < 0.0006
+end
+
 @testset "FieldInputs" begin
     using Luna
     import LinearAlgebra: norm
@@ -116,3 +170,45 @@ end
     @test Eω_single ≈ Eω_tuple
     @test Eω_single ≈ Eω_tuple_of_namedtuples
 end
+
+@testset "Nonlinear coupling" begin
+    using Luna
+    import LinearAlgebra: norm
+    a = 125e-6
+    gas = :Ar
+    pres = 0.167
+    flength = 3
+
+    τfwhm = 10e-15
+    λ0 = 800e-9
+    energy = 150e-6
+
+    # HE11 and LP11 (consisting of HE21 + TM01)
+    modes = (
+        Capillary.MarcatiliMode(a, gas, pres, n=0, m=1, kind=:TM, ϕ=0.0),
+        Capillary.MarcatiliMode(a, gas, pres, n=2, m=1, kind=:HE, ϕ=float(-π/4)),
+        Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=float(0.0)),
+    )
+
+    grid = Grid.RealGrid(flength, λ0, (160e-9, 3000e-9), 0.4e-12)
+
+    densityfun = let dens0=PhysData.density(gas, pres)
+        z -> dens0
+    end
+    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
+
+    field = Fields.GaussField(;λ0, τfwhm, energy=energy/2)
+    inputs = ((mode=1, fields=(field,)), (mode=2, fields=(field,)))
+
+    Eω, transform, FT = Luna.setup(
+        grid, densityfun, responses, inputs, modes, :xy; full=true)
+    nl = similar(Eω)
+
+    transform(nl, Eω, 0.0);
+    # nonlinear polarisation in LP11
+    Inl12 = dropdims(sum(abs2.(nl[:, 1:2]); dims=2); dims=2)
+    # nonlinear polarisation in HE11
+    Inl3 = abs2.(nl[:, 3])
+    # LP11 should not couple nonlinearly to HE11
+    @test norm(Inl3)/norm(Inl12) < 1e-32
+end

test/test_polarisation_field.jl

@@ -60,3 +60,67 @@ import Luna: Output
     @test all(output["stats"]["energy"] .≈ sum(outputp["stats"]["energy"], dims=1))
     @test all(output["stats"]["fwhm_r"] .≈ outputp["stats"]["fwhm_r"])
 end
+
+##
+@testset "Linear, LP11" begin
+    using Luna
+    import LinearAlgebra: norm
+    a = 125e-6
+    gas = :Ar
+    pres = 0.167
+    flength = 0.1
+
+    τ = 10e-15
+    λ0 = 800e-9
+    energy = 150e-6
+    grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12)
+
+    densityfun = let dens0=PhysData.density(gas, pres)
+        z -> dens0
+    end
+    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
+    energyfun, energyfunω = Fields.energyfuncs(grid)
+    modes = (
+        Capillary.MarcatiliMode(a, gas, pres, n=0, m=1, kind=:TM, ϕ=0.0, loss=false),
+        Capillary.MarcatiliMode(a, gas, pres, n=2, m=1, kind=:HE, ϕ=π/4, loss=false),
+    )
+    inf = (Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy/2.0),)
+    # same field in each mode
+    inputs = ((mode=1, fields=inf), (mode=2, fields=inf))
+    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
+                                modes, :xy; full=true)
+    statsfun = Stats.collect_stats(grid, Eω,
+                                Stats.ω0(grid),
+                                Stats.peakintensity(grid, modes),
+                                Stats.fwhm_r(grid, modes),
+                                Stats.energy(grid, energyfunω))
+    output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
+    linop = LinearOps.make_const_linop(grid, modes, λ0)
+    Luna.run(Eω, grid, linop, transform, FT, output, status_period=10)
+
+    modes = (
+        Capillary.MarcatiliMode(a, gas, pres, n=0, m=1, kind=:TM, ϕ=0.0, loss=false),
+        Capillary.MarcatiliMode(a, gas, pres, n=2, m=1, kind=:HE, ϕ=-π/4, loss=false),
+    )
+    inf = (Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy/2.0),)
+    # same field in each mode
+    inputs = ((mode=1, fields=inf), (mode=2, fields=inf))
+    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
+                                modes, :xy; full=true)
+    statsfun = Stats.collect_stats(grid, Eω,
+                                Stats.ω0(grid),
+                                Stats.peakintensity(grid, modes, components=:xy),
+                                Stats.fwhm_r(grid, modes, components=:xy),
+                                Stats.energy(grid, energyfunω))
+    outputp = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
+    linop = LinearOps.make_const_linop(grid, modes, λ0)
+    Luna.run(Eω, grid, linop, transform, FT, outputp, status_period=10)
+
+    Iω = dropdims(sum(abs2.(output.data["Eω"]), dims=2), dims=2)
+    Iωp = dropdims(sum(abs2.(outputp.data["Eω"]), dims=2), dims=2)
+
+    @test norm(Iω - Iωp)/norm(Iω) < 1.07e-12
+    @test all(output["stats"]["peakintensity"] .≈ outputp["stats"]["peakintensity"])
+    @test all(sum(output["stats"]["energy"], dims=1) .≈ sum(outputp["stats"]["energy"], dims=1))
+    @test all(output["stats"]["fwhm_r"] .≈ outputp["stats"]["fwhm_r"])
+end