Add accelerator lattice, starting with quadrupoles (#3063)

* Initial version of accelerator lattice

* Clean up EOL white space

* Small clean up for GPU

* Fixed up consts

* Added hard edge fraction plus other clean ups

* More clean up

* Restructure to work on GPUs

* Now this grabs its own copies of particle info

* Updates, including adding dBdx

* Small cleanup in Quad

* Small fixes for GPU

* More cleanup for GPU

* More GPU cleanup

* Rewrite of the accelerator lattice implementation to better handle GPU

* Fix struct forward definition

* Another forward definition fix

* Bug fix

* Added LatticeElementBase

* Removed zcenters array

* Added CI test case

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

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

* Clean up in CI analysis.py

* Cleanup of coding

* Added CI test hard_edged_quadrupoles_moving

* Added Lorentz transform between boosted frame and lab frame

* Fixes for working in the boosted frame

* Added boosted CI test

* Change input name, adding the prefix "lattice."

* Added plasma lens lattice element

This will replace the external field plasma lens

* Fixed CI analysis script to look for "lattice.quad"

* Added checks of lattice element input

* Added documentation

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

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

* Removed duplicate call to lattice finder UpdateIndices

* Added extensive comments

* Reworked the input to use the MAD like description

This is the same as the method used in ImpactX

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

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

* Remove old lines from inputs_lattice_3d

* Added "lattice" element type

* Fixed some Real and ParticleReals

* [pre-commit.ci] pre-commit autoupdate (#3246)

updates:
- [github.com/hadialqattan/pycln: v2.0.1 → v2.0.3](https://github.com/hadialqattan/pycln/compare/v2.0.1...v2.0.3)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

* ABLASTR: Update Poisson Solver API (#3243)

Update the Poisson Solver API to be more usable. Needed for ImpactX.

* Docs: New OLCF Machine (#3228)

* D-T fusion (#3153)

* 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

* Apply suggestions from code review

* Update PR according to comments

* Update benchmark

* Address additional comments

* Numerical Literals

Co-authored-by: Luca Fedeli <luca.fedeli@cea.fr>
Co-authored-by: Neïl Zaim <49716072+NeilZaim@users.noreply.github.com>

* Docs: gaussian beam `q_tot` is not optional (#3249)

* Fix a bug in GPU version of Hankel Transform (#3253)

amrex::Array4 is a 4D array that can be accessed with three spatial indices
plus an optional component index.  We must always provide all three spatial
indices even in 2D.

* Add Python Callback Call when Checkpointing Signal is Received (#3251)

* CI: Add Missing Regression Analysis (NCI corrector) (#3252)

* Fixes to allow mixed precision, ParticleReal float, Real double (#3239)

* Fixes to allow mixed precision, ParticleReal float, Real double

* Fix for the optical depth

* A different way of fixing QuantumSynchrotronEvolveOpticalDepth

* In the QED code, consistently use ParticleReal

* Use ParticleReal type consistently

* Fix typo Docs/source/usage/parameters.rst

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Fix typo Docs/source/usage/parameters.rst

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* Fix small error in plasma lens lattice documentation

* Small addition to the documentation

* Fix the residence correction to allow short elements

* Updated CI benchmarks

* Added check of lattice to isNoOp

* Updated the hard_edged_quadrupoles CI benchmarks

It is not clear why there was a change, but the difference is
essentially round off in the E field. The important thing is that the
particles are still correct.

* Update Source/AcceleratorLattice/AcceleratorLattice.H

Add include statements

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Update Source/AcceleratorLattice/LatticeElements/LatticeElementBase.H

Add includes

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Renamed to README.rst and updated headers

* Made d_lattice_element_finder optional type

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

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

* Include `<optional>`

* Docs: Developer AccLattice Inclusion

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>
Co-authored-by: Luca Fedeli <luca.fedeli@cea.fr>
Co-authored-by: Neïl Zaim <49716072+NeilZaim@users.noreply.github.com>
Co-authored-by: Weiqun Zhang <WeiqunZhang@lbl.gov>
Co-authored-by: Roelof Groenewald <40245517+roelof-groenewald@users.noreply.github.com>
Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

1 files changed, 116 insertions, 0 deletions

diff --git a/Examples/Tests/AcceleratorLattice/analysis.py b/Examples/Tests/AcceleratorLattice/analysis.py
new file mode 100755
index 000000000..d2fd7f6ff
--- /dev/null
+++ b/Examples/Tests/AcceleratorLattice/analysis.py
@@ -0,0 +1,116 @@
+#!/usr/bin/env python3
+
+# Copyright 2022 David Grote
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+
+"""
+This script tests the quad lattice
+The input file sets up a series of quadrupoles and propagates two particles through them.
+One particle is in the X plane, the other the Y plane.
+The final positions are compared to the analytic solutions.
+The motion is slow enough that relativistic effects are ignored.
+"""
+
+import os
+import sys
+
+import numpy as np
+from scipy.constants import c, e, m_e
+import yt
+
+yt.funcs.mylog.setLevel(0)
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+filename = sys.argv[1]
+ds = yt.load( filename )
+ad = ds.all_data()
+
+gamma_boost = float(ds.parameters.get('warpx.gamma_boost', 1.))
+uz_boost = np.sqrt(gamma_boost*gamma_boost - 1.)*c
+
+# Fetch the final particle position
+xx_sim = ad['electron', 'particle_position_x'].v[0]
+zz_sim = ad['electron', 'particle_position_z'].v[0]
+ux_sim = ad['electron', 'particle_momentum_x'].v[0]/m_e
+
+if gamma_boost > 1.:
+    # The simulation data is in the boosted frame.
+    # Transform the z position to the lab frame.
+    time = ds.current_time.value
+    zz_sim = gamma_boost*zz_sim + uz_boost*time;
+
+# Fetch the quadrupole lattice data
+quad_starts = []
+quad_lengths = []
+quad_strengths_E = []
+z_location = 0.
+def read_lattice(rootname, z_location):
+    lattice_elements = ds.parameters.get(f'{rootname}.elements').split()
+    for element in lattice_elements:
+        element_type = ds.parameters.get(f'{element}.type')
+        if element_type == 'drift':
+            length = float(ds.parameters.get(f'{element}.ds'))
+            z_location += length
+        elif element_type == 'quad':
+            length = float(ds.parameters.get(f'{element}.ds'))
+            quad_starts.append(z_location)
+            quad_lengths.append(length)
+            quad_strengths_E.append(float(ds.parameters.get(f'{element}.dEdx')))
+            z_location += length
+        elif element_type == 'lattice':
+            z_location = read_lattice(element, z_location)
+    return z_location
+
+read_lattice('lattice', z_location)
+
+# Fetch the initial position of the particle
+x0 = [float(x) for x in ds.parameters.get('electron.single_particle_pos').split()]
+ux0 = [float(x)*c for x in ds.parameters.get('electron.single_particle_vel').split()]
+
+xx = x0[0]
+zz = x0[2]
+ux = ux0[0]
+uz = ux0[2]
+
+gamma = np.sqrt(uz**2/c**2 + 1.)
+vz = uz/gamma
+
+def applylens(x0, vx0, vz0, gamma, lens_length, lens_strength):
+    """Use analytic solution of a particle with a transverse dependent field"""
+    kb0 = np.sqrt(e/(m_e*gamma*vz0**2)*abs(lens_strength))
+    if lens_strength >= 0.:
+        x1 = x0*np.cos(kb0*lens_length) + (vx0/vz0)/kb0*np.sin(kb0*lens_length)
+        vx1 = vz0*(-kb0*x0*np.sin(kb0*lens_length) + (vx0/vz0)*np.cos(kb0*lens_length))
+    else:
+        x1 = x0*np.cosh(kb0*lens_length) + (vx0/vz0)/kb0*np.sinh(kb0*lens_length)
+        vx1 = vz0*(+kb0*x0*np.sinh(kb0*lens_length) + (vx0/vz0)*np.cosh(kb0*lens_length))
+    return x1, vx1
+
+# Integrate the particle using the analytic solution
+for i in range(len(quad_starts)):
+    z_lens = quad_starts[i]
+    vx = ux/gamma
+    dt = (z_lens - zz)/vz
+    xx = xx + dt*vx
+    xx, vx = applylens(xx, vx, vz, gamma, quad_lengths[i], quad_strengths_E[i])
+    ux = gamma*vx
+    zz = z_lens + quad_lengths[i]
+
+dt = (zz_sim - zz)/vz
+vx = ux/gamma
+xx = xx + dt*vx
+
+# Compare the analytic to the simulated final values
+print(f'Error in x position is {abs(np.abs((xx - xx_sim)/xx))}, which should be < 0.01')
+print(f'Error in x velocity is {abs(np.abs((ux - ux_sim)/ux))}, which should be < 0.002')
+
+assert abs(np.abs((xx - xx_sim)/xx)) < 0.01, Exception('error in x particle position')
+assert abs(np.abs((ux - ux_sim)/ux)) < 0.002, Exception('error in x particle velocity')
+
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, filename)