CC BY
37
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
Том XXIII 1992 М I
УДК 533.6.011.35 : 629.7.025.1 + 681.3.06
МОДЕЛИРОВАНИЕ ТРАНСЗВУКОВОГО ОБТЕКАНИЯ ИЗОЛИРОВАННОГО КРЫЛА НА МНОГОПРОЦЕССОРНОМ КОМПЛЕКСЕ ИЗОТ-1014, ЕС-2706
А. С. Дебабов, А. В. Смирнов
Рассмотрен пример использования многопроцессорного комплекса для аэродинамического расчета трансзвукового обтекания крыла. Показана возможность уменьшения времени счета в 20 раз и более. Приведено решение практических задач.
Появление вычислительных комплексов с периферийными матричными процессорами дает возможность эффективной реализации многомерных систем моделирования сложных аэродинамических процессов, а также позволяет существенно повысить быстродействие существующих и разрабатываемых пакетов прикладных программ. Тем не менее в силу нетрадиционности архитектуры таких комплексов, для того чтобы потенциальная возможность могла быть реализована, необходим тщательный подход к вопросам организации счета. Это включает в себя оптимизацию ввода-вывода и согласования его- со счетом, распараллеливание алгоритма вычислительной модели, учитывающее архитектурные особенности матричного процессора и особенности вычислительной программы.
В данной статье рассматривается реализация системы моделирования трансзвукового обтекания изолировэнного стреловидного крыла на базе вычислительного комплекса, включающего управляющую ЭВМ (УМ) ИЗОТ-Ю14 и несколько матричных процессоров (МП) ЕС-2706.
Управляющая ЭВМ имеет производительность порядка 1,2 млн. операций в секунду для чисел с плавающей запятой, оперативную память 16 Мбайт и способна передавать данные по каналу со скоростью 2,5 Мбайт/с. Она обеспечивает интерфейс «человек — ЭВМ» и предназначена для разработки. Отладки и диагностирования программного обеспечения на уровне связи с МП. В задачу УМ входит также первоначальная загрузка микропрограмм матричного процессора, запуск программы и подкачка данных в МП. Матричный процессор ЕС-2706 (аналоги: №5-100, РРБ-164, АР-120В, АР-190Ь и другие) представляет собой высокопроизводительный синхронный конвейерно-параллельный микропрограммируемый арифметический процессор. Он содержит (рис. 1) несколько запоминающих устройств, подпитывающих два конвейерных арифметических устройства (умножитель и сумматор) через многочисленные магистрали данных. МП управляется синхронно от одного тактового генератора с периодом 167 не (6 МГц). Разрядность обрабатываемых данных 38 бит в формате с плавающей точкой. Команды в процессоре ЕС-2706 64-разрядные, и каждая команда управляет работой всех устройств машины, причем выполнение одной команды осуществляется с максимальной скоростью 6 млн. команд в секунду и эквивалентно одновременному выполнению до 10 функций. Такой принцип управления известен как горизонтальное микропрограммирование. Максймальная скорость вычислений составляет 12 млн. арифметических операций плюс *§0 млн. команд в секунду (адресные и другие вспомогательные операции).
Положительный опыт использования многопроцессорной системы описан в работе [1] на примере анализа ряда двумерных задач физики плазмы, где достигнуто увеличение быстродействия в 10 раз при использовании одного МП и в 32 раза при работе на четырех МП. Далее приведен пример использования МП для ускорения выполнения алгоритма аэродинамического расчета трансзвукового обтекания крыла, известного как РЬО-22 [2] в модификации [3]. Отличие в реализации поставленной задачи от [1] заключается в способе программирования на МП. Так, при адаптации алгоритма работы [1] на МП использовались стандартные подпрограммы библиотеки АРМАТСН, реализующей более 400 различных действий над векторами.
■I____________
Дисплейная станция
ЕС-7927 | Щ-ті\ [ЕС-7927
Контроль диско!
Г
ЕС-50ВЗ\ \ЕС-50ВЗ\ IЕС-5063 ----Г—1 *—Т-1 ----------------
Контроль iucxot .....~~l---------
—|fg-27^| |fg-27^| [fy-27^ I IEC-2706
Рис. 1. Конфигурация многопроцессорного комплекса ЕС-1037, ЕС-2706
В описанном ниже алгоритме не удается выделить длинные линейные участки программы, в которых выполнялись бы элементарные действия над векторами. Поэтому при программировании был применен ассемблер МП АРАЬ, позволяющий распараллелить на уровне МП алгоритм любой трудности и достичь при этом использования не менее 80% ресурсов матричного процессора.
Программа расчета трансзвукового обтекания крыла основана на решении уравнения полного потенциала, моделирующего течение идеального газа:
(а2 — <й + (а2 — <р») <р„,, + (а2 — ф|) фгг — 2 ф, ф,, фЛ!) — 2 ф, фг ф„, — 2ф, фг ф« ±= 0, (1)
где а — скорость звука; ф — потенциал скорости.
Потенциальная модель не приводит к большим погрешностям в расчете при условии малых значений числа М перед скачком уплотнения, замыкающим местнук)' сверхзвуковую зону. Уравнение (1) решается совместно с граничными условиями:
Зф/дп = 0; ду/дг |г=о = 0;
Ф = х cos а — у sin а,
(2)
где первое — условие непротекания на теле (ft — нормаль к телу), второе — условие в плоскости симметрии; третье — условие в дальнем поле.
Используемые в (1), (2) величины ф, а отнесены к невозмущенной скорости, X, у, г — к некоторому линейному размеру.
Решение (1), (2) осуществляется конечно-разностным методом в преобразованной системе координат х, у, г, полученной с помощью ряда конформных и сдвиговых преобразований [2, 3], переводящих пространство х, у, г в х, у, г. После дискретизации уравнение (1) примет вид (3):
AGxx -|- BGyy + CGzi + DGXy + EGXZ -i- FGyZ R — 0, (3)
где G = ф — x cos a — у sin а; А, В, C, D, E, F, R — коэффициенты, возникающие в процессе преобразования уравнения (1) к его представлению в координатах х, у, г.
Решение уравнения (3) осуществляется методом верхней релаксации по линиям x(i) — const, g(j) = var, z(k) — const. При этом за первый проход происходит расчет невязки, а при обратном ходе — уточнение значения потенциала. Затем осуществляется постановка граничных условий и переход к новому значению i. После расчета нового значения потенциала в сечении k = const алгоритм повторяется в сечении k + 1. При расчете невязки в алгоритме предусмотрен механизм изменения конечно-разностного шаблона. Так, для точек, находящихся в дозвуковой области течения, он имеет форму «крест», а для точек в сверхзвуковой зоне смещен вверх по потоку. При этом происходит выделение в расчете областей с большими градиентами параметров потока. Программная реализация данного механизма приводит к ветвлению алгоритма и его усложнению.
Анализ затрат машинного времени центрального процессора на различных участках программы показал, что 85% времени используется на расчет значения невязки при решении сеточных уравнений на линии. Поэтому очевидно, что для уменьшения времени счета необходимо
передать функцию вычисления невязки матричному процессору. Участок программы расчета невязки содержит около 100 операций умножения и 80 операций сложения (вычитания). 50 операций сложения (вычитания) предназначены для вычисления адресов элементов используемых массивов. В процессе расчета из оперативной памяти необходимо выбирать до 70 элементов массивов. Распараллеливание вычислений в МП позволяет в идеальном случае реализовать эти 300 действий (100 + 80 + 50 + 70) в 100 тактов работы МП с частотой 6 МГц за счет одновременного умножения, суммирования, вычисления адресов и считывания данных из памяти МП. Основная трудность при достижении минимального числа рабочих тактов МП состоит в обеспечении непрерывного конвейера вычислений с точной их синхронизацией. В реализованной на МП программе число тактов, необходимых для расчета невязки, составило 180, что больше идеального прогноза и объясняется значительной сложностью программы.
Изменение быстродействия алгоритма в случае адаптации на МП только блока расчета невязки в одной расчетной точке показано на рис. 2. Если длина вычислительного вектора (число расчетных точек) равна единице, время счета увеличивается в 10 раз. Такой отрицательный результат станет понятен при анализе временных затрат всего алгоритма работы управляющей ЭВМ совместно с матричным процессором. Суммарное время работы алгоритма можно представить как
7" = Т'сис + ^обм + Тмп-
где Tnn — время работы матричного процессора; Говм — время обмена данными; Тси<: — время загрузки в МП микропрограмм.
Возникшая ситуация объясняется тем, что, во-первых, вызов программы в МП происходит очень часто и тратится много времени на многократную загрузку микропрограмм; во-вторых, релаксационный метод является неявным и для обеспечения подпрограммы расчета невязки требуется значительное количество подкачиваемой информации.
Таким образом, имеет место соотношение 7"сис » 7мп- То6и » 7\цП. Чтобы увеличить длину вычислительного вектора, реализовано решение уравнения в МП на всей вертикальной линии. Для этого на ассемблер APAL были переложены алгоритмы расчета поправки к значению потенциала и постановки граничных условий и ряд модулей, обеспечивающих решение во всей расчетной области при k = const. Это привело к значительному увеличению загружаемой в.МП, программы, что требует большого (порядка 1 мс) времени 7’сис, но при этом число вызовов сокращается. На рис. 2 приведена зависимость изменения времени работы программы от длины обрабатываемого вектора £мп. Видно, что при LMп, равном 10 (число точек на прямой д = var, х — const), ускорение практически отсутствует. Это соответствует ситуации, когда Гсис й: ~ Тмп и тобм ~ Тмп- Ускорение достигается при расчете в МП поля потенциала во всем сечении г = const.•что соответствует длине вычислительного вектора около 2000 элементов. Для дальнейшего увеличения длины вектора (до 50 000—все расчетное пространство) требуется использование памяти данных МП до 1 Мбайт, что, в принципе, возможно, но приведет к значительному усложнению при программировании на языке APAL. На рис. 2 показано возможное в данном случае ускорение (примерно в 20 раз), реализованное при тестовом прогоне программы без получения достоверного результата. Страничная организация памяти данных МП позволяет во время решения задачи вести обмен информацией между УМ и не используемыми при счете страницами. Это дает возможность снизить потери времени, связанные с обменом, и достичь максимальной производительности и при длине вычислительного вектора /.мп = 2000.
Дальнейшее увеличение эффективности программы расчета крыла основывается на распараллеливании счета на несколько МП. Так, если в ЭВМ последовательного действия необходимо провести расчет по всем сечениям 2 — const и затем перейти к новой итерации, то наличие нескольких присоединенных МП позволяет переходить к новой итерации, не дожидаясь окончания предыдущей. В этом случае модифицированные в одном МП данные для сечения z(k) = const используются в другом МГ1 для расчета в сечении z(k — 1) = const.
Время одной итерации счета программы в секундах на сетке, содержащей около 50 000 узлов, на различных ЭВМ составляет:
micro VAX-И •...........................70
ИЗОТ-1014...............................55
BASF-7/78...............................6,1
1МПЕС-2706 ......................... 3,2
2МПЕС-2706 ...............1,7
Сравнение показывает, что применение одного или нескольких ускоряющих счет процессоров может значительно повысить эффективность программы по сравнению с ЭВМ средней мощности (рис. 3). Причем необходимо подчеркнуть, что использование процессоров эффективно для пакетов прикладных программ, которые применяются длительное время при массовом расчете и практически не модифицируются. Это объясняется большими затратами программиста на адаптацию программ для матричного процессора и сложностью вносить изменения в отлаженный алгоритм, используемый в МП. Быстрота переложения ФОРТРАН-программы на язык АРА1, оценивается разработкой 15—20 микрокоманд в день.
Разработанный вариант программы расчета трансзвукового обтекания крыла с использованием МП ЕС-2706 может быть применен при проведении в короткие сроки большого параметрического счета. Расчетные модули, адаптированные для МП и реализующие решение уравнения полного потенциала в преобразованной системе координат, могут применяться при создании эффективных быстродействующих алгоритмов, например задачи аэродинамического про-
о
~15 ' 30
30
Рис. 2
Число операций с плаНающей точкой
Рис. 3
ектирования (аналога приведенного в работе [4] ), учет трехмерного пограничного слоя на базе алгоритма работы [3], учет влияния на аэродинамику крыла струи от винта по работе [5].
Целью данной работы являлось исследование адаптации вычислительного алгоритма для матричного процессора и ускорения решения аэродинамической задачи. Поэтому в качестве примера того, как реально можно уменьшить время получения итогового результата, рассмотрены две задачи: расчет зависимости сх(М) при фиксированном значении су и различных геометрических характеристиках крыла и расчет неоднородностей потока, вызванных наличием крыла в зоне расположения винта самолета.
Трудоемкость первой задачи объясняется тем, что для получения, например, пяти расчетных точек при сравнении трех различных геометрий крыла требуется провести не менее 45 расчетов (при фиксированном М для обеспечения требуемого су необходимо три расчета). Реально затрачиваемое время для получения графика, приведенного на рис. 4, составляет около 1 —1,5 недели при средней загрузке ЭВМ. При использовании версии программы, адаптированной на МП, данные кривые были получены менее чем за 4 часа. Схема проведения расчета была следующей. Параллельно на счет запускались две программы с независимым использованием двух МП, по полученным значениям су производилась корректировка угла атаки с целью получения во втором запуске на одном МП результата при требуемом с„. Число запусков пру этом составило 30 вместо 45.
Вторая задача отличается тем, что использование МП позволяет практически в интерактивном режиме, задавая режим обтекания крыла, получать поля течения в зоне расположения винтов и производить оценку изменения угла притекания в сечении лопасти винта Да:
где Й — скорость вращения винта; г — текущий радиус лопасти винта; Уа, Ув — осевая и азимутальная компоненты скорости.
Положение диска винта следующее: хо = — 0,11, уо = — 0,055, го = 0,25; нормаль к диску винта выбрана таким образом, чтобы при крейсерском режиме полета М = 0,75 и о =3,925 неоднородность набегающего потока была минимальной. Результат расчета Да, представленный на рис. 5, соответствует изменению условий работы винта при варьировании М с фиксированным
0,03
о/н
о хш2Г,\-в;,т]-з • 23°, •» * *1,0°
а 18°, »
Рис. 5
V 0,В 0,7 0,8 М„ Рис. 4
значением с„. Расчет производился по схеме, описанной выше, и занял не более одного часа. Результаты, показанные на рис. 5, ввиду трудоемкости их получения практически никогда. не исследуются, хотя учет влияния изменения обтекания крыла приводит к локальному изменению угла притекания потока к лопасти на 0,5—1°. Согласно оценкам работы [5J это может привести к изменению динамической нагрузки в сечении лопасти на 3—5% по сравнению с нагрузкой, полученной в однородном потоке.
ЛИТЕРАТУРА
1. Велихов Е. П., Лазаров В л., Марков Ст., Петров И. В., Р о г а л ь с к и й А. В., С а г д е е в Р. 3., Ч о ч и а Г. А., Шевченко В. И. Решение' нелинейных задач физики на многопроцессорном комплексе ЕС-1037 — ЕС-2706 // Препринт ИКИ. — 1986. № 1169.
2. Jamson A., Caughey D. A.. Numerical calculation of transonic flow past a swept wing.— N-Y: University. ERDA report COO-3077-140, 1972.
3. Карась О.В.Довалев В. E. Применение обратного метода расчета трехмерного пограничного слоя к задаче обтекания крыла с учетом влияния вязкости // Ученые записки ЦАГИ.— 1985. Т. 15, № 4.
4. Raj P. Amulti-grid method for transonic wing analysis and design // J. of Aircraft.— 1984. Vol. 21, N 2.
5. Смирнов А. В. Разработка метода расчета влияния струи винта на трансзвуковое обтекание крыла // Ученые записки ЦАГИ.— 1990. Т. 21, № 2.
Рукопись поступила 20/IX 1990 г.
