Научная статья на тему 'Интегрированная среда моделирования спин-орбитального движения заряженных частиц'

Интегрированная среда моделирования спин-орбитального движения заряженных частиц Текст научной статьи по специальности «Математика»

CC BY
192
60
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
СПИН-ОРБИТАЛЬНОЕ ДВИЖЕНИЕ / МАТРИЧНОЕ ИНТЕГРИРОВАНИЕ / КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ЧАСТИЦ / SPIN-ORBIT MOTION / MATRIX INTEGRATION / BEAM DYNAMICS SIMULATION

Аннотация научной статьи по математике, автор научной работы — Иванов Андрей Николаевич

Охарактеризована интегрированная среда моделирования спин-орбитальной динамики заряженных частиц во внешних электромагнитных полях. Приведены уравнения, оценивающие динамику в обобщенных криволинейных системах координат. Динамика спина описана в этих же координатах на основе уравнения Баргмана-Мишеля-Телегди. Также приведен метод матричного интегрирования обыкновенных дифференциальных уравнений, основанный на построении нелинейного отображения. Разработанная среда предоставляет гибкие и удобные механизмы исследования динамики частиц в произвольно заданных управляющих полях. Имеется возможностьгенерации вычислительного кода на языках MATLAB и C++. Предложенный метод решения дифференциальных уравнений обладает свойством естественного распараллеливания и может бытьреализован до произвольного порядка нелинейности. Созданные программные средства протестированы в сравнении с имеющимися аналогами.

i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «Интегрированная среда моделирования спин-орбитального движения заряженных частиц»

УДК 517.938

Вестник СПбГУ. Сер. 10. 2014. Вып. 2

А. Н. Иванов

ИНТЕГРИРОВАННАЯ СРЕДА МОДЕЛИРОВАНИЯ СПИН-ОРБИТАЛЬНОГО ДВИЖЕНИЯ ЗАРЯЖЕННЫХ ЧАСТИЦ

Санкт-Петербургский государственный университет, 199034, Санкт-Петербург, Российская Федерация

Охарактеризована интегрированная среда моделирования спин-орбитальной динамики заряженных частиц во внешних электромагнитных полях. Приведены уравнения, оценивающие динамику в обобщенных криволинейных системах координат. Динамика спина описана в этих же координатах на основе уравнения Баргмана—Мишеля—Телегди. Также приведен метод матричного интегрирования обыкновенных дифференциальных уравнений, основанный на построении нелинейного отображения. Разработанная среда предоставляет гибкие и удобные механизмы исследования динамики частиц в произвольно заданных управляющих полях. Имеется возможность генерации вычислительного кода на языках MATLAB и C+—+. Предложенный метод решения дифференциальных уравнений обладает свойством естественного распараллеливания и может быть реализован до произвольного порядка нелинейности. Созданные программные средства протестированы в сравнении с имеющимися аналогами. Библиогр. 13 назв. Ил. 4. Табл. 1.

Ключевые слова: спин-орбитальное движение, матричное интегрирование, компьютерное моделирование динамики частиц.

A. N. Ivanov

AN INTEGRATED DEVELOPMENT ENVIRONMENT FOR SPIN-ORBIT MOTION SIMULATION OF CHARGED PARTICLES

St. Petersburg State University, 199034, St. Petersburg, Russia Federation

In the article an integrated development environment for spin-orbit dynamics simulation of charged particles motion in electromagnetic fields is described. The necessity of a new computer modeling tool development is caused by deficiencies in existing software. The research is consists of three part. First of all, the equations of a particle spin-orbit dynamics in generalized coordinate systems are considered. The given mathematical model allows not only reflecting the physical properties of considered system, but also its adequate implementing in program codes. Spin dynamics is described in the classical representation of T—BMT equation. The second part is devoted to the numerical implementation of the nonlinear matrix integration of systems of ordinary differential equations. Numerical method that is described in the article is based on Taylor series expansion in matrix form. All computational operations are implemented via addition and multiplication of numerical matrices, that allows to easy realize the approach in parallel techniques. Finally, the developed computer modeling environment is described. The computational results have been compared with existing simulation packages. The given software provides flexible and useful tool for particle dynamics investigation in arbitrary control fields. The program has scalable mechanism that allows to create new control elements by declarative description of electromagnetic fields. The software provides a user graphic interface. Moreover the possibility for computational code generation in MATLAB and C+—+ programming languages exists. The developed method for ordinary differential equations has a native parallel feature and can be easy implemented up to the necessary order of nonlinearity. The computational simplicity allows also to implement the approach in parallel and distributed computational systems. Bibliogr. 13. Il. 4. Table 1.

Keywords: spin-orbit motion, matrix integration, beam dynamics simulation.

Иванов Андрей Николаевич — аспирант; e-mail: 05x.andrey@gmail.com; a.ivanov@fz-juelich.de Ivanov Audrey Nikolaevich — post-graduate student; e-mail: 05x.andrey@gmail.com; a.ivanov@ fz-juelich.de

Введение. Для исследования динамики пучков заряженных частиц используются различные пакеты моделирования. Каждый из них основывается на специфических математических моделях, а для вычислений применяются разные численные алгоритмы. Среди наиболее распространенных программ можно выделить COSY Infinity [1], MAD (DIMAD, SMAD) [2], OptiM, TRANSPORT. Соответствующие программные пакеты отличаются быстродействием и предоставляемым функционалом. Основной проблемой при использовании указанных программ является тот факт, что ни одна из них не удовлетворяет всем требованиям, возникающим при изучении спин-орбитального движения заряженных частиц.

Цель настоящей работы - построение среды программирования, которая позволит пользователю исследовать динамику заряженных частиц в накопительных и ускорительных кольцах. В основу численного метода закладывается алгоритм матричного интегрирования [3, 4] дифференциальных уравнений в нелинейном приближении. К преимуществам разрабатываемого пакета программ, отличающих его от вышеуказанных сред, является наличие в нем следующих возможностей, предлагаемых конечному пользователю:

• задание произвольного вида управляющих полей;

• учет ошибок задания управляющих элементов;

• графический интерфейс, включающий в себя визуальный конструктор ускорительного кольца, текстовый редактор программного кода с подсветкой синтаксиса;

• интеграция со сторонними продуктами и генерация кода (MATLAB, C+—Ъ).

Исследователям-физикам предлагается интегрированная среда разработки

(IDE), облегчающая процесс программирования, описания задач и позволяющая им акцентировать внимание на анализе предмета.

Алгоритмы, заложенные в программное решение, прошли всестороннюю апробацию при проведении научно-исследовательских работ в рамках международной кол-лаборации JEDI, посвященной измерению электрического дипольного момента элементарных частиц.

Предлагаемая программная среда построена на основе гибкой и масштабируемой архитектуры, позволяющей вносить изменения и модифицировать алгоритмы без полной перекомпиляции исходного кода. Решение позволяет как проводить численные расчеты, так и готовить декларативное описание вычислительных задач, которое может быть выполнено на сторонних платформах. В частности, реализована функция генерирования файлов MATLAB для предоставления возможности проводить вычисления в этой среде моделирования. Также доступна кодогенерация на языке C+—Ъ, которая позволяет запускать расчетные задачи на гетерогенных вычислительных ресурсах.

Основные уравнения. Движение заряженной частицы в электромагнитных полях в релятивистском случае описывается уравнением Ньютона-Лоренца

^ = Q(E+[vxB]),

1

р = m07V, 7 =

л/а-^/с2)

в котором р - вектор импульса частицы, Е и В - вектора напряженности электрического поля и магнитной индукции, с - скорость света, то, Q - масса покоя и заряд частицы, - вектор и модуль скорости частицы, 7 - фактор Лоренца.

Далее будем считать все поля статическими, а также использовать обозначения, введенные Силадьи в [5]. В криволинейной ортогональной системе координат это уравнение запишется в виде

dpi | / qj dhj qi+1 dhi+1 \ ^ f щ dhj qi+2 dhi+2 \ _ dt \hi+1 dqi+1 hi dqi JPt+1 \hi+2 dqi+2 Ы dqi JРг+2 (1) = Q(Ei + hi+i<ji+iBi+2 — hi+2q_i+2Bi+i), i = 1, 2, 3,

где qi - криволинейные координаты; hi - коэффициенты Ламе; индексы меняются циклически (q4 = qi, q5 = q2); оператор « • » обозначает дифференцирование по времени. Соотношения (1) описывают изменение криволинейных координат частицы с течением времени. Наряду с таким подходом применяется запись уравнений, когда в качестве независимой переменной выбирается одна из координат (например, q3), и оператор дифференцирования по времени заменяется на дифференцирование по координате q3:

£ =_v__— (О)

dt (h2q2 + h2q'2 + hi)1/2 dq3' У>

Здесь и далее оператор « '» обозначает дифференцирование по выбранной координате. Учитывая (2), соотношения (1) преобразуются в уравнения для проекций траектории на две взаимно перпендикулярные плоскости

т + <«-- (4 -«t) + + <

- = ^ - + ms* - №>, й + и - + („; - +

- ¥1Г^V^ = -^r-d - v'/C')''HHE,/V + inBl - /,„;в3),

h2 oq2 J moh2V

где H и G - функции криволинейных координат и скоростей:

H =(hlq[ 2+h2q2 2 + hl)i/2, G = -^-(hiq[B2 - h2q'2Bi + НЕз/v) -

moh3

" (i-wrv'LjJlp-lLp) + &

(3)

Hh3 \ \ hi dqi h3 dq3

+ h2' f—^i - si Obi) + Obi)

2 2 V h2 dq2 h3 dq3 J dq3 J

дq3 ) '

Формулы (3) представляют собой релятивистские уравнения траектории в обобщенной ортогональной системе координат. Далее в качестве траектории, вдоль которой рассматривается движение частиц, будем выбирать некоторую кривую, а как независимая переменная интегрирования выступит длина вдоль этой кривой. Координатные оси ql = х, q2 = у будем выбирать как нормаль и бинормаль к касательной в точке. Таким образом, вводится естественная система координат вдоль заданной кривой, когда в качестве фазовых переменных, описывающих траекторное движение

частиц, выступают пары х — х', у — у', Ь — V. Наряду с этим представлением можно применять и другие координаты. Учитывая, что гамильтоновый характер динамики частиц [6] налагает требование симплектичности [7], предъявляемое к конечному решению, перепишем уравнения для пар канонических координат (х,рх), (у,ру), (Ь, Т), где рх и ру являются проекциями импульса на оси х и у, Т соответствует кинетической энергии частицы.

В качестве опорной будем использовать плоскую кривую Н\ = Нх = 1, Н2 = Ну = 1, Нз = На. Это допущение оправдано в виду конструктивных особенностей моделируемых физических полей, в которых опорные кривые являются либо прямыми (свободные промежутки, мультипольные линзы), либо дугами окружностей (поворотные магниты, электростатические конденсаторы). При этом проекции импульса

на координатные оси (£ € {х,у, в}) равны

р* = —£ • (5)

Разрешив равенства (5) относительно х', у', получим

Г= , (6)

У(т 0^)2 — Р2х — Ру Дифференцируя (5) по в и учитывая соотношения (6), имеем

7' V С ( Рх х'' ру у'' НаН'а

/ I I ' I > I гх ~ гу а '»а-а (:-7\

Ра =Р£ — + — +то«7т7 -Р£ -"77 +-"77 + ~пт • 7

^ XV 7 / Н \moVY Н Н Н2 )

В уравнения (7) явно входят значения v,v',J,J', зависящие от фазовых координат в данный момент времени (независимой переменной в). Запишем их как функции V = v(T), V' = v'(T,T'), 7 = 7(Т), 7' = ^'(Т'), используя соотношения

Т2 2 = шо7с — шос ,

Т' = —дФ'(х, у, в) = д(Ехх' + Еуу' + Еа), ( )

где Ф - электростатический потенциал поля. Из уравнений (8) следует

Т + шос2 , Т'

7 =-5—, 7/ =

22 шос2 шос2

Т + шос2

л/Т2 + 2Тто0с2, (9)

гтг 2 4

1 ст°с :(Т2+2Тт0с2)-^

(Т + шос2)2

Уравнения (2), (6)-(8) задают систему обыкновенных дифференциальных уравнений, описывающих орбитальное движение заряженных частиц во внешних полях. Динамика спина, в свою очередь, описывается уравнением Баргмана—Мишеля— Телегди (БМТ) [6], которое в криволинейной системе координат с плоской опорной кривой примет вид

с

V

тт

о = __г_((1+,с)в * _ (е + ' , 1£21?Г

\ 1+7 т^с2 1 + 7 шов2

здесь 8 - трехмерный вектор спина, О = (д — 2)/2, д - аномальный магнитный момент.

Численный метод интегрирования. В предлагаемой среде моделирования для численного анализа динамики частиц используется матричное интегрирование уравнений (2), (6)-(10). В отличие от традиционных пошаговых схем указанный подход позволяет строить нелинейное отображение целиком для заданного интервала интегрирования [3, 8]. Данная концепция удачным образом отражает предметную область. В накопительном кольце ансамбль частиц последовательно проходит физические элементы (линзы, дефлекторы и пр.), для которых и строятся матрицы нелинейного перехода. Последовательная конкатенация этих отображений представляет собой нелинейное преобразование начального фазового момента в конечное, отвечающее полному обороту частиц в накопительном кольце. Подход численного интегрирования на основе построения отображений заметно выигрывает по производительности у пошаговых методов в виду отсутствия необходимости деления интервала интегрирования на шаги и независимости коэффициентов матричного отображения от начальных условий интегрирования.

Матричное интегрирование основано на разложении правых частей дифференциальных уравнений в ряд Тейлора до заданного порядка нелинейности. Для обеспечения условий сходимости разложения функций в ряд введем замену переменных:

(11)

Х1 = Х, Х2 = У, Хз =

Х4 _ Рх Ро' Х5 _ Ру 1 «^6 ро = 5Т,

Х7 = Sx, Х8 = Б у, Хд =

В (11) ро - начальный импульс, 5Т = (Т — То)/То, То - начальное значение кинетической энергии частицы. В введенных обозначениях систему обыкновенных дифференциальных уравнений, описывающих спин-орбитальное движение заряженных частиц во внешних полях, можно представить в виде

±Х = ¥(з,Х), (12)

ав

где X = (х1,х2,хз,х4,х5,хб,х?,х8,хд)т, а вектор-функция Е задается уравнениями (2), (6)-(10) с точностью до замены переменных.

Уравнение (12) может быть представлено как матричное разложение в ряд Тейлора

^х = ]Гр1'£(5)ХМ, (13)

к=о

в котором Х[к является кронекеровской степенью вектора X, а матрицы Р1к могут быть вычислены следующим образом:

1к 1 д^,Х0)

[ ' ~ (к)\ д(х^)т '

(к)! = к\\... кп!, а каждая компонента вектора Xк есть моном, состоящий из степеней переменных, и вычисляется по формуле (х^1,..., хП), к\ + ... + кп = к.

Решение системы (13) может быть представлено в виде ряда

^

X = £ Е1к да*1. (14)

к=0

Элементы матриц Н1к зависят от независимой переменной в и могут быть вычислены в символьном виде [3]. В разработанной версии метода используется численная реализация алгоритма вычисления отображения, основная идея которого заключается в том, что сначала строится дополнительная система дифференциальных уравнений относительно элементов матриц Н1к, после чего эта новая система решается некоторым пошаговым методом интегрирования. В результате отображение представляется числовыми матрицами, отвечающими заданному интервалу интегрирования. Рассмотрим подробнее алгоритм получения матриц К1к.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Продифференцировав (14) и учитывая (13), имеем систему

¿X _ ^¿Н^з) [к]

Л '

к=0

- ш ®

к=0

?1кI

к=1 к=1

\11

частные производные которой по ХО*1 равны

¿в

к=1 (15)

Уравнения (15) определяют систему обыкновенных дифференциальных уравнений, описывающих эволюцию матриц Н1к вдоль интервала интегрирования.

Для решения системы (15) может быть применен любой метод пошагового интегрирования. Заметим, что низкопроизводительное пошаговое интегрирование запускается один раз на этапе построения матриц, после чего они сохраняются в памяти и используются для вычисления образа вектора начального состояния Хо по формуле (14). В текущей версии среды моделирования пошаговым методом служит неявная схема Рунге—Кутта 4-го порядка (см., например, [9]).

Реализация алгоритма в программном коде. Из уравнений (2), (6)-(10) видно, что для описания спин-орбитальной динамики заряженной частицы достаточно задать управляющие внешние электромагнитные поля Ех, Еу,Еа, Вх, Ву,В8 и выбрать плоскую опорную кривую (к8) для ортогональной системы координат.

Рис. 1. Класс-диаграмма моделируемых электромагнитных элементов

На рис. 1 представлена обобщенная диаграмма разработанных классов электромагнитных элементов. Класс EMElement представляет собой абстрактный класс, инкапсулирующий логику дифференциальных уравнений (2), (6), (7), (8)-(10) и их решения в матричном виде (15). Классы StraightElement и ArcElement также относятся к абстрактным, но определяют метод hs(x,y,s) и описывают движение вдоль прямой (hs = 1) и дуги окружности постоянного радиуса R (hs = 1 + x/R) соответственно. Они являются родительскими классами для произвольного элемента, который задается внешними управляющими полями E,B.

Для реализации разложения нелинейных уравнений в ряд Тейлора была разработана подсистема символьных вычислений, основанная на оперировании с символьными полиномами. Все операции выполняются вплоть до заданного порядка нелинейности. Внутренняя логика символьных структур данных применяет в своей работе хэш-таблицы, что позволяет производить операции на приемлемом уровне производительности. При выборе средств разработки учитывались следующие требования:

• скорость разработки;

• производительность вычислений;

• визуализация полученных результатов.

В качестве языка программирования выбран C#, особенности которого [10] дают возможность достигнуть компромисса между удобством написания кода и скоростью вычислений. Кроме того, при реализации программ, по возможности, не использовались особенности синтаксиса языка, специфичные для платформы .NET. Код написан в С-подобной нотации. Последнее замечание позволяет без дополнительных модификаций переносить некоторые из разработанных вычислительных библиотек на C+—+ для проведения сравнительного анализа производительности вычислений. Средой разработки была выбрана MS Visual Studio 2010. Указанная IDE обладает широкими возможностями по организации кода, проведению его тестирования и сборки исходных файлов в исполняемые модули.

В процессе разработки программного обеспечения применялись технологии компонентного и объектно-ориентированного программирования. Все алгоритмические компоненты и численные методы реализованы в виде переносимых dll библиотек. Разработанные графические интерфейсы написаны на технологии Windows Presentation Foundation (WPF) и соответственно работают в операционной системе Windows с установленным .NET Framework 3.5 или выше.

Сравнение результатов численных расчетов с COSY Infinity. COSY Infinity в настоящее время является единственной программой общего пользования, позволяющей осуществлять численное моделирование спин-орбитального взаимодействия заряженных частиц. Несмотря на ряд спорных решений, применяемых в этой программе (например, незамкнутая референс-орбита в общем случае и, как следствие, невозможность изучения ошибок задания управляющих полей на динамику частиц), COSY Infinity зарекомендовала себя как средство, результаты расчетов на котором полностью совпадают с аналитическими выкладками. Последние проводились в таких приближениях как отсутствие краевых полей и идеальные внешние поля физических элементов, когда референс-орбита строится всегда замкнутой, исходя из физических соображений, и вследствие этого делает корректным применение для расчетов COSY Infinity.

Численный метод интегрирования дифференциальных уравнений в COSY Infinity, так же как и вышеописанный матричный подход, является методом, построенным на основе отображений в виде разложения в ряд Тейлора до заданного порядка нелинейности, однако для оценки элементов отображения (в случае с COSY Infinity используется тензорный формализм) применяется дифференциальная алгебра. Фазовые координаты, описанные в данной работе, полностью соответствуют математической модели в COSY Infinity, что делает возможным сравнение результатов вычислений на разных программах в численном виде.

х-Ю-3 х-Ю-3

Рис. 2. Орбитальное движение (а) и спиновая динамика (б) частицы (COSY Infinity (1) и разработанный метод (2))

Результаты расчета сравнивались: визуально - спин-орбитальная динамика (рис. 2) и численно - скорости вращения спина. Последняя величина является

интегральной и усредненной, что делает ее менее зависимой от возможных ошибок округления и несоответствия математических моделей, заложенных в разные программы. Последовательность тестовых сценариев приведена в таблице. Значения в ней соответствуют времени (в секундах), за которое вектор спина частицы повернется на 2п радиана. Начальные значения всех фазовых координат частицы равны нулю, если не указаны явно, а также Ss = 1. Рассматриваются электростатические элементы.

Скорость вращения спина (с)

Тестовый случай COSY Infinity Разработанный метод

ST = 1 • 10 ST = 3 • 10 Ax = 0.003

ST =1 • 10 ST = 3 • 10 Ax = 0.003

Цилиндрический конденсатор

5749.4 635.6

1184.3

Цилиндрический конденсатор X 5705.1 633.9

1188.5

5749.0 635.5 1184.3 16

5704.6 633.8 1188.5

Кольцо из конденсаторов и квадруполей

ST =1 • 10 ST = 3 • 10 Ax = 0.003

0.2008 0.0704 2072.3

0.2008 0.0704 2072.3

Кольцо из конденсаторов и квадруполей с RF полем

ST =1 • 10

ST = 3 • 10

4438.2 492.9

4415.3 491.7

Цилиндрический конденсатор. В данном случае изучается поле отдельного цилиндрического конденсатора. Подобное соответствие имеет место и для мультипо-лей (квадруполи, секступоли).

Цилиндрический конденсатор х 16. Накопительное кольцо строится из серии последовательно примыкающих друг к другу конденсаторов. Референс-частица в таком случае движется по окружности постоянного радиуса. После построения отображения для отдельного конденсатора матрицы конкатенируются 16 раз для получения суммирующей нелинейной матрицы перехода, отвечающей динамике всего кольца. Следует отметить, что в этом случае происходит потеря точности. Так, например, при объединении двух отображений 2-го порядка нелинейности результирующее отображение имеет элементы 4-го порядка, которые необходимо отбрасывать. И COSY Infinity, и разработанная программа показали одинаковую тенденцию изменения спиновой динамики, связанную с ошибками округления порядка нелинейности.

Кольцо из конденсаторов и квадруполей. Взят пример кольца, состоящего из цилиндрических конденсаторов и квадруполей, находящихся между ними и осуществляющих фокусировку частиц.

Кольцо из конденсаторов и квадруполей RF полем. К предыдущему примеру добавлено RF поле, приводящее к колебанию энергии частиц относительно референс-частицы и соответственно к уменьшению вращения спина.

Описание разработанной среды моделирования. Примененные алгоритмы и методы реализованы в виде интегрированной среды разработки (IDE), содержащей такие элементы как редактор кода с подсветкой синтаксиса, графический дизайнер ускорителя (рис. 3), обозреватель проекта, технология автодополнения кода и всплывающих подсказок. Также поддерживается возможность генерации вычислительного кода на языках C++ и MATLAB.

Рис. 3. Визуальный (а) и текстовый (б) редакторы кода

Для отслеживания ошибок применяется специальный выходной поток. Пользователь получает подробную информацию о ходе проведения расчетов и может реагировать на возникающие ошибки в интерактивном режиме. Используемый в среде моделирования язык программирования имеет С-подобный нетипизированный синтаксис. Также предусмотрена возможность описания процессов в графическом виде, подобно тому, как это осуществляется в таких средах визуального программирования как LabVIEW и МАТЬАВ/ЯтиНпк. Поток вычислительных задач в данном случае дается как схемы-графы. На рис. 4 представлен пример графической нотации, в которой в виде диаграмм выражается последовательность производимых вычислительных действий.

Lattice Designer представляет собой графическое окно визуального редактора последовательности элементов.

Блок Fringe Fields позволяет задать и настроить конфигурацию краевых полей. Имеется возможность указать элементы, которые обладают краевым полем рассеивания, длину поля, закон распределения (например, линейные или в соответствии с функцией Энге).

Solver описывает используемый численный метод интегрирования, а также предоставляет возможности по генерации вычислительного кода на сторонних языках программирования и конвертации начальных данных во входные файлы сторонних пакетов моделирования.

Заключение. В работе описан и реализован алгоритм численного решения дифференциальных уравнений на основе построения нелинейного матричного отображения. Данный подход интегрирован в среду моделирования, которая предоставляет физикам-исследователям гибкий и удобный инструмент проведения расчетов. Корректность алгоритма протестирована в сравнении со сторонней программой численного моделирования COSY Infinity. Отметим, что при использовании традиционных пошаговых методов интегрирования (например, схемы Рунге-Кутта) с каждым новым шагом растет глобальная ошибка конечного решения. Кроме того, такие подходы не сохраняют физических свойств системы (симплектичность и сохранение полной энергии). Существующие модификации [11-13] требуют в своей работе разрешения неявных уравнений, что сказывается на времени вычислений. Попытки распараллеливания таких численных методов не реализуемы на практике из-за сильных зависимостей внутри алгоритма.

Построенная среда моделирования автоматически выполняет все необходимые вычислительные операции и по заданным нелинейным правым частям уравнений строит итоговое отображение в виде набора числовых матриц. Также в программе реализован генератор кода, который конвертирует отображение в вычислительный алгоритм в форме функций на языках программирования высокого уровня. В настоящий момент поддерживаются С++ (с расширением OpenMP), MATLAB (с поддержкой Parallel Computing Toolbox). Сгенерированные исходные коды могут быть при необходимости скомпилированы и выполнены на различных системах высокопроизводительных вычислений. В этом случае от исследователя требуется лишь задание начальных точек. На выходе будут получены конечные состояния системы, вычисленные при помощи найденных отображений. Рассмотренный в настоящей работе метод обладает свойством естественного распараллеливания. Матричный подход дает возможность реализовывать алгоритмы решения дифференциальных уравнений на высокопроизводительных вычислительных системах. Наиболее естественным в данном случае является использование графических процессоров (GPU), которые позволяют отобразить вычислительную модель алгоритма на специализированную конвейерную архитектуру.

Автор благодарит Ю. В. Сеничева за постановку задач и обсуждение результатов моделирования, Д. В. Зюзина за помощь в проведении сравнительных расчетов на программе COSY Infinity, а также своего научного руководителя проф. С. Н. Андрианова.

Литература

1. COSY Infinity // URL: http://www.bt.pa.msu.edu/index_cosy.htm.

2. MAD — Methodical Accelerator Design // URL: http://mad.web.cern.ch/mad.

3. Андрианов С. Н. Динамическое моделирование систем управления пучками частиц. СПб.: Изд-во С.-Петерб. ун-та, 2002. 376 c.

4. Ivanov A., Andrianov S. Matrix formalism for long-term evolution of charged particle and spin dynamics in electrostatic fields // Proc. of ICAP 2012. Rostock, Germany, 2012. P. 187-189.

5. Силадьи М. Электронная и ионная оптика / пер. с англ. И. М. Ахмеджанова и др. М.: Мир, 1990. 639 с. (Szilgyi М. Electron and ion optics.)

6. Balandin V. V., Golubeva N. I. Hamiltonian methods for the study of polarized proton beam dynamics in accelerators and storage rings // URL: http://arxiv.org/pdf/physics/9903032.

7. Канаков О. И., Мильченко Н. А. Симплектические методы интегрирования гамильтоновых систем // Труды XII науч. конференции по радиофизике / под ред. А. В. Якимова, С. М. Грача. Нижний Новгород: Изд-во ТАЛАМ, 2008. С. 87-88.

8. Иванов А. Н. Численная реализация матричного формализма // Процессы управления и устойчивость: Труды 43-й междунар. науч. конференции аспирантов и студентов / под ред. А. С. Ерёмина, Н. В. Смирнова. СПб.: Издат. Дом С.-Петерб. гос. ун-та, 2012. С. 347-353.

9. Oevel W., Sofroniou M. Symplectic Runge-Kutta schemes II: classification of symmetric methods // URL: http://citeseerx.ist.psu.edu/viewdoc/summary? doi=10.1.1.46.5060.

10. Gilani F. Harness the Features of C# to Power Your Scientific Computing Projects // URL: http://msdn.microsoft.com/en-us/magazine/cc163995.aspx.

11. Fasma D., Brigida P. Energy-Preserving Runge-Kutta methods (2010) // URL: http://www.dm. uniba.it/ delbuono/sds10/ LecturePaceSDS10.pdf.

12. Aubry A., Chartier P. Pseudo-symplectic Runge-Kutta methods // URL: http://www.irisa.fr/ ipso/fichiers/hambit.pdf.

13. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи / пер. с англ. И. А. Кульчицкой, С. С. Филиппова. М.: Мир, 1990. 512 c. (Hairer E., Norsett S. P., Wanner G. Solving ordinary differential equations. Nonstiff problems.)

References

1. COSY Infinity // URL: http://www.bt.pa.msu.edu/index_cosy.htm.

2. MAD - Methodical Accelerator Design // URL: http://mad.web.cern.ch/mad.

3. Andrianov S. N. Dinamicheskoe modelirovanie sistem upravlenija puchkami chastic (Dynamical modeling of control for particle beam Systems). St.-Petersburg: Izd-vo S.-Peterb. un-ta, 2004, 376 p.

4. Ivanov A., Andrianov S. Matrix formalism for long-term evolution of charged particle and spin dynamics in electrostatic fields. Proc. of ICAP 2012. Rostock, Germany, 2012, pp. 187-189.

5. Szilgyi М. Jelektronnaja i ionnaja optika (Electron and ion optics). Per. s angl. I. M. Ahmedzhanova e. a. Мoscow: Mir, 1990, 639 p.

6. Balandin V. V., Golubeva N. I. Hamiltonian methods for the study of polarized proton beam dynamics in accelerators and storage rings. URL: http://arxiv.org/pdf/physics/9903032.

7. Kanakov O. I., Milchenko N. A. Simplekticheskie metody integrirovanija gamil'tonovyh sistem (Symplectic methods for Hamiltonian systems integration). Proc. of XII conf. on radiophysics. Nizhniy Novgorod, 2008, pp. 87-88.

8. Ivanov A. N. Chislennaja realizacija matrichnogo formalizma (Numerical implementation of the matrix formalism). Proc. of 43th conf. on control process and stability. Pod red. А. S. Еryominа, N. V. Smirnova. St.-Petersburg: Izdat. Dom S.-Peterb. un-ta, 2012, pp. 347-353.

9. Oevel W., Sofroniou M. Symplectic Runge-Kutta schemes II: classification of symmetric methods. URL: http://citeseerx.ist.psu.edu/viewdoc/summary? doi= 10.1.1.46.5060.

10. Gilani F. Harness the Features of C# to Power Your Scientific Computing Projects. URL: http://msdn.microsoft.com/en-us/magazine/cc163995.aspx.

11. Fasma D., Brigida P. Energy-Preserving Runge-Kutta methods (2010). URL: http://www.dm. uniba.it/ delbuono/sds10/ LecturePaceSDS10.pdf.

12. Aubry A., Chartier P. Pseudo-symplectic Runge-Kutta methods. URL: http://www.irisa.fr/ ipso/fichiers/hambit.pdf.

13. Hairer E., Norsett S. P., Wanner G. Reshenie obyknovennyh differencial'nyh uravnenij. Nezhestkie zadachi (Solving ordinary differential equations. Nonstiff problems). Per. s angl. I. A. Kul'chickoj, S. S. Filippova. Мoscow: Mir, 1990, 512 p.

Статья рекомендована к печати проф. Д. А. Овсянниковым. Статья поступила в редакцию 19 декабря 2013 г.

i Надоели баннеры? Вы всегда можете отключить рекламу.