Birefringence and computations on GPU vs CPU

[AndiGaw]

Hello,

Thank you for sharing your code. It is very useful and I am going to cite it once I am done with my paper(s).

I use your code to solve MMGNLSE for a birefringent fiber and I have a question. Namely, If I understand correctly, when considering birefringence and coupling between different spacial modes, your polarization code (sim.scalar=false) does not longer give correct results and cannot be used to calculate SK. For that reason, I run scalar simulations and provide a separate SK tensor for calSRSK.m file along with the SR tensor (this is since SK =/= SR). Further, for sanity check, I compare calculations done on GPU and CPU but I get different results. Since I save the SK/SR used in the computation, the clear difference between the GPU/CPU calculations is the number of SK/SR elements. I think it is linked to the summation S13 (from the supplementary materials). Is the summation S13 rigorous, or is it some sort of a physical approximation? As far as I understood, it is only used for the GPU computations, right?

Thanks in advance and best regards,
Andrzej

[AaHaHaa]

Hi, first thanks for using my code.

There was a bug in multimode polarization computation, but I fixed it last year. If you run into this problem, please re-download the code. I'm afraid you're using the old version. It can run correctly now.

I've checked very carefully the new code last year, so I think the current one should solve the polarization mode correctly. Could you please confirm with me that you're using the latest code?

If you're using the latest code, then please send me the script. I can test it myself for you.

About the old polarization bug, the problem lies in cuda's summation. I forgot to put int(midx & 1) in if-else statements, and the wrong "if (midx & 1)" makes no sense since "midx & 1" is a binary number. I added int(....) to fix it.

Since the old bug lied in cuda, the CPU computation should always be the correct one. (It should be......I cannot check it now; everything is updated to the new correct one.)

[AaHaHaa]

I am thinking about your statement:
"Namely, If I understand correctly, when considering birefringence and coupling between different spatial modes, your polarization code (sim.scalar=false) does not longer give correct results and cannot be used to calculate SK."

My code should be able to solve multimode polarization cases, so I don't understand your statements. The calculation of SK values is independent of birefringence since they're overlap integrals, only related to functions of x and y.

Due to your inconsistency between CPU and GPU results, as I mentioned previously, I kinda, or strongly, speculate that you might be using my old-version buggy code.

[AaHaHaa]

Another point is that
In GPU cuda, I've optimzied the computation with reducing computations with symmetry from real and imaginary parts (see the details in the supplement of my JOSAB multimode-gain paper). This, however, isn't implemented in CPU. CPU uses MATLAB complex number computations only. I kinda don't assume that people will use CPU for multimode scenarios.
That said, you might see inconsistency in "matrix size" in their SRa (or b) and SK arrays, due to implementations of computational reduction in GPU. However, you should be able to find the same numeric values between CPU's and GPU's.

[AndiGaw]

Hello,

Thank you for your swift reply.

Indeed, I have been using the old version of the code, which I downloaded just before the update.

I am now running a couple of tests to check if the problem is gone. I will let you know asap.

Best regards,
Andrzej

[AaHaHaa]

If you would like to run with some very special cases, such as with higher-order HE modes, I would recommend that you take off refine_scalar/polarized_S or other commands and let the code sums over all potential M^4 overlap integrals. Forget about speed and optimization. Make it work first. You might need to modify the cuda, taking off beginning_nonzeros and ending_nonzeros to remove these optimization processes.
CUDA's summing over all overlap integrals (though some might be redundant) is still faster than CPU computation.

To be honest, it might be better for you to use our old-version code, which has no optimization. It's easier for you to modify it. For example, I think its CUDA has no symmetry-related optimization I mentioned.

[AndiGaw]

My suggestion: Try to run with polarized case, and use your modification: SK = mode_info.SK; in line 158 of calc_SRSK().

When calculating beginning_nonzeros and ending_nonzeros, it also runs refine_scalar/polarized_S, which removes redundant SR/SK's. This operation is different for scalar and polarized cases. Please redownload my calc_SRSK.m; I've added more comments.

Thank you. I re-downloaded it. I have one question concerning it. Namely, while trying to provide my SR/SK and build these structs with modified calc_polarized_SRSK function I get the following error:

Error using parallel.gpu.CUDAKernel/feval
Number of input arguments for this kernel is fixed at 13

If I am not mistaken, this is linked to the fact that sim.cuda_SRSK.NumRHSArguments = 22 for polarized case and 13 for scalar simulations. Since, I cannot find sim.cuda_SRSK file I wanted to ask is this error linked to some mistake in calc_SRSK or I have to change something elsewhere in the code?

A curious question:
Are you doing research on spatial modes that are different with different polarizations?

I am doing research on the x/y polarized LP modes of a birefringent fiber. For the moment, I am not planning to run HE, but thank you for the info :).

To be honest, it might be better for you to use our old-version code, which has no optimization.

Do you mean this one: Multimode nonlinear fiber optics: Massively parallel numerical solver, tutorial, and outlook?

[AndiGaw]

Since, I cannot find sim.cuda_SRSK file I wanted to ask is this error linked to some mistake in calc_SRSK or I have to change something elsewhere in the code?

Nvm, I found it in setup_stepping_kernel ;)

[AaHaHaa]

If I am not mistaken, this is linked to the fact that sim.cuda_SRSK.NumRHSArguments = 22 for polarized case and 13 for scalar simulations. Since, I cannot find sim.cuda_SRSK file I wanted to ask is this error linked to some mistake in calc_SRSK or I have to change something elsewhere in the code?

I think you also need to modify my cuda to change the number of input arguments.

Do you mean this one: Multimode nonlinear fiber optics: Massively parallel numerical solver, tutorial, and outlook?

Not really; that's just the paper discussing what MPA is. You need our old-version code. See my introduction of this package in the first page of Github. If I were you who needs a huge change of code, it would be easier to write one from scratch.

I am doing research on the x/y polarized LP modes of a birefringent fiber. For the moment, I am not planning to run HE, but thank you for the info :).

How does the polarization significantly affect the LP spatial mode, which is the case that you're thinking about? But anyway, this is just my curiosity. Just want to remind you that if polarization can be decoupled from spatial mode, then you don't need to modify any of my code to run multimode polarization simulations, in case that you forget this.

[AndiGaw]

By the way, see my readme.pdf: The polarized fields is arranged as [A1x,A1y, A2x,A2y...] for example. If you want to study the polarized multimode case with excited x-polarized multimode fields, you should multiply with [1,1e-3,1,1e-3] so that there are weakly-excited y-polarized components with mainly x-polarized inputs. Since you mentioned `[1,1,0,0]', I'd like to notify you of this in case you wrongly thought it were [A1x,A2x,A1y,A2y].

I understood the [mode1+ mode1_ mode 2+ mode2_...] arrangement. I was referring to which modes serve as the input when num_modes = length(sim.midx)*2 is used. If I use mode1+ mode1_ as my input, I guess, num_modes = length(sim.midx)*2 treats each of them as a different spatial mode, while, in reality, they are differently polarized modes of the same spatial mode.

[AaHaHaa]

num_modes = length(sim.midx)*2 treats each of them as a different spatial mode, while, in reality, they are differently polarized modes of the same spatial mode.

This sounds confusing to me.

num_modes=length(sim.midx)*2 isn't sent into GMMNLSE_propagate(). There is no num_modes as an input argument of the propagate function. If you would like to create an initial field with polarized field in spatial mode 1 (weakly-excited field in spatial mode 2), then multiply a Gaussian with [1,1,1e-3,1e-3] for example if two (spatial mode 1)'s polarization modes are equally excited without phase difference. Perhaps you were simply rephrasing this; if yes, then we reach consensus.