Использование CUDA в расчете динамики пучка С.Б. Ворожцов, В.Л. Смирнов, Е.Е. Перепелкин Дубна, ОИЯИ 6 апреля

Циклотрон Постановка задачи Численные методы Программная реализация на CUDA Результаты CBDA: Cyclotron Beam Dynamic Analysis code

Постановка задачи

Линия инжекции ЭСД Дуант Магнитный сектор Инфлектор Компьютерная модель циклотрона

Области задания карт полей Инфлектор Электрическое поле Аксиальный канал Магнитное поле Линза Магнитное поле

Ресурсоемкое моделирование Необходимость рассмотреть не менее 5 различных конфигураций центральной зоны; Необходимость ускорять различные ионы; Сложная геометрическая структура; Учет пространственного заряда; Одна итерация требует ~ несколько дней расчетов

Уравнения движения

Пространственный заряд PIC методPP метод

Численные методы

Уравнение движения Уравнение движения из постановки задачи можно представить в упрощенном виде, дополнив его вторым уравнением для определения координат частиц

Пример решения ОДУ Рассмотрим решение обыкновенного дифференциального уравнения (ОДУ) методом Рунге -Кутта Задача Коши

Метод Рунге - Кутта

Решение краевой задачи При поиске коэффициентов Фурье используется алгоритм БПФ (Быстрого Преобразования Фурье)

Задание области для краевой задачи Lx Lz Ly Lx Lz Ly Lz Ly Lx

Раздача плотности заряда Ячейка 7 Ячейка 1 Ячейка 8 Ячейка 3 Ячейка 2 Ячейка 5 Ячейка 6 Узел

Потери частиц tntn t n+1 A B C D Если точка D принадлежит треугольнику ABC, тогда Условие пересечения где ε Δ – допустимое отклонение от поверхности

Программная реализация на CUDA

Функции ядра Track ( карты полей, координаты и скорости частиц ) метод Рунге-Кутта Losses ( геометрия установки, координаты частиц ) проверка пересечений с геометрией Rho ( координаты частиц ) раздача заряда в узлы сетки FFT ( функция плотности заряда или потенциал) БПФ по базисным функциям sin(πn/N) PoissonSolver ( Фурье коэффициенты ) решение краевой задачи E_SC ( потенциал электрического поля ) поиск электрического поля

__global__ void Track ( ) Много входных параметров. Использование типа переменной __constant__ для неизменных параметров: __device__ __constant__ float d_float[200]; __device__ __constant__ int d_int[80]; Каждой частице соответствует нить: int n = threadIdx.x+blockIdx.x*blockDim.x; Количество if, goto, for необходимо максимально сократить

Проблема количества if, goto, for Инфлектор Электрическое поле Аксиальный канал Магнитное поле Линза Магнитное поле

__global__ void Losses ( ) Нити одного блока копируют вершины треугольников из global в shared память. Синхронизация нитей после копирования треугольников __syncthreads() Каждой частице соответствует номер нити: int n = threadIdx.x+blockIdx.x*blockDim.x; Проверка условия пересечения частицей c номером n, загруженных в shared память, треугольников Для каждого блока геометрии есть своя функция Losses

__global__ void Rho Каждая частица с номером n = threadIdx.x + blockIdx.x*blockDim.x дает свой вклад, в окружающие ее узлы. Для этого по координатам частицы определяется какой ячейки она принадлежит Одна частица может дать вклад в 8 ближайших узлов. Таким образом, каждая нить заполняет свои 16 ячеек в общем массиве вклада: 8 – номеров узлов и 8 – значений вклада. Далее производится сложение этих вкладов для каждого узла.

__global__ FFT ( ) Действительное БПФ по базисным функциям sin(πn/N); 3D преобразование состоит из трех последовательных 1D БПФ по осям: X, Y, Z соответственно int n = threadIdx.x+blockIdx.x*blockDim.x; k=(int)(n/(NY+1)); j=n-k*(NY+1); m=j*(NX+1)+k*(NX+1)*(NY+1); FFT_X[i+1]=Rho[i+m]; n = j + k*(NY+1) NZ NY Массив данных для функции Rho трех переменных

__global__ PoissonSolver ( ) Номер нити int n = threadIdx.x+blockIdx.x*blockDim.x; Каждая нить находит значение коэффициентов Фурье PhiF потенциала Phi PhiF ind(i,j,k) = -RhoF ind(i,j,k) / ( kx i 2 + ky j 2 + kz k 2 ) В узле с номером: ind(i,j,k)=i+j*(NX+1)+k*(NX+1)*(NY+1), где k=(int)(n/(NX+1)*(NY+1)); j=(int)(n-k*(NX+1)*(NY+1))/(NX+1); i=n-j*(NX+1)-k*(NX+1)*(NY+1); RhoF – коэффициенты Фурье для функции плотности заряда Rho.

__global__ E_SC ( ) Вычисление электрического поля в узле с номером int n = threadIdx.x+blockIdx.x*blockDim.x+st_ind φnφn φ n + 1 φ n - 1 φ n - ( NX + 1 ) φ n + ( NX + 1 ) φ n - ( NX + 1 )( NY + 1 ) φ n + ( NX + 1 )( NY + 1 )

Результаты

Аксиальная инжекция пучка

Процесс банчировки пучка

Ускорение в циклотроне

Анимация

Потери частиц

Ускорение банчей

Оптимизация центральной области φ RF = 13° φ RF = 15° φ RF = 28° φ RF = 10° С постами Дуант Без постов «Земля» F = ZU RF - W GAP

Выбор оптимальной конфигурации S0S1S2 S3S4

Распределение ускоряющего поля

Производительность на 8800GTX Функции* Время, [мс]Ускорение, [раз] CPU ** GPU Track Losses Rho79614 Poisson/FFT35313 E_SC Total *Размер сетки: 2 5 x 2 5 x 2 5. Число частиц: 100,000 треугольников: 2054 **CPU с частотой 2.4 ГГц

Сравнение CPU и GeForce 8800GTX Число частиц Время вычислений Ускорение, [раз] CPU * GPU 1,0003 мин. 19 c.12 c.17 10,00034 мин. 14 с.42 с ,0005 ч. 41 мин.6 мин.56 1,000,0002 дня 8 ч. 53 мин.1 ч.60 *CPU с частотой 2.4 ГГц

Сравнение CPU с Tesla C1060 Число частиц Время вычислений Ускорение, [раз] CPU 2.5ГГц GPU C ,0003 мин. 12 с.11 с.18 10,00032 мин. 24 с.27 с ,0005 ч. 14 мин. 31 с.3 мин. 34 с.88 1,000,0002 дня 4 ч. 25 мин.34 мин. 29 с.91 БЕЗ пространственного заряда

Сравнение CPU с Tesla C1060 Число частиц Время вычислений Ускорение, [раз] CPU 2.5 ГГц GPU C ,00033 мин. 36 с.44 с ,0005 ч. 28 мин. 12 с.5 мин. 4 с.65 1,000,0002 дня 8 ч. 27 мин.50 мин. 17 с.67 С пространственным зарядом

Эффект пространственного заряда I ~ 0 Потери 24% I = 4 мА Потери 94%

Заключение Очень дешевая технология в сравнении с CPU; Увеличение производительности на 1.5 – 2 порядка дает шанс проведения моделирования ресурсоемких физических моделей; Требует аккуратного программирования.