Add cascade collider #1326

whokion · 2024-07-17T22:05:10Z

Add a helper for intra-nucleus nucleon-nucleon collisions

Add TwodGridBuilder and angular distribution data for sampling cos\theta of the intra-nucleus nucleon-nucleon scattering.
Add CascadeCollider and CascadeParticle for sampling the final state of the ntra-nucleus cascade collision and an associated test.
Add namespace lorentz for FourVector and a couple of utility functions

…eon collisions

sethrj · 2024-07-18T10:37:03Z

src/celeritas/phys/FourVector.hh

+/*!
+ * Add two four-vectors.
+ */
+inline CELER_FUNCTION FourVector add(FourVector const& a, FourVector const& b)


Instead of add we might as well make it operator+, but perhaps first we should do operator+=?

Okay, I added operator+ and the assignment operator+=. I also realize that it may be better to template FourVector with the precision type, which I may do in another MR. I did not add other assignment/arithmetic operators as they are not likely used.

sethrj · 2024-07-18T10:42:15Z

src/celeritas/phys/FourVector.hh

@@ -28,6 +28,8 @@ struct FourVector
    real_type energy{0};  //!< Particle energy
 };

+namespace lorentz
+{


I don't think we need a whole namespace for this... and even if we did the entire file (especially the FourVector) should be in it.

Removed - I though that it might help for clarify.

src/celeritas/phys/FourVector.hh

sethrj · 2024-07-18T10:47:08Z

src/celeritas/neutron/model/NeutronInelasticModel.cc

        CELER_ASSERT(channel_data.par.slope > 0);
        xs_params.push_back(channel_data.par);

        GenericGridBuilder build_grid{&data.reals};
        make_builder(&data.nucleon_xs)
-            .push_back(build_grid(bins, make_span(channel_data.xs)));
+            .push_back(build_grid(xs_energy_bins, make_span(channel_data.xs)));


To allow deduplication of cross section grids etc., the build_grid (and the other builders) should be outside the for loop. I'd also replace make_builder(&data.nucleon_xs) inside the loop with CollectionBuilder nucleon_xs{&data.nucleon_xs}; outside of it. (The make_builder command was before we had C++17 CTAD.)

Good to know! Updated as suggested. Thanks.

sethrj · 2024-07-18T10:51:09Z

src/celeritas/neutron/model/NeutronInelasticModel.hh

    {
        StepanovParameters par;
-        ChannelArray xs;
+        Array<double, 13> xs;
+        Array<double, 19 * 6> cdf;  //! c.d.f[cos_theta_bins][cdf_energy_bins]


Since these are fixed dimensions it would be more sensible to have:

using ArrayEnergy = Array<double, 6>; Array<ArrayEnergy, 19> cdf; // [angle][energy]

Also since the constants each have only 4 digits of precision we could probably use float... or store even a 16-bit integer if we first multiply them by 10000 😆 it's fine to keep them as double though.

Okay, I keep them as they are for now.

@whokion Forgetting about the precision type: the documentation is wrong since the array is accessed and built with {energy, angle}

Suggested change

Array<double, 19 * 6> cdf; //! c.d.f[cos_theta_bins][cdf_energy_bins]

Array<double, 6 * 19> cdf; //! [energy][angle]

sethrj · 2024-07-18T13:46:43Z

src/celeritas/neutron/interactor/detail/CascadeParticle.hh

+    FourVector vec4;  //!< Four momentum in natural MevMomentum and
+                      //!< MevEnergy units


I'm not a big fan of the combined text+number; could it just be "momentum"?

Change vec4 to four_vec.

whokion · 2024-07-18T18:09:16Z

@sethrj Thanks a lot for all review comments. Will follow up later today! Regarding to vec4 --> momentum, I will think an alternative as a component of FourVector is mom (three vector), which is confusing with the proposed momentum (four vector). Maybe changing FourVector(mom, energy) --> (p, e) and then vec4 --> momentum or mom. Welcome to the 4-momentum world^^.

…pace lorentz

amandalund

Nice addition @whokion! Thanks! Just a couple of tiny documentation fixes.

src/celeritas/grid/TwodGridBuilder.hh

src/celeritas/neutron/interactor/detail/CascadeCollider.hh

whokion · 2024-07-19T23:55:30Z

@amandalund Thanks for the review and suggestions!

sethrj

Some of the arithmetic opertions in the collider need to be checked/clarified/simplified...

sethrj · 2024-07-20T10:52:25Z

src/celeritas/neutron/interactor/detail/CascadeCollider.hh

+        size_type idx = cos_grid.size() - 2;
+        real_type cdf_upper = 0;
+        real_type cdf_lower = 1;
+
+        do
+        {
+            cdf_upper = cdf_lower;
+            cdf_lower = calc_cdf({kin_energy_, cos_grid[idx]});
+        } while (cdf_lower > cdf && idx-- > 0);
+
+        real_type frac = (cdf - cdf_lower) / (cdf_upper - cdf_lower);
+        cos_theta = fma(frac, cos_grid[idx + 1] - cos_grid[idx], cos_grid[idx]);


Even though this algorithm is only ten lines, it's a little complicated, I think we'll see it a lot in tabulated sampling, and it could use some extra assertions/testing—after this merge request I'll create a helper function that encapsulates it.

Great! I originally tried an array of GenericGridData [energy bin], so that we can swap x an y after find a corresponding cos_grid_data for each energy bin, then use the GenericGridCalculator to sample cos(theta). I believe that two approaches are basically equivalent, but using Two2GridData is much more compact and cleaner. If we can encapsulate the inverse sampling inside Twod(Sub)gridCalculator, it would be even great! Thanks for your help.

sethrj · 2024-07-20T11:06:19Z

src/celeritas/neutron/model/NeutronInelasticModel.hh

    {
        StepanovParameters par;
-        ChannelArray xs;
+        Array<double, 13> xs;
+        Array<double, 19 * 6> cdf;  //! c.d.f[cos_theta_bins][cdf_energy_bins]


@whokion Forgetting about the precision type: the documentation is wrong since the array is accessed and built with {energy, angle}

Suggested change

Array<double, 19 * 6> cdf; //! c.d.f[cos_theta_bins][cdf_energy_bins]

Array<double, 6 * 19> cdf; //! [energy][angle]

src/celeritas/phys/FourVector.hh

test/celeritas/phys/FourVector.test.cc

src/celeritas/neutron/interactor/detail/CascadeCollider.hh

sethrj · 2024-07-20T12:28:16Z

src/celeritas/neutron/interactor/detail/CascadeCollider.hh

+            result[0].four_vec
+                = {{fv.mom[0] * vscm + fv.mom[1] * vxcm + fv.mom[2] * cm_dir},
+                   fv.energy};


This is exactly why I didn't want to make "array operators": even though this reads like a simple algebra expression, C++ is creating an insane number of temporaries here. I created some test kernels replicating this operation: as written, this takes 35 registers and 280 (!) bytes of stack space. Rewriting inside a for loop reduces the usage to 32 registers and 88 bytes of stack space, or even less (depending on how it's written).

Suggested change

result[0].four_vec

= {{fv.mom[0] * vscm + fv.mom[1] * vxcm + fv.mom[2] * cm_dir},

fv.energy};

for (int i = 0; i < 3; ++i)

{

result[0].four_vec[i] = fv.mom[0] * vscm[i]

+ fv.mom[1] * vxcm[i]

+ fv.mom[2] * cm_dir[i];

}

result[0].four_vec.energy = fv.energy;

Remember that 200 bytes of overflow per thread is 200MB of memory for the expected use case.

Good to know! I thought that they are very much equivalent. Anyway, changed as suggested.

Yeah, unfortunately C++ forces the creation of temporaries that all store local copies of the intermediate compilation. The only ways to get around it are to manually unroll like this or to use a complicated template metaprogramming math library like Eigen or Blitz++ (neither of which are CUDA compatible).

src/celeritas/neutron/interactor/detail/CascadeCollider.hh

sethrj · 2024-07-20T12:34:56Z

src/celeritas/neutron/interactor/detail/CascadeCollider.hh

+
+    if (!degenerated)
+    {
+        Real3 vscm = make_unit_vector(cm_velocity_ - vel_parallel * cm_dir);


This should use an axpy intead of array operators.

Do you mean axpy(-vel_parallel, cm_dir, &cm_velocity_) which will change value for cm_velocity?

Never mind. I got it.

sethrj · 2024-07-20T12:43:47Z

src/celeritas/neutron/interactor/detail/CascadeCollider.hh

+    // Booster vector in the center of mass frame
+    Real3 cm_velocity_;
+    // Momentum magnitude in the center of mass frame
+    real_type cm_p_;
+    // Kinetic energy in the target rest frame
+    real_type kin_energy_;


Suggested change

// Booster vector in the center of mass frame

Real3 cm_velocity_;

// Momentum magnitude in the center of mass frame

real_type cm_p_;

// Kinetic energy in the target rest frame

real_type kin_energy_;

// Boost vector in the center of mass frame [1/c]

Real3 cm_velocity_;

// Momentum magnitude in the center of mass frame [MeV / c]

real_type cm_p_;

// Kinetic energy in the target rest frame [MeV]

real_type kin_energy_;

sethrj · 2024-07-20T12:44:56Z

src/celeritas/neutron/interactor/detail/CascadeCollider.hh

+    // Degenerated if velocity perpendicular to the c.m. momentum is small
+    bool degenerated
+        = (dot_product(cm_velocity_, cm_velocity_) - ipow<2>(vel_parallel)
+           < this->epsilon());
+
+    if (!degenerated)


Are you sure this expression is right? It's comparing the difference in the square of speeds? Maybe it would be better to calculate the vscm and then check the magnitude of that...

Also it appears that epsilon here needs to have units of [1/c^2]? Is that right?

I agree that they should be equivalent, i.e., norm(vscm) == sqrt(dot_product(cm_velocity_, cm_velocity_) - ipow<2>(vel_parallel)). epsilon is small arbitrary magic number from Geant4 and used for condition for the degeneracy [1/c^2] and vscm [1/c] as c = 1 in Geant4. Okay, let's use norm(vscm) as the condition for the degeneracy, which should be equivalent for the final logic.

whokion · 2024-07-20T18:22:02Z

@sethrj Thanks again for all the comments and suggestions! I think that I addressed them for now.

sethrj

Thanks Soon! Looks great.

Soon Yung Jun added 4 commits July 17, 2024 16:49

Add TwodGridBuilder

27e9c5b

Add namespace lorentz and utility functions - add, norm

fd84def

Add angular cumulative distribution function data

49946c9

Add a helper, CascadeCollider for sampling intra-nucleus nucleon-nucl…

057140c

…eon collisions

whokion added enhancement New feature or request physics Particles, processes, and stepping algorithms labels Jul 17, 2024

whokion requested review from sethrj and amandalund July 18, 2024 03:12

sethrj reviewed Jul 18, 2024

View reviewed changes

Soon Yung Jun added 7 commits July 18, 2024 20:29

Add assignment and arithmetic operators for addition and remove names…

a789085

…pace lorentz

Remove namespace lerentz

76c3bdd

Add basic tests for FourVector functions

1a5b90b

Add celeritas/phys/FourVector.test

6899222

Use CollectionBuilder and allow deduplication of grids

7aa101e

Change vec4 to four_vec

99edd88

Fix failures from orange-float-g4@11.0

6ce3220

amandalund approved these changes Jul 19, 2024

View reviewed changes

Soon Yung Jun added 2 commits July 19, 2024 18:49

Fix the file path/name

4dbc028

Address review comments from @amandalund

a2f7d9e

sethrj requested changes Jul 20, 2024

View reviewed changes

Soon Yung Jun added 5 commits July 20, 2024 10:53

Fix the documentation for the ChannelData::cdf array

0990abb

Drop unnecessary std::move

c5a31ff

Use EXPECT_VEC_SOFT_EQ instead of EXPECT_EQ

d4448b0

Address review comments and suggestions from @sethrj

4b62cdb

Fix the degeneracy condition

ca28a8c

Merge branch 'develop' into cascade-collider

87baf70

sethrj approved these changes Jul 22, 2024

View reviewed changes

sethrj enabled auto-merge (squash) July 22, 2024 12:22

sethrj merged commit 7de9e78 into celeritas-project:develop Jul 22, 2024
29 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add cascade collider #1326

Add cascade collider #1326

whokion commented Jul 17, 2024 •

edited by sethrj

Loading

sethrj Jul 18, 2024

whokion Jul 19, 2024

sethrj Jul 18, 2024

whokion Jul 19, 2024

sethrj Jul 18, 2024

whokion Jul 19, 2024

sethrj Jul 18, 2024

whokion Jul 19, 2024

sethrj Jul 20, 2024

sethrj Jul 18, 2024

whokion Jul 19, 2024

whokion commented Jul 18, 2024

amandalund left a comment

whokion commented Jul 19, 2024

sethrj left a comment

sethrj Jul 20, 2024

whokion Jul 20, 2024

sethrj Jul 20, 2024

sethrj Jul 20, 2024 •

edited

Loading

whokion Jul 20, 2024

sethrj Jul 20, 2024

sethrj Jul 20, 2024

whokion Jul 20, 2024

whokion Jul 20, 2024 •

edited

Loading

sethrj Jul 20, 2024

sethrj Jul 20, 2024

whokion Jul 20, 2024 •

edited

Loading

whokion commented Jul 20, 2024

sethrj left a comment

	Array<double, 19 * 6> cdf; //! c.d.f[cos_theta_bins][cdf_energy_bins]
	Array<double, 6 * 19> cdf; //! [energy][angle]

		FourVector vec4; //!< Four momentum in natural MevMomentum and
		//!< MevEnergy units

-            result[0].four_vec
-                = {{fv.mom[0] * vscm + fv.mom[1] * vxcm + fv.mom[2] * cm_dir},
-                   fv.energy};
+                for (int i = 0; i < 3; ++i)
+                {
+                    result[0].four_vec[i] = fv.mom[0] * vscm[i]
+                                            + fv.mom[1] * vxcm[i]
+                                            + fv.mom[2] * cm_dir[i];
+                }
+                result[0].four_vec.energy = fv.energy;

Add cascade collider #1326

Add cascade collider #1326

Conversation

whokion commented Jul 17, 2024 • edited by sethrj Loading

Add a helper for intra-nucleus nucleon-nucleon collisions

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

whokion commented Jul 18, 2024

amandalund left a comment

Choose a reason for hiding this comment

whokion commented Jul 19, 2024

sethrj left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sethrj Jul 20, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

whokion Jul 20, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

whokion Jul 20, 2024 • edited Loading

Choose a reason for hiding this comment

whokion commented Jul 20, 2024

sethrj left a comment

Choose a reason for hiding this comment

whokion commented Jul 17, 2024 •

edited by sethrj

Loading

sethrj Jul 20, 2024 •

edited

Loading

whokion Jul 20, 2024 •

edited

Loading

whokion Jul 20, 2024 •

edited

Loading