Correct and test fusion module in RZ geometry (#3255)

* Add particle rotation in NuclearFusionFunc.H

* Minor

* indent

* initial work

* fixed bugs and added species

* update documentation

* delete unused file

* Add properties for neutron, hydrogen isotopes, helium isotopes

* Update code to be more consistent

* Correct typo

* Parse deuterium-tritium fusion

* Start putting in place the files for deuterium-tritium

* Update documentation

* Prepare structures for deuterium tritium

* Fix typo

* Fix compilation

* Add neutron

* Add correct formula for the cross-section

* Correct compilation error

* Fix nuclear fusion

* Reset benchmarks

* Prepare creation functor for 2-product fusion

* First implementation of momentum initialization

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Use utility function for fusion

* Minor modification of variable names

* Fix GPU compilation

* Fix single precision compilation

* Update types

* Use util function in P-B fusion

* Correct compilation errors

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Correct errors

* Update values of mass and charge

* Correct compilation error

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Correct compilation error

* Correct compilation error

* Correct compilation error

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Reset benchmark

* Use helium particle in proton-boron, to avoid resetting benchmark

* Fixed proton-boron test

* Revert "Fixed proton-boron test"

This reverts commit 73c8d9d0be8417d5cd08a23daeebbc322c984808.

* Incorporate Neil's recommendations

* Reset benchmarks

* Correct compilation errors

* Add new deuterium tritium automated test

* Correct formula of cross-section

* Correct cross-section

* Improve analysis script

* Add test of energy conservation

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Add test of conservation of momentum

* Progress in analysis script

* Fix error in the initial energy of the deuterium particles

* Add check of isotropy

* Clean up the test script

* Rewrite p_sq formula in a way to avoids machine-precision negative numbers

* Add checksum

* Clean up code

* Add test for fusion in RZ geometry

* Update code to take into account actual timestep

* Update RZ test

* Update RZ test

* Increase number of particles

* Impart radial memory on DT particles

* Correct RZ momenta

* Remove unused file

* Update test

* Fix definition of theta

* Add new test

* Add checksum

* Update test

* Update tests

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Fix Python analysis script

* Remove CPU and ID from new benchmark

Co-authored-by: Yinjian Zhao <yinjianzhao@lbl.gov>
Co-authored-by: Luca Fedeli <luca.fedeli@cea.fr>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

1 files changed, 31 insertions, 0 deletions

diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
index a807a4262..759462427 100644
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
@@ -196,12 +196,36 @@ public:
                 multiplier_ratio = max_N;
             }
 
+#if (defined WARPX_DIM_RZ)
+            /* This momentum rotation is analogous to the one in ElasticCollisionPerez.H. */
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(WarpX::n_rz_azimuthal_modes==1,
+              "RZ mode `warpx.n_rz_azimuthal_modes` must be 1 when using the binary nuclear fusion module.");
+            amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta];
+            amrex::ParticleReal * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta];
+#endif
+
             for (int k = 0; k < max_N; ++k)
             {
                 // c1k : how many times the current particle of species 1 is paired with a particle
                 // of species 2. Same for c2k.
                 const int c1k = (k%NI1 < max_N%NI1) ? c1 + 1: c1;
                 const int c2k = (k%NI2 < max_N%NI2) ? c2 + 1: c2;
+
+#if (defined WARPX_DIM_RZ)
+                /* In RZ geometry, macroparticles can collide with other macroparticles
+                 * in the same *cylindrical* cell. For this reason, collisions between macroparticles
+                 * are actually not local in space. In this case, the underlying assumption is that
+                 * particles within the same cylindrical cell represent a cylindrically-symmetry
+                 * momentum distribution function. Therefore, here, we temporarily rotate the
+                 * momentum of one of the macroparticles in agreement with this cylindrical symmetry.
+                 * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation;
+                 * there is a corresponding assert statement at initialization.) */
+                amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]];
+                amrex::ParticleReal const u1xbuf = u1x[I1[i1]];
+                u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
+                u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
+#endif
+
                 SingleNuclearFusionEvent(
                     u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ],
                     u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ],
@@ -211,6 +235,13 @@ public:
                     m_probability_threshold,
                     m_probability_target_value,
                     m_fusion_type, engine);
+
+#if (defined WARPX_DIM_RZ)
+                amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]];
+                u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
+                u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
+#endif
+
                 p_pair_indices_1[pair_index] = I1[i1];
                 p_pair_indices_2[pair_index] = I2[i2];
                 ++i1; if ( i1 == static_cast<int>(I1e) ) { i1 = I1s; }