Picongpu: laser produces NaNs if `FOCUS_POS` is at `initPlaneY`

[hightower8083]

Current laser implementation fails if FOCUS_POS = initPlaneY(different units but you've got the idea). This happens in calculation of curvature radii (see include/picongpu/fields/laserProfiles/GaussianBeam.hpp)

float_X const R_y = -focusPos * ( 1.0_X + ( y_R / focusPos )*( y_R / focusPos ) );

since focusPos can actually be =0

This is easy to fix noting that in you only need 1/R_y, which can be introduced in a safer way, for example:

float_X  R_y_inv = 0.0_X;
if (focusPos != 0.0_X)
{
  R_y_inv = -1.0_X / focusPos / ( 1.0_X + ( y_R / focusPos )*( y_R / focusPos ) );
}

and used instead of 1/R_y

[hightower8083]

can try patching it, if you agree with such modif

[sbastrakov]

@hightower8083 it will be very welcome!
I think your suggestion of fixing is perfectly fine. Probably just moving the focusPos multiplier inside the brackets to avoid uncertainty should also work.

[hightower8083]

Well, the second term will stay infinite, which is logical -- at focus R_curve should go to infinity. I think using R_curve_inv is more healthy, cause it never hits any singularity. Will patch it asap

[ax3l]

The fix you suggest will introduce some more register usage but no branching, looks good to me. maybe use constexpr instead, so the compiler can throw out the code that is unused.

One problem I immediately see is, that one cannot compare floats that way easily for a robust fix. Instead, please compare to the "init plane's cell number" and the equivalent of the focus position rounded to a cell, which are both integers. We will assist you during PR review.

[sbastrakov]

@ax3l to our offline discussion. You are right, we do not want to use my suggestion to just transform the expression and have division by zero remaining.

[hightower8083]

@ax3l actually with inverse radius used everywhere there is no need for the conditional at all — its y/(y*y+yR*yR) and is a regular function

[n01r]

this can be closed now, right?

[hightower8083]

yep