Разработка параллельных алгоритмов в прикладной гидрофизической задаче с использованием современных средств профилирования кода

Тюгин Дмитрий Юрьевич кандидат физико-математических наук научный сотрудник, лаборатория моделирования природных и техногенных катастроф, Нижегородский Государственный Технический Университет им. Р. Е. Алексеева 603155, Россия, Нижегородская область, г. Нижний Новгород, ул. Минина, 24, ауд. 1117 Tyugin Dmitry PhD in Physics and Mathematics Researcher at Laboratory for Modelling Natural and Technogenic Catastrophes at the Alekseeva Nizhny Novgorod State Technical University  603155, Russia, Nizhegorodskaya oblast', g. Nizhnii Novgorod, ul. Minina, 24, aud. 1117 dtyugin@gmail.com Другие публикации этого автора

Представленные результаты получены за счет средств гранта Российского научного фонда (проект №17-71-10101)

Введение

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

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

Среди профилировщиков кода для многоядерных систем с общей памятью можно выделить следующие программные продукты: Valgrind ^[1], OpenSpeedShop ^[2], ompP ^[3], TAU ^[4], VTune ^[5], MAQAO ^[6] и др. Некоторые из них требуют предварительной подготовки кода для анализа приложений, путем вставки специальных директив или вызовов функций внешней библиотеки. Кроме того, не все из них имеют пользовательский интерфейс, а выполняют только профилировку, визуализация и анализ результатов может выполняться в другом приложении, поддерживающем формат вывода. После анализа доступных профилировщиков выбор был сделан в пользу программного продукта VTune. Он содержит набор приложений для профилировки, имеет пользовательский интерфейс и не требует вставки директив в исходный код для анализа времени выполнения участков кода.

После анализа кода распараллеливание может быть выполнено с применением различных библиотек и технологий: OpenMP ^[7,8], TBB ^[9], QtConcurrent ^[10]. Стандарт OpenMP базируется на включении директив препроцессора в исходный код, и поддерживается большинством компиляторов. Следует также отметить наличие выстроенных средств параллелизма в новых стандартах языков программирования. Например, в C++ версии 11 появилась поддержка потоков в библиотеке STL (std::threads), которую также поддерживает большинство современных компиляторов.

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

Для решения задачи распараллеливания была выбрана библиотека TBB. Выбор связан с тем, что в библиотеке реализована эффективная поддержка планировщика и распараллеливание на основе задач.

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

Рассмотрим задачу численного моделирования внутренних гравитационных волн в океане в рамках слабонелинейной теории. Данный тип волн возникает в океане вследствие вертикальной неоднородности воды по температуре и солёности. Это приводит к образованию в толще жидкости слоев различной плотности, на границах которых и образуются данные волны. Рассмотрим классическое уравнение для описания распространения внутренних волн – уравнение Гарднера ^[11]

Где n является смещением изопикнической поверхности в максимуме вертикальной моды; c – фазовая скорость внутренних волн; a – коэффициент квадратичной нелинейности; a1 – коэффициент кубической нелинейности; b – коэффициент дисперсии. Параметры c, a, a1 и b зависят от стратификации плотности и глубины океана.

В частности, такие данные могут быть получены путем натурных измерений температуры и солёности в океане с помощью контактных датчиков подводного зондирования, опускаемых с корабля, либо могут быть рассчитаны на основе математических моделей. Эти данные распространяются на свободной основе организациями, которые занимаются исследованием климата планеты. Для расчета будем использовать набор данных RCO полученный на базе модели глобальной циркуляции ^[12].

На основе карт двумерных распределений параметров уравнения Гарднера можно проводить моделирование распространения волн. Таким образом расчет карт данных параметров является необходимым этапом численного эксперимента. Подробное изложение математической теории расчета данных параметров представлено в работах ^{[11, 13]}. В статье ограничимся рассмотрением программной реализации алгоритмов расчета.

Основной алгоритм расчета карт параметров:

1. Для каждой географической точки карты

1.1.Загрузить данные солености, температуры, глубины

1.2.Рассчитать параметры плотности

1.3.Рассчитать коэффициенты

1.3.1. Методом пристрелки найти решение краевой задачи Штурма-Лиувилля

1.4.Сохранить коэффициенты в контейнер данных

Анализ кода и разработка параллельного алгоритма

Рассмотрим пример анализа кода с помощью профилировщика VTune. Для этого необходимо собрать исследуемое приложение в конфигурации Release с включением символов отладки (параметры компилятора /Zi /DEBUG для компилятора MSVC). Затем необходимо выполнить HotSpot анализ путем запуска приложения через профилировщик. Данный тип анализа собирает информацию о времени выполнения кода предоставляя распределение в процентах по функциям и строкам кода, рис. 1.

Рис. 1 Анализ последовательного алгоритма профилировщиком VTune

Интересуемый участок расчётного кода представлен в листинге 1.

Листинг 1.

void calcCoeffs(igwdatamodel::IDataProfileIterator::Ptr tempProfileIt,

igwdatamodel::IDataProfileIterator::Ptr salProfileIt,

igwdatamodel::IDataWriter::Ptr dst)

{

while(!tempProfileIt->atEnd() && !salProfileIt->atEnd())

{

igwdatamodel::IDataProfile::Ptr tempProfile = tempProfileIt->get();

igwdatamodel::IDataProfile::Ptr salProfile = salProfileIt->get();

Coeffs result = coeffs(getZ(tempProfile),

getBvf(tempProfile, salProfile),

getDepth(tempProfileIt));

dst->write(result);

tempProfileIt->next();

salProfileIt->next();

}

}

Как видно из анализа большая часть времени выполнения (99.8%) принадлежит функции deriveddatalib::coeffs. Анализ стека этой функции показывает, что основное время состоит из вызовов методов интегрирования системы дифференциальных уравнений, реализованной в библиотеке boost - integrate_adaptive (61.6%) и (22.9%). Данные методы не могут быть распараллелены так как шаги интегрирования имеют зависимость между предыдущим и следующим. Поэтому распараллеливание целесообразно провести по данным на основе независимости расчетов коэффициентов в каждой отдельно взятой географической точке. Таким образом логично поместить множественные вызовы функции deriveddatalib::coeffs в отдельные задачи.

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

DataHolder, листинг 2.

Листинг 2.

class DataHolder{

igwdatamodel::IDataWriter::Ptr dst;

igwdatamodel::IDataProfileIterator::Ptr tempProfileIt,salProfileIt;

tbb::mutex readMutex;

tbb::mutex writeMutex;

public:

Profile getNextProfile(){

tbb::mutex::scoped_lock lock(readMutex);

Profile profile;

while(!tempProfileIt->atEnd() && !salProfileIt->atEnd())

{

profile = getProfile(tempProfileIt->get(),salProfileIt->get())

tempProfileIt->next();

salProfileIt->next();

return profile;

}

return Profile();

}

void setResult(const Coeffs& result){

tbb::mutex::scoped_lock lock(writeMutex);

dst.write(result);

}

};

Данный класс содержит в себе указатели на контейнеры входных данных температуры и солености (tempProfileIt, salProfileIt), и контейнер для записи данных (dst). Так как контейнеры не имеют атомарных средств доступа то необходимо синхронизировать вызов этих функций так чтобы две задачи не могли одновременно выполнить их вызов. В функциях чтения (getNextProfile) и записи (setResult) реализован вызов синхронизационных примитивов типа мьютекс ( tbb::mutex::scoped_lock lock(writeMutex)). Библиотека TBB содержит специальный класс обертку для реализации блокировки мьютекса (tbb::mutex::scoped_lock), позволяющий заблокировать мьютекс при создании объекта обертки и освободить при разрушении.

Выполнение расчета характеристик в точке перенесено в класс CalcTask, наследуемый от tbb::task. Данный класс реализует разделение работ по задачам, листинг 3.

Листинг 3.

class CalcTask: public tbb::task{

DataHolder* m_holder;

public:

CalcTask(DataHolder* dataHolder):

m_holder(dataHolder){}

tbb::task* execute() {

Profile data = m_holder-> getNextProfile();

if(data.valid){

Coeffs result = coeffs(data.z, data.bvf, data.depth);

m_holder->setResult(result);

tbb::task::recycle_as_continuation();

return this;

}

return nullptr;

}

};

Класс CalcTask принимает на вход указатель на объект доступа к данным и переопределяет виртуальную функцию родительского класса execute, тело которой определяет код, выполняемый параллельно.

Вызов функции tbb::task::recycle_as_continuation() помечает задачу как задачу для повторного выполнения. Таким образом, задача будет перезапущена планировщиком повторно. Данный подход позволяет использовать пул задач и переиспользовать без удаления старых и создания новых задач, тем самым минимизируя накладные расходы, связанные с созданием и удалением объектов. Реализация параллельного запуска задач представлена в листинге 4.

Листинг 4.

void calcCoeffsParallel(igwdatamodel::IDataProfileIterator::Ptr temperature,

igwdatamodel::IDataProfileIterator::Ptr salinity,

igwdatamodel::IDataWriter::Ptr dst)

{

tbb::task_scheduler_init init;

DataHolder holder(temperature,

salinity,

dst);

tbb::task_list list;

for(int i=0;i<init.default_num_threads();++i){

CalcTask& t = *new(tbb::task::allocate_root()) CalcTask(&holder);

list.push_back(t);

}

tbb::task::spawn_root_and_wait(list);

}

Запуск параллельных задач осуществляется в новой функции calcCoeffsParallel. В ней создается объект DataHolder, инициализируется данными. Далее создается объект tbb::task_scheduler_init, который позволяет получить доступное на системе число потоков, как правило равное числу физических потоков поддерживаемое аппаратным обеспечением. Список заполняется количеством задач равным числу потоков. Число задач может быть и больше количества потоков, так как планировщик будет выполнять балансировку нагрузки автоматически. Но в данном случае задачи при старте и окончании расчета выполняют доступ к данным и время доступа мало. Поэтому увеличение задач не приведет к более эффективному распараллеливанию. В том случае если доступ к данным может занимать значительное время, то целесообразно создать отдельную параллельную задачу для наполнения данных. В таком сценарии создание количества задач большее, чем число физических потоков целесообразно, так как заполненные данными для расчета задачи будут ожидать в очереди выполнения и начинать расчет сразу в то время пока данные из выполненных задач будут сохраняться.

Анализ достигнутой производительности

Расчеты проводились на сервере со следующей конфигурацией: 2 процессора Intel Xeon CPU E5-2650 v4 2.20Ггц, 64 Гб оперативной памяти, 512 Гб SSD диск Samsung Evo 850. Каждый процессор имеет 14 ядер и может выполнять до 28 потоков одновременно с помощью технологии Hyper Threading. Суммарно сервер может выполнять до 56 потоков одновременно.

Результаты профилировки распараллеленного приложения в профилировщике VTune представлено на рис. 2.

Рис.2 Результат профилировки параллельного алгоритма в профилировщике VTune

Анализ стека вызовов и времени выполнение показывает, что основное время работы приложения находится в расчетных задачах.

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

Рис.3 График ускорения параллельного алгоритма

На рис. 3 приведен анализ ускорения параллельного алгоритма при расчете на 2348 точках исходного набора данных. Пунктирной линией помечено идеальное ускорение, которое равно количеству запускаемых расчетных задач. В реальности (синяя кривая) практическое ускорение почти никогда не совпадает с идеальным, так как расчеты зачастую связаны с сохранением и загрузкой данных, а также синхронизацией, которые вносят накладные расходы процессорного времени. С увеличением количества потоков возрастает интенсивность вызова синхронизационных примитивов. Кроме того, дополнительные потоки в режиме Hyper Threading имеют более низкую производительность чем потоки запущенные в режиме равном количеству физических ядер процессора.

На рис. 4 представлена карта фазовой скорости, рассчитанная параллельным алгоритмом. Время расчета сократилось в 39 раз и составило около 165 секунд против 6437 секунд последовательного алгоритма. Данный расчетный блок входит в состав программного комплекса для исследования волн в океане ^[14] и позволит проводить вычисления значительно быстрее для новых наборов данных, охватывающих большие территории.

Рис. 4 Карта распределения параметра фазовой скорости

Заключение

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