Использование CUDA в расчете динамики пучка

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

21. __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” необходимо максимально сократить

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

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

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

25. __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) Массив данных для функции Rho трех переменных NY NZ

26. __global__ PoissonSolver ( ) • Номер нити int n = threadIdx.x+blockIdx.x*blockDim.x; • Каждая нить находит значение коэффициентов Фурье PhiF потенциала Phi PhiFind(i,j,k) = -RhoFind(i,j,k) / ( kxi2 + kyj2 + kzk2 ) В узле с номером: 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.

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

28. Результаты

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

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

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

32. Анимация

33. Анимация

34. Анимация

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

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

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

38. S0 S1 S2 S3 S4 Выбор оптимальной конфигурации

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

40. Производительность на 8800GTX *Размер сетки: 25 x 25 x 25. Число частиц: 100,000 треугольников: 2054 **CPU с частотой 2.4 ГГц

41. Сравнение CPU и GeForce 8800GTX *CPU с частотой 2.4 ГГц

42. Сравнение CPU с Tesla C1060 БЕЗ пространственного заряда

43. Сравнение CPU с Tesla C1060 С пространственным зарядом

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

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