Adding EPAC2004 space charge benchmarks

examples/CMakeLists.txt

@@ -639,3 +639,21 @@ add_impactx_test(hvkicker_madx.py
         examples/kicker/analysis_kicker.py
         OFF  # no plot script yet
 )
+
+# FODO + RF EPAC2004 ##################################################
+#
+# with space charge
+add_impactx_test(fodo_rf_sc
+    examples/epac2004_benchmarks/input_fodo_rf_SC.in
+      ON  # ImpactX MPI-parallel
+      OFF  # ImpactX Python interface
+    examples/epac2004_benchmarks/analysis_fodo_rf_SC.py
+    OFF  # no plot script yet
+)
+add_impactx_test(fodo_rf_sc.py
+    examples/epac2004_benchmarks/run_fodo_rf_SC.py
+      OFF  # ImpactX MPI-parallel
+      ON   # ImpactX Python interface
+    examples/epac2004_benchmarks/analysis_fodo_rf_SC.py
+    OFF  # no plot script yet
+)

examples/epac2004_benchmarks/analysis_fodo_rf_SC.py

@@ -0,0 +1,109 @@
+#!/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}")
+
+atol = 0.0  # ignored
+rtol = 3.0 * num_particles**-0.5  # from random sampling of a smooth distribution
+print(f"  rtol={rtol} (ignored: atol~={atol})")
+
+assert np.allclose(
+    [sigx, sigy, sigt],
+    [9.84722273e-4, 6.96967278e-4, 4.486799242214e-03],
+    rtol=rtol,
+    atol=atol,
+)
+
+print(
+    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
+)
+
+atol = 4.0e-8
+print(f"  atol={atol}")
+
+assert np.allclose(
+    [emittance_x, emittance_y, emittance_t],
+    [0.0, 0.0, 0.0],
+    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}")
+
+atol = 0.0  # ignored
+rtol = 3.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],
+    [9.84722273e-4, 6.96967278e-4, 4.486799242214e-03],
+    rtol=rtol,
+    atol=atol,
+)
+
+print(
+    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
+)
+
+atol = 4.0e-8
+print(f"  atol={atol}")
+
+assert np.allclose(
+    [emittance_x, emittance_y, emittance_t],
+    [0.0, 0.0, 0.0],
+    atol=atol,
+)

examples/epac2004_benchmarks/input_fodo_rf_SC.in

@@ -0,0 +1,128 @@
+###############################################################################
+# Particle Beam(s)
+###############################################################################
+beam.npart = 10000
+beam.units = static
+beam.energy = 250.0
+beam.charge = 1.42857142857142865e-10
+beam.particle = proton
+beam.distribution = kurth6d
+beam.sigmaX = 9.84722273e-4
+beam.sigmaY = 6.96967278e-4
+beam.sigmaT = 4.486799242214e-03
+beam.sigmaPx = 0.0
+beam.sigmaPy = 0.0
+beam.sigmaPt = 0.0
+beam.muxpx = 0.0
+beam.muypy = 0.0
+beam.mutpt = 0.0
+
+
+###############################################################################
+# Beamline: lattice elements and segments
+###############################################################################
+lattice.elements = monitor fquad dr gapa1 dr dquad dr gapb1 dr fquad monitor
+
+monitor.type = beam_monitor
+monitor.backend = h5
+
+dr.type = drift
+dr.ds = 0.1
+dr.nslice = 4
+
+fquad.type = quad
+fquad.ds = 0.15
+fquad.k = 2.4669749766168163
+fquad.nslice = 6
+
+dquad.type = quad
+dquad.ds = 0.3
+dquad.k = -2.4669749766168163
+dquad.nslice = 12
+
+gapa1.type = rfcavity
+gapa1.ds = 1.0
+gapa1.escale = 0.042631556991578
+gapa1.freq = 7.0e8
+gapa1.phase = 45.0
+gapa1.mapsteps = 100
+gapa1.nslice = 10
+gapa1.cos_coefficients =           \
+                0.120864178375839  \
+                -0.044057987631337 \
+                -0.209107290958498 \
+                -0.019831226655815 \
+                0.290428111491964  \
+                0.381974267375227  \
+                0.276801212694382  \
+                0.148265085353012  \
+                0.068569351192205  \
+                0.0290155855315696 \
+                0.011281649986680  \
+                0.004108501632832  \
+                0.0014277644197320 \
+                0.000474212125404  \
+                0.000151675768439  \
+                0.000047031436898  \
+                0.000014154595193  \
+                4.154741658e-6     \
+                1.191423909e-6     \
+                3.348293360e-7     \
+                9.203061700e-8     \
+                2.515007200e-8     \
+                6.478108000e-9     \
+                1.912531000e-9     \
+                2.925600000e-10
+gapa1.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
+
+gapb1.type = rfcavity
+gapb1.ds = 1.0
+gapb1.escale = 0.042631556991578
+gapb1.freq = 7.0e8
+gapb1.phase = -1.0
+gapb1.mapsteps = 100
+gapb1.nslice = 10
+gapb1.cos_coefficients =            \
+                0.120864178375839   \
+                -0.044057987631337  \
+                -0.209107290958498  \
+                -0.019831226655815  \
+                0.290428111491964   \
+                0.381974267375227   \
+                0.276801212694382   \
+                0.148265085353012   \
+                0.068569351192205   \
+                0.0290155855315696  \
+                0.011281649986680   \
+                0.004108501632832   \
+                0.0014277644197320  \
+                0.000474212125404   \
+                0.000151675768439   \
+                0.000047031436898   \
+                0.000014154595193   \
+                4.154741658e-6      \
+                1.191423909e-6      \
+                3.348293360e-7      \
+                9.203061700e-8      \
+                2.515007200e-8      \
+                6.478108000e-9      \
+                1.912531000e-9      \
+                2.925600000e-10
+gapb1.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
+
+
+###############################################################################
+# Algorithms
+###############################################################################
+algo.particle_shape = 2
+algo.space_charge = true
+
+amr.n_cell = 56 56 64
+geometry.prob_relative = 4.0
+
+###############################################################################
+# Diagnostics
+###############################################################################
+diag.slice_step_diagnostics = true

examples/epac2004_benchmarks/input_thermal.in

@@ -0,0 +1,35 @@
+###############################################################################
+# Particle Beam(s)
+###############################################################################
+beam.npart = 10000
+#beam.npart = 2
+beam.units = static
+beam.energy = 2.0e3
+beam.charge = 1.0e-9
+beam.particle = proton
+beam.distribution = thermal
+beam.k = 1.0
+beam.kT = 1.0
+beam.kT_halo = 1.0
+beam.halo = 0.0
+
+###############################################################################
+# Beamline: lattice elements and segments
+###############################################################################
+lattice.elements = monitor
+
+monitor.type = beam_monitor
+monitor.backend = h5
+
+constf1.type = constf
+constf1.ds = 2.0
+constf1.kx = 1.0
+constf1.ky = 1.0
+constf1.kt = 1.0
+
+
+###############################################################################
+# Algorithms
+###############################################################################
+algo.particle_shape = 2
+algo.space_charge = false