Add mbtrack2 vs elegant benchmark of TBL

Simulation Files/elegant/mbtrack2_benchmark_transiant_beam_loading/benchmark_ele_save/cleanUp

+#!/bin/bash
+
+\rm track.w* twiss.param* *.twi* track.h* *.rec bunch beam track_h.w* track_h.h*
+

Simulation Files/elegant/mbtrack2_benchmark_transiant_beam_loading/benchmark_ele_save/makeBeam.ele

+&run_setup
+    lattice = "test.lte",
+    p_central_mev = 2.75e3,
+    use_beamline="no_line"
+    output = beam
+&end
+
+
+&run_control
+&end
+
+&sdds_beam
+        input = bunch,
+        track_pages_separately = 0,
+        use_bunched_mode = 1
+        n_duplicates = 405,
+        duplicate_stagger[4] = 2.842361460761627e-09,
+&end
+
+&track &end
+

Simulation Files/elegant/mbtrack2_benchmark_transiant_beam_loading/benchmark_ele_save/makeBunch.ele

+&run_setup
+    lattice = "test.lte",
+    p_central_mev = 2.75e3,
+    use_beamline="no_line"
+    output = bunch
+&end
+
+&twiss_output
+        matched = 1
+&end
+
+&run_control
+&end
+
+&bunched_beam
+        n_particles_per_bunch = 10000,
+        emit_x = 8.252315e-11,
+        beta_x = 1.0,
+        beta_y = 1.0,
+        sigma_dp = 8.823319e-04,
+        sigma_s = 2.7e-3,
+        distribution_type[0] = 3*"gaussian",
+        distribution_cutoff[0] = 3*3
+&end
+
+&track &end
+

Simulation Files/elegant/mbtrack2_benchmark_transiant_beam_loading/benchmark_ele_save/runScript

+#!/bin/bash
+ 
+# Transient beam loading benchmark with mbtrack2
+ 
+elegant twiss.ele
+ 
+elegant makeBunch.ele
+
+elegant makeBeam.ele
+ 
+# Track w/ RFmode
+time mpirun -n 80 Pelegant track.ele
+
+# Plotting of bunch length and energy spread
+sddsplot -layout=1,2 -join=x -col=Pass,Sdelta -col=Pass,St track.w1 -sep=1
+# Plot phase space
+sddsplot -col=dt,p -graph=dot -split=page -sep=page -same track.w2 -title=@Pass
+# Plot energy (delta) and time (dCt) centroid
+sddsplot -layout=1,2 -join=x -col=Pass,Cdelta -col=Pass,dCt track.w1 -sep=1
+# Plot longitudinal histogram
+sddsplot -column=dt,dtFrequency track.h1
+# plot passive cavity voltage
+sddsplot -col=Pass,V track.rec

Simulation Files/elegant/mbtrack2_benchmark_transiant_beam_loading/benchmark_ele_save/test.lte

+rf: rfca, VOLT=1.8e6, PHASE=166.55978566105205, FREQ=351820131.44377375, FIDUCIAL="light"
+rf_h: rfmode, Q=1e8,Rs=90e8,freq=1407420550.0705662,bin_size=1e-12,&
+        n_bins=2000,ramp_passes=0,record="%s.rec",sample_interval=10, PRELOAD=1, PRELOAD_CHARGE=5.912111838384184e-07, BUNCHED_BEAM_MODE=1
+
+ring1: ilmatrix, L=354.48130800001866, nux=0.2, nuy=0.3, BETAX=1, BETAY=1, ALPHAC=0.000105
+sr1: sreffects, JX=1.919076, JDELTA=1.080924, EXREF=8.252315e-11, SDELTAREF=8.823319e-04, DDELTAREF=-0.00015213639999999998, PREF=5380.615754828152
+c : CHARGE, TOTAL=5.912111838384184e-07
+r : RECIRC
+
+ring: line=(ring1,sr1,rf,rf_h)
+
+w1: watch,filename="%s.w1",mode="parameters",flush_interval=10, START_PID=1, END_PID=10000
+h1: histogram, filename="%s.h1", interval=1000, START_PID=1, END_PID=10000
+
+w2: watch,filename="%s.w2",mode="parameters",flush_interval=10, START_PID=1000001, END_PID=1010000
+h2: histogram, filename="%s.h2", interval=1000, START_PID=1000001, END_PID=1010000
+
+w3: watch,filename="%s.w3",mode="parameters",flush_interval=10, START_PID=2000001, END_PID=2010000
+h3: histogram, filename="%s.h3", interval=1000, START_PID=2000001, END_PID=2010000
+
+w4: watch,filename="%s.w4",mode="parameters",flush_interval=10, START_PID=3000001, END_PID=3010000
+h4: histogram, filename="%s.h4", interval=1000, START_PID=3000001, END_PID=3010000
+
+
+w5: watch,filename="%s.w5",mode="parameters",flush_interval=10, START_PID=4000001, END_PID=4010000
+h5: histogram, filename="%s.h5", interval=1000, START_PID=4000001, END_PID=4010000
+
+
+ring_h: line=(c, r, sr1,rf,rf_h,ring1,w1,w2,w3,w4,w5,h1,h2,h3,h4,h5)
+ring_no_h: line=(c, r, sr1,rf,ring1,w1,w2,h1)
+
+d0: drift,l=0
+no_line: line=(d0)

Simulation Files/elegant/mbtrack2_benchmark_transiant_beam_loading/benchmark_ele_save/track.ele

+&run_setup
+    lattice = "test.lte",
+    p_central_mev = 2.75e3,
+    use_beamline="ring_h"
+&end
+
+&twiss_output
+        matched = 1
+&end
+
+&run_control
+        n_passes = 20000
+&end
+
+&sdds_beam
+        input = beam
+	track_pages_separately = 0,
+        use_bunched_mode = 1
+&end
+
+&track &end

Simulation Files/elegant/mbtrack2_benchmark_transiant_beam_loading/benchmark_ele_save/twiss.ele

+! Computes Twiss parameters and related values for the
+! nominal PAR lattice.
+
+&run_setup
+    lattice = "test.lte",
+    p_central_mev = 2.75e3,
+    use_beamline="ring",
+    parameters = %s.param
+&end
+
+&twiss_output filename = "%s.twi",
+        statistics=1, radiation_integrals=1,
+&end
+
+&run_control &end
+&bunched_beam &end
+&track &end
+

Simulation Files/mbtrack2/Elegant_benchmark_transiant_beam_loading/elegant_benchmark_mp1e6.py

+# -*- coding: utf-8 -*-
+"""
+Created on Tue May  4 17:10:31 2021
+
+@author: gamelina
+"""
+import numpy as np
+from mbtrack2 import BeamLoadingEquilibrium
+from mbtrack2.tracking.rf import CavityResonator, RFCavity
+from machine_data import v0356
+from scipy.constants import c
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+
+plt.rcParams['figure.figsize'] = [12, 5]
+plt.rcParams['figure.dpi'] = 200
+
+mpi_switch = True
+
+Vc = 1.8e6
+ring = v0356()
+ring.U0 = 4.183751e-01*1e6
+ring.L = 354.48130800001866
+I0 = 0.5
+
+MC = RFCavity(ring, 1, Vc, np.arccos(ring.U0/Vc))
+
+m = 4
+Rs = 90e8
+Q = 1e8
+QL = 1e8
+detune = 140e3
+HHC = CavityResonator(ring, m, Rs, Q, QL,detune)
+HHC.Vg = 0
+HHC.theta_g = 0
+
+
+import numpy as np
+from machine_data import v0356, load_TDR_1_wf
+from mbtrack2.tracking import Beam, CavityResonator, LongitudinalMap, SynchrotronRadiation
+from mbtrack2.tracking.monitors import  BeamMonitor, CavityMonitor, ProfileMonitor
+from mbtrack2 import WakeField, WakePotential, WakePotentialMonitor
+
+
+long = LongitudinalMap(ring)
+rad = SynchrotronRadiation(ring)
+
+buffer = 100
+tot_turns = 40e3
+
+name = "benchmark_elegant"
+
+HCmon = CavityMonitor("HHC", ring, file_name=name, save_every=100,
+                 buffer_size=buffer, total_size=tot_turns/100, mpi_mode=mpi_switch)
+
+bbmon = BeamMonitor(h=ring.h, file_name=None, save_every=10, buffer_size=buffer,
+                 total_size=tot_turns/10, mpi_mode=mpi_switch)
+
+promon0 = ProfileMonitor(bunch_number=0, save_every=100, buffer_size=buffer, total_size=tot_turns/100, 
+                 dimensions="tau", n_bin=50, file_name=None, mpi_mode=mpi_switch)
+
+promon100 = ProfileMonitor(bunch_number=100, save_every=100, buffer_size=buffer, total_size=tot_turns/100, 
+                 dimensions="tau", n_bin=50, file_name=None, mpi_mode=mpi_switch)
+
+promon200 = ProfileMonitor(bunch_number=200, save_every=100, buffer_size=buffer, total_size=tot_turns/100, 
+                 dimensions="tau", n_bin=50, file_name=None, mpi_mode=mpi_switch)
+
+promon300 = ProfileMonitor(bunch_number=300, save_every=100, buffer_size=buffer, total_size=tot_turns/100, 
+                 dimensions="tau", n_bin=50, file_name=None, mpi_mode=mpi_switch)
+
+promon400 = ProfileMonitor(bunch_number=400, save_every=100, buffer_size=buffer, total_size=tot_turns/100, 
+                 dimensions="tau", n_bin=50, file_name=None, mpi_mode=mpi_switch)
+
+
+bb = Beam(ring)
+filling = np.ones(ring.h)
+filling[-10:] = 0
+filling = filling*I0/(ring.h-10)
+bb.init_beam(filling, mp_per_bunch=1e6, track_alive=False, mpi=mpi_switch)
+
+
+HHC.init_phasor(bb)
+
+for i in range(int(tot_turns+1)):
+    if (i % 1000 == 0):
+        if mpi_switch and (bb.mpi.rank == 0):
+            print(i)
+        else:
+            print(i)
+
+    long.track(bb)
+    rad.track(bb)
+
+    if mpi_switch:
+        bb.mpi.share_distributions(bb)
+
+    MC.track(bb)
+    HHC.track(bb)
+
+    bbmon.track(bb)
+    HCmon.track(bb, HHC)
+
+    promon0.track(bb)
+    promon100.track(bb)
+    promon200.track(bb)
+    promon300.track(bb)
+    promon400.track(bb)
+
+bbmon.close()
+
+