Implement binary Coulomb collisions between plasma species #676

MaxThevenet · 2022-02-25T18:18:51Z

This PR proposes an implementation of binary Coulomb collisions between plasma species. It is largely inspired by the WarpX implementation (see ECP-WarpX/WarpX#539 for the original PR), based on this article by F. Pérez et al. In particular files ComputeTemperature.H, ElasticCollisionPerez.H, ShuffleFisherYates.H, and UpdateMomentumPerez.H were copy-pasted before being modified, and are kept as close as possible to the WarpX implementation on purpose, in view of a potential future merge. The main authors @Yin-YinjianZhao @ax3l @RemiLehe @atmyers are co-authors of this commit.

The module can be used by adding e.g.

plasmas.collisions = collision1
collision1.species = plasma plasma
collision1.CoulombLog = 20

in an input file (assuming the electron plasma species is called plasma).

Main changes

As mentioned above, files ComputeTemperature.H, ElasticCollisionPerez.H, ShuffleFisherYates.H, and UpdateMomentumPerez.H were copy-pasted from WarpX, and modified to work with the specifics of quasi-static PIC:

HiPACE++ does not store the longitudinal momentum of plasma particles, it uses the pseudo-potential psi instead. Therefore, uz is replaced by psi in the argument of a few functions.
The core of the collision module in UpdateMomentumPerez.H needs the actual uz. Conversions were added between psi and uz in this file.

Other changes are:

A new class CoulombCollision was created (see CoulombCollision.H/cpp). For each type of collisions enabled (i.e., between 2 plasma species, potentially twice the same), a new instance of this class is created. The vector of objects is stored in class MultiPlasma.
Minor changes to TileSort were implemented so it works on CPU and GPU.

In the loop over slices in HiPACE::Evolve, after each particle push by 1 longitudinal cell, plasma particles of the relevant species are sorted per-cell and collided pairwise within each cell.

Tests

Linear regime: comparison with WarpX

This HiPACE++ input file and this WarpX input file were run on GPU, giving the following:
Lineout: transverse cut of the fields, where we compare the general shape of the fields and the difference between with and without collisions:

x-z slice of the Ez field:

x-z slice of the Ex - c*By field:

The number of ppc was increased to 8x8x1 to observe a good agreement.

Non-linear regime: comparison with WarpX

This HiPACE++ input file and this WarpX input file were run on GPU, giving the following:
Lineout: transverse cut of the fields, where we compare the general shape of the fields and the difference between with and without collisions:

x-z slice of the Ez field:

x-z slice of the Ex - c*By field:

The agreement is not excellent for the non-linear regime (even without collisions), in particular because these runs are a bit far from convergence. However, changing the parameters turned out to be difficult because of ECP-WarpX/WarpX#2886.

TODO (for future PRs)

Fix the PhysConstSI stuff, that's not great.
Run with new psi+1
Fix dt: each plasma particle has a different dt depending on its psi, which has an impact on the algorithm in the relativistic regime.
Enable tiling: the current implementation has been tested on GPU and on CPU without tiling, but it should be crossed-checked (and maybe improved) when tiling and openMP threading are used.
Inter-species 1: electrons and (heavy, so immobile) protons (species have opposite charges)
Inter-species 2: electrons and (light, so mobile) positrons
Inter-species 3: electrons and heavier ions (e.g. Helium)

The three last tests can be performed in the WarpX vs. HiPACE++ framework.

…SA PIC loop Co-authored-by: Axel Huebl <axelhuebl@lbl.gov> Co-authored-by: Andrew Myers <atmyers@lbl.gov> Co-authored-by: Yinjian Zhao <yinjianzhao@lbl.gov> Co-authored-by: Remi Lehe <rlehe@lbl.gov>

docs/source/run/parameters.rst

…to collision2

SeverinDiederichs

To continue with this work, let's merge it with minor suggestions after @AlexanderSinn reviewed it.

Then, the following things must be fixed:

normalized units
adding the constant of motion (or unifying it with psi, as discussed offline, this needs some changes in the particle pusher)
not changing psi (this is non-physical), but rather the constant of motion
currently, the safety nets are minimal: if a user fails at providing correct names in collision.species = name name it just segfaults. It should be asserted that the names in the collision.species are in plasmas.names

src/particles/CoulombCollision.cpp

src/particles/ElasticCollisionPerez.H

src/particles/ComputeTemperature.H

tests/collisions.SI.1Rank.sh

Co-authored-by: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com>

…to collision2

AlexanderSinn

See review

src/particles/UpdateMomentumPerez.H

src/particles/PlasmaParticleContainer.H

MaxThevenet · 2022-03-16T08:33:02Z

Update: I merged dev into this branch to include the new psi+1 handling. Now it all works as expected (CI passes, as well as the tests above). This is ready for review/merge.

MaxThevenet · 2022-03-29T17:25:44Z

src/particles/UpdateMomentumPerez.H

+    // particle's Lorentz factor
+    const amrex::Real g1_tmp = (1.0_rt+u1x*u1x*inv_c2 + u1y*u1y*inv_c2 + psi1*psi1) / (2.0_rt*psi1);
+    amrex::Real u1z = PhysConstSI::c * (g1_tmp - psi1);
+
+    // particle's Lorentz factor
+    const amrex::Real g2_tmp = (1.0_rt+u2x*u2x*inv_c2 + u2y*u2y*inv_c2 + psi2*psi2) / (2.0_rt*psi2);
+    amrex::Real u2z = PhysConstSI::c * (g2_tmp - psi2);


This is what's different from WarpX

MaxThevenet · 2022-03-29T17:26:03Z

src/particles/UpdateMomentumPerez.H

+        const amrex::Real ga = std::sqrt( T_R(1.0) + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_c2 );
+        psi1 = ga - u1z/PhysConstSI::c;


as well as backward conversion here

MaxThevenet · 2022-03-29T17:26:09Z

src/particles/UpdateMomentumPerez.H

+        const amrex::Real ga = std::sqrt( T_R(1.0) + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_c2 );
+        psi2 = ga - u2z/PhysConstSI::c;


MaxThevenet and others added 3 commits February 25, 2022 18:27

add header files copy-pasted from WarpX, with modifications for the Q…

6991993

…SA PIC loop Co-authored-by: Axel Huebl <axelhuebl@lbl.gov> Co-authored-by: Andrew Myers <atmyers@lbl.gov> Co-authored-by: Yinjian Zhao <yinjianzhao@lbl.gov> Co-authored-by: Remi Lehe <rlehe@lbl.gov>

add rest of the implementation

cd56650

fix error in CMake file

2554dc7

MaxThevenet added documentation Improvements or additions to documentation component: plasma About the plasma species labels Feb 25, 2022

MaxThevenet requested review from AlexanderSinn and SeverinDiederichs February 25, 2022 18:18

MaxThevenet added 2 commits February 25, 2022 19:22

fix doxygen doc

4eaeda7

add PhysConstSI object

3c31cad

MaxThevenet changed the title ~~Implement binary Coulomb collisions between plasma species~~ [WIP] Implement binary Coulomb collisions between plasma species Feb 25, 2022

MaxThevenet added 2 commits February 25, 2022 20:54

update doc

b2839fc

Also expose SI constants as static constexpr

1a0c2e9

MaxThevenet changed the title ~~[WIP] Implement binary Coulomb collisions between plasma species~~ Implement binary Coulomb collisions between plasma species Feb 26, 2022

MaxThevenet mentioned this pull request Feb 28, 2022

[WIP] Collisions #675

Closed

9 tasks

Merge branch 'development' into collision2

5892620

MaxThevenet commented Mar 9, 2022

View reviewed changes

docs/source/run/parameters.rst Show resolved Hide resolved

MaxThevenet and others added 11 commits March 9, 2022 17:00

Update docs/source/run/parameters.rst

0299bea

Merge branch 'development' into collision2

d176d1d

add CI for collisions with SI

dfbe56a

add benchmark

66995fb

update script to reset benchmarks

d599cdd

fix CI

dace608

fix

1cca484

fix few typo and compiler warnings

55ea529

Merge branch 'collision2' of https://github.com/MaxThevenet/hipace in…

7e33807

…to collision2

fix benchmark by @SeverinDiederichs

04a32e3

Merge branch 'collision2' of https://github.com/MaxThevenet/hipace in…

a57ae37

…to collision2

SeverinDiederichs approved these changes Mar 14, 2022

View reviewed changes

MaxThevenet and others added 2 commits March 14, 2022 13:11

Apply suggestions from code review

9f0b120

Co-authored-by: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com>

Merge branch 'development' into collision2

d3558b7

MaxThevenet added 3 commits March 14, 2022 13:22

check user input

eb99fe8

Merge branch 'collision2' of https://github.com/MaxThevenet/hipace in…

3ebb159

…to collision2

compiler warning

3fb77b7

AlexanderSinn requested changes Mar 14, 2022

View reviewed changes

src/particles/UpdateMomentumPerez.H Show resolved Hide resolved

src/particles/UpdateMomentumPerez.H Show resolved Hide resolved

src/particles/PlasmaParticleContainer.H Outdated Show resolved Hide resolved

MaxThevenet commented Mar 14, 2022

View reviewed changes

src/particles/PlasmaParticleContainer.H Outdated Show resolved Hide resolved

MaxThevenet added 5 commits March 14, 2022 16:43

Apply suggestions from code review

c9a6857

Merge branch 'development' into collision2

d971f58

update to the new psi

6bde74e

fix typo

fd5be9c

fix other typo, now code runs, let's see CI

635e57a

MaxThevenet merged commit 95c4f2c into Hi-PACE:development Mar 16, 2022

MaxThevenet deleted the collision2 branch March 16, 2022 10:27

MaxThevenet mentioned this pull request Mar 16, 2022

Collisions #715

Open

6 tasks

MaxThevenet commented Mar 29, 2022

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement binary Coulomb collisions between plasma species #676

Implement binary Coulomb collisions between plasma species #676

MaxThevenet commented Feb 25, 2022 •

edited

Loading

SeverinDiederichs left a comment

AlexanderSinn left a comment

MaxThevenet commented Mar 16, 2022

MaxThevenet Mar 29, 2022

MaxThevenet Mar 29, 2022

MaxThevenet Mar 29, 2022

		const amrex::Real ga = std::sqrt( T_R(1.0) + (u1xu1x+u1yu1y+u1zu1z)inv_c2 );
		psi1 = ga - u1z/PhysConstSI::c;

		const amrex::Real ga = std::sqrt( T_R(1.0) + (u2xu2x+u2yu2y+u2zu2z)inv_c2 );
		psi2 = ga - u2z/PhysConstSI::c;

Implement binary Coulomb collisions between plasma species #676

Implement binary Coulomb collisions between plasma species #676

Conversation

MaxThevenet commented Feb 25, 2022 • edited Loading

Main changes

Tests

Linear regime: comparison with WarpX

Non-linear regime: comparison with WarpX

TODO (for future PRs)

SeverinDiederichs left a comment

Choose a reason for hiding this comment

AlexanderSinn left a comment

Choose a reason for hiding this comment

MaxThevenet commented Mar 16, 2022

MaxThevenet Mar 29, 2022

Choose a reason for hiding this comment

MaxThevenet Mar 29, 2022

Choose a reason for hiding this comment

MaxThevenet Mar 29, 2022

Choose a reason for hiding this comment

MaxThevenet commented Feb 25, 2022 •

edited

Loading