Simulation de tracking d'extraction lente¶

In [1]:
import matplotlib.pyplot as plt
import numpy as np
import xobjects as xo
import xplt
import xtrack as xt
from scipy.constants import c as clight

%config InlineBackend.figure_format = "retina"

Chargement de la maille PIMMS¶

In [2]:
# Chargement de la maille et des réglages optimisés des sextupoles d'extraction
env = xt.load("inputs/pimms.json")
env.vars.load_json("inputs/extraction_sextupoles.json")
ring = env.pimms
ring.configure_bend_model(edge="full", core="adaptive", num_multipole_kicks=5)
twiss = ring.twiss4d()

Création d'un paquet de particules¶

Pour cet exemple nous préparons manuellement les coordonées pour avoir, longitudinalement, un faisceau continu. Le paquet xpart expose des fonctionalités pour créer automatiquement différents types de faisceaux adaptés à la machine.

In [3]:
num_particles = 2500  # macroparticules
beam_intensity = 1e10  # protons

# Distributions Gaussiennes en espace des phases normalisés transverses
x_norm = np.random.normal(size=num_particles)
px_norm = np.random.normal(size=num_particles)
y_norm = np.random.normal(size=num_particles)
py_norm = np.random.normal(size=num_particles)

# Distribution Gaussienne des momenta longitudinaux (rms spread 5e-4)
# et temps d'arrivée des particules réparti sur un tour
delta = 5e-4 * np.random.normal(size=num_particles)
zeta = np.random.uniform(size=num_particles) * ring.get_length()

# Création de l'objet de particules
particles = ring.build_particles(
    x_norm=x_norm,
    px_norm=px_norm,
    y_norm=y_norm,
    py_norm=py_norm,
    delta=delta,
    zeta=zeta,
    method="4d",
    weight=beam_intensity / num_particles,
    nemitt_x=1.5e-6,
    nemitt_y=1e-6,
)

# Sauvegarde de la distribution initiale
p0 = particles.copy()
In [4]:
xplt.PhaseSpacePlot(particles)
plt.tight_layout()
No description has been provided for this image

Définition du comportement temporel des sextupôles d’extraction¶

In [5]:
# Il est possible de définir une fonction temporelle (ici linéaire)
ring.functions["func_sext"] = xt.FunctionPieceWiseLinear(x=[0, 0.1e-3, 0.5e-3], y=[0, 0, 1.0])
In [6]:
# On définit ensuite un lien de l'alimentation des sextupoles à cette fonction
ring.vars["kse1"] *= ring.functions["func_sext"](ring.vars["t_turn_s"])
ring.vars["kse2"] *= ring.functions["func_sext"](ring.vars["t_turn_s"])

# Inspectons l'expression
ring.vars["kse1"]._expr
Out[6]:
(4.9635634700966795 * f['func_sext'](vars['t_turn_s']))

Définition d'une limite d'ouverture réaliste du septum¶

In [7]:
ring["septum_aperture"].max_x = 0.035

Utilisation d'un contexte multi-threaded pour la performance¶

In [8]:
ring.discard_tracker()
context = xo.ContextCpu()
# context = xo.ContextCpu(omp_num_threads="auto")  # pour multi-thread
ring.build_tracker(_context=context)
Out[8]:
<xtrack.tracker.Tracker at 0x7ff96405f610>

Tracking des particules¶

In [9]:
def compute_intensity(line: xt.Line, particles: xt.Particles) -> float:
    """Return the bunch's intensity, in charges."""
    nplike = particles._context.nplike_lib
    particles.hide_lost_particles()
    intensity = particles.q0 * nplike.sum(particles.weight)
    particles.unhide_lost_particles()
    return intensity
In [10]:
# Definition de quantités à enregistrer pendant le tracking
log = xt.Log(
    "kse1",  # alimentation du premier sextupôle d'extraction
    "kse2",  # alimentation du second sextupôle d'extraction
    intensity=compute_intensity,  # fonction définie par l'utilisateur
)
In [11]:
# Activation de la mise à jour temporelle des variables pour le tracking
ring.enable_time_dependent_vars = True
In [12]:
# Track!
ring.track(particles, num_turns=5000, with_progress=25, log=log)

Visualisation des quantités enregistrées¶

In [13]:
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(5, 5))
ax1.plot(ring.log_last_track["kse1"], label="kse1")
ax1.plot(ring.log_last_track["kse2"], label="kse2")
ax1.set_ylabel("Alimentation sextupôle [$m^{-2}$]")

ax2.plot(np.array(ring.log_last_track["intensity"]) * 1e-10, c="C2", label="Intensité")
ax2.set_ylabel("Intensité [$10^{10}$ppb]")
ax2.set_xlabel("Tour")
ax2.set_ylim(bottom=0.0)

for axis in (ax1, ax2):
    axis.legend()
    axis.yaxis.set_major_locator(plt.MaxNLocator(5))
fig.align_ylabels((ax1, ax2))
plt.tight_layout()
plt.show()
No description has been provided for this image

Visualisation des particules après tracking¶

In [14]:
alive = particles.state > 0
lost = ~alive
In [15]:
fig, ax = plt.subplots()
ax.plot(particles.x[alive], particles.px[alive], ".", ms=1, c="green", label="En circulation")
ax.plot(particles.x[lost], particles.px[lost], ".", ms=1, c="red", label="Extraites")
ax.set(xlabel="$x$ [m]", ylabel="$p_x$")
plt.xlim(-0.03, 0.05)
plt.ylim(-3.5e-3, 3.5e-3)
ymin, ymax = ax.get_ylim()
ax.axvline(x=ring["septum_aperture"].max_x, color="k", alpha=0.4, linestyle="--")
ax.text(
    ring["septum_aperture"].max_x * 0.98,
    ymax * 0.95,
    "Septum",
    rotation=90,
    va="top",
    ha="right",
    alpha=0.5,
    c="k",
)
ax.xaxis.set_major_locator(plt.MaxNLocator(5))
ax.yaxis.set_major_locator(plt.MaxNLocator(5))
plt.legend()
plt.tight_layout()
plt.show()
No description has been provided for this image

Introduction d'une excitation transverse pour contrôler l'extraction¶

In [16]:
# Defintion d'un nouveau type d'élément
class SpillExcitation:
    def __init__(self):
        self.amplitude = 0
        self.tune = 0
        self.f_rev = 1 / twiss.T_rev0

    def track(self, particles):
        f_excit = self.tune * self.f_rev
        # Time of arrival corrected for delay within the turn
        t_particle_pass = (
            particles.at_turn[particles.state > 0] / self.f_rev
            - particles.zeta[particles.state > 0] / particles.beta0[0] / clight
        )
        particles.px[particles.state > 0] += self.amplitude * np.sin(2 * np.pi * f_excit * t_particle_pass)
In [17]:
# Installation de l'élement
ring.discard_tracker()
ring.insert_element("spill_exc", SpillExcitation(), at_s=0)
ring.build_tracker(_context=context)
Out[17]:
<xtrack.tracker.Tracker at 0x7ff95dfba350>
In [18]:
# Utilisation de l'amplitude et la fréquence de l'excitation pour contrôler l'extraction
# Definition de fonctions temporelles de ces propriétés
ring.functions["fun_excit"] = xt.FunctionPieceWiseLinear(x=[0, 0.5e-3, 0.7e-3, 4e-3], y=[0, 0, 1e-5, 3.2e-5])
ring.functions["fun_freq"] = xt.FunctionPieceWiseLinear(
    x=[0, 0.5e-3, 1.5e-3, 4e-3], y=[0.6615, 0.6615, 0.6617, 0.6622]
)
In [19]:
# On définit un lien de l'amplitude et la fréquence à ces fonctions
ring.vars["ampl_excit"] = ring.functions["fun_excit"](ring.vars["t_turn_s"])
ring.vars["freq_excit"] = ring.functions["fun_freq"](ring.vars["t_turn_s"])
In [20]:
ring.element_refs["spill_exc"].amplitude = ring.vars["ampl_excit"]
ring.element_refs["spill_exc"].tune = ring.vars["freq_excit"]

Remise à zéro du model temporel¶

In [21]:
ring.vars["t_turn_s"] = 0  # temps de la simulation
particles = p0.copy(_context=context)  # état initial des particules

Nouveau tracking des particules¶

In [22]:
# Definition de quantités à enregistrer pendant le tracking
log = xt.Log(
    "kse1",  # alimentation du premier sextupôle d'extraction
    "kse2",  # alimentation du second sextupôle d'extraction
    "t_turn_s",  # temps de simulation
    "ampl_excit",  # amplitude de l'excitation
    "freq_excit",  # fréquence de l'excitation
    intensity=compute_intensity,  # fonction définie par l'utilisateur
)
In [23]:
# Track!
ring.track(particles, num_turns=5000, with_progress=25, log=log)

Visualisation des quantités enregistrées¶

In [24]:
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, sharex=True, figsize=(5, 7))
t_ms = np.array(ring.log_last_track["t_turn_s"]) * 1e3  # temps en [ms]

ax1.plot(t_ms, ring.log_last_track["kse1"], label="kse1")
ax1.plot(t_ms, ring.log_last_track["kse2"], label="kse2")
ax1.set_ylabel("Alim. sext. [$m^{-2}$]")

ax2.plot(t_ms, np.array(ring.log_last_track["intensity"]) * 1e-10, c="C2", label="Intensité")
ax2.set_ylabel("Intensité [$10^{10}$ppb]")
ax2.set_ylim(bottom=0.0)

ax3.plot(t_ms, ring.log_last_track["ampl_excit"], c="C3", label="Amplitude")
ax3.set_ylim(bottom=0.0)
ax3.set_ylabel("Ampl. d'excitation")

ax4.plot(t_ms, ring.log_last_track["freq_excit"], c="C4", label="Fréquence")
ax4.set_ylabel("Fréq. d'excitation")
ax4.set_xlabel("Temps [ms]")

for axis in (ax1, ax2, ax3, ax4):
    axis.legend()
    axis.yaxis.set_major_locator(plt.MaxNLocator(5))

fig.align_ylabels((ax1, ax2, ax3, ax4))
plt.tight_layout()
plt.show()
No description has been provided for this image

Visualisation des particules après second tracking¶

In [25]:
alive = particles.state > 0
lost = ~alive
In [26]:
fig, ax = plt.subplots()
ax.plot(particles.x[alive], particles.px[alive], ".", ms=1, c="green", label="En circulation")
ax.plot(particles.x[lost], particles.px[lost], ".", ms=1, c="red", label="Extraites")
ax.set(xlabel="$x$ [m]", ylabel="$p_x$")
plt.xlim(-0.03, 0.05)
plt.ylim(-3.5e-3, 3.5e-3)
ymin, ymax = ax.get_ylim()
ax.axvline(x=ring["septum_aperture"].max_x, color="k", alpha=0.4, linestyle="--")
ax.text(
    ring["septum_aperture"].max_x * 0.98,
    ymax * 0.95,
    "Septum",
    rotation=90,
    va="top",
    ha="right",
    alpha=0.5,
    c="k",
)
ax.xaxis.set_major_locator(plt.MaxNLocator(5))
ax.yaxis.set_major_locator(plt.MaxNLocator(5))
plt.legend()
plt.tight_layout()
plt.show()
No description has been provided for this image