Распараллеливание противопотоковой разностной lu-схемы на кластерных системах с использованием MPI Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 9, № 6, 2004

РАСПАРАЛЛЕЛИВАНИЕ ПРОТИВОПОТОКОВОЙ РАЗНОСТНОЙ LU-СХЕМЫ НА КЛАСТЕРНЫХ СИСТЕМАХ С ИСПОЛЬЗОВАНИЕМ MPI

А. Д. Рычков, Н.Ю. ШокинА Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: rych@ict.nsc.ru, shokina@hlrs.de

У. КЮСТЕР

High Performance Computing Center, Stuttgart, Germany e-mail: kuester@hlrs.de

The parallelization of the upwind difference LU-scheme is considered by the example of the problem on three-dimensional non-stationary flow in the airbag combustion chamber. The data on the calculation speed acceleration are provided, obtained using two cluster systems in High Performance Computing Center Stuttgart (Germany).

Введение

В настоящее время практически все ведущие автомобилестроительные фирмы и концерны оснащают свои автомобили среднего класса и выше подушками безопасности (айрбэ-гами), позволяющими значительно снизить число летальных исходов при автомобильных авариях. В качестве высокопроизводительных источников газа для наддува мягкой оболочки айрбэгов широко применяются пиротехнические газогенераторы, унитарные твердотопливные составы которых обеспечивают на выходе экологически безопасные газообразные продукты сгорания с относительно низкой температурой. Количество различных конструкций таких газогенераторов может достигать десятков единиц — от миниатюрных, для натяжения ремней безопасности, до газогенераторов наддува айрбэгов переменной газопроизводительности, которая может изменяться, например, в зависимости от веса пассажира или водителя. О важности и перспективности использования этих средств говорит и тот факт, что уже в течение десяти последних лет регулярно проводится международная конференция АЙРБЭГ-2000 (Fraunhofer Institut Chemische Technologie, Pfinztal (Berghausen)), специально посвященная вопросам разработки, эксплуатации и развитию этих средств безопасности.

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

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

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

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

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

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

В данной работе (см. также [1]) рассмотрена технология распараллеливания проти-вопотоковой разностной ЬИ-схемы на примере задачи о трехмерном нестационарном течении в камере сгорания айрбэга. Приведены данные об ускорении счета. Расчеты были выполнены на двух кластерных системах в Центре высокопроизводительных вычислений (НЬКБ) в г. Штутгарт (Германия).

1. Физическая модель

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

Топливные гранулы, заполняющие остальной объем камеры сгорания, имеют цилин-

щель

Рис. 1. Схема камеры сгорания айрбэга.

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

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

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

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

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

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

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

2. Математическая модель

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

в векторной форме записывается как

dQ 5R dF dG

dt dx dy dz

H,

(1)

R

где Q = (p, pu, pv, pw, pE, pk, pe)T, p = "pT, Ro, " — универсальная газовая постоянная и молекулярный вес газа соответственно;

R = (pu, pu2 +p-rxx, pUV-TXy, pUW Txz ,u(pE +p)-urxx-vrxy WTxz-qxx, puk-Txk, pu£-Txe)T;

F = (pv, puv -Tyx,pv2 +p-Tyy, pVW-Tyz ,v(pE +p)-uTyx -VTvy WTyz -Qyy, pvk -Tyk, pve-Ty£)T; G = (pw, puW-Tzx, pvW-Tzy, pW2+p-Tzz ,w(pE+p)-uTzx-vTzy-WTzz-Qzz, pwk-Tzk, pWE-Tzef;

H =(0, 0, 0, 0, (Ф - pe), (C\ЛФ - C12f2pe)e/k);

(du dv

xy

Tyx Ve ( r> +

Tx

Tz

Ve

\dy dx

du dw dz dx

Ty

yz

Tz

zy

dv dw\ V dz dy) ,

Ty

yy

2 du 2 / du + dv + dw

e dx 3 \ dx dy dz

2 dv 2 / du + dv + dw

e dy 3 \ dx dy dz

2 dw 2 / du + dv + dw

e dz 3 \ dx dy dz

Qx

dT

dT

dT

Xe dx, Qyy = К dy, Qzz = К dz,

Vt\ dk

Txk = {Vt + VJ a£, Tyk

Tx

Vt\ de

Vt + V£) a£, Ty£

Vt\ dk

Vt +--T-, Tzk

(?k J dy

Vt\ de

Vt +--T-, Tz£ -

O£ dy

Vt

Vt +— Ok

Vt

Vt +— O £

dk

dz'

de

dzJ

. du dv dw

=Vt 1 Txx dx + Tyy dy + Tzz +

2 k f du + dv + dw

3 \dx dy dz

du + dv dy dx

+

du + dw dz dx

2

+

dv dz

+

dw dy

C1 = 1.44, C2 = 1.92, fi = 1, f2 = 1 - 0.22exp

-Rt 36

Ok

1,

O£

1.3,

Я = pk2 и

Rt = -, Ve

Ve

Cp = 0.09,

л л Cp ™ p pk

= V + Vt, К = A + VtTT~, Vt = Cfp-,

Prt e

U = [1 - exp (aiRk + «2Ri + a3Rl)]1/2 ,

2

2

= —1.5 ■ 10-4, «2 = —1.0 ■ 10-9, «з = —5 ■ 10-10, Кк = Рк1/2Уп.

ц

Здесь уп — расстояние до ближайшей стенки; ц, Л — молекулярная вязкость и теплопроводность; к, е — турбулентная кинетическая энергия и скорость ее диссипации; Е = и2 + V2 + и)2

е +--^--полная внутренняя энергия на единицу объема; е — внутренняя энергия, причем e = СЕ для совершенного газа постоянного состава; С и Ср — удельные теплоемкости при постоянном объеме и постоянном давлении соответственно.

3. Численное решение

Для численного решения системы уравнений (1) применяется метод конечных объемов с использованием противопотоковой разностной ЬИ-схемы второго или третьего порядков точности [3]. С этой целью в левую часть системы (1) добавляется производная по псев-дЦ

довремени ——, по которой организуется итерационный процесс на каждом временном дт

слое, затем система линеаризуется по Ньютону и записывается в так называемой "дельта"-форме. Это означает, что линеаризация проводится относительно значений переменных на п-м временном слое, причем не учитывается зависимость матриц Якоби от этих переменных, и в левой части полученной системы сохраняются только члены, отвечающие за первый порядок аппроксимации для невязких членов, и члены, аппроксимирующие повторные производные по соответствующим направлениям для вязких членов. Линеаризованная таким способом система записывается в следующем виде:

ТдТ + 2 ТД^ + А-+1/2Л*+1/2 + Л+-1/2Д*-1/2 + А-+1/2А^+1/2 + А+_1/2'1/2 +

+А_+1/2Дк+1/2 + А+_1/2Дк-1/2+] = (2)

2 — + ЦпП_,1 Сп+1Л

2Д-Чзк + [ШЬ^к) ,

где

Дфз+1 = (ц^1)5+1 — (цп+1

8

)

ДЯ^.к = — [ (К ■ Ь)г+1 /2 — (К ' Ь)г_ 1/2 + (Г ' Ь)^'+1/2 — (Г ' Ь)^_1/2 +

+ (С ' Ь)к+1/2 — (С ' Ь)к_1/2

+ Тг.7.к Н.

Ьг±1/2, Ь^±1/2, Ьк±1/2 — площади граней объема ^¿к; (К ' Ь),±1/2, (Г ■ Ь)^±1/2, (С ■ Ь)к±1/2 — соответствующие суммарные (вязкие и невязкие) разностные потоки через эти грани; в и в + 1 — временные шаги по псевдовремени.

После того как итерации по псевдовремени сойдутся, левая часть системы (2) обратится в нуль, поскольку Д^8+1 = 0. Тогда система разностных уравнений (2) будет аппроксимировать систему (1) на (п + 1) -м временном слое со вторым порядком точности по времени и со вторым или третьим в зависимости от способа аппроксимации разностных потоков через грани по пространству.

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

скорость сходимости при этом несколько уменьшается, но зато простота структур матриц позволяет свести решение (2) к последовательным скалярным рекуррентным вычислениям типа бегущего счета, что существенно повышает эффективность решения (2), особенно на ЭВМ параллельного действия. Опуская все несущественные технические детали таких преобразований, запишем ЬИ-алгоритм для решения (2), который выполняется в два этапа:

Ь-обход (последовательное увеличение значений индексов г, ], к)

Л„/,* — Д-1 Г 3 - + Qгj,fc I /отт

Дфг-,к — В [--^^-Уг-,к + + (д)

+А+-1/2Дфг*-1-,к + А+-1/2Дфг*--1,к + А+-1/2Аф*^,к-1 ;

и -обход (последовательное уменьшение значений индексов г, ], к)

Дф- = Дф]к - В-1 [А~+1/2М1] + А4-+1/2Аф*]+1,к + А-+1/2Аф*],к+1^ . (4)

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

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

Дф* — В-1 Г 3 (®],к) - + QП-,к т/ + +

Дфг-,к — В Г---— + \ШЬ—к) + (5)

+Аи,2ЬФи- + А-+-1/2ЬФ*г--1

Дф]к — ДФ— - В-1 [А-1/2Дф*+1- + А-+1/2Афф1+1\. (6)

Областью решения в этом случае является прямоугольник, и тогда для наиболее полной загрузки процессоров за направление счета (по индексам) нужно выбрать направление, по которому имеется наибольшее число счетных ячеек. Пусть, например, это будет индекс г. Тогда решение уравнений (5) можно свести к совокупности решения одномерных (по индексу г) задач, реализуемых "бегущим счетом". При этом значения величин Дф-1 -, Дф*--1 определяются из граничных условий задачи как симметричные или антисимметричные значения в дополнительных фиктивных счетных ячейках. Область интегрирования разрезается на полосы, каждая из которых размещается в одном процессоре. Если каждой такой полосе приписать индекс ], то первым стартует процессор с номером ] — 0, затем, после того как он вычислил значение величины Дф0 0 (значения всех матриц и других величин рассчитываются с нижнего итерационного слоя в, и их расчет проблем не вызывает), запускается процессор с номером ] + 1. Затем, после расчета им величины Дф0 1, запускается процессор с номером ] — 2 и т. д. Общая структура вычислений имеет ступенчатую форму с отставанием на одну счетную ячейку по индексу ] (рис. 2).

Последовательный простой процессоров начинается с того момента, когда первый процессор (] — 0) достиг правой границы, поскольку для реализации и-обхода (6) необходимо

Рис. 2. Использование процессоров.

Рис. 3. Применение топологии типа "кольцо" для более полной загрузки процессоров.

вычислить все величины А—* ■. При реализации такой схемы распараллеливания на MPI важно выбрать оптимальное число процессоров, так как с увеличением их числа, с одной стороны, уменьшается время счета, а с другой — растут затраты на межпроцессорные обмены.

Используя топологию типа "кольцо" [4], можно достичь более полной загрузки процессоров при некотором одновременном снижении числа межпроцессорных обменов (рис. 3). В этом случае область изменения по индексу j делится на несколько частей (секций) и счет начинается с самой первой из них 0 < j < N0. Внутри этой секции процессоры распределяются, как и прежде, по каждому индексу j и счет организуется, как было описано выше.

При достижении первым процессором правой границы он переключается на счет нижней границы второй секции, все необходимые данные для которой уже были подготовлены No-м процессором первой секции. По мере достижения остальными процессорами первой секции правой границы они переключаются на счет своих j-линий во второй секции N0 + 1 < j < N1 и т. д., пока не будут вычислены все величины А—* j. Реализация U-обхода для расчета величины А—1 осуществляется аналогично, отличаясь лишь тем, что начинается с верхней правой счетной ячейки.

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

4. Результаты расчетов

Численные расчеты были выполнены в Центре высокопроизводительных вычислений (НЬКБ) в г. Штутгарт (Германия).

Ниже приведены описания двух использовавшихся кластерных систем: 1а64 и strider. Более подробная документация находится на сайте НЬИ^З:

— ia64: http://www.hlrs.de/hw-access/platforms/zx6000/

— strider: http://www.hlrs.de/hw-access/platforms/strider/

4.1. IA64

The IA64 cluster platform (table 1) consists of one front node for interactive access (ia64.hww.de) and several nodes for execution of parallel programs. The cluster consists of 8 IA64 dual CPU nodes with 900MHz ITANIUM2 CPU's + 8GByte memory on each node and a 4way frontend

T a b l e 1

Short overview of installed nodes of IA64 cluster

Name CPU Memory Disk NodeID PBS feature

HP zx6000 2 x Itanium2 900 MHz 8 GB 36GB ia64-1 - ia64-8 zx6000

HP rx5670 4 x Itanium2 1.0 GHz 48 GB 4x36 GB ia64.hww.de master

node with 4 ITANIUM2 1GHz CPU's + 48GByte memory. Additionally a RAID system with 1.7 TByte is available. The local disks of each node (26 GByte) serves as scratch disks. Features:

— Cluster of 8 dual SMPs nodes HP zx6000 servers with 8GByte memory

— Frontend node is a 4way HP rx5670 server with 48GByte memory

— Node-Node interconnect Myrinet 2000 Network + Gigabit Ethernet

— Disk 1.7 TByte home + 200 GByte scratch

— Batchsystem: OpenPBS, Maui scheduler

— PVFS scratch filesystem

— Operating System: RedHat Advanced Workstation 2.1

— HP Chipset ZX1

— Intel Compilers Main data:

— Peak Performance: 73.6 GFLOP/s

— Processors: 12 Memory: 64 + 48 GB = 112 GB

— Shared Disk: 1.7 TB (Raid5)

— Distributed Disks: 8 x 36 GB = 288 GB

— Number of Nodes: 8

— Node-node data transfer rate: 250 MB/s

4.2. Strider

The CRAY Opteron cluster (strider) platform (table 2) consists of one front node for interactive access (strider.hww.de) 125 nodes for execution of parallel programs and 2 I/O nodes. The cluster consists of 127 AMD dual CPU nodes with 2GHz Opteron CPU's + 4GByte memory on each node and a 2way frontend node with 2 Opteron 2GHz CPU's + 8GByte memory. Additionally a RAID system with 2 TByte is available. The local disks of each node (26 GByte) serves as scratch disks. Features:

T a b l e 2

Short overview of installed nodes of CRAY Opteron cluster

Name CPU Memory Disk NodelD PBS feature

compute node 2 x AMD Opteron 2000 MHz 4 GB 36GB c0-01 - c0-36, c1-02 - c1-28, c2-01 - c2-32, c3-02 - c3-31

master node 2 x AMD Opteron 2000 GHz 8 GB 2x70 GB strider.hww.de

ТаблицаЗ Ускорение счета

Число процессоров 2 6 12

Ускорение (IA64) 1.74 4.53 7.25

Ускорение (strider) 1.69 3.15 5.86

— Cluster of 127 dual SMPs (AMD Opteron's) servers with 4GByte memory

— Frontend node is a 2way (AMD Opteron) server with 8GByte memory

— Node-Node interconnect Myrinet 2000 Network

— Disk 1000 GB home + 1000 GB scratch + 3250 GByte local scratch

— Batchsystem: OpenPBS, Maui scheduler

— PBSPro, CRAY Utils

— Operating System: SuSE SLES 8 (AMD64) (United Linux 1.0)

— PGI Compilers Main data:

— Peak Performance: 1048 GFLOP/s

— Processors: 256

— Memory: 516 GB

— Shared Disk: 2000 GB (Raid5)

— Distributed Disks: 3250 GB

— Number of Nodes: 125 compute, 1 frontend, 2 I/O

— Node-node data transfer rate: 250 MB/s

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

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

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

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

[1] Ryohkoy A.D., Shokina N.Yu., Neutz J., Eisenreich N. Numerical modeling of processes in combustion camber of automotive airbags // 33rd Intern. Annual Conf. of ICT, Karlsruhe, Germany, 2002.

[2] ЗЕльдович Я.Б., ЛЕйпунский О.И., Ливрович В.Б. Теория нестационарного горения пороха. М.: Наука, 1975. 132 с.

[3] Yoon S., Jameson A. An LU-SSOR scheme for Euler and Navier — Stokes equations // AIAA 87-600, January, 1987.

[4] Корнеев В.Д. Параллельное программирование в MPI. М.; Ижевск: Ин-т компьютерных исследований, 2003. 304 с.

Поступила в редакцию 21 октября 2004 г.