LupoLab / Luna.jl

Nonlinear optical pulse propagator
MIT License
65 stars 27 forks source link

Processing multi-mode parameter scans #316

Open phockett opened 1 year ago

phockett commented 1 year ago

Thanks for the great tool @chrisbrahms, @jtravs.

I have what is, hopefully, a basic (==stupid user, and new to Julia) question: how does one process the results from a parameter scan when running with multiple modes?

I've run parameter scans as per the docs (http://lupo-lab.com/Luna.jl/dev/scans.html), and processed the resulting files without issue (per http://lupo-lab.com/Luna.jl/dev/scans.html#Processing-scan-output), but when I add modes to the prop_capillary run I can't work out how to do the processing step... the example code breaks (seems to be due to array dimension changes), and I've tried to fix it in a few inelegant, probably silly, ways, but now find myself stuck.

Is there some equivalent of modes=:sum that I can invoke for the processing step, or a getModes or similar line that needs to be added to the example?

Thanks in advance for any help on this.

chrisbrahms commented 1 year ago

Hi! Very glad you like Luna!

The difference in the output from a multimode simulation is that the field in the output no longer has size (Nω, Nz) but (Nω, Nm, Nz) where Nm is the number of modes. In the example in the docs, that means that also has size (Nλ, Nm, Nz), so the last line inside scanproc, which tries to slice into it with only two indices, fails.

The solution depends on what you want to know:

  1. If you want to have the mode-resolved output, you can simply replace the Iλ[:, end] with Iλ[:, :, end] and the output will retain the extra dimension, for you to do with what you want
  2. If you only care about the fundamental mode, you can pick that out with Iλ[:, 1, end]
  3. If you want the total spectral power density in all modes, you simply need to sum over them and drop the (now-singleton) dimension: dropdims(sum(Iλ[:, :, end]; dims=2); dims=2)

(by the way, I just realised that the example uses a silly way to extract the final slice of the field--you can also pass a zslice value to getIω to only extract one slice to begin with and save some processing)

The same thing applies to the output of getEt, though here you need to be careful to take the abs2 before summing over the modes if you need that (summing the fields makes no sense).

Another thing to note is that the line

max_peakpower = maximum(output["stats"]["peakpower"])

in the example gives the global maximum (i.e. the maximum peak power in any mode) when used on multimode outputs. The size of any statistic arrays which can be mode-resolved changes from (Nz,) to (Nm, Nz) as well, so depending on what you want, you may also need to sum over the modes or pick out the fundamental mode etc.

Hope this helps! If you have any further issues or weird errors, feel free to paste an example code in here and we can have a look.

phockett commented 1 year ago

Excellent, thanks @chrisbrahms. I've done a quick test with case (3), and that is working nicely for me I think.

I do plan to take a closer look at the modes later, but for the initial pass I'm just interested in scanning over somewhat large parameter ranges, and looking at the integrated outputs (we're currently interested in seeing what we can do for deep UV/VUV generation in a capillary, but don't have much experience here as yet, hence the broad sweeps to get started).

Meanwhile, I clearly need to spend some more time getting familiar with Julia!


For reference for anyone else struggling, here's the processing snippet with a sum over modes:

# Case with sum over modes
λ, Iλ, zstat, edens, max_peakpower = Processing.scanproc(outputdir) do output
    λ, Iλ = Processing.getIω(output, :λ)
    zstat = Processing.VarLength(output["stats"]["z"])
    edens = Processing.VarLength(output["stats"]["electrondensity"])
    max_peakpower = maximum(output["stats"]["peakpower"])
#     Processing.Common(λ), Iλ[:, end], zstat, edens, max_peakpower  # Default case (modes=1)
    Processing.Common(λ), dropdims(sum(Iλ[:, :, end]; dims=2); dims=2), zstat, edens, max_peakpower
end