Add combined function bend element.

docs/source/usage/examples.rst

@@ -25,6 +25,7 @@ This section allows you to **download input files** that correspond to different
    examples/quadrupole_softedge/README.rst
    examples/positron_channel/README.rst
    examples/cyclotron/README.rst
+   examples/cfbend/README.rst
 
 
 Unit tests

docs/source/usage/parameters.rst

@@ -235,6 +235,21 @@ Lattice Elements
 * ``<element_name>.type`` (``string``)
     Indicates the element type for this lattice element. This should be one of:
 
+         * ``cfbend`` for a combined function bending magnet. This requires these additional parameters:
+
+            * ``<element_name>.ds`` (``float``, in meters) the segment length
+
+            * ``<element_name>.rc`` (``float``, in meters) the bend radius
+
+            * ``<element_name>.k`` (``float``, in inverse meters squared) the quadrupole strength
+
+                = (magnetic field gradient in T/m) / (magnetic rigidity in T-m)
+
+              * k > 0 horizontal focusing
+              * k < 0 horizontal defocusing
+
+            * ``<element_name>.nslice`` (``integer``) number of slices used for the application of space charge (default: ``1``)
+
         * ``drift`` for a free drift. This requires these additional parameters:
 
             * ``<element_name>.ds`` (``float``, in meters) the segment length

docs/source/usage/python.rst

@@ -410,6 +410,18 @@ This module provides elements for the accelerator lattice.
       :param madx_file: file name to MAD-X file with beamline elements
       :param nslice: number of slices used for the application of space charge
 
+.. py:class:: impactx.elements.CFbend(ds, rc, k, nslice=1)
+
+   A combined function bending magnet.  This is an ideal Sbend with a normal quadrupole field component.
+
+   :param ds: Segment length in m.
+   :param rc: Radius of curvature in m.
+   :param k:  Quadrupole strength in m^(-2) (MADX convention)
+              = (gradient in T/m) / (rigidity in T-m)
+              k > 0 horizontal focusing
+              k < 0 horizontal defocusing
+   :param nslice: number of slices used for the application of space charge
+
 .. py:class:: impactx.elements.ConstF(ds, kx, ky, kt, nslice=1)
 
    A linear Constant Focusing element.

examples/CMakeLists.txt

@@ -565,3 +565,21 @@ add_impactx_test(cyclotron.py
     examples/cyclotron/analysis_cyclotron.py
     OFF  # no plot script yet
 )
+
+# Combined-function bend ########################################################
+#
+# w/o space charge
+add_impactx_test(cfbend
+    examples/cfbend/input_cfbend.in
+      ON  # ImpactX MPI-parallel
+      OFF  # ImpactX Python interface
+    examples/cfbend/analysis_cfbend.py
+    OFF  # no plot script yet
+)
+add_impactx_test(cfbend.py
+    examples/cfbend/run_cfbend.py
+      OFF  # ImpactX MPI-parallel
+      ON   # ImpactX Python interface
+    examples/cfbend/analysis_cfbend.py
+    OFF  # no plot script yet
+)

examples/cfbend/README.rst

@@ -0,0 +1,47 @@
+.. _examples-cfbend:
+
+Combined Function Bend
+======================
+
+A single combined function bending magnet (an ideal sector bend with an upright quadrupole field component added).  The magnet
+parameters are based a single CSBEND element appearing in the ELEGANT input file for the ALS-U lattice.
+
+The beam parameters are based on:
+C. Steier et al, "Status of the Conceptual Design of ALS-U", IPAC2017, WEPAB104, `DOI:10.18429/JACoW-IPAC2017-WEPAB104 <https://doi.org/10.18429/JACoW-IPAC2017-WEPAB104>`__ (2017).
+
+A 2 GeV electron bunch with normalized transverse rms emittance of 50 pm undergoes a 3.76 deg bend.
+
+In this test, the initial and final values of :math:`\sigma_x`, :math:`\sigma_y`, :math:`\sigma_t`, :math:`\epsilon_x`, :math:`\epsilon_y`, and :math:`\epsilon_t` must agree with nominal values.
+
+
+Run
+---
+
+This example can be run as a Python script (``python3 run_cfbend.py``) or with an app with an input file (``impactx input_cfbend.in``).
+Each can also be prefixed with an `MPI executor <https://www.mpi-forum.org>`__, such as ``mpiexec -n 4 ...`` or ``srun -n 4 ...``, depending on the system.
+
+.. tab-set::
+
+   .. tab-item:: Python Script
+
+       .. literalinclude:: run_cfbend.py
+          :language: python3
+          :caption: You can copy this file from ``examples/cfbend/run_cfbend.py``.
+
+   .. tab-item:: App Input File
+
+       .. literalinclude:: input_cfbend.in
+          :language: ini
+          :caption: You can copy this file from ``examples/cfbend/input_cfbend.in``.
+
+
+Analyze
+-------
+
+We run the following script to analyze correctness:
+
+.. dropdown:: Script ``analysis_cfbend.py``
+
+   .. literalinclude:: analysis_cfbend.py
+      :language: python3
+      :caption: You can copy this file from ``examples/cfbend/analysis_cfbend.py``.

examples/cfbend/analysis_cfbend.py

@@ -0,0 +1,104 @@
+#!/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.8 * 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],
+    [
+        5.0e-6,
+        8.0e-6,
+        0.0599584916,
+        1.277171100130637e-14,
+        1.277171100130637e-14,
+        5.396264243e-005,
+    ],
+    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.8 * 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],
+    [
+        1.98301912761436476154e-005,
+        1.92110147814673522465e-006,
+        5.99584916026070271843e-002,
+        3.90166131470905952893e-010,
+        1.27717110013063679910e-014,
+        5.396264244141050920556e-005,
+    ],
+    rtol=rtol,
+    atol=atol,
+)

examples/cfbend/input_cfbend.in

@@ -0,0 +1,47 @@
+###############################################################################
+# Particle Beam(s)
+###############################################################################
+beam.npart = 10000
+beam.units = static
+beam.energy = 2.0e3  #2 GeV
+beam.charge = 1.0e-9
+beam.particle = electron
+beam.distribution = waterbag
+beam.sigmaX = 5.0e-6  #5 um
+beam.sigmaY = 8.0e-6  #8 um
+beam.sigmaT = 0.0599584916  #200 ps
+beam.sigmaPx = 2.5543422003e-9 #exn = 50 pm-rad
+beam.sigmaPy = 1.5964638752e-9 #eyn = 50 pm-rad
+beam.sigmaPt = 9.0e-4  #approximately dE/E
+beam.muxpx = 0.0
+beam.muypy = 0.0
+beam.mutpt = 0.0
+
+
+###############################################################################
+# Beamline: lattice elements and segments
+###############################################################################
+lattice.elements = monitor cfbend1 monitor
+
+lattice.nslice = 25
+
+cfbend1.type = cfbend
+cfbend1.ds = 0.5       # projected length 0.5 m, angle 3.76 deg
+cfbend1.rc = 7.613657587094493   # bending radius [m]
+cfbend1.k = -7.057403   # (upright) quadrupole component [m^(-2)]
+
+monitor.type = beam_monitor
+monitor.backend = h5
+
+
+###############################################################################
+# Algorithms
+###############################################################################
+algo.particle_shape = 2
+algo.space_charge = false
+
+
+###############################################################################
+# Diagnostics
+###############################################################################
+diag.slice_step_diagnostics = false

examples/cfbend/run_cfbend.py

@@ -0,0 +1,67 @@
+#!/usr/bin/env python3
+#
+# Copyright 2022-2023 ImpactX contributors
+# Authors: Chad Mitchell, Axel Huebl
+# License: BSD-3-Clause-LBNL
+#
+# -*- coding: utf-8 -*-
+
+import amrex.space3d as amr
+from impactx import ImpactX, RefPart, distribution, elements
+
+sim = ImpactX()
+
+# set numerical parameters and IO control
+sim.particle_shape = 2  # B-spline order
+sim.space_charge = False
+# sim.diagnostics = False  # benchmarking
+sim.slice_step_diagnostics = True
+
+# domain decomposition & space charge mesh
+sim.init_grids()
+
+# load a 5 GeV electron beam with an initial
+# normalized transverse rms emittance of 1 um
+energy_MeV = 2.0e3  # reference energy
+bunch_charge_C = 1.0e-9  # used with space charge
+npart = 10000  # number of macro particles
+
+#   reference particle
+ref = sim.particle_container().ref_particle()
+ref.set_charge_qe(-1.0).set_mass_MeV(0.510998950).set_energy_MeV(energy_MeV)
+
+#   particle bunch
+distr = distribution.Waterbag(
+    sigmaX=5.0e-6,  # 5 um
+    sigmaY=8.0e-6,  # 8 um
+    sigmaT=0.0599584916,  # 200 ps
+    sigmaPx=2.5543422003e-9,  # exn = 50 pm-rad
+    sigmaPy=1.5964638752e-9,  # eyn = 50 pm-rad
+    sigmaPt=9.0e-4,  # approximately dE/E
+    muxpx=0.0,
+    muypy=0.0,
+    mutpt=0.0,
+)
+sim.add_particles(bunch_charge_C, distr, npart)
+
+# add beam diagnostics
+monitor = elements.BeamMonitor("monitor", backend="h5")
+
+# design the accelerator lattice
+ns = 25  # number of slices per ds in the element
+
+bend = [
+    monitor,
+    elements.CFbend(ds=0.5, rc=7.613657587094493, k=-7.057403, nslice=ns),
+    monitor,
+]
+
+# assign a lattice segment
+sim.lattice.extend(bend)
+
+# run simulation
+sim.evolve()
+
+# clean shutdown
+del sim
+amr.finalize()

examples/cfbend/run_cfbend_madx.py

@@ -0,0 +1,54 @@
+#!/usr/bin/env python3
+#
+# Copyright 2022-2023 ImpactX contributors
+# Authors: Chad Mitchell, Axel Huebl
+# License: BSD-3-Clause-LBNL
+#
+# -*- coding: utf-8 -*-
+
+
+import amrex.space3d as amr
+from impactx import ImpactX, RefPart, distribution, elements
+
+sim = ImpactX()
+
+# set numerical parameters and IO control
+sim.particle_shape = 2  # B-spline order
+sim.space_charge = False
+# sim.diagnostics = False  # benchmarking
+sim.slice_step_diagnostics = True
+
+# domain decomposition & space charge mesh
+sim.init_grids()
+
+# load a 5 GeV electron beam with an initial
+# normalized transverse rms emittance of 1 um
+bunch_charge_C = 1.0e-9  # used with space charge
+npart = 10000  # number of macro particles
+
+#   reference particle
+ref = sim.particle_container().ref_particle().load_file("chicane.madx")
+
+#   particle bunch
+distr = distribution.Waterbag(
+    sigmaX=2.2951017632e-5,
+    sigmaY=1.3084093142e-5,
+    sigmaT=5.5555553e-8,
+    sigmaPx=1.598353425e-6,
+    sigmaPy=2.803697378e-6,
+    sigmaPt=2.000000000e-6,
+    muxpx=0.933345606203060,
+    muypy=0.933345606203060,
+    mutpt=0.999999961419755,
+)
+sim.add_particles(bunch_charge_C, distr, npart)
+
+# design the accelerator lattice
+sim.lattice.load_file("chicane.madx", nslice=25)
+
+# run simulation
+sim.evolve()
+
+# clean shutdown
+del sim
+amr.finalize()