The eternal return of the Chinese boss

Global Sensitivity Analysis on Numerical Solver Parameters of Particle-In-Cell Models in Particle Accelerator Systems Matthias Frey∗, Andreas Adelmann arXiv:2001.10411v4 [physics.comp-ph] 2 Sep 2020 Paul Scherrer Institut, CH-5232 Villigen, Switzerland Abstract Every computer model depends on numerical input parameters that are chosen according to mostly conservative but rigorous numerical or empirical estimates. These parameters could for example be the step size for time integrators, a seed for pseudo-random number generators, a threshold or the number of grid points to discretize a computational domain. In case a numerical model is enhanced with new algorithms and modelling techniques the numerical influence on the quantities of interest, the running time as well as the accuracy is often initially unknown. Usually parameters are chosen on a trial-and-error basis neglecting the computational cost versus accuracy aspects. As a consequence the cost per simulation might be unnecessarily high which wastes computing resources. Hence, it is essential to identify the most critical numerical parameters and to analyze systematically their effect on the result in order to minimize the timeto-solution without losing significantly on accuracy. Relevant parameters are identified by global sensitivity studies where Sobol’ indices are common measures. These sensitivities are obtained by surrogate models based on polynomial chaos expansion. In this paper, we first introduce the general methods for uncertainty quantification. We then demonstrate their use on numerical solver parameters to reduce the computational costs and discuss further model improvements based on the sensitivity analysis. The sensitivities are evaluated for neighbouring bunch simulations of the existing PSI Injector II and PSI Ring as well as the proposed DAEδALUS Injector cyclotron and simulations of the rf electron gun of the Argonne Wakefield Accelerator. Keywords: adaptive mesh refinement, Particle-In-Cell, global sensitivity analysis, polynomial chaos expansion, particle accelerator, high intensity ∗ Corresponding author Email addresses: [email protected] (Matthias Frey), [email protected] (Andreas Adelmann) 1. Introduction Numerical models in scientific research disciplines are usually extremely complex and computationally intensive. A common feature to all is the dependency on numerical model parameters that do not represent an actual property of the underlying scientific problem. They are either prescribed in the source code and, thus, hidden to the user or can be chosen at runtime. Examples are the seed for pseudo random number generators, an error threshold in a convergence criterion, the number of grid points in mesh-based models or the step size in time integrators. Latter two are often chosen to satisfy memory constraints or time limits. When applying new algorithms and modelling techniques the sensitivity of such input values on the response is usually studied by varying only a single parameter. While this captures the main influence of the tested parameter, possible correlations with other parameters are missed. A remedy is the evaluation of Sobol’ indices [1] which are variance-based global sensitivity measures to express both individual and correlated parameter influences. Instead of Monte Carlo estimates, these quantifiers can easily be obtained by surrogate models based on polynomial chaos expansion (PCE). Uncertainty quantification (UQ) based on PCE is generally used in many areas of scientific computing and modelling [2, 3, 4, 5] and several frameworks exist such as [6, 7, 8]. The coefficients of the truncated PCE that are required to determine Sobol’ indices are mostly computed using the projection or regression method [9]. Since some numerical parameters are limited to integers, the former method is not applicable. In this paper we study the sensitivity of adaptive mesh refinement (AMR) and multi bunch [10, 11] parameters in Particle-In-Cell (PIC) simulations of high intensity cyclotrons where we use the new AMR capabilities of OPAL (Object Oriented Parallel Accelerator Library) [12] as presented in [13]. We further explore the sensitivity of a rf electron gun model with respect to the number of macro particles, the energy binning and the time step. The sensivitities are evaluated using ordinary least squares and Bayesian compressive sensing (BCS) [14, 15]. Both methods are part of the uncertainty quantification toolkit UQTk [16, 17] (version 3.0.4). The results are cross-checked using Chaospy [7] (version 3.0.5) together with the orthogonal matching pursuit [18, 19] regression model of the Python machine learning library scikit-learn [20, 21] (version 0.21.2). Although we apply UQ on numerical solver parameters, it is a general method used for example in [5] to evaluate the sensitivities and to predict the quantities of interest due to uncertainties in physical parameters. 2 The paper is organised as follows: In Sec. 2.1 we elaborate the physical applications and their numerical modelling. An introduction to UQ is given in Sec. 3. The experimental setup and its results are presented in Sec. 4 and Sec. 5, respectively. In Sec. 6 are the final conclusions. 2. Applications 2.1. High Intensity Cyclotrons Cyclotrons are circular machines that accelerate charged particles (e.g. protons) or ions (e.g. H2+ ). Depending on the particle species and delivered energy, these machines find different appli- cations ranging from isotope production [22, 23] and neutron spallation [24] to cancer treatment [25, 26]. An example that provides a beam for neutron spallation is the High Intensity Proton Accelerator (HIPA) facility at Paul Scherrer Institut (PSI) consisting of two cyclotrons, i.e. the PSI Injector II and PSI Ring. At a frequency of 50.65 MHz 10 mA DC (direct current) proton bunches at 870 keV are injected into PSI Injector II in which they are collimated to approximately 2.2 mA and accelerated up to 72 MeV (∼ 37 % speed of light), before being transported to the PSI Ring where they reach a kinetic energy of 590 MeV (∼ 79 % speed of light) at extraction. Another example is the planned facility of the DAEδALUS and IsoDAR (Isotope Decay At Rest) experiments for neutrino oscillation and CP violation. It consists of two cyclotrons where the DAEδALUS Injector Cyclotron (DIC) [27, 28] is the first acceleration stage delivering a 60 MeV/amu H2+ beam to the Superconducting Ring Cyclotron (DSRC) with extraction energy of 800 MeV/amu. Since cyclotrons are isochronous, i.e. the magnetic field is increased radially in order to keep the orbital frequency constant, i.e. the revolution time per turn is energy-independent and, thus, the bunches lie radially on axis. A sketch showing five bunches denoted as circles on adjacent turns is depicted in Fig. 1a. As shown in [11], a small turn separation causes interactions between neighbouring bunches which yields to more halo particles (cf. Fig. 1b). In order to resolve these effects, the open source beam dynamics code OPAL [12] got recently enhanced with AMR capabilities [13] which adds more complexity to the numerical model. The influence of the AMR solver parameter settings on the statistical measures of the particle bunches is yet unknown and a too conservative AMR regrid frequency worsens the time-to-solution. Furthermore, the applied energy binnnig technique [10] to fulfill the electrostatic assumption (cf. Sec. 2.1.1) increases the computational costs considerably. Hence, the goal of this study is to quantify the impact of AMR solver and energy binning parameters in order to improve further computational investigations of these bunch interactions. 3 r∼E E3 E4 E5 E2 y E1 x (a) Five bunches evolving radially on axis due to (b) Separation of a particle bunch into core (blue) the isochronism of the cyclotron. The origin of the and halo (red) particles. In the PSI Ring the over- coordinate system denotes the center of the ma- all loss is on the order of 10−4 which corresponds chine. The orbit radius r of each bunch is propor- to a beam intensity of about 0.2 µA, i.e. all par- tional to its energy E1 < E2 < · · · < E5 . ticles outside of approximately 3σ of a Gaussian distribution with standard deviation σ. Figure 1: Sketch of neighbouring bunches (left) in the context of isochronous cyclotrons and characterization of a single bunch (right). 2.1.1. Numerical Model The numerical model of neighbouring bunch simulations in OPAL, as presented in [11], is based on [10]. Due to the energy difference of the particle bunches on neighbouring turns a single transformation into the particle rest frame does not fully satisfy the requirements of the electrostatic assumption to solve Poisson’s equation. Instead, each of the N macro particles is assigned to an energy bin b due to its momentum βγ according to $
% sinh-1 (βγ) − sinh-1 mini={1,N } (βγ)i b= η (1) where the binning parameter η is a measure of the energy spread. In each time step the force on a particle exerted by all others is the sum of the electric field contributions of each energy bin b evaluated in the appropriate rest frame of the particles obtained by a Lorentz transform with the proper relativistic factor γb . The algorithm is summarised in Alg. 1. The computation of the electric field of an energy bin involves only particles of that bin, thus, the charge deposition applies only on a subset of particles M ⊂ {1, . . . , N }. However, the field on the mesh is interpolated and applied to all N particles. 4 Algorithm 1 Electrostatic Particle-In-Cell with B energy bins and N macro particles. 1: En ← 0 ∀n ∈ {1, . . . , N } // Electric field at particle location 2: for b ∈ {1, . . . , B} do // Loop over energy bins 3: x̃n ←LorentzTransform(xn , γb ) ∀n ∈ {1, . . . , N } 4: ρ̃i,j,k ←DepositCharge(x̃m , qm ) ∀m ∈ M ⊂ {1, . . . , N } 5: Ẽbi,j,k 6: Ẽbn ← GatherEfield(Ẽbi,j,k , x̃n ) 7: Ebn ←BackLorentzTransform(Ẽbn , γb ) 8: xn ←BackLorentzTransform(x̃n , γb ) 9: // Interpolate charge onto mesh ←PoissonSolve(ρ̃i,j,k ) En ← En + Ebn ∀n ∈ {1, . . . , N } // Get field at particle location ∀n ∈ {1, . . . , N } ∀n ∈ {1, . . . , N } ∀n ∈ {1, . . . , N } // Add field contribution 10: end for 2.1.2. RF Electron Gun Model To study the effect of energy binning we further use the example of the Argonne Wakefield Accelerator (AWA) [29, 30, 31] facility, an experiment setup for beam physics studies and accelerator technology developments. The facility is equipped with a photocathode rf electron gun that emits high intensity electron beams at high accelerating gradients (≫ 1 MV/m). Due to the high gradients the electrostatic approximation is invalidated and, hence, energy binning is necessary. In OPAL, we model the particle emission by s 2 Ekin ptotal = + 1 −1 mc2 px = ptotal sin (ϕ) cos (θ) py = ptotal sin (ϕ) sin (θ) pz = ptotal |cos (ϕ)| √
x ϕ = 2 cos−1 with x ∈ [0, 1] and θ ∈ [0, π] uniformly randomly sampled [12]. 2.2. Quantities of interest In accelerator physics interesting quantities of interest (QoI) also denoted as observables, of the co-moving frame are the rms beam size σω = p hω 2 i εω = p hω 2 ihp2ω i − hωpω i2 and the projected emittance ∀ω = x, y, z 5 (2) that describes the phase space volume per dimension. The bracket h·i represents the moment. In order to quantify halo (cf. Fig. 1b), i.e. the tails of a particle distribution, we use two statistical definitions for bunched beams by [32, 33], the spatial-profile parameter hω 4 i 15 − σω4 7 (3) √ p ω 15 3 I4 − 2 ε2ω 7 (4) hω = and the phase-space halo parameter Hω = with Eq. (2) and fourth order invariant (cf. also [34]) I4ω = ω 4 p4ω + 3 ω 2 p2ω 2 − 4 ωp3ω ω 3 pω . In case the bunch has uniform density Eq. (3) is zero due to the constant 15/7. An important quantity of interest in the rf electron gun model is the energy spread ∆E ∝ p hp2z i. (5) 3. Non-Intrusive Uncertainty Quantification Simulations of physical phenomena usually rely on measured input data. Depending on the accuracy of the measurement and the model sensitivity, the response may vary significantly. UQ introduces methods to quantify this variability in order to estimate the reliability of the obtained results. UQ distinguishes two approaches which are called intrusive and non-intrusive UQ. In contrast to intrusive UQ, non-intrusive UQ uses the computational model as a black box. In this paper we only give a short overview following the description and notation of [5, 9, 1, 35, 36]. A detailed introduction to UQ in general is found in e.g. [37, 38]. In Sec. 3.1, the surrogate model based on the polynomial chaos expansion (PCE) is explained in general. The sections 3.2 to 3.5 describe methods to obtain the coefficients of the expansion where special focus is given on the methods applied in this paper. The definition of the Sobol’ indices and their computation with the PCE is given in Sec. 3.6. In order to check the error bounds of the estimated sensitivities, we use the bootstrap method explained in Sec. 3.7. 3.1. Surrogate Model based on Polynomial Chaos Expansion The PC-decomposition originates from [39], where a random variable of a Gaussian distribution is represented as a series of multivariate Hermite polynomials of increasing order. And 6 as stated by the Cameron-Martin theorem [40], any functional in L2 space can be represented in a series of Fourier-Hermite functionals. Later, this method was rediscovered and applied by [41] to a stochastic process and then generalized by [42] to other probability measures and their corresponding orthogonal polynomials. Let a multivariate polynomial Ψαi (ξ) of dimension d ∈ N \ {0} and multiindex αi = (α1 , α2 , . . . , αd ) ∈ Nd be defined by Ψαi (ξ) = d Y ψαj (ξj ) j=1 with orthogonal univariate polynomials {ψαj }dj=1 . The response of a model m(x) with random input vector x ∈ Ω1 × · · · × Ωd , where Ωj ∀j = {1, . . . , d} denotes the sample space of the j-th random variable, can then be represented as m(x) = ∞ X cαi Ψαi (T (x)) (6) i=0 where the basis of a random input component is determined by its probability distribution (cf. Tab. 1) and T : x 7→ ξ denotes an isoprobabilistic transform. In case of dependent input components, for example, T represents the Rosenblatt transform [43] that yields independent random variables. Another method for dependent variables presented in [44] applies the GramSchmidt orthogonalization. Under the assumption of only independent input variables, the transform T reduces to a simple linear mapping of every component of x onto the defined interval of the corresponding univariate polynomial, e.g. ξj ∈ [−1, 1] for Legendre polynomials. In numerical computations the sum in Eq. (6) is truncated at some polynomial degree p, hence the expansion is only an approximation of the exact model m(x), i.e. m̂(x) = P −1 X cαi Ψαi (T (x)). (7) i=0 The truncation scheme is not clearly defined. A common rule, which is also used here, is the so-called total order truncation that keeps all multiindices α for which ||α||1 ≤ p. This yields a number of P = (p + d)! p!d! (8) multiindices. Three other schemes are explained in [35]. In the next sections we describe methods to compute the coefficients cαi of Eq. (7). 7 PDF of ξj Polynomial Basis Support {ψαj (ξj )} Ωj Gaussian Hermite ] − ∞, ∞[ Gamma Laguerre [0, ∞[ Uniform Legendre [a, b] with a, b ∈ R Table 1: Examples of the Wiener-Askey polynomial chaos of random variables ξj with appropriate probability density function (PDF) [42]. 3.2. Projection Method The (spectral) projection method computes the coefficients of Eq. (7) making use of the orthogonality of the basis functions, i.e. Ψαi (ξ)Ψαj (ξ) = 0 with ∀i 6= j. Thus, the PC coefficients are given by c αi = m(T −1 (ξ))Ψαi (ξ) . Ψ2αi (ξ) While the denominator is evaluated by analytic formulas (see examples in the appendix of [9]), the numerator is computed by Gaussian quadrature integration where N = (p + 1)d integration points, i.e. high fidelity model m(x) evaluations, are required. 3.3. Linear Regression Method The coefficients of Eq. (7) can also be computed with regression-based methods 1 ĉ = arg min 2 c N −1 X m(xj ) − j=0 P −1 X cαi Ψαi (ξ j ) i=0 ! 2 + λ1 2 ||c||2 + λ2 ||c||1 2 (9) 2 with regularization parameters λ1 , λ2 ≥ 0, and the l1 norm and l2 norm denoted by ||·||1 and ||·||2 , respectively. The minimization problem is called ordinary least squares if λ1 = λ2 = 0, elastic net [45] if λ1 , λ2 > 0, ridge regression [46] (or Tikhonov regularization) if only λ1 > 0 and Lasso [47] if only λ2 > 0. In matrix form the problem reads ĉ = arg min c 1 λ1 2 2 ||y − Ac||2 + ||c||2 + λ2 ||c||1 2 2 (10) with model response y = (m(x0 ), . . . , m(xN −1 )) ∈ RN ×1 , unknown coefficient vector c =
⊺ cα0 , . . . , cαp ∈ RP ×1 and system matrix A ∈ RN ×P . In case of λ2 = 0, the coefficients of ⊺ 8 Eq. (10) are obtained in closed form by ĉ = (A⊺ A + λ1 I) −1 A⊺ y with P × P identity matrix I. In contrast to the projection method (cf. Sec. 3.2), this method does not require a fixed number of samples N . However, [9] gives an empirical optimal training sample size of N = (d − 1)P (11) with P defined in Eq. (8). 3.4. Orthogonal Matching Pursuit The matching pursuit (MP) is a greedy algorithm developed by [18] which was enhanced by [19] to obtain better convergence. This improved method is called orthogonal MP (OMP). In terms of PCE, the algorithm searches a minimal set of non-zero coefficients to represent the model response, i.e. 2 ĉ = arg min ||y − Ac||2 c subject to ||c||0 ≤ Nc where ||c||0 denotes the number of non-zero coefficients in c with a user-defined maximum Nc [48]. The vectors and matrices are defined according to Eq. (10). It is an iterative procedure where in the (i + 1)-th step a new coefficient vector ci+1 is searched that maximizes the inner product to the current residual r i = y − y i . We refer to the given references for details. 3.5. Bayesian Compressive Sensing As stated in [14, 15], the linear regression model Eq. (9) can be interpreted in a Bayesian manner, i.e. p (c|D) = p (D|c) p (c) p (D) with posterior distribution p (c|D), likelihood p (D|c), prior p (c) and evidence p (D) of training N −1 data D = {x, y}j=0 [36]. The likelihood is assigned a Gaussian noise model   N −1 2 X (m (x ) − m̂ (x )) 1 i i  exp − p (D|c) = 2 N/2 2σ (2πσ 2 ) j=0 with variance σ 2 . It is a measure of how well the high fidelity model is represented by the surrogate model Eq. (7). In order to favour a sparse PCE solution, a Laplace prior !  P +1 P X λ |ci | exp −λ p (c) = 2 i=0 9 (12) is chosen. Using Eq. (12) in the maximum a posteriori (MAP) estimate for c, i.e. arg max log [p (D|c) p (c)] , (13) c the Bayesian approach is equivalent to Eq. (9) with λ1 = 0 [49], since Eq. (13) is identical to a minimization of arg min − log [p (D|c) p (c)] . c An iterative algorithm to obtain the coefficients is described in [15]. It requires a user-defined stopping threshold ε that basically controls the number of kept basis terms, with more being skipped the higher the value is. The overall method is known as Bayesian Compressive Sensing (BCS) [14, 15]. 3.6. Sensitivity Analysis Sobol’ indices [1] are good measures of sensitivity since they provide information about single and mixed parameter effects. In addition to these sensitivity measures, there are various other methodologies such as Morris screening [50]. A survey is presented in [51] on the example of a hydrological model. Instead of Monte Carlo, Sobol’ indices are also easily obtained by surrogate models based on PCE as discussed in the following subsections. 3.6.1. Sobol’ Sensitivity Indices In [1] Sobol’ proposed global sensitivity indices that are calculated on an analysis of variance (ANOVA) decomposition (Sobol’ decomposition) of a square integrable function f (x) with x ∈ I d := [0, 1]d , i.e. [52] f (x) = f0 + d X i=1 X fi (xi ) + 1≤i1 <···<is ≤d with mean f0 = and Z fi1 is (xi1 , xis ) + · · · + f12...d (x1 , x2 , . . . , xd ) Z (14) f (x)dx Id 1 fi1 ...is (xi1 , . . . , xis )dxk = 0, (15) 0 for k = i1 , . . . , is and s = 1, . . . , d. Since Eq. (15) holds, the components of Eq. (14) are mutually orthogonal. Therefore, the total variance of Eq. (14) is Z f 2 (x)dx − f02 D= Id 10 that can also be written as D= d X X Di + i=1 1≤i1 <···<is ≤d where Di1 ...is = Z Is Di1 is + · · · + D123...d (16) fi21 ...is (xi1 , . . . , xis )dxi1 · · · dxis , (17) with 1 ≤ i1 < · · · < is ≤ d. Based on Eq. (16) and Eq. (17), the Sobol’ indices are defined as Si1 ...is := with d X Si + i=1 X 1≤i<j≤d Di1 ...is D Sij + · · · + S12...d = 1. (18) The first order indices Si are also known as main sensitivities. They describe the effect of a single input parameter on the model response. The total effect of the i-th design variable on the model response, proposed by [53], is the sum of all Sobol’ indices that include the i-th index, i.e. SiT = X Si i∈I with I = {i = (i1 , . . . , is ) : ∃k, 1 ≤ k ≤ s ≤ d, ik = i}. 3.6.2. Sobol’ Indices using Polynomial Chaos Expansion Instead of Monte Carlo techniques, Sobol’ indices can be estimated using surrogate models based on PCE since the truncated expansion can be rearranged like Eq. (14). The Sobol’ estimates are then given by [9] Ŝi1 ...is = where Ii1 ,...,is = and variance    α: X 1 D̂ c2α Ψ2α α∈Ii1 ,...,is αk > 0 ∀k = 1, . . . , n, αk = 0 ∀k = 1, . . . , n, D̂ = P −1 X c2αi Ψ2αi . i=1 The main and total sensitivities are computed by Ŝi = 1 X D̂ α∈Ii 11 c2α Ψ2α  k ∈ (i1 , . . . , is ) k 6∈ (i , . . . , i ) 1 s with Ii = {α = (α1 , . . . , αd ) : αi > 0 ∧ ∀k 6= i, αk = 0} and ŜiT = 1 X D̂ c2α Ψ2α , α∈Ii with Ii = {α = (α1 , . . . , αd ) : αi > 0}, respectively. 3.7. Confidence Intervals using Bootstrap In this subsection we briefly outline the computation of confidence intervals for the estimates of Sobol’ indices using the bootstrap method [54]. In the context of PCE, the bootstrap method has already been applied in [55], where it is referred to as bootstrap-PCE (or bPCE). The bootstrap method, in general, generates B independent samples each of size N by resampling from the original dataset. Each bootstrap sample, that may contain a point several times, is then considered as a new training sample to compute the coefficients of Eq. (7). In order to calculate the 95 % confidence interval for Sobol’ indices we follow the description of [56], where the bounds are given by Ŝi1 ...is ± 1.96 · s.e.(Ŝi1 ...is ) with 1 ≤ s ≤ d and standard error (s.e.) of B ∈ N>1 bootstrap samples v u B  2 u 1 X (b) s.e.(Ŝi1 ...is ) = t Si1 ...is − Si∗1 ...is B−1 b=1 and bootstrap sample mean Si∗1 ...is = B 1 X (b) Si1 ...is . B b=1 4. Experiment Design 4.1. High Intensity Cyclotrons In order to study the effect of AMR solver parameters and energy binning in neighbouring bunch simulations we perform sensitivity experiments with three different high intensity cyclotrons, the PSI Ring [57], the PSI Injector II [58] and the DAEδALUS Injector Cyclotron (DIC) [27, 28]. We always accelerate 5 particle bunches with 106 particles. The coarsest level grid is kept constant with 243 mesh points which is refined twice. For the PSI Injector II and PSI Ring the particles are integrated in time over one turn using 2880 steps per turn and for the DIC over three turns with 1440 steps per turn. The experimental setup is summarized in Tab. 2. In all 12 no. turns steps/turn no. bunches particles/bunch PIC base grid no. AMR levels 1 or 3 1440 or 2880 5 106 24 × 24 × 24 2 Table 2: Experimental setup of the PSI Ring, PSI Injector II and DAEδALUS Injector Cyclotron model. experiments the initial particle distribution is read from a checkpoint file to guarantee identical conditions for all training and validation points of a UQ sample. A list of the design variables under consideration is given in Tab. 3. While the resolution is basically controlled by the maximum number of AMR levels, the refinement policy affects its location. As described in [13] the OPAL library provides several refinement criteria such as the charge density per grid point, the potential as well as the electric field. Here, we want to analyze the effect of the threshold λ ∈ [0, 1] of the electrostatic potential refinement policy, where a grid cell (i, j, k) on a level l is refined if |φli,j,k | ≥ λ max |φli,j,k | i,j,k holds. Due to the motion of the particles in space the multi-level hierarchy has to be updated regularly to maintain the resolution which is defined by the regrid frequency fr . It should be noted that the regrid frequency defines the number of steps until the AMR hierarchy is updated. Hence, if fr = 1, the AMR levels are updated in each time step. Whenever this happens, the electric self-field needs to be recalculated by solving Poisson’s equation. The number of Poisson solves is controlled by the number of energy bins and therefore by the binning parameter η (cf. Sec. 2.1.1). The lower the value of η, the smaller the bin width and, hence, the more expensive the model is. As an upper limit of the regrid frequency fr we choose 120 integration steps. Since we perform either 1440 or 2880 steps per turn, this corresponds to an azimuthal angle of 30◦ and 15◦ , respectively. The choice of the binning parameter η in Eq. (1) depends mainly on the energy difference between bunches. The upper bound of the sampling range was selected such that we have at most as many energy bins as bunches in simulation. However, the sampling from the range is not straightforward since there exist more states with fewer energy bins (cf. Fig. 2) due to Eq. (1). Instead, we sample the binning parameter in subranges of equal bin count in order to avoid a biased sample set. 13 number of bins 5 PSI Ring PSI Injector II DIC 4 3 2 1 0.005 0.010 0.015 η 0.020 0.025 Figure 2: Number of energy bins in the PSI Ring, PSI Injector II and DAEδALUS Injector Cyclotron (DIC) with respect to the binning parameter η. The shown binning curves are with respect to the initial simulation energies used in this study. A straightforward uniform sampling from the full range yields a biased sample set. symbol design variable sampling range fr regrid frequency [1, 120] λ refinement threshold [0.5, 0.9] binning Ring (10−3 ) [4.7, 5.7] ∪ [5.8, 7.6] ∪ [7.7, 11.5] ∪ [11.6, 23.0] ∪ [23.1, 27.1] η binning Inj-2 (10 −3 ) binning DIC (10−3 ) [4.0, 4.9] ∪ [5.0, 6.5] ∪ [6.6, 9.8] ∪ [9.9, 19.7] ∪ [19.8, 23.8] [4.1, 5.0] ∪ [5.1, 6.6] ∪ [6.7, 10.0] ∪ [10.1, 20.0] ∪ [20.1, 24.1] Table 3: List of design variables and their sampling ranges for the neighbouring bunch simulations. 4.2. RF Electron Gun Model Like in the neighbouring bunch model, the time-to-solution in the rf electron gun model is dominated by the Poisson solver and the time integration. A reduction of the computational effort with regard to the Poisson solver is achieved by smaller PIC meshes and fewer energy bins NE . The costs of the time integrator is cheapened with coarser time steps ∆t and fewer macro particles. Instead of the AMR model, the rf electron gun model uses the Fast Fourier Transform (FFT) Poisson solver of OPAL where we put a Lx × Ly × Lz uniform mesh of Lx = Ly = 64 and Lz = 32 grid points. The final number of emitted macro particles is given by N p = pf L x L y L z where the particle multiplication factor pf is an integer. This parameter basically controls the number of particles per grid cell and, hence, the noise of the PIC model. The design variables and sampling ranges are given in Tab. 4. We model the rf electron gun of the Argonne Wakefield Accelerator (AWA) that has a length of approximately 30 cm. 14 symbol design variable sampling range pf particle factor [1, 5] NE number of bins [2, 10] ∆t time step (0.1 ps) [1, 10] Table 4: List of design variables and their sampling ranges for the rf electron gun model. The time step is the only floating point variable. 4.3. Surrogate Model Selection In order to avoid overfitting we proceed like [35] where the truncation order of the PC expansion and the settings of the regression models are chosen such that the relative l2 error v u PN −1 2 u t i=0P[m(xi ) − m̂(xi )] (19) N −1 2 i=0 m (xi ) between the surrogate m̂(x) and high fidelity model m(x) of the training and validation set is approximately of equal magnitude. As an additional error measure we also compare the relative l1 error PN −1 |m(xi ) − m̂(xi )| . PN −1 i=0 |m(xi )| i=0 (20) The number of samples N in Eq. (19) and Eq. (20) corresponds either to the number of training Nt or validation points Nv . The total number of N = 100 samples was randomly partitioned into disjoint training and validation sets with Nt = 0.5N and Nv = 0.5N , respectively. Since we have d = 3 design variables, we satisfy Eq. (11) with Nt = 50 up to polynomial order p = 3. 5. Results The estimated sensitivities are obtained from PC surrogate models where we use either ordinary least squares (OLS) and Bayesian compressive sensing (BCS) of UQTk [16, 17] or orthogonal matching pursuit (OMP) of scikit-learn [20, 21] together with Chaospy [7] to compute the expansion coefficients. A summary of the PC model setups is given in Tab. 5. In order to study the evolution of the sensitivities we construct the PC surrogate models at equidistant steps of the accelerator models and evaluate their sensitivities. These steps correspond to azimuthal angles in the cyclotrons or longitudinal positions in the rf electron gun model. In the examples below we only show the first order Sobol’ indices since their sum is already almost one which is the maximum per definition (cf. Eq. (18)). 15 model PC order stopping criterion ε (BCS) Nc (OMP) PSI Ring 2 1 × 10−7 5 PSI Injector II 2 1 × 10−4 7 DIC 2 1 × 10−8 6 AWA 2 1 × 10−9 7 Table 5: PC surrogate model settings for all accelerator model examples. The stopping criterion of Bayesian Compressive Sensing (BCS) and Orthogonal Matching Pursuit (OMP) is discussed in Sec. 3.5 and Sec. 3.4, respectively. 5.1. High Intensity Cyclotrons In order to study the effect of the input parameters we evaluate the sensitivities of the halo parameters Eq. (3) and Eq. (4) with respect to the center bunch of the 5 bunches (cf. Fig. 1a). The initial kinetic energy of the center bunch in the different models is approximately 98 MeV, 25 MeV and 17 MeV for the PSI Ring, PSI Injector II and DIC, respectively. As shown in Fig. 4, Fig. 8 and Fig. 12, the relative l1 and l2 errors (cf. Eq. (20) and Eq. (19)) between training and test samples are in good agreement for all cyclotron examples. A similar observation is done at a single angle in Fig. 5, Fig. 9 and Fig. 13. The average errors are given in Tab. 7, Tab. 8 and Tab. 9. The computation methods OLS, BCS and OMP yield similar results. In case of the PSI Injector II, the refinement threshold has more than 80 % impact on the halo. The energy binning parameter η and regrid frequency fr play a negligible role. The increase of the 95 % bootstrap confidence intervals in Fig. 6 correlates with the decrease of the standard deviation in Fig. 3. It is best observed for hx at around 215◦ or Hx between 195◦ and 255◦ . In contrast to the PSI Injector II, the DIC also strongly depends on the regrid frequency. It has an average main sensitivity of approximately 60 % for hx . The parameters also exhibit more correlations as observed between the main and total sensitivities (cf. Tab. 8). The standard deviations for the DIC are one order of magnitude smaller than for the PSI Injector II, causing the confidence intervals to increase as illustrated in Fig. 10. This effect is even stronger in the PSI Ring where Coulomb’s repulsion is less dominant and the halo parameters are smaller (cf. Fig. 11) compared to the PSI Injector II. The standard deviation is in the order of O(10−4 ) denoting no significant influence of the input parameters on the model response, hence, the confidence intervals in Fig. 14 exhibit large 16 ranges. For this reason we can make no reliable statement about the sensitivities for the PSI Ring. Nevertheless, these findings give rise to computational savings. Due to the small deviations, it is sufficient for the PSI Ring to select a cheap model. According to Tab. 6, the cheapest model among all N = 100 samples is 2.47 times faster than the most expensive model. model design variables PSI Ring PSI Injector II DIC time [s] fr λ η 111 0.8272 0.0227 7938 3 0.5022 0.0052 19 613 4 0.6455 0.0224 10 526 3 0.5022 0.0045 21 422 90 0.8584 0.0157 9560 74 0.5958 0.0216 31 709 Table 6: Most expensive and cheapest cyclotron models with respect to runtime among all N = 100 samples. mean hx Hx hy Hy 101 std 10−2 10−4 50 100 150 200 250 azimuth [deg] 300 350 400 Figure 3: Evolution of the mean and standard deviation (std) of the spatial-profile parameters hx , hy and the phase-space halo parameters Hx , Hy as defined in Eq. (3) and Eq. (4), respectively. Based on the mean of Hx , the location of the dipoles in the PSI Injector II can be detected, i.e. at 90◦ , 180◦ , 270◦ and 360◦ . 17 relative l2 error [%] relative l1 error [%] hx Hx hy Hy 10−1 10−3 10−1 10−3 50 100 150 200 250 azimuth [deg] 300 350 400 Figure 4: Evolution of the relative l2 and l1 error between the surrogate and the true model of the PSI Injector II. The full lines are the errors to the surrogate model obtained with the training set and the dashed lines are the errors to the surrogate model computed with the validation set. For each quantity, the dashed and full lines are close to each other, indicating no overfitting of the surrogate model. 18 y=x training points validation points 3.4 Ĥx ĥx 3.6 3.3 3.2 3.4 3.1 3.2 3.4 3.4 hx 3.6 Hx 5.0 8.5 ĥy Ĥy 4.9 4.8 8.0 4.7 4.6 8.0 8.5 4.6 4.8 5.0 Hy hy Figure 5: Comparison between the high fidelity (x-axis) and PC surrogate model (y-axis) at 390◦ of the PSI Injector II simulation. The blue and red dots indicate the training and validation points, respectively. In the best case all points coincide with the dashed black line. 19 Ŝ[hx ] fr λ η 1.0 0.5 Ŝ[Hx ] 0.0 1.0 0.5 Ŝ[hy ] 0.0 1.0 0.5 Ŝ[Hy ] 0.0 1.0 0.5 0.0 50 100 150 200 250 azimuth [deg] 300 350 400 Figure 6: Evolution of the estimated first order Sobol’ indices in the PSI Injector II. The error bars denote the 95 % bootstrapped (B = 100) confidence interval (cf. Sec. 3.7). The main sensitivities are evaluated for the spatial-profile parameters hx , hy and the phase-space halo parameters Hx , Hy as defined in Eq. (3) and Eq. (4), respectively. The bars are coloured with respect to the regrid frequency fr , AMR refinement threshold λ and energy binning parameter η. The refinement threshold has the highest impact on the halo measures. 20 QoI hx Hx hy Hy method l1 error [%] l2 error [%] Sobol’ sensitivity indices train test train test Ŝfr ŜfTr Ŝλ ŜλT Ŝη ŜηT OLS 0.23 0.33 0.33 0.44 0.03 0.06 0.91 0.94 0.02 0.04 BCS 0.26 0.30 0.38 0.44 0.03 0.05 0.92 0.95 0.02 0.04 OMP 0.24 0.32 0.34 0.45 0.03 0.04 0.92 0.95 0.02 0.04 OLS 0.15 0.22 0.21 0.29 0.08 0.16 0.80 0.89 0.01 0.05 BCS 0.16 0.20 0.24 0.28 0.07 0.15 0.81 0.90 0.01 0.05 OMP 0.15 0.21 0.22 0.30 0.07 0.16 0.80 0.90 0.01 0.05 OLS 0.18 0.27 0.25 0.36 0.06 0.10 0.88 0.92 0.01 0.04 BCS 0.20 0.24 0.30 0.35 0.05 0.09 0.88 0.92 0.01 0.04 OMP 0.19 0.26 0.27 0.37 0.05 0.08 0.89 0.92 0.01 0.04 OLS 0.11 0.17 0.15 0.22 0.06 0.10 0.88 0.92 0.01 0.03 BCS 0.12 0.15 0.17 0.21 0.06 0.10 0.88 0.92 0.01 0.03 OMP 0.11 0.16 0.16 0.22 0.05 0.09 0.89 0.93 0.01 0.04 Table 7: Average relative l1 and l2 errors between the high fidelity model and the PC surrogate models for the training and validation sets as well as the average main and total sensitivities for the PSI Injector II. OLS: Ordinary Least Squares; BCS: Bayesian Compressive Sensing; OMP: Orthogonal Matching Pursuit. 21 hx Hx hy Hy mean 6 × 100 4 × 100 3 × 100 std 10−3 10−4 200 400 600 azimuth [deg] 800 1000 Figure 7: Evolution of the mean and standard deviation (std) of the spatial-profile parameters hx , hy and the phase-space halo parameters Hx , Hy as defined in Eq. (3) and Eq. (4), respectively. The mean of all quantities shows a more or less periodic pattern along the three turns of the DAEδALUS Injector Cyclotron. 22 relative l1 error [%] 10−2 relative l2 error [%] hx 10−2 Hx hy Hy 10−3 10−3 200 400 600 azimuth [deg] 800 1000 Figure 8: Evolution of the relative l2 and l1 error between the surrogate and the true model of the DAEδALUS Injector Cyclotron. The full lines are the errors to the surrogate model obtained with the training set and the dashed lines are the errors to the surrogate model computed with the validation set. For each quantity, the dashed and full lines are close to each other, indicating no overfitting of the surrogate model. 23 y=x training points 4.640 7.391 4.639 ĥx Ĥx 7.392 validation points 7.390 4.638 7.389 4.637 7.390 7.392 4.638 hx 4.640 Hx 2.2910 3.917 2.2900 Ĥy ĥy 2.2905 3.916 2.2895 3.915 2.2890 2.289 2.290 2.291 3.915 3.916 3.917 Hy hy Figure 9: Comparison between the high fidelity (x-axis) and PC surrogate model (y-axis) at 1120◦ of the DAEδALUS Injector Cyclotron simulation. The blue and red dots indicate the training and validation points, respectively. In the best case all points coincide with the dashed black line. 24 Ŝ[hx ] fr λ η 1.0 0.5 Ŝ[Hx ] 0.0 1.0 0.5 Ŝ[hy ] 0.0 1.0 0.5 Ŝ[Hy ] 0.0 1.0 0.5 0.0 200 400 600 azimuth [deg] 800 1000 Figure 10: Evolution of the estimated first order Sobol’ indices in the DAEδALUS Injector Cyclotron. The error bars denote the 95 % bootstrapped (B = 100) confidence interval (cf. Sec. 3.7). The main sensitivities are evaluated for the spatial-profile parameters hx , hy and the phase-space halo parameters Hx , Hy as defined in Eq. (3) and Eq. (4), respectively. The bars are coloured with respect to the regrid frequency fr , AMR refinement threshold λ and energy binning parameter η. Beside the refinement threshold, the regrid frequency has also a significant impact on the halo measures. 25 QoI hx Hx hy Hy l1 error [10−2 %] l2 error [10−2 %] train test train test Ŝfr ŜfTr Ŝλ ŜλT Ŝη ŜηT OLS 0.29 0.45 0.38 0.61 0.59 0.71 0.27 0.33 0.02 0.09 BCS 0.30 0.46 0.39 0.61 0.60 0.72 0.25 0.33 0.02 0.09 OMP 0.30 0.45 0.39 0.60 0.62 0.72 0.26 0.32 0.01 0.07 OLS 0.21 0.31 0.28 0.41 0.29 0.40 0.57 0.67 0.02 0.06 BCS 0.23 0.32 0.30 0.42 0.29 0.40 0.57 0.67 0.02 0.05 OMP 0.23 0.31 0.31 0.40 0.25 0.34 0.64 0.73 0.01 0.02 OLS 0.34 0.55 0.44 0.72 0.11 0.24 0.67 0.76 0.07 0.15 BCS 0.35 0.55 0.45 0.71 0.11 0.24 0.67 0.76 0.07 0.15 OMP 0.36 0.54 0.46 0.69 0.10 0.22 0.71 0.77 0.06 0.13 OLS 0.24 0.37 0.31 0.47 0.47 0.58 0.39 0.48 0.02 0.07 BCS 0.25 0.38 0.32 0.49 0.47 0.58 0.39 0.48 0.02 0.07 OMP 0.25 0.37 0.33 0.47 0.44 0.53 0.46 0.53 0.01 0.04 method Sobol’ sensitivity indices Table 8: Average relative l1 and l2 errors between the high fidelity model and the PC surrogate models for the training and validation sets as well as the average main and total sensitivities for the DAEδALUS Injector Cyclotron. OLS: Ordinary Least Squares; BCS: Bayesian Compressive Sensing; OMP: Orthogonal Matching Pursuit. 26 hx Hx hy Hy mean 3 × 100 2 × 100 100 std 10−4 10−5 10−6 10−7 100 150 200 250 300 azimuth [deg] 350 400 450 Figure 11: Evolution of the mean and standard deviation (std) of the spatial-profile parameters hx , hy and the phase-space halo parameters Hx , Hy as defined in Eq. (3) and Eq. (4), respectively. The variability of these quantities in the PSI Ring cyclotron is on the order of O(10−4 ) which is two orders of magnitude smaller than for the PSI Injector II (cf. Fig. 3). 27 relative l2 error [%] relative l1 error [%] hx Hx hy Hy 10−2 10−3 10−4 10−2 10−3 10−4 150 200 250 300 azimuth [deg] 350 400 450 Figure 12: Evolution of the relative l2 and l1 error between the surrogate and the true model of the PSI Ring cyclotron. The full lines are the errors to the surrogate model obtained with the training set and the dashed lines are the errors to the surrogate model computed with the validation set. For each quantity, the dashed and full lines are close to each other, indicating no overfitting of the surrogate model. 28 y=x training points validation points ×10−3 + 1.757 1.8780 1.8779 Ĥx ĥx 0.5 0.0 1.8778 1.8778 1.8779 1.8780 0.0 hx 4 ×10−4 + 1.4 ×10−4 + 2.024 4 2 3 0 Ĥy ĥy 0.5 ×10−3 + 1.757 Hx 2 −2 1 −4 5 10 ×10−4 + 1.399 hy 2 4 ×10−4 + 2.024 Hy Figure 13: Comparison between the high fidelity (x-axis) and PC surrogate model (y-axis) at 471◦ of the PSI Ring simulation. The blue and red dots indicate the training and validation points, respectively. In the best case all points coincide with the dashed black line. 29 Ŝ[hx ] fr λ η 1.0 0.5 Ŝ[Hx ] 0.0 1.0 0.5 Ŝ[hy ] 0.0 1.0 0.5 Ŝ[Hy ] 0.0 1.0 0.5 0.0 150 200 250 300 azimuth [deg] 350 400 450 Figure 14: Evolution of the estimated first order Sobol’ indices in the PSI Ring. The error bars denote the 95 % bootstrapped (B = 100) confidence interval (cf. Sec. 3.7). The main sensitivities are evaluated for the spatial-profile parameters hx , hy and the phase-space halo parameters Hx , Hy as defined in Eq. (3) and Eq. (4), respectively. The bars are coloured with respect to the regrid frequency fr , AMR refinement threshold λ and energy binning parameter η. Due to the high uncertainty of the sensitivities, no reliable conclusion can be drawn for this machine. 30 QoI hx Hx hy Hy l1 error [10−2 %] l2 error [10−2 %] train test train test Ŝfr ŜfTr Ŝλ ŜλT Ŝη ŜηT OLS 0.85 0.78 1.13 1.01 0.15 0.19 0.50 0.54 0.27 0.34 BCS 0.92 0.80 1.24 1.06 0.15 0.19 0.50 0.54 0.27 0.34 OMP 0.86 0.78 1.14 1.01 0.14 0.16 0.50 0.51 0.34 0.36 OLS 0.47 0.45 0.62 0.58 0.11 0.16 0.59 0.62 0.22 0.29 BCS 0.60 0.50 0.79 0.65 0.11 0.16 0.59 0.62 0.22 0.29 OMP 0.48 0.45 0.64 0.57 0.07 0.13 0.63 0.64 0.24 0.30 OLS 0.80 0.76 1.05 0.98 0.10 0.19 0.60 0.63 0.20 0.28 BCS 0.96 0.80 1.28 1.06 0.10 0.19 0.60 0.63 0.20 0.28 OMP 0.81 0.77 1.07 0.98 0.09 0.17 0.59 0.61 0.24 0.31 OLS 0.42 0.39 0.55 0.50 0.13 0.22 0.51 0.55 0.25 0.34 BCS 0.47 0.39 0.62 0.49 0.13 0.23 0.50 0.55 0.25 0.34 OMP 0.43 0.40 0.56 0.50 0.10 0.23 0.52 0.54 0.24 0.37 method Sobol’ sensitivity indices Table 9: Average relative l1 and l2 errors between the high fidelity model and the PC surrogate models for the training and validation sets as well as the average main and total sensitivities for the PSI Ring. OLS: Ordinary Least Squares; BCS: Bayesian Compressive Sensing; OMP: Orthogonal Matching Pursuit. 31 5.2. RF Electron Gun Model In order to approximate the high fidelity model we use PC surrogate models of second order where the BCS method uses a tolerance of ε = 10−9 and the OMP method is stopped once 7 non-zero coefficients are found. In Fig. 16 are the relative l2 and l1 errors evaluated along the rf electron gun model. The mean errors are summarized in Tab. 10. It shows that the l1 and l2 errors on the test and training points match with an absolute difference of O(10−2 ) and O(10−1 ), respectively. An example of a comparison between the PC surrogate and high fidelity model is illustrated in Fig. 17. The first order Sobol’ indices and their 95 % bootstrapped confidence intervals are illustrated in Fig. 18. Except to the sensitivities of the horizontal projected emittance εx , we observe a convergence of the model parameter influences. The energy spread ∆E and the projected emittance εs strongly depend on the time step (Ŝ[εs ], Ŝ[∆E] > 0.90). This high influence is due to the momentum component in their definitions (cf. Eq. (5) and Eq. (2)) and the fact that the smaller the time step, the better the process of acceleration (i.e. the evolution of the momentum) is resolved. The rms beam size in longitudinal direction is dominated by the energy binning (Ŝ[NE ] ≈ 0.45) and the particle multiplication factor (Ŝ[pf ] ≈ 0.41). While a higher pf value improves the statistics of the beam size and reduces the numerical noise of PIC, the energy binning is coupled with Coulomb’s repulsion that affects the beam size. In transverse direction, NE and ∆t are important instead. The convergence of the relative errors is correlated with the convergence of the variances of the quantities of interest as observed in Fig. 15. The model might therefore be improved with an adaptive time stepping scheme that addresses this effect. The cheapest AWA rf electron gun model, i.e. ∆t = 1 ps, NE = 2 and pf = 1, is 15 times faster than the most expensive model which has ∆t = 0.1 ps, NE = 10 and pf = 5. 32 QoI σx εx σs εs ∆E method l1 error [%] l2 error [%] Sobol’ sensitivity indices train test train test Ŝ∆t T Ŝ∆t Ŝpf ŜpTf ŜNE T ŜN E OLS 0.03 0.04 0.04 0.05 0.32 0.32 0.12 0.12 0.56 0.56 BCS 0.04 0.04 0.05 0.05 0.32 0.32 0.12 0.12 0.56 0.56 OMP 0.03 0.04 0.04 0.05 0.32 0.32 0.12 0.12 0.55 0.56 OLS 0.13 0.16 0.16 0.20 0.43 0.44 0.09 0.10 0.46 0.47 BCS 0.14 0.15 0.17 0.19 0.43 0.44 0.09 0.10 0.46 0.47 OMP 0.13 0.16 0.17 0.20 0.44 0.45 0.10 0.10 0.45 0.46 OLS 0.02 0.02 0.03 0.03 0.13 0.14 0.41 0.41 0.45 0.46 BCS 0.02 0.03 0.03 0.03 0.13 0.14 0.41 0.41 0.45 0.46 OMP 0.02 0.02 0.03 0.03 0.13 0.13 0.39 0.39 0.48 0.48 OLS 0.69 0.62 1.04 0.93 0.95 0.96 0.01 0.01 0.04 0.05 BCS 0.93 0.86 1.22 1.12 0.95 0.96 0.01 0.01 0.04 0.04 OMP 0.69 0.63 1.04 0.95 0.95 0.96 0.01 0.01 0.04 0.04 OLS 0.11 0.12 0.16 0.17 0.94 0.94 0.01 0.01 0.05 0.05 BCS 0.14 0.15 0.18 0.19 0.94 0.94 0.01 0.01 0.05 0.05 OMP 0.11 0.12 0.16 0.17 0.94 0.94 0.01 0.01 0.05 0.05 Table 10: Average relative l1 and l2 errors between the high fidelity model and the PC surrogate models for the training and validation sets as well as the average main and total sensitivities for the rf electron gun model of the AWA. OLS: Ordinary Least Squares; BCS: Bayesian Compressive Sensing; OMP: Orthogonal Matching Pursuit. 33 σx εx σs εs ∆E 15 20 longitudinal position [cm] 25 mean 100 10−1 10−2 10−2 std 10−3 10−4 10−5 0 5 10 30 Figure 15: Evolution of the mean and standard deviation (std) of the energy spread ∆E, the projected emittances εx , εs and the rms beam sizes σx , σs for rf electron gun model of the AWA. 34 relative l2 error [%] relative l1 error [%] σx εx σs εs ∆E 15 20 longitudinal position [cm] 25 100 10−1 10−2 100 10−1 0 5 10 30 Figure 16: Evolution of the relative l2 and l1 error between the surrogate and the true model of the AWA rf electron gun. The full lines are the errors to the surrogate model obtained with the training set and the dashed lines are the errors to the surrogate model computed with the validation set. For each quantity, the dashed and full lines are close to each other, indicating no overfitting of the surrogate model. 35 y=x training points 2.052 2.18 σ̂s [mm] 7.4 2.19 ε̂x [mm] σ̂x [mm] ×10 validation points −2 7.3 2.050 2.048 7.2 2.046 2.18 2.19 7.2 σx [mm] 2.5 4.8 2.0475 2.0500 2.0525 σs [mm] ×10−2 2.4 d [MeV] ∆E ε̂s [mm] ×10−2 εx [mm] 7.4 ×10−2 2.3 4.6 2.2 4.6 4.8 ×10−2 εs [mm] 2.2 2.4 ×10−2 ∆E [MeV] Figure 17: Comparison between the high fidelity (x-axis) and PC surrogate model (y-axis) at the exit of the rf electron gun model of the AWA, i.e. s ≈ 30 cm. The blue and red dots indicate the training and validation points, respectively. In the best case all points coincide with the dashed black line. 36 pf Ŝ[σx ] ∆t NE 1.0 0.5 Ŝ[εx ] 0.0 1.0 0.5 Ŝ[σs ] 0.0 1.0 0.5 Ŝ[εs ] 0.0 1.0 0.5 Ŝ[∆E] 0.0 1.0 0.5 0.0 0 5 10 15 20 longitudinal position [cm] 25 30 Figure 18: Evolution of the estimated first order Sobol’ indices in the rf electron gun model of the AWA. The error bars denote the 95 % bootstrapped (B = 100) confidence interval (cf. Sec. 3.7). The main sensitivities are evaluated for the energy spread ∆E, the projected emittances εx , εs and the rms beam sizes σx , σs . The bars are coloured with respect to the time step ∆t, particle multiplication factor pf and the number of energy bins NE . The impact of the parameters on the quantities in longitudinal direction converges. 6. Conclusions In this paper we discussed uncertainty quantification based on polynomial chaos expansion and gave a brief introduction to four numerical methods to compute the polynomial coefficients. The choice of the method depends on the problem and its dimension. While the projection method is the most accurate, the number of high-fidelity evaluations grows exponentially with the dimension which is not the case for the other presented methods. Bayesian compressive sensing and matching orthogonal pursuit favour sparse solutions by the selection of the most important contributions. The least squares method solves a linear system which may fail in case 37 the matrix is ill-conditioned which happens when input dimensions can only take a few discrete values. This can be the case with integer input, for example, when specifying the number of AMR levels or the grid sizes for the Poisson solver. Beside a cheap surrogate model that mimics the high fidelity model, polynomial chaos based uncertainty quantification has the additional benefit to easily evaluate the Sobol’ sensitivity indices. As demonstrated in this paper, this technique is not only suitable to gain knowledge about the sensitivity of physical parameters (e.g. initial beam properties) on the quantities of interest but also numerical parameters of computer codes. Since some tested numerical parameters might be limited to integers, the projection method to obtain the polynomial coefficients is, however, not applicable. Instead, regression-based methods, Bayesian Compressive Sensing or Orthogonal Matching Pursuit and others have to be applied. A further difficulty with numerical parameters is a fair random sampling. In some cases (cf. Fig. 2) a straightforward, uniform sampling of the parameter yields biased input data and, hence, may induce wrong conclusions. To circumvent, we perform a stratified sampling that guarantees a well-balanced distribution. The sensitivity studies of the three high intensity cyclotrons show that the sensitivity results are different among accelerators of the same type. While the AMR threshold is the most important parameter in the PSI Injector II with a sensitivity of about 90 %, the regrid frequency is relevant in the DAEδALUS Injectory Cyclotron (DIC), too. Large bootstrap confidence intervals for the Sobol’ indices indicate a failure of the analysis since the contributed variation of the model response is rather due to noise than the tested input parameters. In such a case no reliable statements based on the sensitivity estimates can be done. In contrast to our intuition the standard deviation of the halo parameters remain pretty constant throughout one turn in the PSI Ring and PSI Injector II and the considered three turns in the DIC. Nevertheless, these findings give rise to computational savings. Without losing significantly on accuracy (cf. Fig. 3, Fig. 7 and Fig. 11), energy binning can be totally switched off for these cyclotrons. In addition, this reduces the amount of AMR hierarchy updates which reduces the time-to-solution even further since the operators of the adaptive multigrid solver to solve Poisson’s equation do not need to be set up in every time step. To illustrate this, we take the benchmark example in [13] that solves Poisson’s equation 100 times using a three level AMR hierarchy with a base level of 5763 grid points. The benchmark running on 14 400 CPU (Central Processing Unit) cores shows that the matrix setup due to AMR regriding takes up 42.15 % computing time. A reduction of the regrid frequency by a factor 10 yields a speedup of 7.10 in the matrix setup timing. In our UQ samples, the speedup between the cheapest and most expensive model is at least 2.0 and at most 3.3. 38 Another interesting case we have studied is the rf electron gun model of the AWA. Relevant parameters for this model are the energy binning NE and the time step ∆t. The particle multiplication factor pf , that basically controls the number of particles per grid cell, is only important for the longitudinal beam size. Although NE and pf have together an average main sensitivity of 86 %, ∆t is the dominating parameter close to the cathode. An adaptive energy binning and time step scheme that is based on Sobol’ sensitivity indices is therefore a possible future enhancement to reduce the time-to-solution for a target accuracy. 7. Acknowledgments The authors thank the developers of UQTk, especially B. Debusschere and K. Sargsyan. The computer resources were provided by the Swiss National Supercomputing Centre (CSCS), Paul Scherrer Institut (PSI) and the Laboratory Computing Resource Center at the Argonne National Laboratory (ANL). This project is funded by the Swiss National Science Foundation (SNSF) under contract number 200021 159936. References [1] Sobol’, I. M., Mathematics and Computers in Simulation 55 (2001) 271, The Second IMACS Seminar on Monte Carlo Methods. [2] Najm, H. N., Annual Review of Fluid Mechanics 41 (2009) 35. [3] Reagan, M. T. et al., Combustion Theory and Modelling 8 (2004) 607. [4] Eck, V. G. et al., International Journal for Numerical Methods in Biomedical Engineering 32 (2016) e02755. [5] Adelmann, A., SIAM/ASA Journal on Uncertainty Quantification 7 (2019) 383. [6] Marelli, S. and Sudret, B., UQLab: A Framework for Uncertainty Quantification in Matlab, in Vulnerability, Uncertainty, and Risk, Second International Conference on Vulnerability and Risk Analysis and Management (ICVRAM) and the Sixth International Symposium on Uncertainty, Modeling, and Analysis (ISUMA), pages 2554–2563, 2014. [7] Feinberg, J. and Langtangen, H. P., Journal of Computational Science 11 (2015) 46 . 39 [8] Adams, B. M. et al., Dakota, a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantification, and sensitivity analysis: Version 6.0 users manual, 2014, Sandia Technical Report SAND2014-4633, Updated November 2015 (Version 6.3). [9] Sudret, B., Reliability Engineering & System Safety 93 (2008) 964 , Bayesian Networks in Dependability. [10] Fubiani, G., Qiang, J., Esarey, E., Leemans, W. P., and Dugan, G., Phys. Rev. ST Accel. Beams 9 (2006) 064402. [11] Yang, J. J., Adelmann, A., Humbel, M., Seidel, M., and Zhang, T. J., Phys. Rev. ST Accel. Beams 13 (2010) 064201. [12] Adelmann, A. et al., arXiv e-prints (2019) arXiv:1905.06654. [13] Frey, M., Adelmann, A., and Locans, U., Computer Physics Communications (2019) 106912. [14] Ji, S., Xue, Y., and Carin, L., IEEE Transactions on Signal Processing 56 (2008) 2346. [15] Babacan, S. D., Molina, R., and Katsaggelos, A. K., IEEE Transactions on Image Processing 19 (2010) 53. [16] Debusschere, B. et al., SIAM Journal on Scientific Computing 26 (2004) 698. [17] Debusschere, B., Sargsyan, K., Safta, C., and Chowdhary, K., Uncertainty Quantification Toolkit (UQTk), pages 1807–1827, Springer International Publishing, Cham, 2017. [18] Mallat, S. G. and Zhifeng Zhang, IEEE Transactions on Signal Processing 41 (1993) 3397. [19] Pati, Y. C., Rezaiifar, R., and Krishnaprasad, P. S., Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition, in Proceedings of 27th Asilomar Conference on Signals, Systems and Computers, pages 40–44 vol.1, 1993. [20] Pedregosa, F. et al., Journal of Machine Learning Research 12 (2011) 2825. [21] Buitinck, L. et al., API design for machine learning software: experiences from the scikitlearn project, in ECML PKDD Workshop: Languages for Data Mining and Machine Learning, pages 108–122, 2013. 40 [22] Morrissey, D. J., Sherrill, B. M., Steiner, M., Stolz, A., and Wiedenhoever, I., Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 204 (2003) 90 , 14th International Conference on Electromagnetic Isotope Separators and Techniques Related to their Applications. [23] Alonso, J. R., Barlow, R., Conrad, J. M., and Waites, L. H., Nature Reviews Physics 1 (2019) 533. [24] Fischer, W. E., Physica B: Condensed Matter 234-236 (1997) 1202 , Proceedings of the First European Conference on Neutron Scattering. [25] McCarthy, D. W. et al., Nuclear Medicine and Biology 24 (1997) 35 . [26] Weber, D. C. et al., International Journal of Radiation Oncology Biology Physics 63 (2005) 401 . [27] Alonso, J. R. and DAEδALUS Collaboration, AIP Conference Proceedings 1525 (2013) 480. [28] Adelmann, A. and DAEδALUS Collaboration, AIP Conference Proceedings 1560 (2013) 709. [29] Gai, W., Power, J. G., and Jing, C., Journal of Plasma Physics 78 (2012) 339345. [30] Jing, C. et al., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 898 (2018) 72 . [31] M. E. Conde et al., Research Program and Recent Results at the Argonne Wakefield Accelerator Facility (AWA), in Proc. of International Particle Accelerator Conference (IPAC’17), Copenhagen, Denmark, 14-19 May, 2017, number 8 in International Particle Accelerator Conference, pages 2885–2887, Geneva, Switzerland, 2017, JACoW. [32] Wangler, T. P. and Crandall, K. R., Beam halo in proton linac beams, Number 20 in International Linac Conference, 2000. [33] Allen, C. K. and Wangler, T. P., Phys. Rev. ST Accel. Beams 5 (2002) 124202. [34] Lysenko, W. P. and Overley, M. S., AIP Conference Proceedings 177 (1988) 323. [35] Sargsyan, K. et al., International Journal for Uncertainty Quantification 4 (2014) 63. 41 [36] Sargsyan, K., Surrogate Models for Uncertainty Propagation and Sensitivity Analysis, pages 673–698, Springer International Publishing, Cham, 2017. [37] Ghanem, R., Higdon, D., and Owhadi, H., editors, Handbook of Uncertainty Quantification, Springer International Publishing, Cham, 2017. [38] Sullivan, T. J., Introduction to Uncertainty Quantification, Springer International Publishing, Cham, 2015. [39] Wiener, N., American Journal of Mathematics 60 (1938) 897. [40] Cameron, R. H. and Martin, W. T., Annals of Mathematics 48 (1947) 385. [41] Spanos, P. D. and Ghanem, R., Journal of Engineering Mechanics 115 (1989) 1035. [42] Xiu, D. and Karniadakis, G., SIAM Journal on Scientific Computing 24 (2002) 619. [43] Rosenblatt, M., The Annals of Mathematical Statistics 23 (1952) 470. [44] Jakeman, J. D., Franzelin, F., Narayan, A., Eldred, M., and Pflüger, D., Computer Methods in Applied Mechanics and Engineering 351 (2019) 643 . [45] Zou, H. and Hastie, T., Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67 (2005) 301. [46] Hoerl, A. E. and Kennard, R. W., Technometrics 12 (1970) 55. [47] Tibshirani, R., Journal of the Royal Statistical Society: Series B (Methodological) 58 (1996) 267. [48] Rubinstein, R., Zibulevsky, M., and Elad, M., Efficient Implementation of the K-SVD Algorithm using Batch Orthogonal Matching Pursuit, Technical report, CS Technion, 2008. [49] Figueiredo, M., Adaptive sparseness using jeffreys prior, in Advances in Neural Information Processing Systems 14, edited by Dietterich, T. G., Becker, S., and Ghahramani, Z., pages 697–704, MIT Press, 2002. [50] Morris, M. D., Technometrics 33 (1991) 161. [51] Gan, Y. et al., Environmental Modelling & Software 51 (2014) 269 . [52] Sobol’, I. M., Matematicheskoe Modelirovanie 2 (1990) 112, (in Russian), English version in: Mathematical Modelling and Computational Experiments, 1(4):407-414, 1993. 42 [53] Homma, T. and Saltelli, A., Reliability Engineering & System Safety 52 (1996) 1 . [54] Efron, B., Ann. Statist. 7 (1979) 1. [55] Marelli, S. and Sudret, B., Structural Safety 75 (2018) 67 . [56] Archer, G. E. B., Saltelli, A., and Sobol’, I. M., Journal of Statistical Computation and Simulation 58 (1997) 99. [57] Frey, M., Snuverink, J., Baumgarten, C., and Adelmann, A., Phys. Rev. Accel. Beams 22 (2019) 064602. [58] Kolano, A., Adelmann, A., Barlow, R., and Baumgarten, C., Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 885 (2018) 54 . 43