Ionization module #260

MaxThevenet · 2019-08-08T23:01:10Z

This PR proposes an implementation for field ionization (using the ADK model) in WarpX. The main points are:

Ionization energies and initialization

The tables for ionization energies were taken from the NIST Atomic Spectra Database. This is handled with:

Utils/atomic_data.txt: the text file generated by NIST website.
Utils/write_atomic_data_cpp.py: a python script that reads atomic_data.txt and extracts ionization levels into a C++ file (see line below).
Utils/IonizationEnergiesTable.H: C++ file with ionization tables data and metadata.
PhysicalParticleContainer::InitIonizationModule handles the init of the ionization module. The mapping between ionizable species and their product species is done at the MultiParticleContainer level in mapSpeciesProduct.

Creation of new electrons

As ionization is a two-species process, it is handled at the MultiParticleContainer in MultiParticleContainer::doFieldIonization, in a 2-step approach. Within each tile:

An ionization mask is built (1 if particle is ionized, 0 otherwise) in PhysicalParticleContainer::buildIonizationMask.
The mapping between indices of ionized ions (source) and the corresponding created electron (product) is done via a prefix sum of the ionization mask (using amrex::Gpu::inclusive_scan for portability). Particle quantities are copied from source to product in createIonizedParticles.

Peripheral changes

Ionization required a number of peripheral changes, listed below:

The particle pusher and current deposition needed to be modified to account for ionization level. This is done with a new argument, ion_lev, that is a pointer to the array of ionization levels when ionization is on (nullptr otherwise).
The ionization level is the first runtime integer particle attribute in WarpX. This required some small additions in WarpXParticleContainer (GetiAttribs is the equivalent of GetAttribs for int attributes, and the mapping between the names and indices of int attributes is done in particle_icomps, similar to particle_comps for Real attributes).

Caveats

In principle, some energy should be extracted from the electromagnetic fields when ionization occurs. This can be done by depositing some current to account for the energy taken by the ionization process, but it is not included in this PR.
Only one single-level ionization can occur per ion per time step.
The ionization is currently done at the beginning of each time step. For more accuracy, it should be done right after the field gather (and after a half-push in position to have all particle data synchronized).

Tests

This implementation was tested for plasma mirror in 2D and 3D (the image below shows the reflection of a laser pulse on a silicon target with 1 electron ionized per ion initially, in 2d). The same simulation was run on GPU, and gave similar image (and ionization does not trigger CPU-GPU exchanges at each iteration). Other tests were done for laser-plasma acceleration, in the lab frame and in the boosted frame.

…nents

…ocal)

…on and definition

MaxThevenet · 2019-08-12T18:06:44Z

Thanks for starting this review @rlehe. I implemented your comments in a bunch of small commits.

I enabled ionization when sub cycling is on, and ran the automated inputs.rt test with:

MR, no sub-cycling (3200 iterations)
MR, sub-cycling (1600 iterations)

This required changes in the input file:

amr.max_level = 1
warpx.fine_tag_lo = -3.e-6    6.e-6 
warpx.fine_tag_hi =  3.e-6   11.e-6
warpx.do_subcycling = 0

In both cases, the patch occupies half of the plasma and the test is successful.

I can add these to the automated tests if you think it is useful, though they take a little bit of time.

RemiLehe

Thanks again for this PR.
Here are a few more comments.
Also: did you test the code on GPU (using the Python automatic tests)?

RemiLehe · 2019-08-19T20:22:02Z

Examples/Modules/ionization/inputs.bf.rt

+ions.physical_element = N
+ions.do_continuous_injection=1
+
+electrons.mass = m_p


Should this be m_e?

I implemented the changes you suggested in the input files.

RemiLehe · 2019-08-19T20:22:26Z

Examples/Modules/ionization/inputs.rt

+ions.ionization_product_species = electrons
+ions.physical_element = N
+
+electrons.mass = m_p


RemiLehe · 2019-08-19T20:26:43Z

Examples/Modules/ionization/inputs.bf.rt

+particles.nspecies = 2
+particles.species_names = electrons ions
+
+ions.mass = m_p


Technically, this should be ions.mass = 14*m_p (approximate value), this the ions are Nitrogen.

RemiLehe · 2019-08-19T20:26:58Z

Examples/Modules/ionization/inputs.rt

+particles.nspecies = 2
+particles.species_names = electrons ions
+
+ions.mass = m_p


RemiLehe · 2019-08-19T20:30:18Z

Examples/Modules/ionization/inputs.2d

@@ -0,0 +1,79 @@
+max_step = 500 


Could this file be removed?
It does not seem to be used in the automated tests.

RemiLehe · 2019-08-19T20:36:14Z

Source/Particles/Deposition/CurrentDeposition.H

@@ -26,6 +26,7 @@ void doDepositionShapeN(const amrex::Real * const xp,
                        const amrex::Real * const uxp,
                        const amrex::Real * const uyp,
                        const amrex::Real * const uzp,
+                        const int * const ion_lev,


Could you add the documentation from ion_lev in the above Doxygen description?
In particular, could you mention that:

This is needed in order to have the charge of each macroparticle (since q is a scalar here).

For non-ionizable species, ion_lev is a null pointer.

Comments added.

RemiLehe · 2019-08-19T20:36:54Z

Source/Particles/Deposition/CurrentDeposition.H

@@ -36,6 +37,7 @@ void doDepositionShapeN(const amrex::Real * const xp,
                        const amrex::Real stagger_shift, 
                        const amrex::Real q)
 {
+    const bool do_ionization = ion_lev;


Could you add a comment saying that this line checks whether ion_lev is a null pointer?

Comments added.

RemiLehe · 2019-08-19T20:58:20Z

Source/Particles/Deposition/CurrentDeposition.H

+            amrex::Real wq = q*wp[ip];
+            if (do_ionization){
+                wq *= ion_lev[ip];
+            }


This is maybe a bit dangerous, because here, the charge q is directly the charge that the user enters in the input script. Now for ionizable species, users might be a bit confused as to what charge should be entered (e.g. they might think that <species>.charge should be 0, if the ionizable species is initially neutral).
For this reason, it might be better if wq is independent of the user-entered charge, in the case of an ionizable species:

if (do_ionization){ wq = wp[ip] * ion_lev[ip] * PhysConst::q_e; } else { wq = wp[ip] * q; }

An alternative would be to explicitly check that the user-entered charge is q_e for an ionizable species.

I implemented the second option: if the user specifies a charge that is not q_e, the charge is over-written and the code raises a Warning.

RemiLehe · 2019-08-19T23:05:11Z

Source/Particles/MultiParticleContainer.cpp

+{
+    // For particle i in mfi, if is_ionized[i]=1, copy particle
+    // particle i from container pc_source into pc_product
+    static void createIonizedParticles (


Why does this need to be static void and in a namespace?
Could this go into a separate file, for ionization?

I removed the static keyword

RemiLehe · 2019-08-20T01:41:50Z

Also: could you add the modified ChargeDeposition (taking into account the ion_lev) to this PR?

RemiLehe · 2019-08-20T15:11:10Z

Source/Particles/Deposition/CurrentDeposition.H

+            amrex::Real wq  = q*wp[ip];
+            if (do_ionization){
+                wq *= ion_lev[ip];
+            }


See remark below, in the Esirkepov kernel.

RemiLehe · 2019-08-20T15:11:30Z

Source/Particles/Deposition/CurrentDeposition.H

@@ -179,6 +184,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp,
                                  const amrex::Real * const uxp,
                                  const amrex::Real * const uyp,
                                  const amrex::Real * const uzp,
+                                  const int * ion_lev,


Same remark as before: could you add the doc for ion_lev?

RemiLehe · 2019-08-20T15:17:23Z

Source/Particles/MultiParticleContainer.cpp

+                const int tile_id = mfi.LocalTileIndex();
+                pc_product->GetParticles(lev)[std::make_pair(grid_id,tile_id)];
+                pc_product->DefineAndReturnParticleTile(lev, grid_id, tile_id);
+            }


@MaxThevenet Could you add a comment in the above code, explaining why you need to touch each tile, and briefly explaining what DefineAndReturnParticleTile does?

RemiLehe · 2019-08-20T15:36:06Z

Source/Particles/PhysicalParticleContainer.cpp

            [=] AMREX_GPU_DEVICE (long i) {
+                Real qp = q;
+                if (ion_lev){ qp *= ion_lev[i]; }


Same remark as for the current deposition:
maybe it is safer to do

if (ion_lev) { qp = ion_lev[i] * PhysConst::q_e; } else { qp = q; }

RemiLehe · 2019-08-20T15:36:19Z

Source/Particles/PhysicalParticleContainer.cpp

            [=] AMREX_GPU_DEVICE (long i) {
+                Real qp = q;
+                if (ion_lev){ qp *= ion_lev[i]; }


Same remark

RemiLehe · 2019-08-20T15:42:01Z

Source/Particles/PhysicalParticleContainer.cpp

+    pp.query("ionization_initial_level", ionization_initial_level);
+    pp.get("ionization_product_species", ionization_product_name);
+    pp.get("physical_element", physical_element);
+    // Add Real component for ionization level


Should this comment be Add integer component for ionization level ?

RemiLehe · 2019-08-20T15:43:05Z

Source/Particles/PhysicalParticleContainer.cpp

+    pp.get("physical_element", physical_element);
+    // Add Real component for ionization level
+    AddIntComp("ionization_level");
+    plot_flags.resize(PIdx::nattribs + 1, 1);


Is the above line really needed? It seems to me that there is already a similar line in the ParticleIO.cpp file...

RemiLehe · 2019-08-20T15:48:28Z

Source/Utils/write_atomic_data_cpp.py

@@ -0,0 +1,67 @@
+import re, os


Could you add a string at the beginning of this file, explaining its purpose?
(e.g. this could be copy/pasted from this PR's description)

RemiLehe · 2019-08-20T16:59:32Z

Source/Particles/MultiParticleContainer.cpp

+        // The advantage of inclusive_scan is that the sum of is_ionized
+        // is in the last element, so no other reduction is required to get
+        // number of particles.
+        amrex::Gpu::inclusive_scan(is_ionized.begin(), is_ionized.end(), i_product.begin());


Could you add a comment here, saying that Gpu::inclusive_scan runs on the current GPU stream, and synchronizes with the CPU, so that the next line (executed by the CPU) has the updated values of i_product?

Comment added.

MaxThevenet · 2019-08-21T19:56:47Z

@RemiLehe Thanks for your remarks, I implemented the changes. A few comments:

For now, let's keep createIonizesParticles in MultiParticleContainer.cpp, and we'll update it after @atmyers makes changes to get rid of runtime components.
The code throws a warning if the user does not use charge=q_e for ionizable species, and the charge is overwritten
an error is raised if ionization is on when compiled for GPU, due to conflict between recent Redistribute and runtime particle components.

MaxThevenet added 30 commits August 3, 2019 17:27

add ionization_level component

7b9aa6d

add ionization initialization

58d728b

add ionization table data

43c0660

add target species handling

a0f579f

minor cleaning

d76aabe

typo

1617bfa

get dt from WarpX instance, to initialize ionization prefactors

dd12326

cleaning, and add capability for boosted frame runtime particle compo…

77df972

…nents

escape if no ionized particle

2245a69

use c++ deposition by default. includes old attribs

dab4a00

remove copy of old attribs, use current values instead for simplicity

f2d1b0f

cleaning, remote outdated functions

4be45ed

add example

e3e6577

add profile data, remove typo in example

30fa100

use amrex::inclusive_scan for portability

a27dbd9

merge conflicts

73a98af

exclusive_scan requires ManagedDeviceVector

b2754b6

Add doc

97882ad

Merge branch 'dev' into ion_lev

3596c60

Merge branch 'dev' into ion_lev

7ec19bd

Merge branch 'dev' of https://github.com/ECP-WarpX/WarpX into ion_lev

e84c1f7

change ionization level from Real to Int

c7a56a0

Merge branch 'dev' into ion_lev

b39198d

add some const, and remove unnecessary includes

8dcbb28

catch errors in inputs

e3bb554

add few comments

33c7e58

minor cleaning, remove commented lines

6123a53

remove Langmuir 2d psatd tests as it fails randomly (on Travis, not l…

95a0cd7

…ocal)

use Weiqun's convention: space before brackets for function declarati…

74f85ac

…on and definition

final minor cleaning

c87ee4b

MaxThevenet added 3 commits August 12, 2019 10:00

remove unnecessary comments in input files

74cf515

cleaning and allow ionization with subcycling

5131a6c

typo

fab45a5

add plotting capability in ionization rt analysis script

8469c54

MaxThevenet mentioned this pull request Aug 19, 2019

[WIP] Target species for ionization and QED #189

Closed

RemiLehe requested changes Aug 19, 2019

View reviewed changes

Merge branch 'dev' into ion_lev

930140e

RemiLehe requested changes Aug 20, 2019

View reviewed changes

RemiLehe reviewed Aug 20, 2019

View reviewed changes

MaxThevenet added 11 commits August 20, 2019 14:25

copy member variables to tmp local variables

83f9298

cleaning ionization examples

3d97f28

add comments and make sure charge=q_e for ionizable species

9795e23

Merge branch 'dev' into ion_lev

5ec7a45

charge deposition depends on ionization level

16dd4d0

add comments for charge deposition

d871e6f

add back element to plot flags

fe26d39

fix bug in particle output

b675832

had deleted call to DepositCharge by mistake. This fixes it

751c811

typo

2f5d6b1

error if ionization on GPU

e71f1cc

RemiLehe approved these changes Aug 21, 2019

View reviewed changes

RemiLehe merged commit 4dd5eff into dev Aug 21, 2019

MaxThevenet mentioned this pull request Aug 26, 2019

Implement ionization #153

Closed

RemiLehe deleted the ion_lev branch August 30, 2019 18:00

lucafedeli88 mentioned this pull request Sep 16, 2019

QED: now particles initialize a real component (optical depth) to be used for QED modules #344

Merged

ax3l mentioned this pull request Sep 23, 2019

Modernize SplitParticles for GPU #391

Closed

ax3l mentioned this pull request Jan 20, 2021

SpeciesProperties: C, N, O, Cu #1638

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Ionization module #260

Ionization module #260

MaxThevenet commented Aug 8, 2019 •

edited

Loading

MaxThevenet commented Aug 12, 2019 •

edited

Loading

RemiLehe left a comment

RemiLehe Aug 19, 2019

MaxThevenet Aug 21, 2019

RemiLehe Aug 19, 2019

RemiLehe Aug 19, 2019

RemiLehe Aug 19, 2019

RemiLehe Aug 19, 2019

RemiLehe Aug 19, 2019

MaxThevenet Aug 21, 2019

RemiLehe Aug 19, 2019

MaxThevenet Aug 21, 2019

RemiLehe Aug 19, 2019

MaxThevenet Aug 21, 2019

RemiLehe Aug 19, 2019

MaxThevenet Aug 21, 2019

RemiLehe commented Aug 20, 2019

RemiLehe Aug 20, 2019

RemiLehe Aug 20, 2019

RemiLehe Aug 20, 2019

RemiLehe Aug 20, 2019

RemiLehe Aug 20, 2019

RemiLehe Aug 20, 2019

RemiLehe Aug 20, 2019

RemiLehe Aug 20, 2019

RemiLehe Aug 20, 2019

MaxThevenet Aug 21, 2019

MaxThevenet commented Aug 21, 2019

Ionization module #260

Ionization module #260

Conversation

MaxThevenet commented Aug 8, 2019 • edited Loading

Ionization energies and initialization

Creation of new electrons

Peripheral changes

Caveats

Tests

MaxThevenet commented Aug 12, 2019 • edited Loading

RemiLehe 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

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

RemiLehe commented Aug 20, 2019

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

MaxThevenet commented Aug 21, 2019

MaxThevenet commented Aug 8, 2019 •

edited

Loading

MaxThevenet commented Aug 12, 2019 •

edited

Loading