CC BY
11
УДК 681.326
О МЕТОДАХ ОПТИМИЗАЦИИ И РАСПАРАЛЛЕЛИВАНИЯ ЗАДАЧ ВЫЧИСЛИТЕЛЬНОЙ ГИДРОДИНАМИКИ
ФИЛАТОВ Е.Ю., асп.
Решается задача моделирования обтекания воздухом одиночного здания в двумерном случае. Для решения данной задачи предлагаются методы оптимизации, позволяющие увеличить точность и сократить время расчета, и методы распараллеливания.
Ключевые слова: гидродинамика, метод сгущающихся сеток, метод асинхронного интегрирования по времени, распараллеливание.
OPTIMIZATION METHODS AND PROBLEM PARALLELING OF COMPUTATIONAL HYDRODYNAMICS
FILATOV E.Yu., postgraduate
The article concerns the problem of airflow around separate building simulation in two dimension case. The optimization methods are suggested to solve this problem which enables to increase the accuracy and to reduce the calculation time and paralleling methods.
Key words: hydrodynamics, thickened networks method, asynchronous time integration method, paralleling.
Постановка задачи. Дано одиночное здание в форме параллелепипеда и дана скорость ветра, направление которой параллельно одной из сторон здания. Требуется рассчитать поля скорости и давления вокруг данного здания.
Задача решается в двумерном случае, т.е. в плоскости, параллельной поверхности земли. Ниже будет предложено решение данной задачи, а также методы оптимизации и распараллеливания для предложенного решения. Распространение изложенных методов на трехмерный случай не представляет большого труда.
Решение задачи. Для решения данной задачи используется система «давление - скорость». Для
нахождения поля скорости и решаются уравнения Навье-Стокса:
5Ux г*. .. .. 1 5P
— = -UVU + vAUv - -
dt
dU,
= -Uvuv dt v
p dx
AM 1 dP
vAUv —X". V P dy
(1)
(2)
Для нахождения поля давления Р решается уравнение Пуассона:
AP = -р
dD
"дГ
+ UVD -vAD + Q
2 2 g(j. dU:
Q = 22-^ -zj =
i=1j=1 dXj dx,
gUx I2+2 dUx dUy
gx J gy gx
fgUy \ gy
(3)
(4)
й = ^и. (5)
В качестве начальных условий для поля скорости на всей рассчитываем ой области берется некоторое
начальное значение ио , а для поля давления - некоторое начальное значение Р0.
В качестве граничных условий для расчета поля скорости на выходной границе (рис. 1, правая граница) используется граничное условие II рода, при котором производная по нормали к границе равна о, на остальных свободных границах используются граничные условия I рода, при которых скорость равна
ио , а на здании используется граничное условие
I рода, при котором скорость равна 0. При решении уравнения Пуассона для нахождения поля давления везде (на свободных границах и на здании) используются граничные условия II рода, при которых производная по нормали к границе равна 0.
Поскольку рассматриваемое течение турбулентно (т.к. число Рейнольдса велико), то требуется также рассчитывать поле турбулентной вязкости. Для этого используется модель турбулентности А.Н. Секундова. Турбулентная вязкость в данной модели рассчитывается из следующего уравнения:
gv,
turb
gt
. у и gvturb
H j gxj
2
= H
j=1 —j
g
=1 gx;
gv
+v
turb
чат
turb
Murb J
turb
gxj
(6)
8vi
lam
D T^min (vlam + Pvturb)v
turb’
где
(
D =
ygUi
U gxj
gu, gu
gxj gx,
1/2
f(z) = 0,2
z2 + 1,47z + 0,2.
(7)
(8)
22 -1,472 +1
а = 2; р = 0,06; у = 1; (9)
Ьт1п - кратчайшее расстояние до ближайшей неподвижной стенки; V|am - ламинарная вязкость.
Таким образом, вязкость V, входящая в уравнения Навье-Стокса и Пуассона, состоит из ламинарной и турбулентной составляющих:
у = у1ат +"^игЬ. (10)
Для того чтобы поле турбулентной вязкости сформировалось, требуется задать его начальное значение отличным от 0, например,
у°игЬ = Ю00У|ат. (11)
В качестве граничных условий на входной границе (рис. 1, левая граница) и на здании используются граничные условия I рода, при которых турбу-
лентная вязкость равна 0, а на остальных свободных границах используются граничные условия II рода, при которых производная по нормали к границе равна 0.
Методы оптимизации. Изменение искомых полей в рассматриваемой области происходит неравномерно. Наиболее динамично изменение происходит вблизи здания. В удалении же от здания поля меняются довольно медленно. Кроме того, вблизи здания желательно иметь большую точность расчета, так как поля вблизи здания в значительной степени влияют на расчет полей во всей области. Исходя из этого, расчет полей был бы более эффективным, если бы шаг сетки и шаг по времени вблизи здания были меньше, чем в удалении от здания. Этого позволяют добиться методы сгущающихся сеток и асинхронного интегрирования по времени.
Метод сгущающихся сеток состоит в том, что в разных частях рассматриваемой области выбирается разный шаг по расстоянию. Так, разделим всю область на 3 зоны (рис. 1). Малая зона будет расположена вокруг здания, средняя зона - вокруг малой зоны и большая - вокруг средней. В большой зоне возьмем шаг сетки И, в средней - И/2, а в малой - И/4.
Рис. 1. Метод сгущающихся сеток и метод асинхронного интегрирования по времени
Кроме того, при уменьшении шага по расстоянию для сохранения устойчивости требуется также уменьшать шаг по времени, причем примерно в той же степени, что и шаг по расстоянию. Так, если в большой зоне шаг по времени равен т, то в средней зоне возьмем шаг по времени т/2, а в малой зоне -т/4. Метод, при котором в разных частях расчетной области используются разные шаги по времени, будем называть методом асинхронного интегрирования по времени. Шаг по времени т будем называть большим шагом по времени.
Упрощенно алгоритм расчета полей в условиях сгущающихся сеток выглядит следующим образом:
1. Расчет большого шага по времени в каждой зоне.
2. Обмен значениями полей на границах между зонами.
Под расчетом большого шага по времени понимается расчет 4-х шагов по времени т/4 в малой зоне и расчет 2-х шагов по времени т/2 в средней зоне. При расчете полей на границах зон используется граничное условие I рода и значения, полученные от соседней зоны.
Существуют два подхода в реализации метода сгущающихся сеток. Они различаются распределением узлов на границе между зонами. При 1-м подходе (рис. 2) граничные (0-е) узлы меньшей (я+1)-й зоны не совпадают с узлами большей р-й зоны. При 2-м подходе (рис. 3) граничные (0-е) узлы меньшей (Я+1)-й зоны совпадают с (А-1)-ми узлами большей
р-й зоны. Здесь граничными узлами большей р-й зоны являются А-е узлы.
Таким образом, обмен значениями полей на границах схематично можно представить следующим образом. Для 1-го подхода:
Ч +1 ^ ч:
О[Ч]А - С[ч+1]1'
-1:
Ч ^ Ч + ': ^Й+1]0
Для 2-го подхода:
Ч +1 ^ ч : О,
Ч ^ Ч +1:
_ О[ч]А-1
Э[Ч]А
[Ч]А _ °[Ч+1]2’
О[Ч+1]0 _ О[Ч]А-1.
(12)
(13)
(14)
(15)
Здесь О - данное поле; в квадратных скобках указывается номер зоны; последний индекс - номер узла по одной из координат.
А+2-А+1 ■
А А-1 •
5
4
2
0
д+1
Рис. 2. Метод сгущающихся сеток: 1-й подход
А+1 А А-1
А-2
4
3 О
1
0
Рис. 3. Метод сгущающихся сеток: 2-й подход
Если сравнивать эти два подхода, то реализация 2-го подхода будет заметно проще. Это связано с тем, что при 1-м подходе при передаче значений узлов на границе с р-й зоны на (р+1)-ю (13) используется дополнительная линейная интерполяция. Кроме того, при 2-м подходе в трехмерном случае расчет полей в приземной области около здания будет точнее, чем при 1-м. Однако отметим, что поскольку при 1-м подходе сетки перекрываются по одному узлу меньшей зоны, а при 2-м подходе - по двум узлам, то при 2-м подходе число узлов сеток будет больше, чем при 1-м, и следовательно, время расчета также будет больше.
Данные методы оптимизации - метод сгущающихся сеток и метод асинхронного интегрирования по времени - позволяют увеличить точность вычислений в тех местах рассчитываемой области, где она актуальна, и уменьшить время расчета вследствие уменьшения шагов по расстоянию и по времени в остальной части рассматриваемой области. Результаты расчета полей показаны на рис. 4.
Метод распараллеливания. Опишем метод распараллеливания в условиях изложенных методов оптимизации. Реализуем распараллеливание по пространству. Заметим, что структура сгущающихся сеток заметно усложняет алгоритм распараллеливания.
2
Рис. 4. Результаты расчета поля скорости и линий тока: стрелками изображено поле скорости, а серыми жирными линиями - поле функции тока, вычисленное из поля скоростей
Прежде всего отметим, что при распараллеливании используется клиент-серверная архитектура. Так, среди всех процессоров выбирается один процессор - сервер, который выполняет организационную работу, в частности, разделяет область расчета на участки между другими процессорами. Остальные процессоры - клиенты - производят расчет искомых полей, каждый на своем участке.
Предлагаемый метод распараллеливания включает в себя три алгоритма:
1) алгоритм разделения области расчета на участки между процессорами;
2) алгоритм регистрации границ;
3) алгоритм обмена значениями граничных узлов.
Остановимся подробнее на первом алгоритме.
Сначала все процессоры распределяются между зонами в зависимости от количества узлов в зонах. Число узлов в зонах подсчитывается с учетом асинхронного интегрирования по времени, т.е. в малой зоне число узлов умножается на 4, а в средней - на 2. Это связано с предположением, что время расчета пропорционально числу узлов в области расчета и числу выполняемых шагов по времени.
Процессоры между зонами следует распределять, по возможности, равномерно. Для этого каждой зоне выделяется сначала по одному процессору. Затем для каждой р-й зоны вычисляется величина Ы / рр, где Ы - число узлов в р-й зоне, а рр - число процессоров, ей принадлежащих. Той зоне, где данное отношение наибольшее, выделяется еще один процессор. Эти действия выполняются до тех пор, пока все процессоры не будут распределены.
Далее каждая зона разделяется на участки между принадлежащими ей процессорами. Для этого она разделяется сначала на вертикальные полосы, а затем каждая полоса разделяется горизонтальными линиями на участки (рис. 5).
Объясним, как это происходит. Очевидно, распараллеливание будет наиболее эффективным, если все участки будут содержать примерно равное число узлов и будут иметь форму квадратов. Квадратная форма участка сводит к минимуму периметр
участка, и следовательно, уменьшает время на обмены граничными узлами.
Рис. 5. Разделение зоны на участки между процессорами
Найдем, чему равна сторона такого квадратного участка. Пусть N - число узлов в зоне, а Р - число процессоров в зоне. Тогда т - среднее число узлов в участке:
т = Ы ; (16)
а - средняя сторона участка:
а = [\/т]. (17)
Отрежем вертикальную полосу шириной а. Данная ширина - предварительная. Впоследствии она будет уточняться. Подсчитаем п - число узлов в данной предварительной полосе. Заметим, что полоса может пересекать другую зону или собственно здание. Затем найдем р - число процессоров в данной полосе:
Р =
(18)
Здесь округление следует производить до ближайшего целого. Причем заметим, что полоса при пересечении с другой зоной или зданием может разделиться на две несвязные половины. В этом случае предпочтительнее, если число процессоров на ней р будет четным.
Затем следует уточнить ширину полосы для обеспечения равномерности разделения зоны на участки, поскольку при нахождении величины р про-
исходит достаточно сильное округление (если число процессоров, принадлежащих зоне, мало). Для этого находится уточненное число узлов на данной полосе п', равное произведению числа процессоров в полосе р и среднему числу узлов в участке т: п'= рт. (19)
Затем находится уточненная ширина полосы а', при которой полоса содержит примерно п' узлов.
Заметим также, что эффективность распараллеливания повышается, если граница полосы совпадает с вертикальной границей следующей зоны или здания, поскольку таким образом уменьшается общая протяженность границ.
После того как была найдена очередная полоса, ее следует разделить горизонтальными линиями на участки. Алгоритм такого разделения достаточно прост. Здесь лишь необходимо следить за тем, чтобы участки содержали примерно равное число узлов. Эффективность распараллеливания при этом также повышается, если горизонтальные линии совпадают с горизонтальными границами следующей зоны или здания.
Перед тем как начинать расчет полей, необходимо зарегистрировать все границы между соседними процессорами и сохранить их в реестре границ.
Алгоритмы регистрации границ и обмена границами идейно просты, однако технически довольно сложны. В частности, это связано с тем, что для каждого процессора существует 3 вида границ: граница с процессором из той же зоны, граница с про-
цессором из меньшей зоны и граница с процессором из большей зоны.
Алгоритм расчета полей при распараллеливании напоминает алгоритм расчета полей в однопроцессорном режиме в условиях сгущающихся сеток:
1. Расчет большого шага по времени каждым процессором.
2. Обмен значениями полей на границах между процессорами.
Данная задача была запущена на вычислительной системе МВС-100 на 6 процессорах (сервер и 5 клиентов). Эффективность распараллеливания получилась равной 76 %:
(20)
р
Е = -! Р рт,
Здесь Т1 - время расчета в однопроцессорном режиме; Тр - время расчета в многопроцессорном режиме на р клиентских процессорах; р - число клиентских процессоров.
Это показывает высокую эффективность предложенного алгоритма распараллеливания.
Список литературы
1. Роуч П. Вычислительная гидродинамика. - М.: Мир, 1980.
2. Численные методы и параллельные вычисления для задач механики жидкости, газа и плазмы: Учеб. пособие / Э.Ф. Балаев, Н.В. Нуждин, В.В. Пекунов и др.; Иван. гос. энерг. ун-т. - Иваново, 2003.
Филатов Евгений Юрьевич, ГОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», аспирант кафедры программного обеспечения компьютерных систем, телефон (4932) 26-98-29.
