Picongpu: CFL condition in LWFA example

[StevE-Ong]

HI, I run the LWFA simulation with the following resolution:

    /** Duration of one timestep
     *  unit: seconds */
    constexpr float_64 DELTA_T_SI = 0.70312e-16;

    /** equals X
     *  unit: meter */
    constexpr float_64 CELL_WIDTH_SI = 0.234375e-6;
    /** equals Y - the laser & moving window propagation direction
     *  unit: meter */
    constexpr float_64 CELL_HEIGHT_SI = 4.39453125e-8;
    /** equals Z
     *  unit: meter */
    constexpr float_64 CELL_DEPTH_SI = CELL_WIDTH_SI;

The grid size is "256 2048 256". In the output, CFL condition gives

   PIConGPUVerbose PHYSICS(1) | Sliding Window is ON
   PIConGPUVerbose PHYSICS(1) | used Random Number Generator: RNGProvider3XorMin seed: 42
   PIConGPUVerbose PHYSICS(1) | Courant c*dt <= 2.01515 ? 1
   PIConGPUVerbose PHYSICS(1) | species e: omega_p * dt <= 0.1 ? 0.0153627

The simulation is running good until step 20000. The laser is doing good too. There is something wrong with the number displayed. I selected the time step to be cdt~0.5dy.

[steindev]

Hey, according to the standard reference for the Yee scheme by Taflove A, Hagness SC. Computational electrodynamics: the finite-difference time-domain method. Artech house; 2005, the CFL condition on the timestep for numerical stability of the Yee field solver is

c * dt <= 1 / sqrt( (1/dx)**2 + (1/dy)**2 + (1/dz)**2 ).

Note that this is different from the Wikipedia: Courant-Friedrichs-Lewy condition since the relation above is specific for the Yee field solver.

According to the above condition, your timestep should be smaller than 1.41690e-16.
Your choice dt = 0.5 dy / c = 7.329289201131271e-17 is therefore correct, but you chose a value which can be almost a factor of two larger in order to reduce the accumulation of numerical errors from timestepping.

You need to read the above output as

c dt [<=](needs to be smaller than or equal to) 2.01515 [?](and its value is) 1

Therefore you are ready to go. Actually, if you do not fulfill the CFL condition, the simulation will not start.

I admit, the writing in the output is somewhat confusing. Which originates in the fact that the calculation is done in scaled units. All variables with a dimension of length are measured in a multiple of c dt, which you gave as input, all velocities are measured in multiples of c and time is measured in multiples of dt. Therefore the calculation of the right hand side in the above formula uses dx -> dx/(c*dt) (same for y and z) which results in 2.01515 for your case, and the value of c dt is always 1.

By the way, we usually use

dt = 0.995 / ( c * sqrt( (1/dx)**2 + (1/dy)**2 + (1/dz)**2 ) )

which can be easily implemented as a formula in grid.param.

[sbastrakov]

Thanks, @steindev, for a detailed explanation. I have also been confused by this output more than once. Maybe something like that would be more readable:
Courant [or CFL]: c*dt <= 1 / sqrt( (1/dx)**2 + (1/dy)**2 + (1/dz)**2 ) ? - Fulfilled [or True]: 1 <= 2.01515

[StevE-Ong]

Thanks @steindev for the explanations. Could it be possible if the timestep is automatically set from the Courant condition instead of manually input it and then check if it is correct. Maybe something like

c*dt <= Factor / sqrt( (1/dx)**2 + (1/dy)**2 + (1/dz)**2 )

where Factor is something less than 1 where we can set.

[sbastrakov]

@StevE-Ong I think this could be a good feature to have. However, this might not be that simple technically due to our compile-time parameters. For now, you can manually implement this behavior in a grid.param file like this.

[StevE-Ong]

Thanks @sbastrakov for the sample.

[sbastrakov]

@StevE-Ong is this issue resolved? If so, please feel free to close it. Our confusing output of condition checks is noted.

[steindev]

@sbastrakov: I think, something like Courant: c*dt = 1 <= 2.01515 ? True or Courant: c*dt <= 2.01515 ? c*dt = 1! would clear things up but still keeps the comprehensive form, with the latter being my personal favorite.

@StevE-Ong: Including a function by default is not useful. It depends on your setup whether you want to calculate the timestep from a given spatial resolution or vice versa. Or, as you did it, define a step independently of the grid spacing. I think, providing an understandable output is best to comply with the multitude of use cases. I will put forward a draft for a new output of condition check.

[ax3l]

Thanks to @steindev, yes the original publication of Yee did not report the proper CFL and Taflove A, Hagness SC fixed the convergence criteria.

Please be aware that the current output is useful, since you simply multiply your dt with the factor on the right to go close to the numerically recommended value of the CFL as possible ;) It gives you a ratio to the "ideal-CFL" which you will often find in dispersion analysis, e.g. https://zenodo.org/record/15924 (2.3.1).

For now, you can manually implement this behavior in a grid.param file like this.

Try not to go too close, with single precision, from my tests.