Update example files.

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/run_rfcavity_multiple.py → examples/rfcavity/run_rfcavity_distinct.py

@@ -14,7 +14,7 @@
 sim.particle_shape = 2  # B-spline order
 sim.space_charge = False
 # sim.diagnostics = False  # benchmarking
-sim.slice_step_diagnostics = False
+sim.slice_step_diagnostics = True
 
 # domain decomposition & space charge mesh
 sim.init_grids()
@@ -46,6 +46,18 @@
 
 # design the accelerator lattice
 
+RF_cos_coefs = [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_coefs = [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]
+
 #   Drift elements
 dr1 = elements.Drift(name="dr1", ds=0.4, nslice=1)
 dr2 = elements.Drift(name="dr2", ds=0.032997, nslice=1)
@@ -56,64 +68,24 @@
     escale=62.0,
     freq=1.3e9,
     phase=85.5,
-    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,
-    ],
-    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,
-    ],
-    mapsteps=100,
-    nslice=4,
+    cos_coefficients=RF_cos_coefs,
+    sin_coefficients=RF_sin_coefs,
+    mapsteps=4,
+    nslice=100,
+)
+rf2 = elements.RFCavity(
+    name="rf2",
+    ds=1.31879807,
+    escale=62.0,
+    freq=1.3e9,
+    phase=85.5,
+    cos_coefficients=RF_sin_coefs,
+    sin_coefficients=RF_cos_coefs,
+    mapsteps=4, 
+    nslice=100,
 )
 
+
 # add beam diagnostics
 monitor = elements.BeamMonitor("monitor", backend="h5")
 
@@ -125,14 +97,7 @@
         rf,
         dr2,
         dr2,
-        rf,
-        dr2,
-        dr2,
-        rf,
-        dr2,
-        dr2,
-        rf,
-        dr2,
+        rf2,
         monitor,
     ]
 )