Fix bug in unique ID specification for RF cavities.

examples/rfcavity/README.rst

@@ -52,3 +52,54 @@ We run the following script to analyze correctness:
    .. literalinclude:: analysis_rfcavity.py
       :language: python3
       :caption: You can copy this file from ``examples/rfcavity/analysis_rfcavity.py``.
+
+
+.. _examples-rfcavity-distinct:
+
+Acceleration by Multiple Distinct RF Cavities
+===============================================
+
+Beam accelerated through a pair of 2 distinct RF cavities (without space charge).
+
+We use a 230 MeV electron beam with initial normalized rms emittance of 1 um.
+
+The first RF cavity coincides with the RF cavity design used in ``examples-rfcavity``.
+
+In this test, the initial and final values of :math:`\lambda_x`, :math:`\lambda_y`, :math:`\lambda_t`, :math:`\epsilon_x`, :math:`\epsilon_y`, and :math:`\epsilon_>
+
+
+Run
+---
+
+This example can be run **either** as:
+
+* **Python** script: ``python3 run_rfcavity_distinct.py`` or
+* ImpactX **executable** using an input file: ``impactx input_rfcavity_distinct.in``
+
+For `MPI-parallel <https://www.mpi-forum.org>`__ runs, prefix these lines with ``mpiexec -n 4 ...`` or ``srun -n 4 ...``, depending on the system.
+
+.. tab-set::
+
+   .. tab-item:: Python: Script
+
+       .. literalinclude:: run_rfcavity_distinct.py
+          :language: python3
+          :caption: You can copy this file from ``examples/rfcavity/run_rfcavity_distinct.py``.
+
+   .. tab-item:: Executable: Input File
+
+       .. literalinclude:: input_rfcavity_distinct.in
+          :language: ini
+          :caption: You can copy this file from ``examples/rfcavity/input_rfcavity_distinct.in``.
+
+
+Analyze
+-------
+
+We run the following script to analyze correctness:
+
+.. dropdown:: Script ``analysis_rfcavity_distinct.py``
+
+   .. literalinclude:: analysis_rfcavity_distinct.py
+      :language: python3
+      :caption: You can copy this file from ``examples/rfcavity/analysis_rfcavity_distinct.py``.

examples/rfcavity/analysis_rfcavity_distinct.py

@@ -0,0 +1,97 @@
+#!/usr/bin/env python3
+#
+# Copyright 2022-2023 ImpactX contributors
+# Authors: Axel Huebl, Chad Mitchell
+# License: BSD-3-Clause-LBNL
+#
+
+import numpy as np
+import openpmd_api as io
+from scipy.stats import moment
+
+
+def get_moments(beam):
+    """Calculate standard deviations of beam position & momenta
+    and emittance values
+
+    Returns
+    -------
+    sigx, sigy, sigt, emittance_x, emittance_y, emittance_t
+    """
+    sigx = moment(beam["position_x"], moment=2) ** 0.5  # variance -> std dev.
+    sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
+    sigy = moment(beam["position_y"], moment=2) ** 0.5
+    sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
+    sigt = moment(beam["position_t"], moment=2) ** 0.5
+    sigpt = moment(beam["momentum_t"], moment=2) ** 0.5
+
+    epstrms = beam.cov(ddof=0)
+    emittance_x = (sigx**2 * sigpx**2 - epstrms["position_x"]["momentum_x"] ** 2) ** 0.5
+    emittance_y = (sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2) ** 0.5
+    emittance_t = (sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2) ** 0.5
+
+    return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)
+
+
+# initial/final beam
+series = io.Series("diags/openPMD/monitor.h5", io.Access.read_only)
+last_step = list(series.iterations)[-1]
+initial = series.iterations[1].particles["beam"].to_df()
+final = series.iterations[last_step].particles["beam"].to_df()
+
+# compare number of particles
+num_particles = 10000
+assert num_particles == len(initial)
+assert num_particles == len(final)
+
+print("Initial Beam:")
+sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(initial)
+print(f"  sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
+print(
+    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
+)
+
+atol = 0.0  # ignored
+rtol = 1.5 * num_particles**-0.5  # from random sampling of a smooth distribution
+print(f"  rtol={rtol} (ignored: atol~={atol})")
+
+assert np.allclose(
+    [sigx, sigy, sigt, emittance_x, emittance_y, emittance_t],
+    [
+        4.29466150443e-4,
+        2.41918588389e-4,
+        7.0399951912e-5,
+        2.21684103818e-9,
+        2.21684103818e-9,
+        1.83412186547e-8,
+    ],
+    rtol=rtol,
+    atol=atol,
+)
+
+
+print("")
+print("Final Beam:")
+sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(final)
+print(f"  sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
+print(
+    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
+)
+
+atol = 0.0  # ignored
+rtol = 1.5 * num_particles**-0.5  # from random sampling of a smooth distribution
+print(f"  rtol={rtol} (ignored: atol~={atol})")
+
+assert np.allclose(
+    [sigx, sigy, sigt, emittance_x, emittance_y, emittance_t],
+    [
+        3.52596000000e-4,
+        2.41775000000e-4,
+        7.0417917357e-5,
+        1.70893497973e-9,
+        1.70893497973e-9,
+        1.413901564889e-8,
+    ],
+    rtol=rtol,
+    atol=atol,
+)

examples/rfcavity/input_rfcavity_distinct.in

@@ -0,0 +1,120 @@
+###############################################################################
+# Particle Beam(s)
+###############################################################################
+beam.npart = 10000  # outside tests, use 1e5 or more
+beam.units = static
+beam.kin_energy = 230
+beam.charge = 1.0e-10
+beam.particle = electron
+beam.distribution = waterbag
+beam.lambdaX = 0.352498964601e-3
+beam.lambdaY = 0.207443478142e-3
+beam.lambdaT = 0.70399950746e-4
+beam.lambdaPx = 5.161852770e-6
+beam.lambdaPy = 9.163582894e-6
+beam.lambdaPt = 0.260528852031e-3
+beam.muxpx = 0.5712386101751441
+beam.muypy = -0.514495755427526
+beam.mutpt = -5.05773e-10
+
+
+###############################################################################
+# Beamline: lattice elements and segments
+###############################################################################
+lattice.elements = monitor dr1 dr2 rf dr2 dr2 rf2 monitor
+
+monitor.type = beam_monitor
+monitor.backend = h5
+
+dr1.type = drift
+dr1.ds = 0.4
+dr1.nslice = 1
+
+dr2.type = drift
+dr2.ds = 0.032997
+dr1.nslice = 1
+
+rf.type = rfcavity
+rf.ds = 1.31879807
+rf.escale = 62.0
+rf.freq = 1.3e9
+rf.phase = 85.5
+rf.mapsteps = 4
+rf.nslice = 100
+rf.cos_coefficients =                    \
+                0.1644024074311037       \
+                -0.1324009958969339      \
+                4.3443060026047219e-002  \
+                8.5602654094946495e-002  \
+                -0.2433578169042885      \
+                0.5297150596779437       \
+                0.7164884680963959       \
+                -5.2579522442877296e-003 \
+                -5.5025369142193678e-002 \
+                4.6845673335028933e-002  \
+                -2.3279346335638568e-002 \
+                4.0800777539657775e-003  \
+                4.1378326533752169e-003  \
+                -2.5040533340490805e-003 \
+                -4.0654981400000964e-003 \
+                9.6630592067498289e-003  \
+                -8.5275895985990214e-003 \
+                -5.8078747006425020e-002 \
+                -2.4044337836660403e-002 \
+                1.0968240064697212e-002  \
+                -3.4461179858301418e-003 \
+                -8.1201564869443749e-004 \
+                2.1438992904959380e-003  \
+                -1.4997753525697276e-003 \
+                1.8685171825676386e-004
+rf.sin_coefficients = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  \
+                0 0 0 0 0 0 0
+
+rf2.type = rfcavity
+rf2.ds = 1.31879807
+rf2.escale = 62.0
+rf2.freq = 1.3e9
+rf2.phase = 85.5
+rf2.mapsteps = 4
+rf2.nslice = 100
+rf2.cos_coefficients = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  \
+                0 0 0 0 0 0 0
+rf2.sin_coefficients =                    \
+                0.1644024074311037       \
+                -0.1324009958969339      \
+                4.3443060026047219e-002  \
+                8.5602654094946495e-002  \
+                -0.2433578169042885      \
+                0.5297150596779437       \
+                0.7164884680963959       \
+                -5.2579522442877296e-003 \
+                -5.5025369142193678e-002 \
+                4.6845673335028933e-002  \
+                -2.3279346335638568e-002 \
+                4.0800777539657775e-003  \
+                4.1378326533752169e-003  \
+                -2.5040533340490805e-003 \
+                -4.0654981400000964e-003 \
+                9.6630592067498289e-003  \
+                -8.5275895985990214e-003 \
+                -5.8078747006425020e-002 \
+                -2.4044337836660403e-002 \
+                1.0968240064697212e-002  \
+                -3.4461179858301418e-003 \
+                -8.1201564869443749e-004 \
+                2.1438992904959380e-003  \
+                -1.4997753525697276e-003 \
+                1.8685171825676386e-004
+
+
+###############################################################################
+# Algorithms
+###############################################################################
+algo.particle_shape = 2
+algo.space_charge = false
+
+
+###############################################################################
+# Diagnostics
+###############################################################################
+diag.slice_step_diagnostics = true