Исследование производительности высокопроизводительных вычислительных систем тема диссертации и автореферата по ВАК РФ 05.13.15, доктор наук Снытников Алексей Владимирович

• Снытников Алексей Владимирович
• доктор наук
• 2019

Оглавление диссертации доктор наук Снытников Алексей Владимирович

Введение

Глава 1. Описание реализации метода частиц в ячейках для

высокопроизводительных ВС

1.1 Краткое описание метода частиц в ячейках

1.2 Модель высокотемпературной бесстолкновительной плазмы

1.3 Решение уравнений Максвелла и Власова

1.4 Параллельная реализация

1.4.1 Декомпозиция расчетной области

1.5 Ход вычислений в программе

1.6 Программная реализация вычислительных методов

1.6.1 Вычисление электромагнитного поля во всей расчетной области

1.6.2 Расчет движения модельной частицы

1.6.3 Общее замечание о методике подсчета количества операций

1.7 Список входных и выходных данных программы

1.7.1 Список входных данных

1.7.2 Список выходных данных

1.7.3 Выходные данные программы-теста

Глава 2. Физико-математические задачи и вычислительные методы в исследованиях, проводимых с

использованием высокопроизводительных ВС

2.1 Масштабируемость и увеличение числа ПЭ

2.2 Адаптация вычислительных методов к архитектуре ВВС

2.3 Использование ускорителей вычислений

Глава 3. Комплексная оценка производительности ВС

3.1 О влиянии организации данных на результат измерения производительности процессоров

3.2 Расчет пропускной способности системы памяти

3.3 Расчет производительности процессорных элементов

3.4 Расчет производительности процессорных элементов на основе

движения модельных частиц

3.4.1 Использование результатов измерения

производительности для улучшения

технико-эксплуатационных показателей ВС

3.5 Расчет пропускной способности коммуникационной сети

3.6 Оценка возможности выполнения крупномасштабных трехмерных расчетов

3.7 Формула для комплексной оценки ВС

3.8 Сравнение с известными тестами производительности

3.8.1 ТестЫРЬ

3.8.2 Актуальность создания нового теста для определения производительности ВВС

3.8.3 Тест ЫРСО

Глава 4. Анализ масштабируемости, параллельной

эффективности и ускорения параллельной ВС

4.1 Формулы для анализа данных о масштабируемости

4.2 Структура временного шага метода частиц в ячейках

4.3 Вычисление характеристик коммуникационного оборудования ВС на основе измеренной масштабируемости метода частиц в ячейках

4.4 Анализ масштабируемости как интегральной характеристики ВС

4.4.1 Описание пересылаемых данных

4.4.2 Вспомогательные величины для анализа масштабируемости

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

4.5 Оценки параллельной масштабируемости на основе измерений времени прохождения сообщений

4.5.1 Краткое описание проведенных тестов

4.5.2 Использование факторного анализа для интерпретации результатов теста

4.6 Опредедение зависимости масштабируемости от наличия и типа ускорителей вычислений

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

основе данных о пересылке модельных частиц

4.8 Сравнение с известными тестами производительности

4.8.1 Краткое описание теста IMB

4.8.2 Сравнение результатов теста IMB и PIC-MANAS

Глава 5. Анализ производительности узлов

мультиархитектурной ВС

5.1 Анализ производительности узлов с графическими ускорителями

5.2 Реализация программы теста на различных типах мультиархитектурных ВВС

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

5.2.2 Универсальная процедура запуска

5.2.3 Унифицированная сигнатура расчетных процедур

5.2.4 Механизм передачи параметров расчетных процедур

5.3 Анализ производительности узлов с многоядерными процессорами и ускорителями вычислений

Заключение

Список сокращений и условных обозначений

Список литературы

Список рисунков

Список таблиц

Приложение А. Основные характеристики процессоров и GPU,

использованных в тестовых расчетах

А.1 Процессоры семейства Intel Xeon

А.1.1 Устаревшие процессоры семейства Intel Xeon

А.1.2 Современные процессоры семейства Intel Xeon

А.2 Процессор Phenom II X6 1055T

А.3 Графические ускорители (GPU)

Приложение Б. Краткое описание основных свойств сетевого

оборудования и сетевых протоколов, использованных в тестовых расчетах

Б.1 Gigabit Ethernet

Б.2 Infiniband

Б.3 Omnipath

Введение

Актуальность темы исследования и степень ее разработанности.

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

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

Это означает необходимость создания специализированных программ-тестов для определения быстродействия конкретной ВС. В настоящее время существует большое количество подобных программ: Linpack, HPCG, NAS Parallel benchmarks, SPEC Benchmarks, ECP Exascale Proxy Applications, Mantevo miniapps benchmarks. Кроме того, существуют оценки параллелльной масшабируемости и эффективности высокопроизводительных ВС, например, в книге [179]. Все это многообразие показывает, что задача создания теста производительности ВС остается актуальной.

Желательно, чтобы тесты создавались на базе программ, используемых для решения реальных задач для того, чтобы результатом работы теста были бы не только числовые показатели, но и конкретное подтверждение способности тестируемой ВС решать сложные задачи, актуальные с научной или производственной точки зрения. Так, тест High Performance Lapack создан на основе методов решения системы линейных уравнений с помощью LU-факторизации,

тест HPCG - на основе метода сопряженных градиентов, NAS parallel benchmark ориентирован на задачи гидро и газодинамики.

Кроме тестирования скорости выполнения арифметических операций очень важное значение имеет производительность системы ввода-вывода. Это проблема, имеющая самостоятельное научное значение, о чем свидетельствует, в частности, доклад «Do you know what your I/O is doing» («Знаете ли вы, что делает ваша система ввода-вывода») [55], в котором была рассмотрена работа системы ввода-вывода на суперЭВМ Blue Waters при решении больших задач, и было показано что во многих случаях именно ввод-вывод является основным фактором, замедляющим решение задач и работу суперЭВМ в целом. В приложении к конкретным задачам математического моделирования время собственно решения задачи на ВВС ( высокопроизводительной вычислительной системе) иногда оказывается меньше, чем время записи результатов на жесткие диски. Таким образом анализ производительности системы ввода-вывода должен быть включен в тест производительности. Одним из немногих тестовых пакетов, решающих эту задачу, является Intel Cluster Checker - набор программ, производящий всестороннее тестирование кластерной системы, и каждой отдельной подсистемы, и вопросы совместного функционирования подсистем.

Для сокращения временных и ресурсных затрат на проведение тестирования необходимо создавать тесты на базе программ, реализующих такие математические методы, которые предъявляют высокие требования сразу к нескольким, в идеале - ко всем подсистемам тестируемой ВС.

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

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

— процессорные ядра,

— графические ускорители,

— контроллеры памяти,

— межпроцессорные коммуникационные интерфейсы,

— контроллеры коммуникационной сети,

— сетевые коммутаторы,

— средства доступа к энергонезависимой подсистеме хранения данных (дисковой подсистеме).

Необходимо, чтобы тест допускал возможность определения реального времени вычислений с учетом того, что данные находятся не в регистрах процессора и должны быть туда загружены. Таким образом, актуальной задачей является определяется реального быстродействия системы на широком спектре задач математического моделирования (sustained performance).

Целью данной работы является разработка методики и создание реализующей ее программы для комплексной оценки производительности ВС с учетом функционирования всех подсистем и выявления скорости решения реальных задач. Таким образом, предполагается представить альтернативу традиционной методике оценки производительности ВВС, основанной только лишь на количестве операций в секунду.

Для достижения поставленной цели необходимо было решить следующие задачи:

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

2. Реализовать метод исследования коммуникационной структуры параллельных распределенных ВС для выработки рекомендаций по более оптимальному распределению процессов в приложении на узлах распределенной ВС.

3. Разработать и обосновать методику расчета абсолютной, не основанной на производительности процессоров, оценки пригодности ВС для решения реальных задач.

4. Разработать метод комплексного анализа производительности узлов мультиархитектурной ВС, позволяющий делать прогнозы эффективности данной мультиархитектурной ВС для решения конкретных задач, более достоверные по сравнению с синтетическими тестами.

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

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

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

Преимущества по сравнению с другими вариантами тестирования заключаются в:

— выявлении производительности ВС на реальном приложении, а не на специально сконструированном тесте - это важно с точки зрения динамики нагрузки на ВС;

— оценке интегральной производительности ВС, складывающейся из производительности процессора, ускорителей вычислений, кэш-памяти разных уровней, оперативной памяти и системы хранения данных;

— возможности выявления в ходе одного теста конкретной подсистемы ВС, наиболее заметно снижающей скорость счета реальных приложений.

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

Впервые проводится анализ производительности узлов мультиархитектур-ной ВС, оснащенных графическими ускорителями, ускорителями Intel Xeon Phi, различными типами многоядерных процессоров с точки зрения выявления их производительности, реально достижимой при решении конкретных научно-технических задач.

На основе разделения метода частиц в ячейках на этапы различной сложности, предъявляющих высокие требования к определенным подсистемам тестируемой ВС создана новая оригинальная методика, позволяющая предсказывать производительность при решении различных научных и народнохозяйственных задач, в частности задач гидро- и газодинамики, методов Монте-Карло, алгоритмов динамической балансировки загрузки, алгоритмов, основанных на интегральных преобразованиях, задач анализа генетической информации, параллельную эффективность и масштабируемость приложений, использующих парные и коллективные коммуникационные операции. Также разработанный тест может использоваться для получения рекомендации по настройке системы параллельных коммуникаций, например, задания оптимальных значений таймаутов MPI.

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

Созданный в рамках диссертационной работы комплексный тест позволяет проводить экспериментальное исследование функционирования вычислительных машин, комплексов и компьютерных сетей в широком диапазоне от больших высокопроизводительных ВС и метакомпьютеров, образованных со-

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

Существует также возможность использования результатов тестирования для разработки методов и алгоритмов организации параллельной и распределенной обработки информации, многомашинных и специальных вычислительных систем. Реализация метода частиц в ячейках в представленной программе включает в себя также и возможность счета на графических ускорителях, и ускорителях Intel Xeon Phi. Таким образом проведенные тесты покажут эффективность вычислений с использованием специализированных вычислителей в рамках анализируемой ВС, а также эффективность работы интерфейса, соединяющего графические и другие ускорители с основной частью ВС.

В силу того, что скорость расчета по методу частиц в целом и эффективность распараллеливания в частности принципиально зависит от топологии внутренней сети высокопроизводительной ВС, а также от параметров коммуникационного оборудования и системного ПО, то результаты тестирования могут быть использованы для разработки структур и топологий компьютерных сетей, сетевых протоколов и служб передачи данных в компьютерных сетях.

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

Основные этапы исследования выполнены при поддержке грантов Российского фонда фундаментальных исследований №№ 13-07-00589,14-07-00241 (руководитель), , 16-07-00434, 16-07-00916, 18-07-00364 (руководитель), а также при финансовой поддержке Министерства образования и науки РФ (уникальный идентификатор работ (проекта) RFMEFI57417X0145, соглашение

№14.574.21.0145 от 26.09.2017г.) в рамках федеральной целевой программы «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2014—2020 годы».

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

Параллельная реализация метода частиц в ячейках может быть выполнена с использованием трех основных вариантов декомпозиции расчетной области (эйлерова декомпозиция, лагранжева декомпозиция и смешанная эйлерово-лагранжева декомпозиция), а также с использованием или без использования мелкозернистого распараллеливания т.е. расчета движения каждой отдельной частицы в отдельной нити CUDA (англ. CUDA thread, поток с ограниченным объемом оперативной памяти и сокращенным набором доступных операций, выполняемые на GPU) или потока OpenMP.

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

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

Подсчет количества операций с плавающей точкой в секунду выполняется в том числе и для реализации вычислительной части алгоритма на мульти-архитектурных ВС, оснащенных многоядерными процессорами, графическими ускорителями, и ускорителями Intel Xeon Phi с использованием нескольких технологий распараллеливания, таких как OpenMP, CUDA, OpenACC, что обеспечивает возможность комплексной оценки производительности многоядерного процессора или вычислительного узла в рамках ВС.

Оценка быстродействия коммуникационной подсистемы ВС производится путем анализа продолжительности парных и коллективных операций обмена данными, используемых при различных вариантах распараллеливания метода частиц в ячейках (т.е. эйлеровом, эйлерово-лагранжевом и лагранжевом вариантах), а также путем анализа скорости прохождения сообщений между различными узлами многопроцессорной ВС.

В диссертационной работе создана оригинальная методика комплексного тестирования производительности ВС.

Основные положения, выносимые на защиту:

1. Создан программный комплекс, основанный на моделировании динамики плазмы методом частиц в ячейках на высокопроизводительных ВС для всестороннего исследования производительности ВС, позволяющий в рамках одного запуска программы определить конкретную подсистему, наиболее заметно снижающую скорость счета реальных приложений [52; 93; 165; 167; 177].

2. Реализован и протестирован метод исследования коммуникационной структуры высокопроизводительных ВС, позволяющий давать рекомендации по более оптимальному распределению процессов в приложении на узлах высокопроизводительной ВС, а также вычислять экстраполяцию реально полученной производительности на аналогичные системы с большим количеством узлов и процессоров [173; 177].

3. Разработана и обоснована методика расчета абсолютной оценки пригодности данной ВС для решения реальных задач, основанной на сбалансированности производительности различных подсистем конкретной ВС, в частности оперативной памяти, коммуникационной сети, дисковой подсистемы, процессоров и ускорителей вычислений и поз-

воляющей сравнивать ВС безотносительно используемых программ и решаемых задач [30; 134; 136; 168; 176; 177].

4. Предложен и протестирован метод комплексного анализа производительности узлов мультиархитектурной ВС, оснащенной многоядерными процессорами и графическими ускорителями или ускорителями Intel Xeon Phi, основанный на программе для моделирования динамики плазмы методом частиц в ячейках и позволяющий делать прогнозы эффективности данной мультиархитектурной ВС для решения конкретных задач, более достоверные по сравнению с синтетическими тестами [128; 161; 165; 171].

Реализованная в диссертационной работе тестовая программа получила условное наименование PIC-MANAS (Particle-In-Cell Multi-Archtitecture Numerical Analysis & Simulation).

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

Введение диссертации (часть автореферата) на тему «Исследование производительности высокопроизводительных вычислительных систем»

Апробация работы.

Результаты диссертационной работы известны научной общественности. Основные результаты диссертационной работы докладывались и обсуждались на Международных научных конференциях в России и за рубежом: на международной научных конференциях серии Parallel Computing Technologies (Нижний Новгород, 2003, Красноярск, 2005, Новосибирск, 2009), международной конференции International Conference on Computational Science (Амстердам, 2009), международной конференции Open Magnetic Systems for Plasma Confinement (Новосибирск, 2010), международной конференция «Параллельные вычисления и задачи управления» (Москва, 2010), международной суперкомпьютерной конференции «Научный сервис в сети Интернет» (Новороссийск, 2009, 2011, 2014), международной научной конференции Russian Supercomputing Days (Москва, 2015, 2016), международной конференции «Супервычисления и математическое моделирование» (Саров, 2016), обсуждались на семинарах в Институте Вычислительной Математики и Математической Геофизики СО РАН, Институте Ядерной Физики СО РАН, Институте Вычислительных Технологий СО РАН, Институте Теоретической и Прикладной Механики СО РАН, Институте Прикладной Математики РАН, Научно-Исследовательском Вычислительном Центре МГУ.

Личный вклад. Выносимые на защиту результаты получены соискателем лично. В совместных работах постановки задач и разработка методов их решения осуществлялись при непосредственном участии соискателя.

Публикации. Основные результаты по теме диссертации изложены в 16 печатных изданиях, 10 из которых изданы в журналах, рекомендованных ВАК; 5 - в международных изданиях, индексируемых Scopus и Web of Science; раздел в монографии.

Основные результаты диссертации опубликованы в работах [30; 52; 93; 128; 134; 136; 159—161; 165; 167; 168; 172; 173; 176; 178].

Объем и структура работы. Диссертация состоит из введения, пяти глав, заключения и двух приложений. Полный объём диссертации составляет 174 страницы, включая 32 рисунка и 50 таблиц. Список литературы содержит 182 наименования.

Глава 1. Описание реализации метода частиц в ячейках для

высокопроизводительных ВС

1.1 Краткое описание метода частиц в ячейках

Основываясь на [163], приведем описание основной идеи метода. Пусть задача записана в абстрактной операторной форме:

I + ^ = 0. (1.1)

Решение задачи (1.1) представляется в виде следующей интерполяционной формулы:

N

Ч = £ (1)). (1.2)

Этот переход называют разбиением среды на модельные частицы. Функция Д(и,у) называется ядром модельной частицы. Представление (1.2) позволяет свести решение задачи (1.1) к интегрированию динамической системы частиц.

Система уравнений движения модельных частиц получается, как указано в книге [163] из уравнения Власова подстановкой функции распределения в виде

т

/(1 ,х,и) = ^^ Я(х — хз (^)Ь(и — из (£)) (1.3)

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

Для удобства будем рассматривать наши уравнения в безразмерных величинах

х ( )

= Щ (*), (1.4)

dхj (£)

~1Г = иг

(И ~

= Е (х,), (1.5)

] = 1,... ,т

- уравнения движения частиц (совпадают с уравнениями характеристик кинетического уравнения Власова).

1.2 Модель высокотемпературной бесстолкновительной плазмы

Математическая модель высокотемпературной бесстолкновительной плазмы представляется кинетическим уравнением Власова и системой уравнений Максвелла [25; 163], которые в безразмерной форме имеют следующий вид:

. г 1-)-|\^/г,е п (л г\

+ + + = 0, (1-6)

(1.7)

(1.8) (1.9)

(1.10)

Здесь индексами i и е помечены величины, относящиеся к ионам и электронам, соответственно; qe = —1, qi = те/т,-; Де(£,г,р) - функция распределения частиц; р^ е, г^е - масса, импульс, положение иона или электрона; E, B - напряженности электрического и магнитного полей. Для перехода к безразмерному виду в качестве единиц используются следующие базовые ве-личины:(скорость света с, масса электрона,

— плотность плазмы п0 = 1014 см —3,

— время t = Ш-—1 где плазменная электронная частота шре = 5,6 • 1011c—1) В начальный момент времени в трехмерной области решения, имеющей

форму прямоугольного параллелепипеда

ж G [0,L], y,z G [0,LJ,

находится плазма, состоящая из электронов и ионов водорода, и пучок электронов. Заданы плотности электронов пучка щ и электронов плазмы пе = 1 — щ. Плотность ионов плазмы в рамках используемой в данной работе модели, подробнее описанной в [167] равна сумме плотностей электронов пучка и электронов плазмы. Температура электронов плазмы Те и пучка Т^ задаются в соответствии с решаемой физической задачей. Температура ионов считается нулевой Ti = 0. Начальное распределение частиц по скоростям максвелловское с плотностью распределения

Ж = — j fB = —rotE,

divE = p,

divB = 0.

= ^ ехр(—,

где Ау - разброс частиц по скоростям (Тъ = Ау2), у0 - средняя скорость пучка. Средняя скорость ионов и электронов фона нулевая. Все частицы распределены по области равномерно, начальная средняя скорость пучка направлена по х и равна у0 = 0.2. Граничные условия периодические.

1.3 Решение уравнений Максвелла и Власова

Решение уравнения Власова проводится методом частиц-в-ячейках [25; 163; 182]. Плазма представляется набором модельных частиц, траекториями движения которых являются характеристики уравнения Власова

% = — (Е + [у,Б]),

дл,

^ = к(Е + [У,Б]),

(1.11)

дп

= V

Для решения системы уравнений (1.11) используется схема с перешаги-

ванием

ш+1/2 то—1/2

Рг 1 — Рг 1

(

= И Е то +

то+1/2 , то—1/2 V + V-

ч 1 ч

Б то

,.то+1 то Г г — Г г = v то+1/2

V г .

Плотность заряда вычисляется по положениям частиц га = (ха,уа,га) с использованием Р1С-ядра

р(г, г) = ^ ^-Я(г, =1

Щх)=1 К (1 — Т) , |Х < К (1.12)

[ 0, |х| > К.

Для решения системы уравнений Максвелла используется метод Лэнг-дона-Лазински, описанный в работе [74], в котором поля определяются из разностных аналогов законов Фарадея и Ампера:

2

-_-= _ гои Е ,

Ет+1_Ет _ го^ Вт+1/2 _ ^+1/2 I1-10;

Значения компонент напряженностей электрического и магнитного полей вычисляются в узлах сеток, смещенных относительно друг друга по времени и пространству:

Вт+1/2 _ (д 1 д2 д3 )т+1/2

В _ (Дг+1/2, 7,/, 7+1/2,/, Пг. п./+1/2^ ,

+1/2, , , +1/2, , , +1/2 Е т _(Е 1 Е2 Е3 )т

Е _ (^г,?+1/2,/+1/2,^г+1/2Л/+1/2, ^г+1/2,^+1/2,/) , •т+1/2 _ ( ■ 1 ■ 2 • 3 )т+1

1 _ (Зг,з+1/211+1/2,-)г+1/21з1 /+1/2,^г+1/2^+1/2,/) ,

рт _ „т

Р _ р г +1/2,.7+1/2,/+1/2°

(1.14)

Разностные операторы го^ и на такой сетке имеют следующий вид (в индексах опустим 1/2)

го^ Н _

/ Нк

1

Н1з-1,1 Н1з,1 _Н1з,1-1

1у Н1з,1

Н1,3,1-1

"Н-1 ,3,1 ни и 1 1 Нг,з-1,1

К ку

\ К К /

(1.15)

и 1 _ ^ 1 7^2 _ 2 и3 _ т^3

Н _ Дг+1,.7',/ Дг,.7_1,/ + ^¿¿+1,/ Дг,.7,/ + +1 (116)

Эта схема имеет порядок аппроксимации 0(т2 + К2). Если в начальный момент времени В0 _ 0, а токи вычисляются таким образом, что выполнено уравнение неразрывности

рт+1 _ рт

р-^ + ^ г+1/2 _ о, (1.17)

Т

то разностные аналоги уравнений (1.9) и (1.10) выполняются автоматически. Для того, чтобы удовлетворить этим условиям, токи вычисляются по схеме, приведенной в

работах [158; 163], В0 _ 0.

Компоненты плотности тока ]

г+1/2 х,1,1—1/2,к—1/2

=

Ах

(1 — 6у )(1 — ) +

■ТО+1/2 = Ах

•1 х,%,1 — 1/2,к+1/2 У т

■ТО+1/2 = пАх

•1 х,%,1+1/2,к—1/2 У т

АуА г 12КУ Кг}

(1 — 6У )6г —

6У (1 — 6,) —

■ то+1/2 = пАх

•1 х,Ц+1/2,к+1/2 У

т

АуАх

12КуКг АуА г"

12КуКг

6„6 7 +--——

у г + 12К„К

,то+1/2 = пАу

<>у;1—1/2,1,к—1/2 У т

(1 — 6х)(1 — 6,) +

АхА г

12КхК ,

А'то+1/2 = пАу

3у,г—1/2,1,к+1/2 Ч т

г+1/2 у,г+1/2,1, к—1/2

А

= Яат

(1 — 6х)6, —

6х(1 — 6г) —

АхА 2'

12 КхК% АхА х

12К Т.*

,то+1/2 = пАу

3у,г+1/2,1,к+1/2 У т

6х6г +

АхАг

12К хК

,то+1/2 = пАг

<1 г,г — 1/2,1—1/2,к У т

(1 — 6х)(1 — 6у) +

,то+1/2 = пАг

<1 г,г—1/2,1 +1/2,к У т

■то+1/2 =

г,¡+1/2,1—1/2,.к У

т

АхАу

12 КхКу

(1 — 6х)6у —

6х(1 — 6у) —

■то+1/2 = пАг

г,г+1/2,1+1/2, к У т

АхАу

12К хК АхАу

12 К х Ку

6х6ь +

АхАу

х 12К хК

х

(1.18)

где Ах = х2 — х1; Ау = у2 — у1; А;г = г2 — г1 - приращение координат, индексом 1 обозначена координата частицы до сдвига, индексом 2 - после; Более подробно эти вычисления описаны в [93; 165; 167]

т

1.4 Параллельная реализация

1.4.1 Декомпозиция расчетной области

Были рассмотрены и реализованы следующие три варианта декомпозиции (рис. 1.1):

Рисунок 1.1: Различные варианты декомпозиции расчетной области. На рисунке показано разделение расчетной области на подобласти вдоль одной из координат. Различные символы (круг, квадрат, звезда, треугольник) означают модельные частицы, находящиеся на разных процессорных элементах. Стрелки означают межпроцессорные коммуникации: прямые стрелки - МР1_Бе^ и МР1_Кееу, изогнутые - МР1_А11ге^ее.

— Лагранжева декомпозиция: разделение расчетной области (пространства моделирования) между процессорными элементами (ПЭ) на части по количеству подвижных элементов (разделяются частицы, независимо от координаты) , схематически лагранжева декомпозиция показана на рис. 1.1.

— Эйлерова декомпозиция: разделение расчетной области (пространства моделирования) между процессорными элементами (ПЭ) на части по

координате (ячейки сетки и находящиеся в них частицы, схематически эйлерова декомпозиция показана на рис. 1.1.

— Смешанная Эйлерово-лагранжева декомпозиция [161]. Расчетная область разделяется на подобласти для решения уравнений Максвелла, и на каждую подобласть назначается группа из М ПЭ, далее частицы каждой подобласти разделяются дополнительно между этими М ПЭ. Этот вариант декомпозиции показан на рис. 1.1.

Эффективность распараллеливания, достигаемая при использовании этих трех вариантов декомпозиции расчетной области показана в таблице 1.

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

Декомпозиция тип коммуникаций эффективность (МСЦ РАН) на 100 ПЭ

эйлерова «точка-точка» 77%

лагранжева коллективные по всем ПЭ 15 %

смешанная «точка-точка» и коллективные по небольшим группам ПЭ 92 %

Для параллельной реализации метода частиц в ячейках в данной работе используется эйлерово-лагранжева декомпозиция расчетной области. Такой вариант выбран из соображений наилучшей масштабируемости вычислительного алгоритма. При этом расчетная область разделяется на подобласти для решения уравнений Максвелла, и на каждую подобласть назначается группа из М ПЭ, далее частицы каждой подобласти разделяются дополнительно между этими М ПЭ. Такой вариант декомпозиции одновременно обеспечивает минимизацию фрагмента вычислительного алгоритма и также соответствует условию линейности алгоритма, необходимому для достижения высокой масштабируемости. Этот подход близок к виртуальным слоям [72].

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

— Для эйлеровой декомпозиции минимальным фрагментом (Понятие минимального фрагмент было введено В.Э.Малышкиным, 2001) является слой сетки вместе с частицами.

— Для лагранжевой декомпозиции минимальный фрагмент - одна частица.

— Тем не менее, для лагранжевой декомпозиции (из-за необходимости осуществлять коллективные коммуникации по всем ПЭ) нарушается критерий линейности алгоритма (В.Э.Малышкин, 1991).

— Эйлерово-лагранжева декомпозиция позволяет одновременно минимизировать размер фрагмента и обеспечить линейность.

1.5 Ход вычислений в программе

Приведем общую схему вычислений в программе:

1. Создание начального распределения модельных частиц (выполняется независимо всеми MPI-процессами) (координаты и импульсы частиц генерируются в оперативной памяти ).

2. Далее на каждом временном шаге:

а) Вычисление полей. Обмен граничными значениями полей между соседними узлами многопроцессорной ВС.

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

в) Выполнение диагностических процедур, вывод результатов на диск. Вывод выполняется стандартными средствами C/C++.

Приведем также описание выходных данных программы. Результатами работы являются (опционально, каждая выдача может быть отключена):

— Списки модельных частиц (координаты, импульсы)

— Плотности электронов, ионов и электронов пучка

— Электрическое и магнитное поля (все три компоненты)

— Ток (три компоненты)

— Одномерная функция распределения электронов (всех, и пучка, и плазмы) по энергии

С точки зрения тестирования ВС основное значение имеют следующие выдачи (отдельно для каждого MPI-процесса):

— Средняя продолжительность операций MPI, причем по отдельности для разных объемов пересылаемых данных и разных стадий вычислительного алгоритма:

— MPI_Send/Recv - пересылка граничных значений поля и пересылка частиц.

— MPI_Allreduce - сложение значений тока в подобласти при лагранжевой декомпозиции.

— Время пересылки частиц, граничных значений поля, время сборки данных при выполнении коллективных операций.

— Время выполнения каждой отдельной стадии вычислительного алгоритма:

— расчет поля,

— пересылка граничных значений,

— расчет движения частиц,

— пересылка частиц в соседние процессоры,

— сборка значений тока при лагранжевой декомпозии

— Время записи файлов на диск

— Время копирования данных на GPU и обратно

Расчет производительности ВС в целом и отдельных ее элементов производится после выполнения программы, на основе полученных времен.

1.6 Программная реализация вычислительных методов

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

Таким образом, цель демонстрации всех перечисленных фрагментов программного кода, заключается в том, чтобы можно было подтвердить оценку количества операций с плавающей точкой, затрачиваемых на одну модельную частицу в течение одного временного шага: 440 операций. Это число будет в дальнейшем использоваться для определения производительности вычислительных систем и в частности, производительности их процессорных ядер во флопсах (FLOPS).

1.6.1 Вычисление электромагнитного поля во всей расчетной

области

На листинге 1.1 показана функция emeElement, выполняющая расчет одной из компонент электрического поля в отдельном узле сетки. Как видно на листинге, для расчета используются: соответствующая компонента тока J и две компоненты магнитного поля H1 и Н2. Вспомогательная функция getGlobalCellNumber используется для того, чтобы переходить от трехмерной нумерации узлов сетки к одномерной (все сеточные величины хранятся в виде одномерных массивов).

Листинг 1.1: Функция, выполняющая расчет одной из компонент электрического поля в отдельном узле сетки.

void emeElement(Cell *c,int3 i,double *E,double *H1, double *H2, double * J,double c1,double c2, double tau, int3 d1,int3 d2 )

5

10

{

int n = c->getGlobalCellNumber(i.x,i.y,i.z); int n1 = c->getGlobalCellNumber(i.x+d1.x,i.y+d1.y,i.z+d1.z); int n2 = c-> getGlobalCellNumber(i.x + d2.x,i.y + d2.y,i.z + d2.z) ;

E [n] += c1*(H1 [n] - H1 [n1]) - c2*(H2[n] - H2 [n2]) - tau*J[n ] ;

}

Непосредственно на листинге 1.1 можно увидеть 7 операций с плавающей точкой. Эта функция вызывается Nx х Ny х Nz раз для каждой из трех компонент электрического поля, всего 21Nx х Ny х Nz операций.

Листинг 1.2: Функция, выполняющая первый этап расчета одной из компонент

магнитного поля в отдельном узле сетки.

void emh1_Element( Cell *c , int3 i ,

double *Q,double *H,double *E1 , double *E2 , double c1,double c2 , int3 d1 , int3 d2)

{

15

20

int n = c->getGlobalCellNumber(i.x,i.y,i.z); int n1 = c->getGlobalCellNumber(i.x+d1.x,i.y+d1.y,i.z+d1.z); int n2 = c->getGlobalCellNumber(i.x+d2.x,i.y+d2.y,i.z+d2.z);

double e1_n1 = E1[n1]; double e1_n = E1[n]; double e2_n2 = E2[n2]; double e2_n = E2[n];

double t = 0.5*(c1*(e1_n1 - e1_n)- c2*(e2_n2 - e2_n)); Q[n] = t; H [n] += Q [n] ;

}

Аналогично, на листинге 1.2 видно 7 операций с плавающей точкой. Эта функция также вызывается Nx х Ny х Nz раз для каждой из трех компонент магнитного поля, всего 21 Nx х Ny х Nz операций.

Листинг 1.3: Функция, выполняющая второй этап расчета одной из компонент

магнитного поля в отдельном узле сетки.

void emh2_Element( Cell *c ,

int i,int l,int k , double *Q,double *H)

5

{

int n = c->getGlobalCellNumber(i,l,k); H [n] += Q [n] ;

}

И листинг 1.3 добавляет еще одну операцию. Эта функция также вызывается Nx х Жу х Nz раз для каждой из трех компонент магнитного поля, всего 3Nx х Жу х Nz операций.

В итоге вычислительная нагрузка при вычислении электромагнитного поля составляет 45^х х ^ х Nz операций.

1.6.2 Расчет движения модельной частицы

Функция расчета движения модельной частицы

В этом разделе на листинге 1.4 показана функция, рассчитывающая движение модельной частицы, то что в зарубежной литературе именуется particle pusher. Следует отметить, что это наиболее важная с точки зрения производительности и наиболее вычислительно нагруженная часть программы. Подробный листинг в данном случае приводится для того, чтобы можно было проверить количество операций с плавающей точкой, выполняемое в при сдвиге одной модельной частицы, которое в дальнейшем будет использоваться для подсчета производительности процессоров

Листинг 1.4: Интегрирование уравнений движения модельной частицы

10

15

void Move(Field fd,double tau)

double tau1,u,v,w,ps; double pu1,pv1,pw1; double3 sx3;

//вычисление изменения импульса частицы //1-я стадия: 21 операция

ElectricMove(fd.E,tau,q_m,&tau1,&pu,&pv,&pw,&ps);

// вычисление поворота вектора импульса // под действием магнитного поля: 51 операция MagneticMove(fd.H,ps ,&pu1 ,&pv1 ,&pw1);

//умножение электрического поля на временной шаг

//3 операции

sx3 = mult(tau1 ,fd.E) ;

//вычисление изменения импульса частицы

//2-я стадия: 3 о перации

add(sx3,&pu,&pv,&pw , pu1 , pv1 , pw1) ;

5

30

35

40

45

50

55

//вычисление модуля 4 - вектора импульса частицы

//10 о пер аций

ps = impulse(pu,pv,pw) ;

//вычисление скорости частицы //3 операции

mult(&u,&v,&w,ps,pu,pv, pw) ;

//вычисление изменения координаты частицы //6 о пер аций

x1 = x + tau * u;

y1 = y + tau * v;

z1 = z + tau * w;

//всего 97 операций

}

void ElectricMove(double3 E,double tau, double q_m,double *tau1 , double *pu,double *pv,double *pw,double *ps)

{

//импульс рассчитывается в два этапа //поэтому уменьшаем шаг вдвое: 2 операции *tau1=q_m*tau*0.5;

//вычисление изменения импульса частицы //6 о пер аций *pu += *tau1*E.x; *pv += *tau1*E.y; *pw += *tau1*E.z;

//вычисление модуля 4 - вектора импульса частицы //8 операций и извлечение квадратного корня (5 операций) *ps = (*tau1) * pow(((*pu) * (*pu) + (*pv) * (*pv) + (*pw) * (* pw)) * 1. + 1.0,-0.5);

//всего 21 операция

}

void MagneticMove(double3 H,double ps,double *pu1,double *pv1, double *pw1)

{

65 double bx , by , bz , s1 , s2 , s3 , s4 , s5 , s6 , s , su , sv , sw ;

//вычисление вспомогательных величин //15 о пер аций bx = ps * H . x; by = ps * H.y; bz = ps * H . z ;

su = pu + pv * bz - pw * by ; sv = pv + pw * bx - pu * bz ; sw = pw + pu * by - pv * bx ;

//матрица на основе компо нент магнитного поля //9 о пер аций s 1 = bx * bx ; s2 = by * by; s3 = bz * bz ; s4 = bx * by ; s5 = by * bz ; s6 = bz * bx ; s = s1 + 1. +s2+s3;

// вычисление поворота вектора импульса // 27 о пер аций

* pu 1 = ((s 1 + 1.) * su + (s4 + bz) * sv + (s 6 - by) * sw) / s;

* pv 1 = ((s4 - bz) * su + (s2 + 1.) * sv + (s5 + bx) * sw) / s;

* p w 1 = ((s 6 + by) * su + (s5 - bx) * sv + (s3 + 1.) * sw) / s;

//всего 51 операция }

95

void add(double3 sx3,double *pu,double *pv,double *pw,double pu1 ,double pv1,double pw1)

{

//сложение веторов: 3 операции *pu = pu1 + sx3.x; 100 *pv = pv 1 + sx3 . y ; *pw = pw1 + sx3.z ;

}

double3 mult(double t,double3 t3)

//умножение вектора на число: 3 операции

70

75

80

85

90

115

t3.x *= t ; t3.y *= t; t3.z *= t ;

return t3 ;

}

//вычисление модуля 4 - вектора импульса //10 о пер аций

double impulse(double pu,double pv,double pw) {

return pow (((pu * pu + pv * pv + pw * pw) + 1.0),-0.5);

}

На листинге 1.4, в конце функции Move приведено суммарное количество всех операций с плавающей точкой, выполняемых в этой функции и в вызываемых из нее: 97 операций. Это значительно меньше, чем приведенная в начале раздела цифра в 440 операций на модельную частицу. Но здесь необходимо учесть еще вычисление полей, величин fd.E и fd.H, используемых в 10-15 строках листинга 1.4.

Функция вычисления электромагнитного поля в точке нахождения

модельной частицы

Расчет поля в точке, в которой находится модельная частица, представляет собой линейную интерполяцию, подробности этого процесса показаны на листинге 1.5

Как видно из листинга 1.5, поля и частицы хранятся внутри ячейки, а не в виде одного глобального массива, содержащего значения поля во всей области. Подход с хранением полей и частиц в виде глобальных массивов также был реализован, и важно отметить, что при этом возникает существенная разница по производительности. Эта разница подробно исследуется в 3-й главе (раздел 3.1) и обычно составляет до 200 %. Таким образом для программы-теста можно использовать как тот, так и другой подход, внося соответствующую поправку в полученные результаты.

Листинг 1.5: Расчет электромагнитного поля в точке, в которой находится модельная частица

Field GetField(Particle *p,CellTotalField cf)

{

int3 i , i 1 ;

double s1 , s2 , s3 , s4,s5 , s6 , s11 , s21 , s31 , s41 , s51 , s61 ; 5 double3 x = p->GetX(); Field fd;

//вычисление интерполяционных коэффициентов //27 о пераций

InverseKernel(x, i , i1 , s1 , s2 , s3 , s4 , s5 , s6 ,

s11 , s21 , s31 , s41 , s51 , s61 ,p) ;

15

20

25

30

//вычисление электрического поля //63 о пер ации

fd.E = GetElectricField(i,i1, s1,s2,s3,s4,s5,s6,

s11 , s21 , s31 , s41 , s51 ,s61 ,p , cf . Ex , cf . Ey , cf . Ez) ; //вычисление электрического поля //63 о пер ации fd.H = GetMagneticField(i,i1, s1,s2,s3,s4,s5,s6,

s11 , s21 , s31 , s41 , s51 , s61 ,p, cf . Hx , cf . Hy , cf .Hz) ; return fd;

//всего по функции GetField //153 о пер ации

}

void InverseKernel(double3 x,int3 & i,int3 & i1,

double& s1,double& s2,double& s3,double& s4,double& s5 ,double& s6,

double& s11,double& s21,double& s31,double& s41,double & s51,double& s61 , Particle *p

)

//вычисление интерполяционных коэффициентов //по трем координатам:

45

50

55

60

65

70

//2 о пер ации

I.x = getCellNumber(х.х,х0,Ьх); //2 о пер ации

И.х = getCellNumberCenter(х.х,х0,Ьх); //1 операция

э1 = s1_interpolate(x.x) ;

//2 о пер ации s2 = s2_interpolate(х . х) ; //всего 7 операций

//V

//7 о пер аций

1.у = getCellNumber(х.у,у0,Ьу); + 1.);

И.у = getCellNumberCenter (х . у , у0 , Ьу) ; ;

s3 = s3_interpolate(х.у); //г

s4 = s4_interpolate(х.у); //г1

//2

//7 о пер аций = getCellNumber (х . z , z0 , Ьz) ;

II.z = getCellNumberCenter(х.z,z0,Ьz); ;

s5 = s5_interpolate(х . z) ; s6 = s6_interpolate(х . z) ;

//(гиг) (в2 + 1.);

//(гиг) ( в2 + 1.5)

//г - я2;

//г1 - 0.5 - э2;

//(гиг) (в2

//(гиг) (э2 + 1.5)

в 2;

■ 0.5 - в2;

//(гиг) (в2 + 1.); //(гиг) (в2 + 1.5)

//г - в2; //г1 - 0.5 - в2;

s11 s21 s31 s41 s51 s61

//6 о пер аций

- s 1

- s2

- s3

- s4

- s5 _

- s6 ;

//всего по функции ТиьегвеКвтивЪ: 27 операций

}

double3 GetElectricField( int3 1,int3 И ,

double& s1гdouble& s2гdouble& s3,double& s4гdouble& s5 ,double& s6,

double& s11гdouble& s21гdouble& s31гdouble& s41гdouble & s51,double& s61 ,

SO

S5

90

95

lOO

lO5

Particle *p,CellDouble *Ex1,CellDouble *Ey1,CellDouble

* Ez 1

double3 E; int3 cell_num;

cell_num.x = i.x; cell_num.y = i1.y; cell_num.z = i1.z; //21 операция

E.x = Interpolate3D(Ex1 , & cell_num , s1 , s 11 , s4,s41, s6 , s61 ,p,0)

cell_num.x = i1.x; cell_num.y = i.y; cell_num.z = i1.z;

E.y = Interpolate3D(Ey1 ,&cell_num , s2 , s21 , s3 , s31 , s6 , s61 ,p , 1)

cell_num.x = i 1.x;//i1; cell_num.y = i1.y;//11; cell_num.z = i.z; //k;

E.z = Interpolate3D(Ez1 , & cell_num,s2,s21,s4,s41, s5 , s51 , p , 2)

return E;

//всего 63 операции

}

double3 GetMagneticField( int3 i,int3 i1 ,

double& s1,double& s2,double& s3,double& s4,double& s5 ,double& s6,

double& s11,double& s21,double& s31,double& s41,double & s51,double& s61 ,

Particle *p,CellDouble *Hx1,CellDouble *Hy1,CellDouble

* Hz 1

double3 H; int3 cn;

cn . x = i 1 . x ; cn . y = i . y ;

cn . z = i . z ;

i . z ;

Interpolate3D(Hx1 ,&cn,s2 , s21 , s3 , s31 ,s5,s51 ,p,3) ;

H . x

cn . x = i . x ;

120

cn . y = i1 . y ;

Таким образом, расчет электромагнитного поля в точке, в которой находится модельная частица выполняется с использованием 153 операций с плавающей точкой.

Необходимо отметить, что есть еще одна не охваченная здесь фаза вычислений, а именно расчет токов по вычисленным координатам модельной частицы

Функция расчета тока по координатам модельной частицы

Расчет тока представляет собой взаимодействие частицы с сеткой, только на этот раз выполняется обратная интерполяция, как показано на листинге 1.6. Здесь важно обратить внимание на то, что количество операций, выполняемых в функции CurrentToMesh, может быть несколько (на 5-10 операций) больше , чем указано на листинге: это зависит от того, перемещается ли частица в соседнюю ячейку, в тексте программы это определяет переменная т. С учетом этого можно считать, что количество операций при вычислении тока по сдвигу частицы составляет 190.

Листинг 1.6: Расчет тока по координатам модельной частицы

//вычисление вклада модельной частицы "p" ток по новым и старым

кооординатам частицы void CurrentToMesh(double tau,int * cells,DoubleCurrentTensor *dt ,Particle *p,int nt,int index)

{

double3 x2;

double s ;

int3 i2,i1 ;

int m , i , l , k ;

double3 x = p->GetX () ;

double3 x1 = p->GetX1();

15

20

25

30

35

40

45

double mass = p->m; double q_m = p->q_m; // Doubl eCurrentTensor dt;

CurrentTensor t;

CurrentTensor *t1 = &(dt->t1); CurrentTensor *t2 = &(dt->t2);

*cells = 1;

//вычисление интерполяционных коэффициентов //9 операций

i1.x=getCellNumberCenterCurrent(x.x,x0,hx); i1.y=getCellNumberCenterCurrent(x.y,y0,hy); i1.z=getCellNumberCenterCurrent(x.z,z0,hz);

//9 операций

i2.x=getCellNumberCenterCurrent(x1.x,x0,hx); i2.y=getCellNumberCenterCurrent(x1.y,y0,hy); i2.z=getCellNumberCenterCurrent(x1.z,z0,hz);

//определение перелета частицы

//из ячейки в ячейку

//3 операции

i=abs(i2.x-i1.x);

l=abs(i2.y-i1.y);

k=abs(i2.z-i1.z);

m=4* i +2*l + k;

//далее в зависимости от того, как

//движется частица, вычисления могут быть разными //рассматриваем один из наиболее //затратных вариантов (m не равен 0) switch (m) {

case 1

case 2

case 3

case 4

case 5

case 6

case 7

goto L1 goto L2 goto L3 goto L4 goto L5 goto L6 goto L7;

}

60

65

70

75

80

85

pqr(i1,x,x1 , mass ,tau,11, 0 , p) ; goto L18; L1 :

/ /4 о перации x2.z = getCellTransitAverage(hz,i1.z,i2.z,z0);

//hz * ( ( i 1. z + i2.z) * .5 - 1.); //3 операции

s = getCellTransitRatio(x1.z,x.z,x2.z) ; // (z2 - z__) / (z 1 - z__); //3 операции x2.x = getRatioBasedX(x1.x,x.x,s); // x + (x1 - x) * s; //3 операции x2.y = getRatioBasedX(x1.y,x.y,s); //y + (y 1 - y) * s;

goto L11 ; L2 :

x2.y = getCellTransitAverage(hy,i1.y,i2.y,y0);

// d_1.h2 * ((11 + 12) * .5 - 1.) + y0; s = getCellTransitRatio(x1.y,x.y,x2.y); // (y2 - y) / (y1 - y); x2.x = getRatioBasedX(x1.x,x . x , s) ;

//x + (x1 - x) * s; x2.z = getRatioBasedX(x1.z,x.z,s); goto L11 ; L3 :

x2.y = getCellTransitAverage(hy,i1.y,i2.y,y0);

//d_1.h2 * ((11 + 12) * .5 - 1.) + y0; x2.z = getCellTransitAverage(hz,i1.z,i2.z,z0); //d_1.h3 * ((k1 + k2) * .5 - 1.);

s = (getCellTransitProduct(x1.z,x.z,x2.z) + getCellTransitProduct(x1.y,x.y,x2.y))

/ (pow(x1.z - x.z,2.0) + pow(x1.y - x.y,2.0)); x2.x = getRatioBasedX(x1.x,x.x,s);

goto L11 ; L4 :

x2.x = getCellTransitAverage(hx,i1.x,i2.x,x0);

90

95

100

105

s = getCellTransitRatio(x1.x,x.x,x2.x);

//s = (x2 - x) / (x1 - x); x2.y = getRatioBasedX(x1.y,x.y,s);

//y2 = y + (y1 - y) * s; x2.z = getRatioBasedX(x1.z,x.z,s);

//z2 = z__ + (z1 - z__) * s;

goto L11 ;

x2.x = getCellTransitAverage(hx,i1.x,i2.x,x0);

//x2 = d_1.h,1 * ((i1 + i2) * .5 - 1.); x2.z = getCellTransitAverage(hz,i1.z,i2.z,z0); //z2 = d_1.h3 * ((k1 + k2) * .5 - 1.);

s = (getCellTransitProduct(x1.z,x.z,x2.z) + getCellTransitProduct(x1.x,x.x,x2.x)) / (pow(x1.z - x.z,2.0) + pow(x1.x - x.x,2.0));

x2.y = getRatioBasedX(x1.y,x.y,s); // y2 = y + (y1 - y) * s; goto L11 ;

x2.x = getCellTransitAverage(hx,i1.x,i2.x,x0); //x2 = d_1.h1

* ((i1 + i2) * .5 - 1.);

x2.y = getCellTransitAverage(hy,i1.y,i2.y,y0); //y2 = d_1.h2

* ((11 + 12) * .5 - 1.) + y0;

s = (getCellTransitProduct(x1.y,x.y,x2.y) + getCellTransitProduct(x1.x,x.x , x2 . x)) / (pow(x1.y - x.y ,2.0) + pow(x1.x - x.x,2.0)); x2.z = getRatioBasedX(x1.z,x.z,s); // z2 = z__ + (z1 - z__) * s;

goto L11 ;

x2.x = getCellTransitAverage(hx,i1.x,i2.x,x0); //x2 = d_1.h1

* ((i1 + i2) * .5 - 1.);

x2.y = getCellTransitAverage(hy,i1.y,i2.y,y0); //y2 = d_1.h2

* ((11 + 12) * .5 - 1.) + y0;

x2.z = getCellTransitAverage(hz,i1.z,i2.z,z0); //z2 = d_1.h3

* ((k1 + k2) * .5 - 1.); L11 :

//собственно вычисление тока: //2 функции по 75 операций pqr(i1, x, x2, mass ,tau,t1 , 0 , p) ; pqr(i2, x2, x1, mass ,tau,t2 , 1,p) ;

120

* cells = 2; L18: p->x = p->x1 ;

p->y = p->y1;

p->z = p->z1;

Reflect(p) ;

Таким образом общее количество операций с плавающей точкой, выполняемых при действиях с одной модельной частицей:

— вычисление поля в точке нахождения модельной частиц: 153 операции

— расчет движения модельной частицы: 97 операций

— расчета тока по координатам модельной частицы: 190 операций или всего 440 операций.

1.6.3 Общее замечание о методике подсчета количества операций

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

Список литературы диссертационного исследования доктор наук Снытников Алексей Владимирович, 2019 год

//г - я2;

//г1 - 0.5 - э2;

//(гиг) (в2

//(гиг) (э2 + 1.5)

в 2;

■ 0.5 - в2;

//(гиг) (в2 + 1.); //(гиг) (в2 + 1.5)

//г - в2; //г1 - 0.5 - в2;

s11 s21 s31 s41 s51 s61

//6 о пер аций

- s 1

- s2

- s3

- s4

- s5 _

- s6 ;

//всего по функции ТиьегвеКвтивЪ: 27 операций

}

double3 GetElectricField( int3 1,int3 И ,

double& s1гdouble& s2гdouble& s3,double& s4гdouble& s5 ,double& s6,

double& s11гdouble& s21гdouble& s31гdouble& s41гdouble & s51,double& s61 ,

SO

S5

90

95

lOO

lO5

Particle *p,CellDouble *Ex1,CellDouble *Ey1,CellDouble

* Ez 1

double3 E; int3 cell_num;

cell_num.x = i.x; cell_num.y = i1.y; cell_num.z = i1.z; //21 операция

E.x = Interpolate3D(Ex1 , & cell_num , s1 , s 11 , s4,s41, s6 , s61 ,p,0)

cell_num.x = i1.x; cell_num.y = i.y; cell_num.z = i1.z;

E.y = Interpolate3D(Ey1 ,&cell_num , s2 , s21 , s3 , s31 , s6 , s61 ,p , 1)

cell_num.x = i 1.x;//i1; cell_num.y = i1.y;//11; cell_num.z = i.z; //k;

E.z = Interpolate3D(Ez1 , & cell_num,s2,s21,s4,s41, s5 , s51 , p , 2)

return E;

//всего 63 операции

}

double3 GetMagneticField( int3 i,int3 i1 ,

double& s1,double& s2,double& s3,double& s4,double& s5 ,double& s6,

double& s11,double& s21,double& s31,double& s41,double & s51,double& s61 ,

Particle *p,CellDouble *Hx1,CellDouble *Hy1,CellDouble

* Hz 1

double3 H; int3 cn;

cn . x = i 1 . x ; cn . y = i . y ;

cn . z = i . z ;

i . z ;

Interpolate3D(Hx1 ,&cn,s2 , s21 , s3 , s31 ,s5,s51 ,p,3) ;

H . x

cn . x = i . x ;

120

cn . y = i1 . y ;

Таким образом, расчет электромагнитного поля в точке, в которой находится модельная частица выполняется с использованием 153 операций с плавающей точкой.

Необходимо отметить, что есть еще одна не охваченная здесь фаза вычислений, а именно расчет токов по вычисленным координатам модельной частицы

Функция расчета тока по координатам модельной частицы

Расчет тока представляет собой взаимодействие частицы с сеткой, только на этот раз выполняется обратная интерполяция, как показано на листинге 1.6. Здесь важно обратить внимание на то, что количество операций, выполняемых в функции CurrentToMesh, может быть несколько (на 5-10 операций) больше , чем указано на листинге: это зависит от того, перемещается ли частица в соседнюю ячейку, в тексте программы это определяет переменная т. С учетом этого можно считать, что количество операций при вычислении тока по сдвигу частицы составляет 190.

Листинг 1.6: Расчет тока по координатам модельной частицы

//вычисление вклада модельной частицы "p" ток по новым и старым

кооординатам частицы void CurrentToMesh(double tau,int * cells,DoubleCurrentTensor *dt ,Particle *p,int nt,int index)

{

double3 x2;

double s ;

int3 i2,i1 ;

int m , i , l , k ;

double3 x = p->GetX () ;

double3 x1 = p->GetX1();

15

20

25

30

35

40

45

double mass = p->m; double q_m = p->q_m; // Doubl eCurrentTensor dt;

CurrentTensor t;

CurrentTensor *t1 = &(dt->t1); CurrentTensor *t2 = &(dt->t2);

*cells = 1;

//вычисление интерполяционных коэффициентов //9 операций

i1.x=getCellNumberCenterCurrent(x.x,x0,hx); i1.y=getCellNumberCenterCurrent(x.y,y0,hy); i1.z=getCellNumberCenterCurrent(x.z,z0,hz);

//9 операций

i2.x=getCellNumberCenterCurrent(x1.x,x0,hx); i2.y=getCellNumberCenterCurrent(x1.y,y0,hy); i2.z=getCellNumberCenterCurrent(x1.z,z0,hz);

//определение перелета частицы

//из ячейки в ячейку

//3 операции

i=abs(i2.x-i1.x);

l=abs(i2.y-i1.y);

k=abs(i2.z-i1.z);

m=4* i +2*l + k;

//далее в зависимости от того, как

//движется частица, вычисления могут быть разными //рассматриваем один из наиболее //затратных вариантов (m не равен 0) switch (m) {

case 1

case 2

case 3

case 4

case 5

case 6

case 7

goto L1 goto L2 goto L3 goto L4 goto L5 goto L6 goto L7;

}

60

65

70

75

80

85

pqr(i1,x,x1 , mass ,tau,11, 0 , p) ; goto L18; L1 :

/ /4 о перации x2.z = getCellTransitAverage(hz,i1.z,i2.z,z0);

//hz * ( ( i 1. z + i2.z) * .5 - 1.); //3 операции

s = getCellTransitRatio(x1.z,x.z,x2.z) ; // (z2 - z__) / (z 1 - z__); //3 операции x2.x = getRatioBasedX(x1.x,x.x,s); // x + (x1 - x) * s; //3 операции x2.y = getRatioBasedX(x1.y,x.y,s); //y + (y 1 - y) * s;

goto L11 ; L2 :

x2.y = getCellTransitAverage(hy,i1.y,i2.y,y0);

// d_1.h2 * ((11 + 12) * .5 - 1.) + y0; s = getCellTransitRatio(x1.y,x.y,x2.y); // (y2 - y) / (y1 - y); x2.x = getRatioBasedX(x1.x,x . x , s) ;

//x + (x1 - x) * s; x2.z = getRatioBasedX(x1.z,x.z,s); goto L11 ; L3 :

x2.y = getCellTransitAverage(hy,i1.y,i2.y,y0);

//d_1.h2 * ((11 + 12) * .5 - 1.) + y0; x2.z = getCellTransitAverage(hz,i1.z,i2.z,z0); //d_1.h3 * ((k1 + k2) * .5 - 1.);

s = (getCellTransitProduct(x1.z,x.z,x2.z) + getCellTransitProduct(x1.y,x.y,x2.y))

/ (pow(x1.z - x.z,2.0) + pow(x1.y - x.y,2.0)); x2.x = getRatioBasedX(x1.x,x.x,s);

goto L11 ; L4 :

x2.x = getCellTransitAverage(hx,i1.x,i2.x,x0);

90

95

100

105

s = getCellTransitRatio(x1.x,x.x,x2.x);

//s = (x2 - x) / (x1 - x); x2.y = getRatioBasedX(x1.y,x.y,s);

//y2 = y + (y1 - y) * s; x2.z = getRatioBasedX(x1.z,x.z,s);

//z2 = z__ + (z1 - z__) * s;

goto L11 ;

x2.x = getCellTransitAverage(hx,i1.x,i2.x,x0);

//x2 = d_1.h,1 * ((i1 + i2) * .5 - 1.); x2.z = getCellTransitAverage(hz,i1.z,i2.z,z0); //z2 = d_1.h3 * ((k1 + k2) * .5 - 1.);

s = (getCellTransitProduct(x1.z,x.z,x2.z) + getCellTransitProduct(x1.x,x.x,x2.x)) / (pow(x1.z - x.z,2.0) + pow(x1.x - x.x,2.0));

x2.y = getRatioBasedX(x1.y,x.y,s); // y2 = y + (y1 - y) * s; goto L11 ;

x2.x = getCellTransitAverage(hx,i1.x,i2.x,x0); //x2 = d_1.h1

* ((i1 + i2) * .5 - 1.);

x2.y = getCellTransitAverage(hy,i1.y,i2.y,y0); //y2 = d_1.h2

* ((11 + 12) * .5 - 1.) + y0;

s = (getCellTransitProduct(x1.y,x.y,x2.y) + getCellTransitProduct(x1.x,x.x , x2 . x)) / (pow(x1.y - x.y ,2.0) + pow(x1.x - x.x,2.0)); x2.z = getRatioBasedX(x1.z,x.z,s); // z2 = z__ + (z1 - z__) * s;

goto L11 ;

x2.x = getCellTransitAverage(hx,i1.x,i2.x,x0); //x2 = d_1.h1

* ((i1 + i2) * .5 - 1.);

x2.y = getCellTransitAverage(hy,i1.y,i2.y,y0); //y2 = d_1.h2

* ((11 + 12) * .5 - 1.) + y0;

x2.z = getCellTransitAverage(hz,i1.z,i2.z,z0); //z2 = d_1.h3

* ((k1 + k2) * .5 - 1.); L11 :

//собственно вычисление тока: //2 функции по 75 операций pqr(i1, x, x2, mass ,tau,t1 , 0 , p) ; pqr(i2, x2, x1, mass ,tau,t2 , 1,p) ;

120

* cells = 2; L18: p->x = p->x1 ;

p->y = p->y1;

p->z = p->z1;

Reflect(p) ;

Таким образом общее количество операций с плавающей точкой, выполняемых при действиях с одной модельной частицей:

— вычисление поля в точке нахождения модельной частиц: 153 операции

— расчет движения модельной частицы: 97 операций

— расчета тока по координатам модельной частицы: 190 операций или всего 440 операций.

1.6.3 Общее замечание о методике подсчета количества операций

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

В том, что касается оптимизации средствами компилятора, в рамках данной диссертационной работы было проведено большое количество экспериментов с различными опциями компилятора и вариантами автоматической оптимизации, как-то протяжка констант, раскрутка циклов, автоматическая векторизация и пр., опубликованных в [128; 161; 165]. Эти эксперименты дают полную уверенность, что автоматическая оптимизация программы применительно к методу частиц не дает почти ничего, к сокращению времени счета может привести только переработка самого вычислительного алгоритма. Таким образом, уменьшится в результате работы компилятора названное число не может, что же касается дополнительных операций, добавленных компилятором, то их учет позволит формально увеличить выполняемую работу и повысить

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

1.7 Список входных и выходных данных программы 1.7.1 Список входных данных

Входные данные:

— Размер сетки

— Количество процессоров

— Количество частиц в ячейке

— Размер области

— Размер части области, занятой пучком и плазмой

— Внешнее магнитное поле

— Энергия частиц пучка

— Поперечная температура частиц пучка

— Температура плазмы (X,Y,Z)

1.7.2 Список выходных данных

Выходные данные (опционально, каждая выдача может быть отключена):

— Списки модельных частиц (координаты, импульсы)

— Плотности электронов, ионов и электронов пучка

— Электрическое и магнитное поля (все три компоненты)

— Ток (три компоненты)

— Одномерная функция распределения электронов (всех, и пучка, и плазмы) по энергии

1.7.3 Выходные данные программы-теста

С точки зрения тестирования ВС основное значение имеют следующие выдачи (отдельно для каждого MPI-процесса):

1. Средняя продолжительность операций MPI, причем по отдельности для разных объемов пересылаемых данных и разных стадий вычислительного алгоритма:

— MPI_Send/Recv - пересылка граничных значений поля и пересылка частиц;

— MPI_Allreduce - сложение значений тока в подобласти при лагранжевой декомпозиции;

2. Время пересылки частиц, граничных значений поля, время сборки данных при выполнении коллективных операций;

3. Время выполнения каждой отдельной стадии вычислительного алгоритма:

— расчет поля;

— пересылка граничных значение;

— расчет движения частиц;

— пересылка частиц в соседние процессоры;

— сборка значений тока при лагранжевой декомпозиции;

4. Время записи файлов на диск;

5. Время копирования данных на GPU и обратно.

Глава 2. Физико-математические задачи и вычислительные методы в исследованиях, проводимых с использованием высокопроизводительных ВС

В последние годы (с 2010 по 2016 год) было опубликовано много статей (более 300) в различных международных журналах по тематике «экзафлопсные вычисления» (т.е. вычисления на перспективных ВВС производительностью порядка 1018 операций с плавающей точкой в секунду), из них около трети посвящены решению различных прикладных задач. Данные обзора отдельными частями опубликованы в [30; 52; 93; 128; 134; 136; 159—161; 165; 167; 168; 171; 173; 176; 177].

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

Таблица 2: Распределение по приложениям статей относящихся к тематике «экзафлопсные вычисления»

вычислительная гидродинамика 27

ядерные технологии 17

физика плазмы 11

разработка новых материалов 10

предсказание погоды 10

биомедицинские приложения 9

астрофизика и космология 6

молекулярная динамика 6

мультифизика 6

геофизика 6

финансы 2

Как видно из таблицы 2 особенно активно математическое моделирование используется в следующих областях науки и техники.

— Гидродинамика (космическая газодинамика, расчет обтекания кораблей и самолетов, например, из числа статей, относящихся к пета- и экзафлопсным вычислениям, [99],[80; 108], предсказание наводнений в прибрежных городах Северной Европы и другие экологические вопросы, моделирование океанических течений, [7; 61; 112; 114], расчет гидродинамической турбулентности [1; 5; 41; 95; 143] моделирование многофазных течений [119; 153], построение сеток для гидродинамического моделирования. [100; 124; 150; 151]),

— Ядерные технологии, точнее моделирование различных процессов, протекающих в ядерных реакторах. [6; 13; 19; 20; 34; 37; 54; 84; 94; 101; 104; 106; 115; 120; 129; 138; 141].

— Разработка новых материалов (квантовохимические расчеты структуры и поведения молекул, вычисление супер-многомерных интегралов в различных приближениях [140]), методы молекулярной динамики [46; 97; 107; 121; 152]

— Моделирование ускорителей заряженных частиц [125].

— Плазменные технологии, например[24].

— Плазменные процессы в установках управляемого термоядерного синтеза [8; 39; 60; 83; 126; 147].

— Разработка новых лекарств и другие биомедицинские приложения [18; 22; 26; 67; 87; 110; 132] (наибольшее количество расчетов с использованием большого количества процессоров, производилось именно в этой области).

— Также математические модели для пета-и экзафлопсных ВВС создаются в геофизике [2; 29; 50; 62; 88; 113].

С учетом того, что во многих областях приходится иметь дело с несколькими различными и одновременно протекающими физическими процессами при этом эти различные процессы моделируются совершенно различными способами, существует необходимость разработки высокомасштабируемого программного обеспечения для таких приложений, которые в зарубежной литературе называются мультифизическими (шиНлрЬузюз). Это делается в работах [11; 36; 48; 77; 103; 137].

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

видно из таблицы 2, но это очень обширная и разнообразная область). В основном разрабатываются масштабируемые версии методов Монте-Карло для переноса нейтронов [20; 34; 37; 54; 101; 104; 115; 138], проводится анализ коммуникационных потерь в методе Монте-Карло [6; 13; 84]. Также рассматриваются вопросы взаимодействия пучков частиц с материалами [19] и применение решеточных уравнений Больцмана с использованием квантовой хромодинамики к задачам ядерной физики [94; 120]. Следует особенно отметить работы по созданию интегральной среды для разработки ядерных реакторов [106] и моделирование новых видов реакторного топлива [129].

В частности, одной из важных областей приложения созданных методов оказывается сейсмология, при использовании метода частиц [163; 182] как альтернативы решению уравнения эйконала [47].

Плазменные задачи представляют большой интерес как с точки зрения выбора и обоснования модели (равновесная плазма или неравновесная, устойчивая или неустойчивая, магнитогидродинамическое описание или кинетика), так и с точки зрения вычислительных методов (решение уравнения Пуассона или уравнений Максвелла, выбор одного из многих вариантов решения кинетического уравнения или уравнений МГД), и в особенности с точки зрения разработки и программной реализации вычислительных алгоритмов (различные варианты декомпозиции расчетной области, организации межпроцессорных обменов, достижения оптимальной производительности), т.е. плазменные задачи требуют поиска новых решений на всех уровнях триады Академика А.А.Самарского: «модель, алгоритм, программа» [174].

Все перечисленные вопросы в том или ином виде встречаются и в других областях математического моделирования, тем не менее в физике плазмы есть проблемы, не имеющие аналога по степени сложности (и соответственно, по количеству необходимых для решения вычислительных ресурсов). Эти проблемы возникают при моделировании плазменных неустойчивостей и турбулентностей в установках управляемого термоядерного синтеза (УТС) при высоких температурах (1-5 КэВ, что соответствует 10 млн.градусов), когда плазма не является даже приближенно равновесной в сколь угодно малой окрестности. В этом случае плазма имеет очень большое количество степеней свободы [162], необходимость учета которых приводит к использованию соответствующих больших объемов памяти.

Особенная важность решения задачи о динамике плазмы в установках УТС связана с тем, что управляемый термоядерный синтез является единственным вариантом решения энергетической проблемы. Все существующие альтернативные энергетические концепции имеют очень серьезные недостатки, не позволяющие рассматривать их как решение энергетической проблемы в мировом масштабе: углеводородная энергетика имеет ограниченные запасы топлива, атомная энергетика представляется экологически опасной, электростанции на солнечной энергии не могут строиться из-за (возможно, временной) дороговизны и экологических проблем, возникающих при производстве фотоэлементов в больших количествах, а различного рода «экологически чистые» электростанции, такие как ветряные, приливные, геотермальные не могут строиться в необходимых больших количествах, также, как и крупные ГЭС.

Другая важная область приложения моделирования плазмы - это промышленные плазменные технологии, используемые, в частности, для производства чистого кремния для нужд микроэлектронной промышленности и для производства микропроцессоров (установки молекулярно-лучевой эпитаксии). Также моделирование высокотемпературной плазмы может быть полезно для исследования природы солнечных вспышек и исследования их влияния на магнитосферу Земли и арктическую навигацию [117].

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

На конференциях по высокопроизводительным вычислениям, например, «Научный сервис в сети Интернет: экзафлопсное будущее», Абау-Дюрсо, 2011 год широко обсуждались вопросы, связанные с построением экзафлопс-компью-теров и трудностями создания программного обеспечения для них. В первую очередь это чрезвычайно большое ожидаемое энергопотребление - порядка 1 ГВт. Таким образом, на первый план выходят вопросы эффективной эксплуатации такой исключительно дорогой машины, а энергоэффективность алгоритма становится определяющей чертой для допуска этого алгоритма к расчетам.

Ограничивая мир высокопроизводительных вычислений первой десяткой наиболее мощных ВС мира, а это вполне оправдано в том случае, если речь идет о возможной архитектуре экзафлопс-компьютера, можно сказать, что

этот мир весьма неоднороден. Серьезно отличаются базовые элементы, на основе которых строятся ВВС, как процессоры (Intel Xeon, AMD Opteron, Sun UltraSPARC), так и ускорители вычислений (Intel Xeon Phi, Nvidia Kepler).

С учетом того, что наиболее мощные компьютеры строятся (вероятнее всего и дальше будут строиться [80]) на основе ускорителей вычислений [180], для эффективного использования компьютеров большой мощности с целью решения актуальных физических задач принципиально необходимо значительно упростить освоение гибридных (т.е. построенных с использованием как процессоров общего назначения, так и ускорителей вычислений) ВВС с точки зрения специалиста по моделированию в конкретной предметной области.

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

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

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

Одним из наиболее эффективных принципиальных подходов к созданию эффективных параллельных алгоритмов является сборочная технология параллельного программирования [72; 82; 157; 166], согласно которой параллельная программа собирается из множества готовых атомарных фрагментов, содержащих как данные, так и алгоритмы.

Существует множество подходов, используемых для достижения более высокой производительности:

— декомпозиция расчетной области[27; 72; 75; 130], с целью ее разделения на более мелкие подобласти,

— использование большого количества процессорных элементов (разработка системного программного обеспечения, позволяющего увеличить число ПЭ сверх определенного предела [56] или создание численных методов с меньшим количеством необходимых коммуникаций),

— применение многоядерных процессоров или различных ускорителей вычислений [181], позволяющих увеличить скорость счета без возрастания нагрузки на коммуникационную сеть, оптимизация под архитектуру процессора или ускорителя,

— сокращение количества отправляемых сообщений,

— переход к численным методам, более подходящим для многопроцессорных ЭВМ (от неявных схем к явным, к методам на основе ортогональных преобразований; от методов, рассматривающих объект моделирования как сплошную среду к методам частиц и пр.) или численных методов специально подобранных или специально оптимизированных для решения конкретной задачи.

Однако все перечисленные подходы можно разделить на три принципиально разных непересекающихся направления:

— увеличение числа ПЭ (процессорных элементов);

— адаптация вычислительных методов к архитектуре ВВС;

— использование ускорителей вычислений.

С этой точки зрения был проведен обзор недавних (2010-2016 годы) публикаций, относящихся к тематике «экзафлопсные вычисления» в журналах Future Generation Computer Systems, Procedia Computers Science, Journal of Parallel and Distributed Computing, Parallel Computing, Journal of Computational Physics, Computer Physics Communications и др. По трем указанным выше направлениям они распределились, как показано на рисунке 2.1 и в таблице 3.

0масштабируемость и увеличение числа ПЭ ■ использование ускорителей вычислений □ адаптация вычислительных методов к архитектуре ВВС

Рисунок 2.1: Распределение по направлениям статей относящихся к тематике «экзафлопсные вычисления»

Таблица 3: Распределение по направлениям статей относящихся к тематике «экзафлопсные вычисления»

Условное наименование направления Число статей В % от общего числа

масштабируемость и увеличение числа ПЭ 147 49

адаптация вычислительных методов к архитектуре ВС 88 30

использование ускорителей вычислений 63 21

2.1 Масштабируемость и увеличение числа ПЭ

Первое направление, названное условно «масштабируемость и увеличение числа ПЭ» было, в свою очередь, разделено на несколько тем, как показано в таблице 4.

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

Таблица 4: Распределение по темам статей относящихся к направлению «масштабируемость и увеличение числа ПЭ»

Тема Количество В %

статей в обзоре от общего числа (от 298 статей)

Масштабируемые 27 7.04

архитектуры

Программные средства

для повышения

масштабируемости 18 6.04

Ускорители и

специальное оборудование

для повышения

масштабируемости 14 1.37

Моделирование

функционирования

экзафлопсных систем 19 6.38

Упрощение работы с

экзафлопсными системами 11 3.69

Декомпозиция

расчетной области 10 3.36

Динамическая

балансировка загрузки 24 8.05

Отказоустойчивость 20 6.71

Уменьшение

количества коммуникаций 10 3.36

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

При этом значение имеет только лишь сам факт наличия публикаций по некоторой теме, количество их не имеет принципиального значения, если только оно не очень велико, в силу того что статистическая погрешность приведенных показателей в процента равна 1Д/3ДО « 6%. Таким образом лишь в том случае можно говорить, что по некоторой теме опубликовано больше статей, чем по какой-то другой, когда разница по количеству статей не менее чем 18-20 работ.

Это же относится и к другим таблицам 2-й главы: важен во-первых, сам факт упоминания некоторой тематики в исследованиях по достижению экзафлопсной производительности, и во вторых, имеет значение лишь значительная разница по количеству статей, написанных на ту или другую тему (более чем на 6 %).

К теме «масштабируемые архитектуры» были отнесены такие работы, как, например, алгоритм привязки процессов к узлам на коммуникационной сети ВВС в виде трехмерного тора [123], метаматематические основы проектирования бесконечно расширяемых машин на основе идеальной модельной архитектуры [14], опыт эксплуатации и решения задач на машине K computer (№ 4 в списке Top500 за ноябрь 2015), имеющей хорошие характеристики масштабируемости и, в отличие от № 1 и № 2, построенной не на основе ускорителей вычислений [137]. Также к этой теме отнесена работа, посвященная использованию большого количества ядер в программах, написанных для старых компьютерных архитектур: несмотря на то, что работа посвящена программированию, основные вопросы, которые там решаются - это архитектурные вопросы, точнее преодоление архитектурных различий программными средствами [79].

К теме «программные средства для повышения масштабируемости» отнесены разработки коммуникационного программного обеспечения, например, MPI, специально предназначенного для работы с очень большим количеством процессов [154], или механизмы снижения задержки сообщений в MPI [71]. Также к этой теме относится создание высокомасштабируемых пакетов программ для решения конкретных задач [97; 122].

Следующая тема, «Ускорители и специальное оборудование для повышения масштабируемости», образована такими работами, как исследование масштабируемого ввода-вывода на кластере из ПЛИС (программируемых логических интегральных схем) [68]. Следует отметить, что к этой теме не относятся вопросы конструирования кластеров на ПЛИС и соответствующие проблемы программирования, а только лишь влияние необычного оборудования (в

частности ПЛИС) на масштабируемость. Кроме того, к этой теме отнесены исследование поведения закона Амдала [12] на многоядерных системах [149] и разработка специализированного сетевого оборудования типа «сеть на кристалле» [9], а также создание сопроцессоров для диспетчеризации потоков [51].

В связи с отсутствием реально функционирующих ВС экзафлопсного класса, а также в связи с предполагаемой исключительной дороговизной вычислительного времени на этих системах, особенную важность приобретает тема «Моделирование функционирования экзафлопсных систем», включающая в себя также и моделирование работы прикладных программ на этих системах. К этой теме относятся анализ расхода времени на коммуникации [13] и подбор с помощью моделирования оптимального (по коммуникациям) варианта декомпозиции расчетной области [144], а также моделирование и предсказание аварийных ситуаций [28].

В докладе профессора Н.Н.Миренкова на конференции Parallel Computing Technologies-2009 (PaCT-2009) [102] была высказана точка зрения о том, что на новые поколения ВВС (гигафлопсные, терафлопсные, петафлопсные,...) раз за разом возлагались очень серьезные надежды на прогресс в решении реальных физических задач, тем не менее реальные успехи могут быть охарактеризованы лишь как очень скромные. Причина этого, по мнению Н.Н.Миренкова [102] заключается в том, что специалисты в конкретной предметной области не понимают, как можно использовать ВВС, и не видят смысла в их использовании. Возможный выход заключается в создании удобных, интуитивно понятных инструментов для создания программ на ВВС. Такие работы активно ведутся, и в рамках настоящего обзора они выделены в тему «Упрощение работы с экзафлопсными системами». Сюда относятся создание учебников по конструированию и эксплуатации кластеров [76], разработка средств отладки параллельных приложений [64; 70], а также создание новых средств администрирования [15].

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

распределении различных видов данных, с различными вариантами доступа к ним, например, [78; 127].

В отличие от предыдущего пункта, к теме «Динамическая балансировка загрузки» отнесены вопросы оптимизации загрузки процессорных элементов, проводимой во время работы программы, например, исследование воздействия несбалансированности загрузки на производительность методов Монте-Карло [135] или реализация метода частиц в ячейках с динамической балансировкой загрузки на языке parallel C [66], а также вопросы, связанные с балансировкой загрузки не только прикладной программы, но и вычислительной системы в целом, что особенно важно с точки зрения перспективы экзафлопсных вычислений [4].

Следующая тема, «отказоустойчивость», исключительно важна для подготовки к вычислениям на экзафлопсных машинах, с учетом того, что различного сорта ошибки будут происходит на системах, состоящих из миллионов процессоров практически каждую секунду. Кроме того, стоимость проведения таких расчетов, как ожидается, будет очень высокой (сейчас расчет продолжительностью в 10 млн. процессоро-часов с использованием метода частиц в ячейках оценивается в 0.5 млн. евро [145]). К этой теме относятся обработка контрольных точек, оптимизация и ускорение их чтения и записи [90; 98], диспетчеризация процессов с учетом наличия поврежденных или вышедших из строя процессорных элементов [40], а также анализ производительности реализации MPI, работающей с учетом наличия в системе сбоев [65].

К последней теме («уменьшение количества коммуникаций») в рамках направления «масштабируемость и увеличение числа ПЭ» отнесены статьи, посвященные созданию параллельных алгоритмов, стремящихся к минимуму коммуникаций (англ. communication-avoiding [58]). Важность этой темы очевидна для экзафлопсных вычислений: при наличии миллионов процессорных элементов, обменивающихся сообщениями, количество сообщений оценивается (минимально), как 0(N), где N - число процессоров, поэтому очень важна разработка алгоритмов, способных или вовсе обойтись без коммуникаций или снизить их количество до 0(1). К таким алгоритмам относятся: метод решёточных уравнений Больцмана [32; 119], а также, разработка модели свободной поверхности океана с минимизацией объема коммуникаций[3] или решение уравнений типа уравнения Пуассона алгебраическим многосеточным методом,

масштабируемым в том смысле, что время решения пропорционально количеству неизвестных (а не количеству ПЭ) [92].

2.2 Адаптация вычислительных методов к архитектуре ВВС

Второе направление работ над экзафлопсными вычислениями, выделенное в рамках обзора, названо «адаптация вычислительных методов к архитектуре ВВС».

Сюда относятся численные методы, разработанные специально для эк-зафлопсных машин (методика локального сгущения шага при решении уравнений реакции-диффузии [73], решение дифференциальных уравнений в частных производных с пониженным объемом коммуникаций [91]). Другая тема - адаптация имеющихся численных методов для работы на очень большом количестве ПЭ (например, вариант итерационного метода Якоби с ускорением Андерсона [109] или масштабируемые алгоритмы Монте-Карло для финансового моделирования [10]). Здесь необходимо ответить на вопрос: где проходит грань между алгоритмами, специально разработанными для экзафлопсных машин, и алгоритмами, адаптированными для экзафлопса. Казалось бы, разница незначительна, и более того, поскольку то, что разрабатывается специально для экзафлопсных ВВС, также разрабатывается не на ровном месте, то можно сказать, что разницы нет вовсе, что это просто одно и то же. Тем не менее разница здесь фундаментальна, разница та же, что между свойством и определением. Алгоритмы, разработанные специально для экзафлопсных вычислений - это те алгоритмы, разработка которых изначально проводилась с целью минимизировать или свести к нулю коммуникации и обеспечить высокую энергоэффективность (с достижением физической корректности результатов). С другой стороны, алгоритмы к экзафлопсным вычислениям адаптированные - это те алгоритмы, которые разрабатывались для параллельных вычислений на десятках или сотнях ПЭ, и которые случайно оказались достаточно хорошо масштабируемыми для их эксплуатации на пета- и экзафлопсных системах. Фундаментальный характер имеющихся различий заключается еще и в том, что энергоэффективными эти алгоритмы почти никогда оказаться не могут по той причине, что разрабатывались они для процессоров общего назначения, и изменить это в процессе адаптации невозможно без принципиальных изменений самого алгоритма.

Также к направлению «адаптация вычислительных методов к архитектуре ВВС» отнесены статьи, посвященные работе с данными большого размера.

Это один из важнейших вопросов при подготовке к счету на экзафлопс-ных ВВС, вызванный очевидной необходимостью сохранить результаты счета объемом в сотни петабайт [23] и передавать их для обработки. Разделение направления по темам показано в таблице 5.

Под названием темы «со-дизайн и комплексная разработка программ» подразумеваются подходы к разработке параллельных вычислительных приложений с учетом особенностей вычислительной системы [16], многоуровневый параллелизм [69; 77], реализация параллельного алгоритма с учетом особенностей решаемой физической задачи [131], сочетанием различных подходов к параллельной реализации [31; 59] и др. В работе [45] представлен подход, названный «объединяй и властвуй», который заключается в следующем: параллельная программа, реализующая некоторую математическую модель, разбивается на три элемента, а именно вычисления, работу с данными и коммуникации. При этом все три элемента тесно взаимодействуют между собой, но разрабатываются различным образом и могут быть представлены каждый несколькими компонентами, что обеспечивает адаптацию к структуре ВВС и выбор оптимальных численных методов. В статье [49] излагаются различные возможные подходы аппаратной точки зрения к построению экзафлопсной ВВС на примере нескольких действующих машин-прототипов с аналогичной архитектурой, но меньшей мощности. Моделирование расчетов на экзафлопсных системах проводится на машинах-прототипах с помощью так называемых мини-приложений. В работе [111] разработана методология упрощения создания и улучшения переносимости приложений для экзафлопсных ВВС, которая называется «разработка сверху вниз» (top down methodology). Данная методология основывается на следующем: разработка кода проводится так, чтобы он не был привязан к типу параллелизма или используемому оборудованию. Доработка под реально имеющееся оборудование проводится на следующей стадии разработки и таким образом код может быть перенесен и запущен на любом типе ВВС.

Таблица 5: Распределение по темам статей относящихся к направлению «адаптация вычислительных методов к архитектуре ВВС»

Тема Количество статей в обзоре В% от общего числа

адаптация вычислительных методов к экзафлопсу 19 6.38

Вычислительные методы, специально разработанные для экзафлопса 17 5.70

Обработка данных большого размера 30 10.07

Со-дизайн и комплексная разработка программы 13 4.36

2.3 Использование ускорителей вычислений

Таблица 6: Распределение по темам статей, относящихся к направлению «использование ускорителей вычислений»

Тема Количество статей в обзоре В% от общего числа

ускорение расчетов с помощью ускорителей 16 5.37

упрощение использования ускорителей 8 2.68

новые типы ускорителей 8 2.68

энергоэффективность 31 10.40

Третье направление работ над экзафлопсными вычислениями, названное условно «использование ускорителей вычислений», было также разделено на несколько тем, как показано в таблице 6.

Наиболее очевидное применение ускорителей - ускорение расчетов (первая строка в таблице 6), например, вопрос о том, как получить на методе частиц в ячейках заявленную для ускорителя Intel Xeon Phi производительность в 1 Teraflops [89], или расчеты в области микроволновой томографии для обнаружения рака груди с использованием процессоров IBM Cell [148]. Также рассматриваются затруднения, возникающие при работе с памятью при реализации методов Монте-Карло на многоядерных системах [141], и проблемы эффективной реализации метода частиц в ячейках на GPU [104; 105; 118].

Несмотря на то, что использование ускорителей вычислений, как графических, так и Intel Xeon Phi и ПЛИС, обещает очень серьезные преимущества как по производительности, так и по энергоэффективности, реализация вычислительных алгоритмов на ускорителях остается очень сложным вопросом по сравнению с реализацией их на процессорах общего назначения. В связи с этим актуальной является тема «упрощение использования ускорителей». В частности, предложен набор программных инструментов для упрощения разработки

[63]. Кроме того, рассматриваются вопросы переноса программ между различными ускорителями вычислений [111] и создание настраиваемых алгоритмов для GPU на примере метода частиц [38].

Следующая тема названа «новые типы ускорителей». При формулировке такой темы этом необходимо ответить на вопрос, какие ускорители вычислений, и почему считаются новыми. В настоящее время для построения ВВС используются в основном два типа ускорителей: ускорители Intel Xeon Phi и графические ускорители (GPU) производства компании Nvidia. Несмотря на то, что и ПЛИС, и графические ускорители корпорации AMD (например, AMD Firestream), и многоядерный процессор IBM Cell существуют достаточно давно, они используются значительно реже. Это дает определенные основания считать все эти ускорители, кроме Intel Xeon Phi, а также Nvidia Kepler и Nvidia Tesla новыми (или нетрадиционными) типами ускорителей вычислений. К этой теме отнесено сравнение высокоуровневых инструментов реализации численных алгоритмов на ПЛИС [146] и опыт создания и эксплуатации кластера на процессорах ARM с пониженным энергопотреблением [139]. Следует особенно отметить, что последняя работа посвящена именно вычислениям на таком кластере и проблемам достижения вычислительной производительности.

Естественным образом отсюда можно перейти к следующей теме: «энергоэффективность». Эта тема особенно важна для экзафлопсных вычислений, в силу того что если ВВС будут строится также, как сейчас, то экзафлопсная ВВС будет иметь слишком высокое энергопотребление. Это можно проиллюстрировать следующим образом: производительность 2-го по мощности компьютера в мире (Tianhe-2) в 30 раз меньше, чем 1 экзафлопс, однако если его энергопотребление (17.8 МВт) умножить на 30, получится 534 МВт, что превышает мощность большинства ядерных реакторов. Аналогично для суперкомпьютера № 1, Sunway TaihuLight, эти показатели составляют 125 Петафлопс и 15 МВт соответственно. Разница с экзафлопсом в 8 раз, оценочное энергопотребление такой системы с экзафлопсной производительностью 120 МВт, что значительно меньше, но тем не менее, очень много, при том, что базовые элементы суперкомпьютера Sunway TaihuLight пока не используются массово при производстве ВВС - поэтому ориентироваться на них, т.е. на процессоры Sunway SW26010 260C при прогнозировании перспективных архитектур ВВС нельзя.

Так или иначе это означает необходимость радикального понижения энергопотребления ВВС. К теме «энергоэффективность» в рамках настоя-

щего обзора отнесены следующие работы: исследование производительности методов Монте-Карло с учетом энергопотребления [17], методы измерения, моделирования и управления энергопотреблением [53] и стратегия создания энергоэффективных программ [142].

Таким образом, по материалам обзора можно сказать, что наиболее актуальными вопросами с точки зрения экзафлопсных вычислений являются энергоэффективность (10.4 % статей в обзоре, таблица 6), обработка данных большого размера (10.07 % статей, таблица 5) и разработка масштабируемых архитектур ВВС (9.73 % статей, таблица 4).

С другой стороны недостаточное внимание по результатам обзора в научной литературе, по крайней мере в рассмотренных журналах, уделяется созданию вычислительных методов специально предназначенных для экзафлопсных ВВС (5.7 % в таблице 5), и даже вопросам адаптации вычислительных методов к экзафлопсным вычислительным системам (там же, 6.38 %). Более того, в основном рассматриваются варианты методов Монте-Карло и методы решения больших СЛАУ, преимущественно на основе методов подпространства Крылова. Остальные методы и технологии математического моделирования в применении к экзафлопсным ВВС изучаются значительно ре-

Также сравнительно небольшое внимание (4.36 %, таблица 5) уделяется такому важному вопросу как со-дизайн и комплексная разработка программ. При этом важно отметить, что лишь две работы из представленного списка ([45; 111]) используют подход со-дизайна полностью, остальные ограничиваются лишь отдельными его элементами, в частности, не производится переработка самого численного метода под архитектуру ВВС (хотя возможность такая рассматривается).

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

Глава 3. Комплексная оценка производительности ВС

В настоящий момент списки Тор50 и Тор500 выстроены в порядке убывания пиковой производительности и производительности на тесте ЫпРаек что, разумеется, дает определенную информацию о сравнительной скорости работы представленных там машин. Но очень многие факторы, такие как скорость работы и объем дисков, пропускная способность шины памяти и коммуникационной сети, неоднородность оборудования и т.д. - остаются за пределами рассмотрения. А это именно те проблемы, с которыми придется столкнуться при попытке посчитать на кластере большую задачу. По этой причине тестирование проводится с помощью измерения времени, затрачиваемого на различные этапы программы, решающей реальную физическую задачу.

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

Кроме того, на основе проведенных расчетов измерена скорость счета и скорость перемещения данных для нескольких протестированных ВС. Результаты, показанные в этой главе, опубликованы в [30; 136; 168; 177]

3.1 О влиянии организации данных на результат измерения производительности процессоров

Модельные частицы расположены внутри расчетной области случайным образом. Если даже модельные частицы расположены рядом в массиве, где хранятся их координаты, то сами значения координат будут близкими только вначале. Частицы с близкими значениями координат находятся в соседних ячейках, или даже просто в одной и той же ячейке. В дальнейшем модельные частицы движутся в рамках расчетной области непредсказуемым образом,

и таким образом, перестают находится в соседних ячейках. Это означает, что обращения к трехмерным массивам, содержащим электрическое и магнитное поля, являются неупорядоченными. Если для частицы с номером ] требуется обратиться к ячейке трехмерной сетки (трехмерного массива) с индексами 1,1,к, то для частицы с номером ]+1, с большой вероятностью, потребуется обращение не к ячейке с индексами (1,1,к-1) или (1,1,к+1) и т.д., то есть к соседней ячейке, а наоборот, потребуется обращение к ячейке, значительно удаленной в рамках расчетной области. Из этого можно сделать вывод, что использование кэш-памяти в данном случае не позволяет сократить время счета.

Здесь важно отметить, что предсказать заранее движение модельных частиц нельзя: знание траекторий движения модельных частиц составляет сущность решения задачи о динамике плазмы, если такое знание есть, то по сути дела задача решена, и нет необходимости проводить расчет.

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

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

Модельные частицы, принадлежащие некоторой ячейке сетки, можно хранить в виде связанного списка или в виде массива. Преимущества списка очевидны: нет ограничения на число частиц в ячейке, простота добавления/удаления, но есть и недостатки, а именно большее по сравнению с массивом время доступа. Если же частицы каждой ячейки хранятся в виде массива (статического массива), то применительно к трехмерной задаче для двух сортов частиц это даст 5-мерный массив для одной только координаты X всех частиц (напомним, что модельная частица в данном случае характеризуется шестью признаками).

Но основная проблема в случае статического массива - это максимальное число частиц в ячейке. Это означает, что заранее неизвестно, какой длины массив потребуется для хранения всех частиц в каждой ячейке. Из проведенных расчетов известно, что максимальное значение плотности электронов было равно 5 (в единицах начальной плотности). Таким образом, можно было бы задать

длину массива частиц в каждой ячейке 5Ж, где N - число частиц в ячейке в начальный момент времени. Но в этом случае размер массива частиц увеличится в 5 раз, а его размер составляет 70 Гб, например, для сетки 512x64x64 и N = 150.

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

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

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

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

Таким образом, имеются следующие варианты организации частиц в программе:

1. исходный неоптимизированный вариант;

2. хранение значений поля в 4-мерном массиве;

3. упорядочивание частиц с помощью массивов номеров;

4. упорядочивание частиц с помощью связанного списка.

Далее рассмотрим результаты тестов, показывающих эффективность выполненной оптимизации. Тестовые расчеты проводились на рабочей станции с процессором AMD Phenom (Phenom II X6 1055T), производительность (раздел А.2 в приложении А) 41.9 GFLOPS, и на кластере Новосибирского Государственного Университета, оснащенного процессорами процессор Intel Xeon E5540 производительность (раздел А.1.1 в приложении А) на LU-разложении 2.5 GFLOPS). В обоих случаях была выбрана сетка такого размера, что даже один трехмерный массив, содержащий, например, одну из компонент поля, заведомо не помещается в кэш. Размер сетки для рабочей станции на базе процессора Phenom 64 х 32 х 32 узла, 150 частиц в ячейке, всего 9.8 млн. частиц. Размер сетки для кластера на базе процессора Intel Xeon E5540 в расчете на одно ядро, 512 х 2 х 64 узла, 150 частиц в ячейке, те же 9.8 млн. частиц на ядро. С учетом определенного в главе 1 количества операций в расчете на одну частицу (440, см. подсчет в разделе 1.6.2) можно рассчитать производительность в единицах GFLOPS.

Таблица 7: Время вычислений, в секундах, и производительность, GFLOPS.

№ Вариант оптимизации Время Производительность

Phenom Xeon E5540 Phenom Xeon E5540

1 Исходный неоптимизированный вариант 13.25 7.22 0.369 0.67

2 Хранение значений поля в 4-мерном массиве 8.8 6.72 0.55 0.72

3 Упорядочивание частиц с помощью массивов номеров 12.51 5.67 0.39 0.86

4 Упорядочивание частиц с помощью связанного списка 10.5 10.3 0.46 0.475

5 Сочетание вариантов 2 и 3 10.92 3.67 0.44 1.33

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

частиц и хранения значений поля в одном 4-мерном массиве приводит к значительному сокращению времени вычислений с частицами, в данном случае, почти в 2 раза. Хранение частиц в виде связанных списков оказалось неэффективным для процессора Xeon E5540, но, возможно, окажется полезным на других архитектурах (так позволяют думать измерения времени на процессоре Phenom) или в тех случаях, когда установить максимальное число частиц в ячейке невозможно.

Относительно небольшая доля от пиковой производительности, полученная в этих тестах, показывает значительную разницу между расчетом по методу частиц в ячейках, и например, методами решения линейных уравнений, которые являются наиболее распространенным тестом производительности ВС. Это подчеркивает необходиость отдельного тестирования также и на методе частиц в ячейках, и кроме того, подчеркивает важность также и проверки в рамках тестов не только производительности, но и скорости доступа к оперативной памяти, что будет предметом рассмотрения в одном из следующих разделов (раздел 3.2).

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

3.2 Расчет пропускной способности системы памяти

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

Таблица 8: Основные характеристики кластеров, на которых производились рас-четы,рассмотренные в разделах 3.2,3.4 и 3.5, Номера даны по списку Тор50 от 2009 года, В колонке ТРЮРБ указана пиковая производительность в терафлоп-

№ Название Узлы Сеть TFLOPS

1 МВС-100К 4xXeon E5450, 3 GHz 8.192 GB RAM Infiniband 4x DDR/ 2xGigabit Ethernet/ Gigabit Ethernet 95.04

2 СКИФ-МГУ 2xXeon E5472, 3 GHz 8.192 GB RAM InfiniBand/ Gigabit Ethernet/ СКИФ-ServNet + IPMI 60

14 СКИФ-Cyberia 2xXeon 5150, 2.667 GHz, 4.096 GB RAM QLogic InfiniPath/ Gigabit Ethernet/ СКИФ-ServNet 12.002

20 Кластер НГУ 2xXeon 5355, 2.66 GHz Infiniband 4x DDR/ Gigabit Ethernet/ Gigabit Ethernet 5.4

Переход от фактически измеренной величины (времени выполнения расчета движения модельных частиц) к пропускной способности памяти выполнялся из следующих соображений: на каждое из используемых ядер приходится 2.5 млн. модельных частиц, каждая частица занимает 48 байт, кроме того, для расчета движения частицы необходимы значения электрического и магнитного полей в той ячейке сетки, где находится частица. Это означает, что для каждой из 6 компонент электромагнитного поля загружается 8 значений, соответствующих вершинам параллелепипеда, то есть ячейки сетки.

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

Таким образом для каждой модельной частицы загружается из памяти 432 байта и сохраняется 144 байта, общий поток данных составляет 576 байт на одну частицу.

В итоге, для вычисления пропускной способности памяти (в СБ/эее) при расчете движения модельных частиц WpJc,GB/sec использовалась следующая формула:

^ПС&В/зес —

Wp х Ыр х Ра

~ДГ

(3.1)

здесь:

Wp - количество байт на одну модельную частицу, Wp — 576; Ыр - количество модельных частиц на одно процессорное ядро (в рассмотренном случае 2.5 х 106); Рсоге - количество процессорных ядер; ДЪ - длительность временного шага, сек.

Рисунок 3.1: Пропускная способность памяти на этапе расчета движения модельных частиц на некоторых кластерах. Количество модельных частиц: 2.5 млн. на каждое процессорное ядро. Основные параметры ВС, использованных в эксперименте, показаны в таблице 8.

Следует отметить, что вопрос о сравнении чисел на рис. 3.1 с заявленной максимальной пропускной способностью является второстепенным, тем не менее, сравнение показано в таблице 9. Основной вопрос в данном случае -это измерение пропускной способности памяти, фактически доступной для расчетного приложения.

Вывод, который можно сделать по таблице 9: для реального расчетного приложения доступно не более половины заявленной максимальной пропускной способности оперативной памяти, причем даже этот показатель получен только лишь на ВВС СКИФ-МГУ «Чебышев», в отношении которой проводилась специальная оптимизация взаимодействия процессоров с оперативной памятью.

Таблица 9: Сравнение пропускной способности памяти, измеренной с помощью теста на основе метода частиц с характеристиками процессора.

Название ВС Процессор Пропускная способность памяти, GB/sec

Данные теста PIC-MANAS Максимум

СКИФ-Cyberia Xeon 5150 1.53 10.6

МВС-100К Xeon E5450 6.56 21

СКИФ-МГУ Xeon E5472 13.62 21

Кластер НГУ Xeon 5355 4.81 21

3.3 Расчет производительности процессорных элементов

В этом разделе описана методика измерения производительности процессорных элементов. Для того, чтобы отделить время счета от времени обращения к оперативной памяти было рассмотрено время работы процедуры, «««< HEAD реализующей одномерное преобразование Фурье, которая является частью физической диагностики, используемой при моделировании динамики плазмы. Измеренное время с учетом известного размера данных и количества операций в БПФ [169], переводится во флопсы: размер массива для одномерного фурье-преобразования 64, количество операций, в соответствии с приведенным в статье [169] 1190. Сравнительная производительность процессорных элементов некоторых из рассмотренных в диссертационной работе ВС выглядит как показано на рисунке 3.2: реализующей одномерное преобразование Фурье, которая является частью физической диагностики, используемой при моделировании динамики плазмы. Измеренное время с учетом известного размера данных и количества операций в БПФ, переводится во флопсы (1 флопс - 1 операция с плавающей точкой в секунду). Сравнительная производительность процессорных элемен-

тов некоторых из рассмотренных в диссертационной работе ВС показана на рисунке 3.2:

Производительность процессорных ядер

Преобразование Фурье (Ю)

ЕЗСКИФ-СуЬепа

□ МВС-ЮОК ■ СКИФ-МГУ

□ Кластер H ГУ

Рисунок 3.2: Производительность процессоров Intel Xeon, измеренная в ходе выполнения одномерного преобразования Фурье на некоторых кластерах. Размерность преобразования N = 64.

3.4 Расчет производительности процессорных элементов на основе

движения модельных частиц

Расчет количества операций с плавающей точкой в секунду (FLOPS) при расчете движения модельных частиц Npjc,flops производился следующим образом:

ы Fp Х Np Х Рсоге (ъ о\ WPIC,FL0PS = -Д--(3-2)

здесь: