J-PARC Main Ring ⼤強度運転時における磁場リップルの影響~GPU を使った粒⼦トラッキングシミュレーション〜
⾼エネルギー加速器研究機構
J-PARC Main Ring 電磁⽯・電源グループ
栗本 佳典

⽬次
• ⼤強度陽⼦加速器J-PARC MRの性能と電磁⽯電源• GPUを使った粒⼦トラッキングシミュレーション• 結果と考察Introduction

J-PARC 30 GeV シンクロトロン (Main Ring)
J-PARC MR••

周⻑

1568 m

⼊射エネルギー

3 GeV

取出しエネルギー

30 GeV

バンチ数

8

ビーム取出し周期2.48 s（速取）, 5.20 s（遅取）ビーム強度

500 kW (速取), 50 kW（遅取）2.48s 速い取り出し（ニュートリノ）
5.2s 遅い取り出し（原⼦核ハドロン）• さらなる⼤強度化、ビームロス低減（どちらも）• 取り出しビーム平坦化（遅い取り出し）4電源電流リップルと遅い取り出し
Phys. Rev. Accel. Beams 22, 072802 (2019)取り出しスピルモニターの信号
チューン変動＠30 GeV (20-200 Hz)

内藤⼤地さんらの論⽂

⾚：電流値による予測
⿊：コヒーレント振動
検出による測定

⾚：電流値による補正あり
⿊：補正無し• Spill Dutyが 3 %から31 %に改善

チューン（ベータトロンチューン）：四極磁場できまる横⽅向振動数/1ターン• チューン変動が電源電流から予測される値と⼀致した（200 Hz）• 電流値によるチューン変動をFPGAを使ってリアルタイムで予測しチューンを補正した。
• これにより遅い取り出しビームの平坦性が改善することも分かった。
より⼤強度の速い取り出し運転での電源リップルの影響は？ 空間電荷効果の⼤きい
⼊射エネルギー付近での粒⼦トラッキングシミュレーションが必要GPUを使った粒⼦トラッキングシミュ
レーション

GPUの構成
Block 1

Block 2Block M

Thread1

Thread1

Thread1

SharedMemoryThread2…
ThreadN

Thread2

SharedMemory…

…
ThreadN

SharedMemory…

Thread2

ThreadN

Global Memory••••ブロックとスレッドの多層構造各スレッドに演算が割り当てられる。並列計算が可能
各ブロックにはブロック内スレッドがアクセスできる共有メモリがある。⾼速だが⼩容量48〜96KB
CUDA C/C++ というC++ベースのプログラミング⾔語を使⽤する。

粒⼦トラッキングシミュレーションの概要
粒⼦⽣成

Drift, Magnets, Cavityなど

空間電荷効果単粒⼦⼒学：ある点からある点への写像は並列に計算可能

PIC シミュレーション : 写像には電荷の空間分布が必要→完全に並列にはできない。

単粒⼦⼒学
ハミルトニアン𝐻 𝑥, 𝑝 , 𝑦, 𝑝 , 𝜎, 𝑝 ; 𝑠𝑝𝑝, 𝜎𝑠

1

𝛽 𝑡,ℎ𝑥𝛿1

𝛿ℎ𝑥ℎ𝑥𝑝𝑝𝑝𝑝,𝑒𝑝,𝑒𝑝 ,ハミルトン⽅程式が解けない場合は 𝐻 𝐻
𝑉 （ 𝐻 および𝑉単体では解ける）とハミルトニアンを分割して、以下のシンプレクティック変換を⾏う。（ 𝐿はコンポーネントの⻑さ）

𝑒/𝑒𝑒/または𝑒𝑒𝑒𝑒𝑒𝑎

実装済み：偏向、四極、六極、薄い多極成分、ドリフト

112

13,𝑏1,𝑐213

空間電荷効果
𝜕𝜕𝑥

𝜕
𝜕𝑦

𝜕
𝜑 𝑥, 𝑦, 𝑧𝜕𝑧𝜌 𝑥, 𝑦, 𝑧𝜀1. 電荷分布 𝑓 𝑥, 𝑦 𝑎𝑛𝑑 𝑔 𝑧 の作成

𝑔 𝑧
2D 近似
𝜎 ,𝜎 ≪ 𝜎𝜕𝜕𝑥𝜕𝜕𝑦𝑢 𝑥, 𝑦𝑔 𝑧 𝑓 𝑥, 𝑦𝜀（⻑いビーム）

もっとも時間コストが⾼い演算（後述）2. 境界条件aを含む⼆次元ポアソン⽅程式で𝑢 𝑥, 𝑦 を得る
3. 勾配,,𝐴𝑔 𝑧 𝑢 𝑥, 𝑦

𝐿を計算しキック𝐴𝐿

移動距離

a. ダクトの形に応じて極座標と直交座標⼆種類のソルバーを実装。ソルバーはFFTと⾏列演算中⼼なので、GPUの並列計算が有効。

共有メモリを⽤いた電荷分布の作成
Block NBlock 1𝑥 ,𝑦Thread

𝑥 ,𝑦Thread

𝑥 ,𝑦Thread
AtomicAdd

Shared Memory…

Subhistogram 1

Thread𝑥,𝑦Thread𝑥,𝑦Thread
AtomicAdd

Shared Memory
Subhistogram N

…

…

𝑥 ,𝑦

Global Memory=Subhistogram 1+…

Total
Histogram+

Subhistogram N• 各ブロックの共有メモリにSub-histogramを割り当て、複数のヒストグラムを並列に計算し、あとで和をとる。
• 各スレッドが平⾏に共有メモリに電荷を⾜していくので、他をブロックする和AtomicAdd 関数を使う。（CUDA LIBRARYにある。）

J-PARC MR の場合のシミュレーション

シミュレーション条件
速い取出しモードの運転条件
••
•
•
•
•••••

横⽅向チューン：(x, y) = (21.35, 21.45)
バンチ内粒⼦数：3x1013 ppbマクロ粒⼦数：200000
加速電圧： 基本波 160 kV ⼆倍⾼調波 110 kV
⼊射時バンチファクタ Iaverage/Ipeak：0.2
⼊射時横⽅向2 エミッタンス ： 16 mm・mrad
アパーチャー：60 mm・mrad
コンポーネント数：約1600 / turn
空間電荷効果によるキック：1000ヵ所/turn
エラーなしと電源リップルのみで⽐較

2700 秒/10000 turns using NVIDIA TESLA V100
MRで使⽤しているのシミュレーション(CPU)の10倍以上の速度電流リップル（時間領域）• 磁場リップルとして、実際に測定された電源電流のリップルを使⽤した。
• ダクトで⾼周波磁場がシールドされる効果（森⽥, 2015年会）をフィル電流リップル（周波数領域） タで模擬

シミュレーション結果
チューン分布

ビーム⽣存率
磁⽯による設定チューン（21.35,21.45）〜1 %チューンの広がり = 0.3-0.4〜0.1 %

バンチファクタIaverage/IpeakSimulation

実測⼊射時0.2 1000ターン⽬以降~0.3（⼆倍⾼調波電圧の効果）

電源リップルの効果はビーム⽣存率で0.1 %程度だが⼀⽅でエラー無しです
でに1 %程度のビームロスがある。
こちらの調査を優先すべき

損失粒⼦のチューン分布
損失直前のチューン分布（10 turn 平均）• 損失したマクロ粒⼦のIDを記録しておいて、もう⼀度同じプログラムを⾛らせ、軌道をターン毎に記録した。（エラーなしのとき）損失する粒⼦は、⾮構造共鳴2x-4y=-43と構造共鳴付近2x-2y = 0 にいることが分かった。
• シミュレーションに8極や12極の外部磁場は⼊れていないため、空間電荷分布がつくる共鳴である• ⾮構造共鳴があるので、3回対称性を破るダクトの配置（菱形近似しているUFOダクト）をやめてみたがが、結果は不変。空間電荷効果の計算箇所も3回対称。

損失粒⼦の詳細な運動

チューンx,yの時間変化

作⽤Jx,Jyの時間変化

x-2y = 0• 作⽤の時間変化に明確な結合共鳴の兆候がある。Jx とJyが同じ量だけ、交換している。
• 構造共鳴のため、共鳴幅は広いと考えられる。
• 設定チューンを（21.35,21.45）を右側（21.45,21.45）にずらすと、損失は⼤幅に増え、ほとんどがこの共鳴で損失するようになる。

チューンダイアグラム上の軌跡

縦⽅向位置の時間変化

損失粒⼦の詳細な運動

チューンx,yの時間変化

作⽤Jx,Jyの時間変化

2x-4y=-43• 作⽤の時間変化に明確な共鳴の兆候はな
く、じわじわエミッタンスが増えている。
• すべての粒⼦について似たような傾向。
しかし、図はJx増、Jy減だが、逆の場合やどちらも増⼤の場合もある。
• ⾮構造共鳴のため、共鳴幅が狭く、短期間に起こる少量のエミッタンス変化を繰り返しているのか？
• なぜエラー無し⾮構造共鳴による損失が起きているのか、さらに調査が必要

チューンダイアグラム上の軌跡

縦⽅向位置の時間変化

まとめ• 電源による磁場リップルがJ-PARC MR⼤強度運転にどの程度影響するかをシミュレーションで調べた。
• GPUで動く空間電荷効果を含んだトラッキングシミュレーションコードを開発し、現⾏MRが使っているのCPUのものより、10倍以上速くなった。
• シミュレーションの結果、磁場リップルの存在はビーム損失にして0.1 %しか寄与しないことが分かった。
• しかし、バンチ内陽⼦の電荷分布⾃⾝がつくる共鳴2x-4y=-43および1 %の損失 が存在する。

2x-2y = 0による

ただし、現時点ではあくまで⼀つのシミュレーションの結果に過ぎない。今後、より詳細な解析や、現存するトラッキングシミュレーションとの⽐較で、本結果の信頼度を⾼めていきたい。
Thanks to KEK⼤⾒さん、原研：發知さん、⽥村さん （アドバイス）, NVIDIA (GPU贈与)

Backupダクト径のチューニング
今回の空間電荷によるポテンシャル計算
𝑢/𝛾𝑢を 𝑢𝑢/𝛾で近似

コヒーレントチューンシフト測定
A. Tobayashi IPAC 2019 MOPGW036

虚電荷の効果を⼩さく⾒積もっている
本来は、ダクト境界の解𝑢両⽅が必要。→
無限解 𝑢
本シミュレーション（1 bunch）によるコヒーレントチューンシフト

コヒーレントチューンシフト測定による131 mmチューニング100 mmBMダクト⾼さ100mm

BMダクト⾼さ60mm

𝑢と今後実装予定。

BMダクト60 mm
今回は、コヒーレントチューンシフトの測定がシミュレーションと合うように、偏向磁⽯⽤ダクト⾼さを変更した。

Particle Tracking with Space
Charge Effect using GPU
Yoshi Kurimoto
J-PARC MR Power Converter GroupMotivation (Questions)• What is the source of the beam loss ?
• How do rippling magnetic fields affect the beam loss ?
• Can we simulate the beam faster using GPU?

GPU
Block 1

Block 2Block M

Thread1

Thread1

Thread1

SharedMemoryThread2…
ThreadN
Global Memory••••

Each operation can be assigned to each thread.
Execution of each thread can be parallelly done.
Threads in a common block can access the shared memory.
Shared memory is limited (12288 float words) but very fastThread2…

…
ThreadN

SharedMemory…

Thread2

ThreadN

SharedMemory

Example𝑎𝑥

𝑎 𝑥𝑎 𝑥⋯𝑎 𝑥Matrix multiplication, Fourier transformation …
1. Each term can be calculated in parallel
Parallela1

x1

×a2×x2aN…

2. Summation also can be partially parallel.
ParallelNot Parallel

a1x1a2x2

a3x3

a4x4

a5x5

a6x6

a7x7

a8x8++++a1x1 +a2x2

a3x3 +a4x4

a5x5 +a6x6

a7x7 +a8x8++×xN

Overview of the particle tracking simulation
Particle Generation

Drift, Magnets, Cavity

Space Charge EffectSingle Particle Mechanics : Maps between two locations can be calculated in parallel

PIC simulation : Maps requires spatial distributions of particles

Single Particle Mechanics

Hamiltonian
𝐻 𝑥, 𝑝 , 𝑦, 𝑝 , 𝜎, 𝑝 ; 𝑠𝑝1ℎ𝑥

1ℎ𝑥𝑝, 𝜎𝑠𝛽 𝑡,𝛿𝑝𝛿𝑝ℎ𝑥𝑝𝑝,𝑒𝑝𝑒,𝑝 ,Then just have to solve the Hamilton equations. But not always solvableSymplectic Map
If a Hamilton equation about 𝐻 can not be solved, The Hamiltonian is divided
𝑉 so that 𝑒and 𝑒 can be exactly
into two components like 𝐻 𝐻
calculated (Hamilton Equations for 𝐻 and 𝑉 can be solved ) Then, a map
between two components (Length : L) is approximately calculated as𝑒/𝑒/𝑒or𝑎𝑒𝑒𝑒𝑒𝑒11213,𝑏1,𝑐213

Bending Magnets𝑒𝐴 𝑥

dispersion Weak focus𝐻ℎ𝑥𝑝ℎ 𝑥𝑝hs
𝑉• 𝑒• 𝑒𝑝 𝑠1ℎ 𝑥2slippage

𝐻𝑥 𝑠ℎ𝑥𝑝 0𝑝 0sin ℎ𝑠
1 cos ℎ𝑠ℎℎ𝑝 0 cos ℎ𝑠 ℎ𝑥 0 sin ℎ𝑠 𝑝 0 sin ℎ𝑠𝑥 0 cos ℎ𝑠

𝑝 0𝑝 𝑠𝑝 0 s 𝑥 0𝑦 𝑠
𝑝 0
𝜎 0𝑥 0 sin ℎ𝑠
cos ℎ𝑠𝜎 𝑠ℎ𝑝 0sin ℎ𝑠
𝑝 0 𝑠
ℎ𝑝 𝑠
𝑝 01𝑥 𝑠𝑝 0 𝑝 0 𝑠𝑥 0𝑦 𝑠

𝑦 0𝑝 0𝑝 𝑠

𝑝 0 𝑝 0 𝑠𝑝 0
𝜎 0𝛾𝑝 0𝑝 𝑠𝑝 0𝑝 𝑠𝑝 0𝜎 𝑠𝑝 02𝑠

Quad Magnets
𝐻𝑘 𝑥𝑒𝐴 𝑥𝑦𝑥 𝑠𝑥 0 cos 𝑘 𝑠

𝑝 𝑠𝑝 0 cos 𝑘 𝑠

𝑝 0𝑘𝑉
• 𝑒

sin 𝑘 𝑠𝑝 𝑠

𝑘 𝑥 0 sin 𝑘 𝑠
𝑝 0𝑦 0 cosh 𝑘 𝑠sinh 𝑘 𝑠
𝑘𝑝 0 cosh 𝑘 𝑠
𝑘 𝑦 0 sinh 𝑘 𝑠

𝜎 𝑠

𝜎 0𝑝 𝑠𝑝 0𝑦 𝑠𝑦h𝑝

𝐻• 𝑒

1𝑘 𝑥2𝑥 𝑠

𝑝 0 𝑝 0 𝑠

𝑥 0𝑦 𝑠𝑦 0𝑝 0𝑝 𝑠

𝑝 0 𝑝 0 𝑠𝑝 0
𝜎 0𝛾𝑝 0𝑝 𝑠𝑝 0𝑝 𝑠𝑝 0𝜎 𝑠𝑝 02𝑠

Drift Space𝑒𝐴 𝑥

𝐻𝑝𝑝2𝑝𝑝𝑝20𝑝2𝛾
𝐻• 𝑒• 𝑒𝑝 𝑠

1𝑘 𝑥 02𝑦 0𝑝 𝑠

𝑘 𝑥 0 𝑦 0 𝑠

𝑝 0𝑝 𝑠𝑝 0𝑝 0𝑥 𝑠𝑥 0𝑝 0𝑦 𝑠

𝑦 0𝜎 𝑠𝜎 0𝑥 𝑠𝑥 01

𝑝 0 𝑝 0 𝑠𝑦 𝑠

𝑦 0𝜎 𝑠𝜎 01 𝑝 0 𝑝 0 𝑠
𝑝 0𝑝 0𝑝 02𝛾𝑝 𝑠

𝑝 0𝑝 𝑠𝑝 𝑠𝑠𝑠

𝑝 0Sext Magnet
𝐻1
𝑘 𝑥6𝑒𝐴 𝑥

𝑝𝑝𝑝

𝑝22𝑝2𝛾𝑝1
𝑘 𝑥6𝐻• 𝑒𝑥 01𝑝 0 𝑝 0 𝑠𝑦 𝑠

𝑦 0𝜎 𝑠𝜎 0𝑝 𝑠𝑝 01 𝑝 0 𝑝 0 𝑠
𝑝 0𝑝 0𝑝 02𝛾𝑝 𝑠

𝑝 0𝑝 𝑠

𝑝 0𝑠

3𝑥𝑦
𝑉

• 𝑒𝑥 𝑠3𝑥𝑦𝑝 𝑠1
𝑘 𝑥 02𝑦 0𝑝 𝑠

𝑘 𝑥 0 𝑦 0 𝑠

𝑝 0𝑝 𝑠𝑝 0𝑥 𝑠𝑥 0𝑦 𝑠𝑦 0𝜎 𝑠

𝜎 0𝑠

𝑝 0

Acceleration• Cavity

𝑉 : Amplitude of RF voltage, 𝑓Synchrotron phase𝑥 𝑠: RF frequency, 𝑛 : n-th order harmonics, 𝜑 :𝑥 0𝑦 0 𝜎 𝑠
𝜎 0𝑦 𝑠𝑝 𝑠
𝑝 0 𝑝 𝑠𝑝 0𝑝 𝑠• Adiabatic Damping (turn by turn)
𝑝 ,𝑝 , 1𝑝 ,𝑝 , 1∆𝐸𝑚𝛾 𝛽∆𝐸𝑚𝛾 𝛽∆𝐸 1𝑝 ,𝑝 , 1𝑚𝛾 𝛽1𝛾𝑝 0

sin

𝜎 0𝜑

sin 𝜑

Space Charge Effect••••Overview
Histogram Making (Charge weighting)
Poisson Solver
Electric fieldOverview
𝜕𝜕𝑥𝜕𝜕𝑦𝜕
𝜑 𝑥, 𝑦, 𝑧𝜕𝑧𝜕𝜕𝑥𝜌 𝑥, 𝑦, 𝑧𝑔 𝑧
𝜀2D approximation

𝜕𝜕𝑦𝑔 𝑧 𝑓 𝑥, 𝑦𝜀𝑢 𝑥, 𝑦𝜎 ,𝜎 ≪ 𝜎

1. Make histogram 𝑓 𝑥, 𝑦 𝑎𝑛𝑑 𝑔 𝑧

The most time-consuming part

2. Solve 2D Poisson equation 𝑢 𝑥, 𝑦 with boundary conditions*
3. Calculate gradient (kick),,𝐴𝑔 𝑧 𝑢 𝑥, 𝑦𝐿

𝐴,𝐿𝑀𝑜𝑣𝑖𝑛𝑔 𝐷𝑖𝑠𝑡𝑎𝑛𝑐𝑒* Using “polar” or “rectangular” boundary condition depending on the duct shapeSpace Charge Effect••••Overview
Histogram Making (Charge weighting)
Poisson SolverElectric fieldMaking Histograms (1)Polar Coordinate𝑟

𝑄 𝑟,𝜃,𝜃

𝑄 𝑟𝑟, 𝜃𝑟,𝜃𝑟𝑟,𝜃𝑞,𝜃

,𝜃

𝑄 𝑟,𝜃

𝑄 𝑟,𝜃𝑟𝑟𝑟𝑟𝑞𝑟𝜃𝜃𝑟 𝛿𝜃𝑟

𝜃 𝜃𝑟 𝛿𝜃𝑟

𝜃𝑟𝑟𝑞𝜃𝜃𝑟 𝛿𝜃𝑟𝑞𝑟𝑟𝑟𝜃𝛿𝜃Making Histograms (2)Cartesian Coordinate
𝑥 ,𝑦𝑥,𝑦𝑄 𝑥 ,𝑦𝑄 𝑥,𝑦𝑄 𝑥 ,𝑦𝑞𝑞𝑥𝑥 𝑦𝛿𝑥𝛿𝑦𝑥𝑞𝑞𝑥𝑥, 𝑦𝑄 𝑥𝑥 ,𝑦𝑥,𝑦,𝑦𝑞𝑥

𝑦𝑥 𝑦
𝛿𝑥𝛿𝑦

𝑦𝑥
𝑦𝛿𝑥𝛿𝑦

𝑦𝑥
𝑦𝛿𝑥𝛿𝑦𝑦Sub-histograms using shared memory
Block N

Block 1𝑥 ,𝑦Thread

𝑥 ,𝑦Thread

𝑥 ,𝑦ThreadAtomicAdd

Shared Memory…

Subhistogram 1

Thread𝑥,𝑦Thread𝑥,𝑦Thread
AtomicAdd

Shared Memory
Subhistogram N…

…

𝑥 ,𝑦

Global Memory=Subhistogram 1+…

TotalHistogram+• Sub-histogram in each shared memory can be made in parallel.
• Each thread in a common block adds entries in shared sub-histogram. But use“exclusive add” operation (called atomicAdd function)

Subhistogram N

Space Charge Effect••••

Overview
Histogram Making (Charge weighting)
Poisson Solver
Electric fieldPoisson Solver (Polar Coordinates)
2D Poisson Equation (Polar Coordinates)𝜕𝑢
1 𝜕𝑟𝑟 𝜕𝑟 𝜕𝑟

1𝜕 𝑢𝑟 𝜕𝜃𝑓 𝑟, 𝜃Boundary Condition : 𝑓 1, 𝜃Discretization with 𝜕𝑦𝜕𝑥
1𝑟 ∆𝑟𝑟∆𝑟 𝑈
2where, 𝑟𝑖𝑈,,∆𝑟∆𝑟, 𝜃𝑟

𝑗∆𝜃

𝑦𝑦
∆𝑥

𝑈
∆𝑟 𝑈 ,2∆𝑟

∆𝑟

0

1𝑁

1 𝑈,𝑟,

12,∆𝜃2𝑈 ,
∆𝜃
2𝜋
𝑁

𝑈,𝑓,

Poisson Solver (Polar Coordinates)∆Using 𝜇∆

,𝛾,𝛿

∆𝜇𝑈,

∆𝛾 𝑃

∆,𝜏

∆

𝛿𝐼 𝑈 ,

𝜏𝑈and,

21𝑃

12⋮0101⋯⋱⋯1010⋮2112𝑓,

Here, a matrix Q is introduced1𝑁𝑄

1𝑁

2𝜋 12cos𝑁

1 2𝑁

12𝜋 22cos
𝑁

1 2𝑁

1

⋯

2𝜋 12cos𝑁
2𝜋 22cos𝑁

1 𝑁/2𝑁
1 𝑁/2𝑁

1𝑁1𝑁2𝜋 𝑁2cos
𝑁

1 2𝑁which satisfies 𝑄𝑄

1⋯2𝜋 𝑁2cos
𝑁𝐼, 𝑄𝑃𝑄𝑇

1 𝑁/2𝑁𝑑𝑖𝑎𝑔 𝜎

1112𝜋 212𝜋 𝑁2sin
𝑁

1𝑁, 𝜎

2𝜋 12sin
𝑁2sin
𝑁𝑁
⋮

⋮
1

14 sin𝑁2𝑁𝑁
2𝑁

𝑁
12𝑁2121⋯

2𝜋 12sin𝑁

1 𝑁
𝑁

12𝜋 22sin𝑁

1 𝑁𝑁

11 𝑁𝑁1⋮21⋯2𝜋 𝑁
2sin𝑁

Poisson Solver (Polar Coordinates)
Multiplying QT from the left, the equation becomes
𝜇 𝑄𝑇𝑈
𝑉,𝑉,
⋮𝑉,Using𝛾 𝜎

𝑄𝑈,𝑈,
⋮𝑈,

𝛿𝜇③𝛾 𝜎⋮00𝑓,0𝑈,

𝜏𝑄 𝑈

𝑄𝑇𝑓 ,,

𝛿𝜏

00⋮⋯⋱
⋯

⋯ 𝑓 ,

⋯ 𝑈,

𝛾𝜇𝜎0

𝛿

𝜇𝜏𝛾𝜎

𝛿

②𝑉,𝑉,
⋮𝑉 ,

①𝑄

𝑓,𝑓,
⋮𝑓 ,

𝑇 using FFT② Solve N linear equations about 𝑉 ,③ Calculate 𝑈 ,𝛿𝐼 𝑄 𝑈 ,, we obtain below

𝜏① Calculate 𝑄𝑇 𝑓 ,𝛾 𝑄 𝑃𝑄,𝑇𝑉,

𝑄 𝑉,⋯ 𝑉 ,𝑉,𝑇 𝑗⋯ 𝑉,1⋯𝑁𝑇𝑖

1 ⋯ 𝑁 using FFT

Using advantagesof GPU

Poisson Solver (Rectangular Coordinates)
2D Poisson Equation (Rectangular Coordinates)
𝜕 𝑢
𝜕𝑦

𝜕 𝑢
𝜕𝑥

𝑢,𝑓 𝑥, 𝑦Boundary Condition : 𝑓 𝑥, 0

Discretization with 𝜕𝑦
𝜕𝑥

𝑦𝑢,𝑢,

2𝑢 ,∆𝑥𝑢

2𝑢 ,
∆𝑦

𝑦
∆𝑥
𝑓,0, 𝑢 , , 𝑢 , , ⋯ , 𝑢 , , 0, 𝑢 , , 𝑢𝑙0, ⋯ ⋯ , 2 𝑚

1

1

𝑓 𝑦, 0𝑖, 𝑗

Odd Extension for x𝑈,

𝑓 𝑥, 𝐿, ,⋯,𝑢 ,𝑓 𝑦, 𝐿1, ⋯ ⋯ , 𝑚

0

Poisson Solver (Rectangular Coordinates)
𝐷𝐹𝑇𝑈2𝑈 ,∆𝑥,0𝑈,𝑢 ,1 2𝑒2𝑒𝑒𝑒

2𝑒𝑉,2𝑉 ,∆𝑥𝑉,𝑒𝑢 ,𝑢 ,2𝑢 ,𝑢 , ⋯0𝑒𝑢 ,𝑒𝑢 ,⋯𝑈,0, 𝑈 , , 𝑈 , , ⋯ , 𝑈 , , 0, 𝑈 , , 𝑈 ,𝑙2𝑢 ,𝑢 ,𝑒𝑝𝜋
4 𝑠𝑖𝑛
2 𝑚 1
Odd Extension for y

𝐷𝐹𝑇 𝐷𝐹𝑇𝑢 ,𝑒𝑒𝑉,𝑒0, ⋯ ⋯ , 2 𝑚,⋯, 𝑈,
14𝑝𝜋𝑠𝑖𝑛
∆𝑥
2 𝑚 1

1𝑒𝑒

𝑉,

Poisson Solver (Rectangular Coordinates)
𝐷𝐹𝑇 𝐷𝐹𝑇𝜕 𝑢𝜕𝑥4

𝜕 𝑢
𝜕𝑦→41𝑝𝜋𝑠𝑖𝑛
∆𝑥
2 𝑚 1

1𝑝𝜋𝑠𝑖𝑛∆𝑥
2 𝑚 1𝑞𝜋1𝑠𝑖𝑛
2 𝑚 1
∆𝑦𝑞𝜋
1
𝑠𝑖𝑛
2 𝑚 1
∆𝑦𝑒

𝐷𝐹𝑇 𝐷𝐹𝑇 𝑉

Odd Extension of 𝑓 ,𝐹,0, 0,0, ⋯ , 0, 0,0,0, ⋯ , 0

𝐹,0, 𝑓 , , 𝑓 , , ⋯ , 𝑓 , , 0, 𝑓 , , 𝑓, ,⋯,𝑓,𝐹,0, 𝑓 , , 𝑓 , , ⋯ , 𝑓 , , 0, 𝑓 , , 𝑓, ,⋯,𝑓,𝐹,𝐹,𝐹,0, 𝑓 , , 𝑓 , , ⋯ , 𝑓 , , 0, 𝑓 , , 𝑓 , , ⋯ , 𝑓 ,0, 0,0, ⋯ , 0, 0,0,0, ⋯ , 0
0, 𝑓 , , -𝑓 , , ⋯ , 𝑓 , , 0, 𝑓 , , 𝑓 , , ⋯ , 𝑓 ,⋮𝐷𝐹𝑇 𝐷𝐹𝑇 𝑓→𝑒𝑒𝐹,𝑒𝑉,

Poisson Solver (Rectangular Coordinates)
𝐷𝐹𝑇 𝐷𝐹𝑇41𝑝𝜋𝑠𝑖𝑛∆𝑥2 𝑚 1𝜕 𝑢𝜕𝑥

𝜕 𝑢𝜕𝑦𝑞𝜋1𝑠𝑖𝑛
2 𝑚 1
∆𝑦

𝐷𝐹𝑇 𝐷𝐹𝑇 𝑓

𝐷𝐹𝑇 𝐷𝐹𝑇 𝑉

𝐷𝐹𝑇 𝐷𝐹𝑇 𝑓

𝑝0 𝑜𝑟 𝑞0

𝐷𝐹𝑇 𝐷𝐹𝑇 𝑉

𝐷𝐹𝑇 𝐷𝐹𝑇 𝑓𝑝𝜋𝑞𝜋114𝑠𝑖𝑛𝑠𝑖𝑛
2 𝑚 12 𝑚 1∆𝑦∆𝑥

𝑝0 𝑎𝑛𝑑 𝑞0

𝐷𝐹𝑇 𝐷𝐹𝑇 𝑉

𝑉,𝑉,𝑖𝐷𝐹𝑇 𝑖𝐷𝐹𝑇 𝐷𝐹𝑇 𝐷𝐹𝑇 𝑓𝑝𝜋𝑞𝜋114𝑠𝑖𝑛𝑠𝑖𝑛
2 𝑚 12 𝑚 1∆𝑦∆𝑥
𝑉, → 𝑢 ,0(𝑉 , is odd extension of 𝑢 ,4 FFTs using GPU• Potential @ a collimatorr

Arbitrary Unit

• Potential @ a circular duct (R=0.065 m)

Arbitrary Unit

ExampleSpace Charge Effect••••Overview
Histogram Making (Charge weighting)
Poisson Solver
Electric fieldElectric Field
Polar Coordinate𝑟𝑟,𝜃

𝐸 : 𝐸 𝑟 ,𝜃

,𝜃

𝐸𝑟𝑟,𝜃

𝐸𝑟

𝐸𝑟,𝜃

,𝜃

,𝜃
𝑟𝑟, 𝜃𝐸𝑟,𝜃

𝐸𝑟,𝜃𝑟,𝜃

,𝜃,𝜃𝐸 : 𝐸 𝑟 ,𝜃𝜑 𝑟

𝐸 :𝜑 𝑟𝐸 𝑟,𝜃

𝐸 : 𝐸 𝑟 ,𝜃𝐸 , 𝑟, 𝜃,𝑟

𝜑 𝑟

𝜑 𝑟
2𝛿𝑟,𝜃

,𝜃

𝜑 𝑟 ,𝜃

𝐸 𝑟

2𝛿𝑟

𝜑 𝑟,𝜃

𝜑 𝑟
2𝛿𝑟,𝜃

,𝜃

𝜑 𝑟 ,𝜃

,,,𝜃𝐸

𝐸 cos 𝜃𝐸

𝐸 sin 𝜃𝑟,𝜃,

1𝐸 sin 𝜃𝑟
1𝐸 cos 𝜃𝑟

𝜑 𝑟 ,𝜃
2𝛿𝜃
,𝜃

𝜑 𝑟
2𝛿𝜃

𝜑 𝑟 ,𝜃

𝐸 𝑟 ,𝜃

,,

𝜑 𝑟,𝜃

𝐸 𝑟

2𝛿𝑟,

𝜑 𝑟,𝜃

𝐸 𝑟 ,𝜃

,𝜃
𝜑 𝑟,𝜃

2𝛿𝜃
𝜑 𝑟,𝜃

,𝜃

𝜑 𝑟
2𝛿𝑟,,𝜃Electric Field
Rectangular Coordinate

𝑥,𝑦Interpolate potential using Bezier Surface

𝑥 ,𝑦𝑥,𝑦𝑥 ,𝑦𝑥,𝑦𝑥,𝑦𝑥 ,𝑦𝑥,𝑦𝑥,𝑦𝑥 ,𝑦𝑥,𝑦𝜑 𝑥, 𝑦3!
𝑥 𝑥
𝑚! 3 𝑚 ! 𝛿𝑥𝑦 𝑦𝑦𝑦3!
𝛿𝑦𝑛! 3 𝑛 ! 𝛿𝑦

𝜑 𝑥,𝑦-1𝑥,𝑦𝐸 𝑥, 𝑦

𝑔𝑟𝑎𝑑 𝜑 𝑥, 𝑦𝑥𝑥

𝜕 𝜕,𝜑 𝑥, 𝑦𝜕𝑥 𝜕𝑦𝛿𝑥