Содержание к диссертации
Введение
Глава 1. Современное состояние проблем, связанных с тематикой диссертации 9
1.1 Обзор основных тенденций в области параллельных вычислений 10
1.2 Основные численные методы решения уравнений эллиптического типа с точки зрения параллельности 18
1.3 Основные численные методы линейной алгебры и их параллельность 21
1.4 Численные методы решения жестких систем обыкновенных дифференциальных уравнений и их внутренний параллелизм 23
Глава 2. Расчет электрического поля установки RS-20 с использованием кластера из персональных компьютеров 27
2.1 Первоначальная постановка задачи 27
2.2 Математическая модель и выбор численного метода 29
2.3 Реализация последовательного варианта расчетов и внутренний параллелизм алгоритма 34
2.4 Модели организации параллельных вычислений для комплексов с распределенной памятью 37
2.5 Выбор модели организации параллельных вычислений 43
2.6 Численные эксперименты 54
Глава 3. Численное моделирование динамики «медленного» Z-пинча с использованием кластера из персональных компьютером 59
3.1 Постановка задачи 59
3.2 Математическая модель и использованные численные методы 60
3.3 Анализ на наличие внутреннего параллелизма 66
3.4 Выбор модели организации параллельных вычислений 69
3.5 Численные эксперименты 73
Заключение 81
Список литерачуры 82
- Основные численные методы решения уравнений эллиптического типа с точки зрения параллельности
- Численные методы решения жестких систем обыкновенных дифференциальных уравнений и их внутренний параллелизм
- Реализация последовательного варианта расчетов и внутренний параллелизм алгоритма
- Математическая модель и использованные численные методы
Введение к работе
Актуальность темы. Начиная с 70-х годов XX века, во всем мире активно ведутся исследования в области управляемого термоядерного синтеза. Для этих исследований актуальными являются вопросы создания установок, способных генерировать мощные энергетические импульсы, и вопросы поведения вещества в электромагнитном поле при сверхвысоких температурах в таких установках — вопросы динамики плазмы.
В последнее время для обострения мощности в экспериментальных комплексах часто используются плазменные прерыватели тока (ППТ) [1, 2], плазменные размыкатели тока, установки типа z-пинча и др. Необходимость получения все более мощных импульсов неотвратимо приводит к усложнению конструкций создаваемых установок и увеличению их стоимости. Актуальной становится задача математического моделирования элементов их конструкций еще на этапе проектирования.
Для детального исследования динамики плазмы в подобных установках не достаточно только натурного эксперимента — он не дает полной и подробной информации. Здесь сказываются не только очень высокие скорости процессов, протекающих в плазме (их характерные времена составляют десятки и сотни наносекунд). В силу принципиальных причин диагностические части установок позволяют измерять характеристики плазмы лишь в отдельные дискретные моменты времени. Для снятия этих характеристик использу-
єтся лазерное излучение, длина волны которого сравнима с масштабами неоднородности в плазме. Таким образом, для получения полной картины происходящего актуальным является комплексный подход, сочетающий проведение натурных и численных экспериментов.
Наиболее полное соответствие с натурными экспериментами дает численное моделирование на основе сложных математических моделей, включающих в себя описание многих тонких эффектов (например, квантовомеха-ническос описание динамики заселенности уровней). При этом резко возрастает объем выполняемых расчетов. Лет 10-12 назад использовать такие модели не позволяли возможности вычислительной техники, сейчас проведение подобных численных экспериментов требует применения параллельных вычислительных систем.
Для проведения параллельных расчетов активно разрабатываются специальные параллельные численные методы. Однако существует целый ряд ранее апробированных методов и комплексов программ, хорошо зарекомендовавших себя при использовании на обычных последовательных компьютерах. Актуальным становится анализ их внутреннего параллелизма и адаптация для выполнения вычислений на параллельных ЭВМ.
Цели и задачи исследования. Основными целями и задачами настоящей
диссертационной работы являются:
исследование внутреннего параллелизма существующих численных методов и комплекса программ, применяемых при МаТеМаТИ-
ческом моделировании конструкции установок типа 1І11Г и динамики поведения плазмы;
отладка технологии расчета сложных задач математической физики на параллельном кластере из персональных компьютеров (учебный класс МФТИ);
расчет электрических полей для определения конструкции модернизированной установки РС-20 (чипа ППТ);
моделирование динамики плазмы па примере медленного 7--пипча.
Методы исследования. Для анализа численных методов и комплексов программ па наличие внутреннего параллелизма, а также для выбора наиболее эффективной модели оріанизации параллельных вычислений на системе с распределенной памятью (кластер из персональных компьютеров) и оценки эффективности работы параллельного нарианта использовалась система анализа программ на параллельность BERT 77 [29] (в разработке которой принимал участие и автор данной работы). Дли решения уравнения Лапласа при расчете электрических полей в модифицированной установке РС-20 был применен метод Якоби. Для моделирования динамики медленного Z-пинча использована модификация имеющихся последовательных программ (основанных па использовании вариационных методов).
Научная новизна. В результате проведенной работы впервые исследо-нана многомерная динамика медленного z-пинча с учетом подробной модели заселенности электронных уровней. Результаты расчетов дают хорошее соответствие с экспериментальными данными, полученными в Израиле [67 —
6CJ"[- Разработана технология наиболее эффективного распараллеливания программ решения задач математической физики для конкретного параллельного комплекса (анализ соотношения времен выполнения отдельных частей программ и передачи данных с помощью BERT 77) бел предварительной профилировки последовательных вариантов.
Научная и практическая значимость. На основании эскизного проекта автора, полученного в результате детального изучения математической модели, осуществлена модернизация экспериментальной установки РС-20. успешно эксплуатирующейся в PIІЦ «Курчатовский институт».
Результаты, полученные в работе, расширяют представления о физике медленных Z-пинчей. Они позволяю! объяснять ряд экспериментально наблюдаемых явлений (в частности, образование горячих точек на оси пинча) и механизм их возникновения.
Отлажена технология расчетов па кластерных вычислительных комплексах, показана эффективность использования учебных классов в качестве высокопроизводительной параллельной вычислительной системы.
Па основании выполненных исследований открывается возможность применения современных вычислительных технологий при моделировании других физических объектов (сильноточные диоды, плазменные переключатели чока и т. п.).
Апробация работы. Основные результаты работы докладывались на 3 Международных конференциях (30th l?PS Conference on Controlled Fusion
and Plasma Physics, SanebPeierburg, 2003; 8-я Международная конференция
«Математика, компьютер, образование» - Пущи но, 2001 г: Interstate Sympo
sium on Plasma Electronics and New Acceleration Methods - Kharkov, Ukraine,
July 2000) и одной Всероссийской (28 Звенигородская конференция по физи
ке плазмы и управляемому синтезу - 2001 г.). Результаты также докладыва-
I лись на семинарах Отделения прикладной физики РНЦ КИ (2001-2003 г.г.),
на научных семинарах кафедр нычислительной математики и информатики МФТИ (2001-2003 г.г.) и на секции Научно-Технического Совета МинАтома
РФ (декабрь 2001 г.).
I Публикации. По теме диссертации опубликовано 6 работ: 3 статьи в цен-
тральных журналах [74-76] и 3 тезиса выступлений на конференциях [77-79].
Благодарности, гранты. Проведенная работа выполнена при поддержке РФФИ (проект № 98-01-00225) и 1NTAS (проект № 97-0021). Автор выражает огромную благодарность А, С. Кингсепу, А. В.Бартову, Г. И. Долгачеву за постановку задачи и полезные обсуждения, А, В. Богатову за полезные рекомендации при анализе параллельности программы, И, В. Коваленко тл техническое содействие. Особую благодарность автор выражает Ю. Н, Субботину Ї
и А. В. Зырянову за предоставленную возможность использования для расчетов учебного компьютерного класса.
Основные численные методы решения уравнений эллиптического типа с точки зрения параллельности
При расчете электрофизических параметров установок возникнет необходимость численного решения уравнений Лапласа или Пуассона в области сложной геометрии со многими несвязными границами. Кроме того, коэффяциенты в уравнениях могут претерпевать разрывы на границе раздела сред (вакуум- диэлектрик), подробнее см. ниже (глава 2),
Численному решению задач эллиптического типа посвящена обширная литература. Укажем, в частности, монографии [30—32j. В цитируемых работах подробно исследованы классы разностных схем, доказаны устойчивость (и сходимость) о соответствующих нормах для различных задач.
При наличии разрывов в коэффициентах уравнения одним из ключевых понятий является понятие консервативности схемы. Понятие консервативности было впервые введено в [33]. Как покачано в цитируемой работе, консервативность является необходимым и достаточным условием сходимости разностных схем в классе разрывных коэффициентов.
Для построения консервативных разностных схем разработаны эффективные методики, в частности, интегро-интерполяционный метод [34]. В последнее время большую популярность среди вычислителей приобрели вариационные и проекционные методы конечных элементов [35—371- В отличие от конечно-разностных методик, решение, получаемое методом конечных элементов (МЮ) принадлежит функциональному пространству W1 и. следовательно, является обобщенным решением. Полученные методом конечных элементов решения удовлетворяют свойству консервативности, так как строятся на основе использования интегральных тождеств.
Основным недостатком подхода МКЭ является проблема удачной триангуляции расчетной области, что фактически эквивалентно построению не которой адаптивной расчетной сетки. Исследования возможности распараллеливания МКЭ на машинах с распределенной памятью показали относительно низкую эффективность их работы [61].
В последние годы появились работы, связанные с построением эффективных монотонных (мажоритарных) разностных схем на хаотических сетках [38, 39], При построении монотонных разностных схем результирующая система сеточных уравнений имеет оператор с диагональным преобладанием, что приводит к быстрой сходимости итерационных методов решения сеточных уравнений. Методики [38, 39] позволяют учесть наличие разрывов в коэффициентах и получить сгущение сетки к особенностям решения в окрестности входящих углов 40]. При реализации на параллельной вычислительной системе такие методы могут быть эффективно реализованы только на устройстве с общей памятью, так как в противном случае потребуется большое количество межпроцессорных обменов информацией на этапе построения адаптивной сетки. Эти обмены могут свести к нулю выигрыш, достигнутый за счет разделения вычислительной работы между несколькими процессорами.
Для решения эллиптических систем уравнений хороню зарекомендовали себя мпогосеточные методы Р. П. Федоренко [41, 42]. Хотя данный класс методов и допускает эффективную параллельную реализацию, одной из основных идей этих методов является быстрое подавление высокочастотных гармоник на крупной сетке, что делает неэффективным применение мпогосеточных методик в областях сложной геометрии со многими несвязными границами и разрывными коэффициентами уравнений на границе сред. В таких областях могут существовать большие градиенты решений, связанные как с разрывом коэффициентов, так и со свойствами геометрии расчетной области. К другому классу методов относятся методы граничных элементов 43] позволяющие решать эллиптические уравнения с разрывными коэффициентами в областях со сложной геометрией, в том числе и с несвязными границами. Однако использование таких методик, в конечном счете, сводится к получению решений систем линейных уравнений большой размерности с прямоугольными матрицами. Эти решения являются решениями в обобщенном смысле, в частности, в смысле наименьших квадратов. Алгоритмы метода наименьших квадратов или квазиобращения прямоугольных матриц также не допускают эффективного распараллеливания ни вычислительных комплексах с распределенной памятью.
Численные методы решения жестких систем обыкновенных дифференциальных уравнений и их внутренний параллелизм
При моделировании физических процессов нередко возникает необходимость решать жесткие системы обыкновенных дифференциальных уравнений {ОДУ). Причина появления этой «жесткости» заключается в наличии значительной разницы в собственных числах матриц Якоби систем ОДУ, описывающих физические процессы с существенно отличающимися характерными временами. Подооные системы часто встречаются при моделировании процессов в ядерных реакторах, при решении задач радиофизики, астрофизики, физики плазмы, биофизики, химической кинетики- Последние задачи обычно могут быть записаны в виде: где ик - - концентрации веществ, участвующих в химических реакциях, ско-рости протекания которых характеризуются коэффициентами а [42], Другим примером физических процессов, приводящих к жестким системам ОДУ, служат задачи кинетики ионизации плазмы. Уравнения, описывающие динамику ионизации, могут быть также записаны в виде (І). В этом случае ик характеризуют числа заполнения соответствующих электронных уровней, вероятности переходов между которыми характеризуются коэффициентами а .
Для решения жестких систем ОДУ разработаны различные семейства численных методов. Явные методы зачастую оказываются неустойчивыми, либо требуют для устойчивости специального выбора шагов по времени (см. далее), либо приводят к огромному количеству матричных вычислений [49-51,56,57], поэтому наибольшее распространение получили семейства неявных методов.
Среди одношаговых неявных методов для решения жестких систем ОДУ исторически первыми являются методы Рунге—Кутты 49, 50, 52, 531. Эти методы при надлежащем выборе временного шага являются устойчивыми, но для вычисления вспомогательных векторов, входящих в функцию приращения, необходимо решать большие нелинейные системы алгебраических уравнений, что значительно увеличивает время счета.
Более экономичными с точки зрения вычислений являются многозначные методы (методы Гира) f50, 54, 55j. При использовании семейства этих методов требуется решать лишь одно нелинейное уравнение. Однако многозначные методы обладают существенным недостатком — при повышении порядка аппроксимации они не А-устойчивы, а лишь А(а)-устойчивы. При этом с ростом порядка аппроксимации значение а уменьшается [54, 55]. В силу А(а)-устойчивости при применении методов Гира для жестких систем ОДУ, описывающих физические процессы с быстрыми осцилляциями, возможно проявление численной неустойчивости. Так как в задачах кинетики ионизации плазмы имеет место релаксация к равновесному состоянию, то в них нет колебаний,, и наличие лишь А(а)-устойчивости не является ограничением для использования методов Гира.
Все перечисленные методы очень плохо распараллеливаются на вычислительных системах с распределенной памятью. Более экономичными методами, обладающими высоким внутренним параллелизмом с малым обменом информацией между частями программы, работающими на разных процессорах, являются явные методы решения жестких систем ОДУ, рассмотренные в [56, 57J, Однако, как уже упоминалось ранее, они устойчивы не при произвольных значениях временною шага, а лишь на упорядоченном наборе из нулей полинома Чебышева. В том случае, когда решение жестких систем ОДУ является лишь составной частью решения общей задачи, такое требование к временному шагу зачастую оказывается принципиально невыполнимым. При исследовании поведения Z-пинча шаг по времени выбирается исходя из выполнения условия Куранта для гидродинамической части системы уравнений и никак не может быть отображен на множество нулей полипома Чебышева.
В последнее время среди неявных методов широкое распространение получили одиоитерационные методы Розенброка [50]. В -утих методах, в отличие от меюдов Рупге-Купы, требуется решать не нелинейную систему алгебраических уравнений, а линейную того же порядка. Правда, коэффициенты этой системы получаются весьма громоздкими и трудоемкими для вычислений, а сами методы плохо распараллеливаются на системах с распределенной памятью. Методы Розенброка оказываются наиболее эффективными для систем ОДУ с бесконечно большой жесткостью, к которым не относятся системы уравнений, описывающих динамику ионизации плазмы.
Исходя из того, что решение систем ОДУ при моделировании поведения плазмы будет происходить в последовательном режиме работы, для его реализации следует выбирать самые экономичные с точки зрения вычислений методы. Из з тих соображений предпочтительнее всего оказываются методы Гира в представлении Нордсика [50, 551
Реализация последовательного варианта расчетов и внутренний параллелизм алгоритма
Распараллеливание хорошо зарекомендовавших себя последовательных алгоритмов в качестве неотъемлемой составной части работ включает в себя реализацию последовательного варианта расчетов. Этот последовательный комплекс программ в дальнейшем исследуется на параллельность и служит для получения эталонных результатов расчетов для отладки соответствующих параллельных программ. От эффективности реализации последовательного варианта во многом зависит и эффективность работы параллельного варианта.
Для реализации численного алгоритма вся область интегрирования была разделена на 2 части по высоте (на подобласти / и // па рис. 3).
Сначала рассчитывались значения на границе между обеими частями, а затем расчеты велись в каждой подобласти отдельно, с учетом рассчитанных граничных значений.
Подобласть / обладала наиболее сложной структурой. Прямые вычисления потенциала в ней потребовали бы либо проверки множества условий на каждой итерации для выбора правильных значений диэлектрической проницаемости и измененных значений шага сетки в районе компенсаторных колец, либо заведения дополнительных массивов информации больших размеров (- 24 Mb) для хранения этих значений. И то, и другое было признано неприемлемым. Лишние проверки увеличивают и без того оольшое время расчета отдельной итерации, а значительные массивы информации приводят к увеличению объема свопинга на машинах с относительно небольшим размером оперативной памяти, что негативно сказывается на общей производи тельности вычислительной системы.
Ввиду того, что особенности геометрии области / сконцентрированы в узкой полосе вокруг диэлектрической вставки, было принято решение раз бить подобласть / на три подобласти la, ih7 lc и вести расчеты в них по отдельности. Вне згой вставки значения диэлектрической проницаемости равны единице, а область имеет простую форму. Поэтому для учета разрыва в коэффициенте диэлектрической проницаемости в области fh были введены дополнительные массивы информации для хранения значений коэффициента диэлектрической проницаемости и сеточных параметров общим объемом 2,5 Mb.
Отказ от использования дополнительных массивов в областях 1а и 1с и то обстоятельство, что при расчете нового шага итерационного процесса в них постоянны диэлектрическая проницаемость и сеточные параметры, позволили ускори гь время вычисления в этих областях примерно в 1,8 раза.
Помимо разрывного коэффициента диэлектрической проницаемости, в подобласти lb имеются дополнительные (внутренние) границы расчетной области, связанные с присутствием н ней компенсаторных колец. Так как при выбранной схеме вычислений значение потенциала в узле сетки зависит лишь от значений потенциала в его окрестности на предыдущей итерации, то потенциал на текущей итерации в подобласти lb сначала расечиїьівалея без учета компенсаторных колец, а затем корректировался в узких областях вокруг каждого такого кольца (показаны на рис. 3). Это потребовало повторно-10 пересчета потенциала ь 0.,6% узлов сетки подобласти 1Ь по позволило упростить алгоритм вычислений и избежать утомительной проверки выполнения условий сходимости для всех узлов сетки. Время вычислений в подобласти lb вследствие этого незначительно уменьшилось (примерно в 1,05 раза).
Математическая модель и использованные численные методы
Рассмотрим математическую модель и численные методы, примененные для расчетов в последовательном программном комплексе.
В рамках подхода магнитной гидродинамики берется следующая матеро матичсская модель движения одпожидкостнои двухтемпературнои плазмы в магнитном поле, записанная (в одномерном случае) в безразмерных пере менных:
где Q — потери на излучение (в приближении тормозного излучения), О, — обмен энергиями между электронной и ионной компонентами за счет столкновений. В качестве рабочего газа рассматривается чистый углерод. В приведенных выше формулах индекс е относится к электронам, а і — к ионам. Большинство обозначений традиционные Р без индекса означает полное давление - - сумму парциальных давлений электронов и ионов, функция J(z,f) учитывает потери энергии на ионизацию. Основмые масштабы, примененные при обезразмеривании, были выбраны следующими: время [t] = 10 с, длина [і ]- Ї СМ, плотность [р] = 10 г/смЗ, температура (энергия) [Т] - 1 КэВ. В качестве уравнений состояния в модели используются уравнения состояния идеальной плазмы: здесь и далее: А — число нуклонов в ядре (атомный нес). Уравнения модели записаны в безразмерном виде. Для обмена энергиями за счет столкновений а для потерь на излучение имеем:
Для коэффициентов электронной и ионной теплопроводности приняты следующие зависимости: Коэффициент магнитной вязкости выражается в виде Здесь и далее z ;, - эффективный заряд иона плазмы.
На внешней границе короны ставятся условия границы с вакуумом. При этом азимутальное магнитное поле на границе связано с полным током во внешней электрической цепи известным соотношением, записанным в безразмерных переменных: где / — ток во внешней цепи, R — текущий радиус соответствующей точки короны.
Приведенная выше система уравнений в программном комплексе решается методом растепления по физическим процессам.
На первом этапе расчета учитывается движение плазмы, для чего численно решаются уравнения невязкой магнитной гидродинамики без учета диссипативных эффектов. При построении разностной схемы используются лаі ранжевы массовые координаты. Разностная схема во внутренних гонках сетки имеет второй порядок аппроксимации по пространственным переменным. Для уменьшения осцилляции, возникающих при расчете разрывных решений (ударных волн, контактных разрывов или тангенциальных разрывов), к давлению добавляется искусственная вязкость. В расчетах применяется комбинация линейной и квадратичной вязкости J 34-]. Реализованные явные полностью консервативные разностные схемы аналогичны схемам, описанным в [34]. В программе расчета заложены ограничения на минимальный размер ячейки разностной сетки. В том случае, когда пространственный размер какой-либо ячейки становится меньше допустимого, проводится консервативный пересчет всех неличин на новую сетку. При этом иересчитываются не только «традиционные» величины — температуры, количество движения, масса плазмы, но и средний заряд иона в ячейке, потери на ионизацию. Алгоритм такого пересчета достаточно очевиден и мы не будем останавливаться на его описании. Численные методы, применяемые для расчетов, вполне аналогичны методам, применявшимся в работах [70, 71].
Па следующих этапах решаются уравнения электронной и ионной теплопроводности, а также диффузии магнитного поля. Для этого применяются схемы потоковой прогонки [42].
На заключительном этапе учитываются обмены энергиями между электронами и ионами и потери на излучение. Жесткие системы ОДУ, возникающие на данном этапе расщепления, для каждой ячейки решаются методом трапеций.
На каждом гидродинамическом шаге (или через несколько гидродинамических шагов по времени) для каждой лаграпжевоЙ ячейки сетки численно интегрируется (в собственном кинетическом масштабе времени) система 200-300 ОДУ, определяющих заселенности квантовых состояний ачомов и ионов плазмы. Для решения системы ОДУ разработана программа, использующая для численного интегрирования метод Гира [55] и применяющая соответствующий модуль из распространенной библиотеки NAG. Метод Гира является самостоятельно стартующим и имеет переменный порядок аппроксимации (сходимости). В частности, у метода имеется разгонный участок, где порядок сходимости понижается Г50].
