Refactor laser implementation #5318

steindev · 2025-03-18T09:05:34Z

This refactoring straightens the implementation of laser profiles and removes confusing terms being added and later subtracted.

Major changes are:

redefinition of the internal coordinate system: The focus becomes the actual 'origin' of the coordinate system, while the origin variable refers to the position on the Huygens surface where the field needs to be evaluated in order to feed it into the simulation
=> This allows to implement expressions for laser fields as they are printed in publications and text books w/o additional phase and envelope terms
the GaussianPulse can have different waist sizes along the transverse directions now
the pulse-front tilt feature of the GaussianPulse is removed, because it is unknown how pulse-front tilted Laguerre-Gaussian pulses propagate (our old implementation fed some field in the simulation which propagated physically toward the focus, but this (most probably) did not correspond to the expectation for the in-focus field)
=> For pulse-front tilted standard Gaussian pulses the DispersivePulse can and should be used
dispersions in the DispersivePulse are applied along the polarization axis, while the other transverse axis is dispersion-free

Draft status of this PR will be removed when the following action items are performed

GaussianPulse:
- Compile test simulation w/ and w/o Laguerre modes and for various propagation directions and polarizations
- Compare results of test sim to old implementation
- resolve possible issues
DispersivePulse:
- Compile test simulation for various propagation directions and dispersion values
- Compare results of test sim to old implementation
- resolve possible issues

Close #5269

steindev · 2025-04-23T10:17:40Z

@psychocoderHPC @PrometheusPi @chillenzer @ikbuibui PR is ready for review now.

It is intentionally split in three commits each of which refactors a specific part of the old implementation.

Test results have been discussed offline already and are documented in #5269.
Main takeaway from the tests: Results are equal between old and refactored implementation in various cases. Additionally, stray fields are reduced in the new implementation.

steindev · 2025-06-02T12:57:18Z

ping @PrometheusPi @finnolec @psychocoderHPC @ikbuibui @chillenzer. No review within 3 weeks is actually quite bad... 😞

chillenzer

You are right, here's my first batch of review. I've yet circumvented the two main files but you might want to start with what I've got here. It's mostly questions for getting an overview.

include/picongpu/fields/incidentField/Functors.hpp

include/picongpu/fields/incidentField/profiles/BaseParam.hpp

include/picongpu/fields/incidentField/profiles/DispersivePulse.def

include/picongpu/fields/incidentField/profiles/ExpRampWithPrepulse.def

include/picongpu/fields/incidentField/profiles/Wavepacket.hpp

chillenzer

Can't say much about the physical correctness but I've checked with the formulae from the given reference and it kinda makes sense to me. Some questions for my understanding and a few minor refactoring ideas are in the comments but they are mostly up to personal style I'd guess, so not blocking.

Wouldn't wanna merge myself in order to have a physicist have a look at this as well, so would you, @PrometheusPi, do the honours please?

chillenzer · 2025-06-03T12:29:43Z

include/picongpu/fields/incidentField/profiles/DispersivePulse.hpp

+                    float_X center = Unitless::SD * (Omega - Unitless::OMEGA0)
+                                     - sim.pic.getSpeedOfLight() * alpha * pos[0] / (Unitless::W0 * Unitless::OMEGA0);


Could be extracted into an own function.

Yes, could be. But I use it only once and its length is comparable to the other definitions. i.e. not as long as expandedWaveVectorX().

include/picongpu/fields/incidentField/profiles/DispersivePulse.hpp

include/picongpu/fields/incidentField/profiles/GaussianPulse.hpp

chillenzer · 2025-06-03T13:16:39Z

PS: Generally very nice documentation. With the given references and the additional hints in the comments it was quite straightforward to follow the formulae (mostly).

PrometheusPi · 2025-06-06T14:50:26Z

@steindev Will you first implement @chillenzer's suggestions are do you want to have more reviews at the current state?

psychocoderHPC · 2025-06-10T07:10:22Z

@PrometheusPi please already start the review. @chillenzer mostly ask for renamings.

chillenzer · 2025-06-10T09:23:47Z

Yes, definitely! Please review. My comments were at most low-priority colouring ideas concerning our new bike shed. Did you know that research has shown that bikes in a green shed are 12.657498% less likely to get stolen by people wearing a blue or yellow shirt? (If the shirt is blue and yellow, for example in stripes, the results are inconclusive and further research is needed.)

chillenzer · 2025-06-16T11:01:59Z

Any movement here?

PrometheusPi · 2025-06-16T11:04:10Z

@chillenzer soon (tomorrow)

steindev · 2025-06-17T09:44:34Z

@chillenzer Thanks for your review 🙏
Your questions are answered. I have had issues with one of your requests. Let's talk offline about this.

Also, it seems I do not see any need for changes from your comments, sorry 😅
But if a certain change is important to you, I will do it. Please tell me.

PrometheusPi

only minor things and questions

PrometheusPi · 2025-06-17T08:49:38Z

include/picongpu/fields/incidentField/Functors.hpp

@@ -1,4 +1,4 @@
-/* Copyright 2020-2024 Sergei Bastrakov, Finn-Ole Carstens
+/* Copyright 2020-2024 Sergei Bastrakov, Finn-Ole Carstens, Klaus Steiniger


Year will be updated next year - fine with me to keep it

PrometheusPi · 2025-06-17T08:58:28Z

include/picongpu/fields/incidentField/profiles/PlaneWave.hpp

-                        // unit: sim.unit.time()
-                        static constexpr float_X INIT_TIME = static_cast<float_X>(
-                            (Params::RAMP_INIT * Params::PULSE_DURATION_SI + Params::LASER_NOFOCUS_CONSTANT_SI)
-                            / sim.unit.time());


Just for my understanding: where does this get defined now?

Nowhere, because it has not been used. But RAMP_INIT is used in the code further below again.

PrometheusPi · 2025-06-17T09:01:03Z

include/picongpu/fields/incidentField/profiles/PlaneWave.def

@@ -38,8 +38,9 @@ namespace picongpu
                         *  unit: seconds */
                        static constexpr float_64 LASER_NOFOCUS_CONSTANT_SI = 13.34e-15;

-                        /** The laser pulse will be initialized after TIME_DELAY_SI half of PULSE_INIT times of the
-                         * PULSE_DURATION before and after the plateau unit: none */
+                        /** The laser pulse will be initialized after TIME_DELAY_SI half of RAMP_INIT times of the


I know I introduced this, but the factor half of ramp-init is confusing since the variable name already states init as the up-ramp in my opinion..

🤔 I see what you mean...will change it

PrometheusPi · 2025-06-17T09:19:07Z

include/picongpu/fields/incidentField/profiles/DispersivePulse.def

+             *
+             * unit: none
+             */
+            static constexpr float_X SPECTRAL_SUPPORT = 6._X;


@steindev Why did the spectral support increase from 4 to 6?

Results of numerical experiments. See #5269.

PrometheusPi · 2025-06-17T09:28:49Z

include/picongpu/fields/incidentField/profiles/DispersivePulse.def

+             *
+             *  unit: none
+             */
+            static constexpr float_64 PULSE_INIT = 17.0;


Why did this change from 20 to 17?

Results of numerical experiments. See #5269.

PrometheusPi · 2025-06-17T11:20:14Z

include/picongpu/fields/incidentField/profiles/DispersivePulse.hpp

+                /** Calculate incident field E value for the given position
+                 * The implementation is based on a definition of the electric field E given in frequency-space domain,
+                 * which is printed in eq. (6) of ref. [1]. Compared to the reference, higher order derivatives of the
+                 * frequency dependent propagation angle are set to zero, i.e. theta'' = theta''' = 0.


What are the ' standing for here?

Read ' as 'prime' and it means derivative with respect to frequency. Standard notation. Obviously... 😒 😉

PrometheusPi · 2025-06-17T11:25:56Z

include/picongpu/fields/incidentField/profiles/DispersivePulse.hpp

+                           * (Unitless::OMEGA0 * Unitless::AD * (Omega - Unitless::OMEGA0)
+                              + Unitless::AD * (Omega - Unitless::OMEGA0) * (Omega - Unitless::OMEGA0)
+                              - Unitless::OMEGA0 / 6.0_X * Unitless::AD * Unitless::AD * Unitless::AD
+                                    * (Omega - Unitless::OMEGA0) * (Omega - Unitless::OMEGA0)


for this equation, I would suggest storing (Omega - Unitless::OMEGA0) in a const variable before and use it here. This avoids 6-1=5 subtractions.

Also did it in the other functions.

Hey, 1970 just called: They are super-excited about their latest idea of Common Subexpression Elimination. They say that maybe in a distant future compilers might be capable of applying that to real-world code. I wonder how long we're gonna have to wait for this.

PrometheusPi · 2025-06-17T11:26:58Z

include/picongpu/fields/incidentField/profiles/DispersivePulse.hpp

+                                     - sim.pic.getSpeedOfLight() * alpha * pos[0] / (Unitless::W0 * Unitless::OMEGA0);
+
+                    // gaussian envelope in frequency domain
+                    float_X const envFreqExp = -(Omega - Unitless::OMEGA0) * (Omega - Unitless::OMEGA0)


here as well

PrometheusPi · 2025-06-17T11:27:04Z

include/picongpu/fields/incidentField/profiles/DispersivePulse.hpp

+                    float_X alpha = expandedWaveVectorX(Omega);
+
+                    // Center of a frequency's spatial distribution
+                    float_X center = Unitless::SD * (Omega - Unitless::OMEGA0)


here as well

PrometheusPi · 2025-06-17T11:27:37Z

include/picongpu/fields/incidentField/profiles/DispersivePulse.hpp

+                    float_X phase = Omega * pos[0] / sim.pic.getSpeedOfLight()
+                                    + 0.5_X * Unitless::GDD * (Omega - Unitless::OMEGA0) * (Omega - Unitless::OMEGA0)
+                                    + Unitless::TOD / 6.0_X * (Omega - Unitless::OMEGA0) * (Omega - Unitless::OMEGA0)
+                                          * (Omega - Unitless::OMEGA0)


here as well

steindev · 2025-06-17T19:20:23Z

Thanks for your review @PrometheusPi 🙏
I made changes to the code accordingly and directly squashed.

* getCurrentTime() -> getTminusXoverC() for separable profiles * make focus the origin of internal coordinate system * fix language typos * rename Unitless::w -> Unitless::OMEGA0 * delete PulseFrontTilt laser * Perform parameter range checks in ExpRampWithPrepulse

* update, expand, and fix language typos in doc-strings * use currentTimeOrigin & evaluate for all times * take Unitless::TIME_DELAY into account * apply dispersions only along polarization axis * move doubly performed computations into common place and compute only once * rename Unitless::w -> Unitless::OMEGA0 * decrease indentation level by namespace nesting * Add metadata

* update, expand, and fix language typos in doc-strings * allow for elliptic transverse shape by having different W0 per axis * remove pulse-front tilt (use via DispersivePulse for TEM_00, for all other modes the propagation is not known anyways) * adapt to new definition of internalCoordinates and time * straighten expressions for field computation (implementation follows Pausch et al. paper now) * rename Unitless::w -> Unitless::OMEGA0

steindev · 2025-06-18T20:00:22Z

@PrometheusPi: The latest changes are force-pushed again and run through the CI now.

Second, here is the comparison between 32bit and 64bit regarding the sqrt(sqrt(...)) discussion.

32bit:

64bit:

(no error in labels, really)

PrometheusPi · 2025-06-19T08:15:32Z

@steindev Thanks for providing the comparison between 32bit and 64bit. This look excellent.
Could you please also plot $|E_{32bit} - E_{64bit}| / E_{max}$?
Do you know, what causes the spike at x=[700, 800] and y~30 ?

steindev · 2025-06-20T10:53:04Z

@PrometheusPi

The difference plot
The spike: It lies directly on the incident field plane and does not propagate inside.

PrometheusPi · 2025-06-20T10:54:54Z

Thanks @steindev Looks good.

steindev added this to the 0.9.0 / next stable milestone Mar 18, 2025

steindev requested review from PrometheusPi, finnolec and fabidie March 18, 2025 09:06

steindev force-pushed the refactor-2025-02_Laser-Impl branch 2 times, most recently from 36092a3 to 83ef24e Compare March 19, 2025 09:48

steindev linked an issue Mar 21, 2025 that may be closed by this pull request

incident laser profiles need improvements #4718

Closed

steindev mentioned this pull request Mar 21, 2025

incident laser profiles need improvements #4718

Closed

steindev force-pushed the refactor-2025-02_Laser-Impl branch from 83ef24e to 4568194 Compare March 21, 2025 18:38

steindev force-pushed the refactor-2025-02_Laser-Impl branch from ac4d155 to b01cbd0 Compare April 23, 2025 10:09

steindev marked this pull request as ready for review April 23, 2025 10:11

steindev force-pushed the refactor-2025-02_Laser-Impl branch from b01cbd0 to 3e9c742 Compare April 23, 2025 15:30

steindev force-pushed the refactor-2025-02_Laser-Impl branch 4 times, most recently from d1c5a2f to e7ff6a8 Compare May 11, 2025 19:21

chillenzer reviewed Jun 3, 2025

View reviewed changes

chillenzer previously approved these changes Jun 3, 2025

View reviewed changes

PrometheusPi requested changes Jun 17, 2025

View reviewed changes

steindev dismissed chillenzer’s stale review via 65ec0a4 June 17, 2025 19:18

steindev force-pushed the refactor-2025-02_Laser-Impl branch from e7ff6a8 to 65ec0a4 Compare June 17, 2025 19:18

steindev added 3 commits June 18, 2025 20:32

steindev force-pushed the refactor-2025-02_Laser-Impl branch from 9fd531d to 4115299 Compare June 18, 2025 18:32

PrometheusPi approved these changes Jun 20, 2025

View reviewed changes

PrometheusPi merged commit 70a6709 into ComputationalRadiationPhysics:dev Jun 20, 2025
10 checks passed

steindev deleted the refactor-2025-02_Laser-Impl branch June 20, 2025 10:58

		float_X center = Unitless::SD * (Omega - Unitless::OMEGA0)
		- sim.pic.getSpeedOfLight() * alpha * pos[0] / (Unitless::W0 * Unitless::OMEGA0);

		@@ -1,4 +1,4 @@
		/* Copyright 2020-2024 Sergei Bastrakov, Finn-Ole Carstens
		/* Copyright 2020-2024 Sergei Bastrakov, Finn-Ole Carstens, Klaus Steiniger

Refactor laser implementation #5318

Refactor laser implementation #5318

Uh oh!

Conversation

steindev commented Mar 18, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

steindev commented Apr 23, 2025

Uh oh!

steindev commented Jun 2, 2025

Uh oh!

chillenzer left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

chillenzer left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

chillenzer commented Jun 3, 2025

Uh oh!

PrometheusPi commented Jun 6, 2025

Uh oh!

psychocoderHPC commented Jun 10, 2025

Uh oh!

chillenzer commented Jun 10, 2025

Uh oh!

chillenzer commented Jun 16, 2025

Uh oh!

PrometheusPi commented Jun 16, 2025

Uh oh!

steindev commented Jun 17, 2025

Uh oh!

PrometheusPi left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

steindev Jun 17, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

steindev commented Mar 18, 2025 •

edited

Loading

chillenzer left a comment •

edited

Loading

steindev Jun 17, 2025 •

edited

Loading

PrometheusPi commented Jun 19, 2025 •

edited

Loading