CC BY
84
Вестник ПНИПУ. Аэрокосмическая техника. 2014. № 39
УДК 532.3:519.6
К.С. Кузьмина, И.К. Марчевский
Московский государственный технический университет им. Н.Э. Баумана,
Москва, Россия
ОБ УСКОРЕНИИ ВЫЧИСЛЕНИЙ ПРИ РЕШЕНИИ ДВУМЕРНЫХ СОПРЯЖЕННЫХ ЗАДАЧ ГИДРОУПРУГОСТИ ВИХРЕВЫМИ МЕТОДАМИ
Описано применение вихревых методов к решению задач о численном моделировании динамики подвижного или деформируемого профиля в потоке несжимаемой среды; при этом профиль может двигаться по заданному закону либо перемещаться и деформироваться под действием гидродинамических нагрузок. Приведено краткое описание вихревых методов, указано на альтернативные подходы к учету граничных условий на поверхности профиля: с помощью равенства нормальных (N-схема) либо касательных (Г-схема) компонент скоростей среды и профиля. В обоих случаях задача сводится к численному решению интегрального уравнения относительно интенсивности вихревого слоя, генерируемого на профиле, и системы обыкновенных дифференциальных уравнений, описывающих движение вихревых элементов в области течения. Описаны преимущества использования Г-схемы по сравнению с N-схемой. Разработан программный комплекс POLARA-TVEM, позволяющий проводить расчеты на многопроцессорных вычислительных комплексах с использованием технологии MPI. Выявлены наиболее затратные с вычислительной точки зрения операции вычислительного алгоритма; приведены оценки эффективности их распараллеливания при использовании технологий MPI, OpenMP и CUDA. Также в программном комплексе POLARA-TVEM реализован алгоритм быстрого метода приближенного решения задачи N тел, использование которого позволяет существенно сократить временные затраты на проведение вычислений. Представлена оценка вычислительной сложности быстрого метода и всего алгоритма вихревого метода с его использованием. Суммарное ускорение расчетов, полученное при решении вычислительно сложных задач в результате использования распараллеленного алгоритма быстрого метода, достигает нескольких сотен и даже тысяч раз при проведении расчетов на кластерах, содержащих несколько десятков вычислительных ядер.
Ключевые слова: вихревой метод, вязкая несжимаемая жидкость, сопряженная задача гидроупругости, профиль, вихревой след, отрывное обтекание, интегральное уравнение, задача N тел, параллельный алгоритм, вычислительный кластер.
K.S. Kuzmina, I.K. Marchevskiy
Bauman Moscow State Technical University, Moscow, Russian Federation
ON COMPUTATION SPEEDING UP WHEN SOLVING TWO-DIMENSIONAL HYDROELASTIC COUPLED PROBLEMS BY USING VORTEX METHODS
The application of vortex methods to the problem of numerical simulation of the movable or de-formable airfoil in incompressible flow is investigated. The law of the airfoil's motion can be either given or it can move and deform under hydrodynamic loads. A brief description of vortex methods is given and two alternative approaches to the satisfaction of the boundary conditions on the surface of the airfoil are pointed out. These approaches are based on the equality of normal (W-scheme) or tangential (T-scheme) component of the flow and airfoil velocity. In both approaches the problem is reduced to numerical solution of the integral equation with respect to the intensity of the vortex layer generated on the airfoil, and the system of ordinary differential equations describing the vortex elements motion of in the flow. The advantages of using the T-scheme in comparison with the W-schemes are described. The program package POLARA-TVEM, which allows to carry out computations using multiple-processor computers, using MPI-based parallel algorithm is developed. The computational costs of the basic operations in the algorithm of vortex element method are calculated and estimates for the parallelization efficiency are given for MPI, OpenMP and CUDA technologies. The algorithm for the fast method of approximate solution of 'N bodies' problem is also implemented in POLARA-TVEM software package. Its usage allows reduce significantly time of computing. The estimations of the computational cost of the fast method algorithm and the entire vortex method algorithm are given. The total acceleration of computations obtained when solving computationally difficult problem resulting from the usage of parallelized fast algorithm reaches several hundred or even thousands of times in the calculations on clusters containing several tens of cores.
Keywords: vortex element method, viscid incompressible flow, coupled aeroelastic problem, airfoil, vortex wake, separated flow, integral equation, N-body problem, parallel computational algorithm, multi-processor cluster.
Введение
Задачи вычислительной аэрогидродинамики, в особенности сопряженные задачи аэрогидроупругости, традиционно относятся к числу наиболее сложных; необходимость их решения служит мощным стимулом развития как численных методов, так и (в некоторой мере) вычислительной техники. Среди множества разработанных к настоящему времени методов и подходов к численному моделированию аэ-рогидроупругих явлений следует выделить вихревые методы [1-3], в частности метод вязких вихревых доменов [1], которые относятся к классу лагранжевых методов и не требуют построения сетки в области течения.
Вихревые методы имеют достаточно узкую область применения, ограниченную малыми дозвуковыми скоростями течений, поскольку
они основаны на модели несжимаемой среды, однако в указанных задачах они могут применяться весьма успешно и заметно превосходить по эффективности другие (сеточные) методы. Наибольшее развитие и обоснование к настоящему времени получили вихревые методы решения двумерных задач, однако даже они на сегодняшний день не реализованы в коммерческих программных комплексах.
В связи с этим разработка эффективных алгоритмов решения сопряженных задач аэрогидроупругости, учитывающих современные достижения в области вихревых методов, а также их программная реализация, предполагающая возможность распараллеливания вычислений и проведения расчетов на многопроцессорных вычислительных комплексах, является актуальным направлением исследований.
Краткое описание метода вихревых элементов
Метод вихревых элементов (МВЭ) - это бессеточный лагранжев метод численного моделирования нестационарного обтекания профилей, основанный на моделировании эволюции завихренности. Распределение завихренности моделируется набором вихревых элементов (ВЭ), для каждого из которых задается положение в пространстве и интенсивность. Скорость среды в любой точке течения представляет собой суперпозицию влияний отдельных ВЭ, вычисляемых по закону Био - Савара, а давление определяется при помощи аналога интеграла Коши - Лагранжа [1] или иными способами [2].
Эффективность МВЭ связана с тем, что при решении многих задач область с отличной от нуля завихренностью сравнительно невелика; она расположена вблизи профиля и позади него. Использование завихренности в качестве первичной расчетной величины позволяет сосредоточить вычислительные ресурсы именно в этой области, а не «распылять» их на всю расчетную область. Особенно эффективным МВЭ является при решении внешних задач обтекания, поскольку на каждом шаге расчета по времени выполнение уравнения неразрывности и граничного условия на бесконечности происходит автоматически.
Граничное условие непротекания или прилипания обеспечивается генерацией ВЭ на поверхности обтекаемых профилей, а течение идеальной либо вязкой среды, описываемое уравнениями Эйлера либо Навье - Стокса, моделируется движением имеющихся ВЭ. Интенсивности генерируемых ВЭ находятся из решения соответствующей сис-
темы линейных алгебраических уравнений, аппроксимирующей граничное условие на профиле, а движение ВЭ описывается системой обыкновенных дифференциальных уравнений, правые части которой представляют собой скорости движения ВЭ. Движение и/или деформация обтекаемого профиля, которые могут либо происходить по заданному закону, либо определяться нестационарными гидродинамическими нагрузками в случае решения сопряженной задачи гидроупругости, моделируется присоединенным вихревым слоем и присоединенным слоем источников, интенсивности которых определяются текущими величинами скоростей точек профиля.
Таким образом, алгоритм решения задачи по расчету обтекания профиля содержит несколько основных этапов, представленных на рис. 1 и выполняемых на каждом шаге расчета по времени.
Движение ВЭ, образующих вихревой след
А
Определение нового положения точек профиля и их скоростей
тг
Вычисление нагрузок, действующих на профиль
Рис. 1. Основные этапы расчета обтекания профиля методом вихревых элементов
Наиболее затратным с вычислительной точки зрения является этап расчета движения ВЭ. При этом происходит операция вычисления скоростей всех вихревых элементов, требующая учета взаимного влияния всех пар ВЭ друг на друга, а также несколько менее трудоемкая, но также достаточно затратная операция реструктуризации вихревой пелены, выполняемая на этапе моделирования движения вихревого следа. Основным направлением ускорения расчетов представляется применение параллельных алгоритмов, а также приближенных быстрых методов типа быстрого метода в задаче N тел. Аналог задачи N тел в методе вихревых элементов возникает при вычислении скоростей движения ВЭ, индуцируемых другими вихревыми элементами.
Расчет интенсивностей присоединенных слоев вихрей и источников
Д
Генерация свободных ВЭ на обтекаемом профиле
Д
Пополнение вихревого следа сгенерированными ВЭ
Математическая модель
Движение вязкой несжимаемой среды описывается уравнениями Навье - Стокса:
V-V = 0, —+(V -V)V = vAV-Vp, dt ' р
где V - скорость среды, V = V(r,t); р - давление, p = p(r,t); р -плотность, принимаемая постоянной, р = const; v - коэффициент кинематической вязкости.
Граничные условия на бесконечности (условия затухания возмущений) имеют вид
V(r, tV, p(r, tpx при | r|
где V„ и p^ - постоянные скорость и давление в набегающем потоке; на профиле K выполняется условие прилипания
V(r, t) = Vk (r, t), r e K.
Используя вектор завихренности П (r, t) = VxV(r, t), уравнение Навье - Стокса можно записать в форме Гельмгольца [1]:
^ + Vx(ilxU) = 0, (1)
dt
где U (r, t ) = V (r, t) + W (r, t), W (r, t) - диффузионная скорость:
(Vxn jxfl
W (r, t ) = v-
Inl2
Уравнение (1) означает, что генерации новой завихренности в области течения £ не происходит, а осуществляется только движение существующей вдоль векторных линий поля и; новая завихренность генерируется только на профиле К.
Влияние профиля заменяется суперпозицией влияний присоединенных (связанных с профилем) вихревого слоя уаи (г, I) и слоя источников ^ (г, ^) и свободного вихревого слоя у(г, ^). Все эти слои расположены на профиле, при этом
1ап (м) = V (r, г) • Т (r, г) даа (г г) = V (r, г) • п (r, г) г е к,
где т (г, г) и п (г, г) - орты касательной и внешней нормали к профилю в соответствующей точке, причем п(г, г)хт(г, г) = к , где к - орт нормали к плоскости течения.
Скорость течения может быть вычислена по обобщенному закону Био - Савара:
, ч 1 г "(£, г )х(г - £) 1 (£, г )х(г - £)
V(г, г)=V. + — / ^ > £ 2 ^ ^+— § ^ ^ £ 2 ^ а + 2п 4) |г- £| 2п К() |г- £|
+ ± § У а, (£, г )х(г - £ ) + ^ X да, (£, г )(г - £ ) .
2п 4) |г - £|2 5 2п |) |г - £|2 5
Здесь £() и К() - область течения и контур профиля в момент времени г; у аи =уайк и у = ук - векторы интенсивностей присоединенного и свободного вихревых слоев на профиле соответственно.
Если рассматривать поле скоростей V(г, г) на всей плоскости Я2,
то на обтекаемом контуре оно будет иметь разрыв. Можно показать [3], что предельное значение вектора скорости среды со стороны профиля на его границе
V- (г, г) = V (г, г)-у(1% г)-2У -(г г) т (г, г)+п (г, г). (3)
Интенсивность свободного вихревого слоя у(£, г) может быть
найдена из условия прилипания на поверхности тела, которое с учетом представлений (2), (3) имеет вид
V- (г, г) = V (г, г), г е к.
Для решения задачи могут быть использованы различные расчетные схемы: классическая ^схема [1, 2] и модифицированная Г-схема [3, 4].
Традиционные подходы к обеспечению выполнения указанного выше граничного условия, которые обычно используются в расчетных схемах метода вихревых элементов, основаны на том, что неизвестная
функция ) находится из условия равенства нормальных составляющих предельного значения скорости среды V_(r ) и скорости обтекаемого профиля VK (r), т.е. фактически из условия непротекания:
V_(r )• n (r ) = Vk (r )• n (r ), r G K.
Данное равенство приводит к сингулярному интегральному уравнению с ядром Гильберта, являющемуся уравнением первого рода, интеграл в левой части которого следует понимать в смысле главного значения по Коши.
Из работы [3] известен альтернативный подход к обеспечению выполнения граничного условия: показано, что оно с математической точки зрения эквивалентно условию равенства касательных компонент предельного значения скорости среды V_ (r) и скорости профиля VK (r) :
V_(r )• т (r ) = Vk (r). т (r ), r g K.
Это уравнение, в отличие от предыдущего, приводит к интегральному уравнению второго рода. Однако принципиально важно то, что в случае гладкого обтекаемого профиля K ядро уравнения является ограниченным, т.е. по сути оно представляет собой уравнение Фред-гольма 2-го рода.
Данные подходы получили название N-схемы и Г-схемы соответственно. В работе [4] показано, что Г-схема позволяет получать решение для более широкого класса задач по сравнению с N-схемой, она также обладает более высокой точностью. К важным особенностям разработанной Г-схемы следует отнести ее «симметрию»: в случае поступательно движущегося в потоке жесткого профиля она позволяет получать одинаковое решение в прямом и обращенном движении.
Программный комплекс POLARA-TVEM
Для численного решения описанного класса задач разработан параллельный программный комплекс POLARA-TVEM, позволяющий проводить расчеты на многопроцессорных вычислительных комплексах с использованием технологии MPI. Проведенные оценки показали, что наиболее затратными с вычислительной точки зрения являются следующие операции алгоритма МВЭ:
Операция 1. Вычисление скоростей ВЭ как результата их взаимного влияния.
Операция 2. Реструктуризация вихревой пелены: объединение близкорасположенных ВЭ и исключение ВЭ, удалившихся от профиля на значительное расстояние.
Операция 3. Вычисление нагрузок, которые действуют на профиль со стороны потока, по рассчитанному полю скоростей и известному распределению завихренности.
На рис. 2 представлены оценки вычислительной трудоемкости различных этапов последовательного алгоритма МВЭ.
Наиболее трудоемким является вычисление скоростей ВЭ. Другие две операции менее трудоемки, в то время как затратами машинного времени на прочие операции (например, решение системы линейных уравнений) можно пренебречь. Согласно закону Ам-дала для вычислительной системы с £ процессорами максимально возможное ускорение программы с долей Р параллельного кода 1
Рис. 2. Трудоемкость операций в методе вихревых элементов
а =
р+(1 - Р)
Результаты расчета максимального ускорения, которое теоретически достижимо при распараллеливании только операций 1, операций 1, 2 или операций 1, 2, 3 алгоритма, например для вычислительных систем с 4, 16 и 64 ядрами, представлены в табл. 1.
Таблица 1
Максимальное ускорение параллельного кода
Распараллеленные операции Доля параллельного кода, % Максимальное ускорение
4 ядра 16 ядер 64 ядра
1 95,3 3,5 9,4 16,2
1, 2 98,2 3,8 12,6 30,0
1, 2, 3 99,5 3,9 14,9 48,7
Для вычислений с использованием небольшого числа ядер максимальные значения ускорений близки, но при использовании системы с десятками вычислительных ядер распараллеливание операций 2 и 3 метода вихревых элементов дает значительный вклад в ускорение, несмотря на то, что трудоемкость этих операций составляет вместе лишь малую долю от трудоемкости операции 1 (приблизительно 4,4 %).
Реальное ускорение, наблюдаемое в расчетах, несколько меньше, поскольку оценки по закону Амдала получены для «идеального» распараллеливания, т.е. ситуации, когда параллельный код выполняется в £ раз быстрее на ^-ядерной вычислительной системе.
На рис. 3 для программного комплекса РОЬАЯА-ТУЕМ представлена зависимость ускорения от числа задействованных вычислительных ядер при распараллеливании операций 1, 2 и 3 алгоритма, полученная при решении одной, весьма сложной с вычислительной точки зрения задачи расчета обтекания профиля. Для сравнения показаны максимальные значения ускорения из оценок по закону Амдала.
О 8 16 24 32 40 48 56 64 72 80
Рис. 3. Реальное ускорение и оценки по закону Амдала
Реальное ускорение при распараллеливании операций 1, 2, 3 близко к оценке по закону Амдала для распараллеливания только операций 1, 2 (сказывается наличие пересылок данных, затрат на синхронизацию и прочих «накладных расходов»), поэтому можно сделать вывод о том, что распараллеливание операции 3 важно для достижения высокой эффективности, хотя ее трудоемкость крайне мала - всего около 1,3 %.
Отметим, что интенсивность обменов данными между задействованными в расчете ядрами является весьма существенной, поэтому для повышения эффективности работы программного комплекса была исследована возможность применения для этих целей технологии
OpenMP. С ее помощью были распараллелены те же три наиболее тру-дозатратные операции, что и при использовании MPI.
При проведении расчетов на вычислительных машинах с общей памятью на 2-8 вычислительных ядрах было установлено, что выигрыш по времени, достигаемый за счет использования OpenMP, не превышает нескольких процентов по сравнению с применением MPI. В связи с этим данная возможность практически не использовалась, поскольку применение MPI позволяет сделать программу универсальной - она может работать на системах как с общей, так и с распределенной памятью (кластерах).
Также в вопросе повышения скорости расчетов представляет большой интерес использование вычислительных возможностей современных графических ускорителей путем применения технологии NVidia CUDA. Ее использование требует существенной переработки вычислительных алгоритмов МВЭ, однако данные затраты могут быть вполне оправданными: в работе [5] на примере решения тестовой задачи о расчете обтекания неподвижного кругового профиля показано, что время выполнения одного шага расчета на персональной ЭВМ с графической картой NVidia GTX670 оказывается в 2,8 раза меньше, чем при проведении расчета на 16-ядерном вычислительном кластере [6] (1,4 с против 4 с; при проведении расчетов на кластере использовалась описанная выше схема распараллеливания).
Таким образом, при решении методом вихревых элементов вычислительно сложных задач, когда вихревой след состоит из сотен тысяч ВЭ, наиболее подходящим представляется гибридный алгоритм, способный работать на кластерных системах с распределенной памятью, когда взаимодействие между узлами кластера осуществляется с помощью технологии MPI, а вычисления внутри каждого узла выполняются на графических ускорителях. Эффективная реализация такого алгоритма, однако, представляет собой нетривиальную задачу в силу особенностей алгоритма и имеющихся в нем информационных зависимостей. Данное направление является предметом дальнейших исследований.
Применение быстрого метода расчета вихревого влияния
Помимо применения параллельных алгоритмов, в программном комплексе POLARA-TVEM предусмотрена возможность использования так называемого быстрого метода в задаче N тел, позволяющего сократить временные затраты на выполнение наиболее трудоемкой
операции алгоритма - вычисление скоростей вихревых элементов, обусловленных их взаимным влиянием друг на друга. Быстрые методы решения задачи N тел основаны на том, что при расчете влияния удаленных ВЭ их коллективное воздействие рассчитывается приближенно. Первоначально такой подход был предложен в работе [7] для решения гравитационной задачи; в работе [8] он был адаптирован применительно к вихревым методам.
В основе алгоритма быстрого метода лежит построение дерева -иерархической структуры прямоугольных областей (рис. 4), при этом прямоугольник нулевого уровня содержит все ВЭ. Он делится по длинной стороне на два одинаковых прямоугольника первого уровня. Путем перебора вихревых элементов определяется их принадлежность к одному из прямоугольников. Далее аналогичным образом эти прямоугольники делятся пополам, образуя области второго уровня. Деление прекращается при выполнении заданного критерия по размеру прямоугольника, количеству ВЭ в нем и (или) номеру уровня - глубине дерева. На втором этапе непосредственно вычисляются суммарные циркуляции и «центры тяжести» положительных и отрицательных ВЭ в каждом прямоугольнике нижнего уровня, а аналогичные параметры областей более высокого уровня находятся по характеристикам дочерних областей. Далее находятся скорости ВЭ в областях нижнего уровня: влияние ВЭ, находящихся в той же области, а также близкорасположенных ВЭ из других областей рассчитывается непосредственно по закону Био - Савара, а далее осуществляется обход дерева, и влияние ВЭ, расположенных в достаточно удаленных областях, учитывается при помощи четырех коэффициентов, которые одинаковы для всех ВЭ из данной области.
Рис. 4. Структура дерева, имеющего глубину четыре уровня, и направление его обхода
Точность вычисления скоростей определяется значением параметра 0, который определяет, насколько далеко по сравнению с их собственными размерами должны находиться друг от друга области, соответствующие двум ячейкам дерева, чтобы их влияния можно было вычислять по приближенным формулам. Как показывает практика, его целесообразно выбирать из промежутка 0 <0< 1, при невысоких требования к точности значение 0 может быть увеличено до 2, при этом меньшие значения параметра обеспечивают большую точность и, соответственно, требуют больших затрат времени на проведение расчета.
Приближенная оценка вычислительной трудоемкости алгоритма быстрого метода, которая полагается равной количеству выполняемых операций умножения и деления при вычислении скоростей ВЭ, получена в работе [9] и имеет вид
б =
(
0(72)к
1 -а
1+1п Г ^
(У2)к -1 4Ш
лЛЛ
\\
+
JJ
(4)
V
4 + 0
+ 4п.
JJ
Здесь п - количество ВЭ в области течения; к - максимальная глубина дерева; а ив - некоторые параметры, значения которых подобраны в работе [9] экспериментально, исходя из результатов решения тестовых задач (а = 0,70, р = 1,30).
Оценка (4) является более точной, чем, например, оценка, приведенная в работе [10]. Точность оценки (4) хорошо иллюстрируется графиком (рис. 5), на котором нанесена сплошная кривая, отражающая зависимость (4), и результаты вычислительного эксперимента, в котором моделировалась эволюция вихревого следа, состоящего из п = 157 674 ВЭ, при различных значениях параметра 0.
Отметим, что процедура построения дерева имеет исключительно низкую вычислительную сложность (десятые-сотые доли процента от общего времени выполнения расчета), поэтому при оценке сложности всего алгоритма ею можно пренебречь. Более подробно данный вопрос исследован в работе [11] при моделировании пространственных течений вихревыми методами, однако выводы указанной работы полностью справедливы и в двумерном случае.
Q
5хЮ9 5х109 5х109 5хЮ8
5хЮ8 5хЮ8
0,5 1,0 1,5 2,0 9
Рис. 5. Общее количество операций умножения/деления, требуемое для выполнения одного шага расчета в зависимости от величины параметра 0: линия - теоретическая оценка (4); точки -результат вычислительного эксперимента
При фиксированном значении параметра 0 в алгоритме быстрого метода остается, по существу, единственный параметр, значение которого не влияет на точность вычисления скоростей ВЭ, но влияет на время проведения расчета. При использовании формулы (4) путем непосредственных расчетов легко убедиться в том, что оптимальная глу-*
бина дерева k , обеспечивающая наименьшую вычислительную сложность алгоритма, практически не зависит от величины 0 - критерия близости/дальности ячеек, однако существенно зависит от значения n -количества ВЭ в расчетной схеме. При оптимальном выборе максимальной глубины дерева вычислительная сложность быстрого метода пропорциональна n log n, тогда как вычислительная сложность «прямого» алгоритма определяется формулой Q0 =6n2.
В табл. 2 приведено отношение вычислительных трудоемкостей «прямого» алгоритма МВЭ и алгоритма, в котором используется быстрый метод. Расчеты проведены по формуле (4) для случаев, когда количество ВЭ в вихревом следе равно 40 000, 160 000 и 600 000 при трех различных значениях параметра 0.
Видно, что использование быстрого метода позволяет сократить время выполнения расчетов в десятки-сотни раз в зависимости от решаемой задачи. На практике именно использование быстрого метода делает возможным решение нестационарных задач при n > 50 000 за разумное время порядка нескольких часов-суток.
Таблица 2
Ускорение алгоритма МВЭ за счет использования быстрого метода
Отношение вычислительной сложности QsJQ п = 40 000 к* = 13 п = 160 000 к* = 15 п = 600 000 к* = 17
0 = 0,2 5,9 20,2 67,7
0 = 0,6 50,0 179,3 613,1
0 = 1,0 139,4 503,5 1727,8
Анализ показывает, что ошибка в определении максимальной глубины дерева даже на 1-2-м уровнях ведет к существенному (до 1,5-2 раз) росту времени вычислений. В табл. 3 приведены полученные в расчетах зависимости среднего времени выполнения алгоритма быстрого метода от максимальной глубины дерева (время счета при оптимальной глубине дерева принято за единицу). Значение параметра 0 было принято равным 0,5. Все расчеты проводились в последовательном режиме.
Таблица 3
Зависимость времени счета от глубины дерева
Число уровней дерева Время счета tk^tk.
п = 40 000 к* = 13 п = 160 000 к* = 15 п = 600 000 к* = 17
12 1,37 5,95 -
13 1,00 2,52 9,92
14 1,43 1,49 4,84
15 2,70 1,00 2,33
16 - 1,43 1,26
17 - 2,86 1,00
18 - - 1,49
При п = 40 000 расчеты проводились лишь при к < 15, а при п = 160000 - только при к < 17, поскольку задание большей максимальной глубины дерева будет означать, что процесс его построения практически всегда будет останавливаться, когда в ячейке будет находиться один вихревой элемент, а сама ячейка, следовательно, имеет нулевой размер. В то же время формула (4) получена в предположении, что ячейка имеет конечные размеры и в ней находятся несколько ВЭ.
Результаты проведенных расчетов хорошо согласуются с теоретической оценкой (4); различие обычно не превышает 10-12 %; в редких случаях наблюдается расхождение в 15-20 %, которое можно объяснить, по-видимому, особенностью формы вихревого следа за профилем. Впрочем, даже если принять 20%-ную погрешность формулы (4), во всех случаях можно верно установить оптимальную глубину дерева.
Соотношение трудоемкостей различных этапов алгоритма МВЭ при использовании быстрого метода, полученное при решении той же задачи, что была рассмотрена выше (см. рис. 2), имеет следующий вид (рис. 6).
В сравнении с предыдущими оценками видно, что доля затрат на вычисление скоростей ВЭ сокращается, хотя и осталась доминирующей.
Использование имеющейся иерархической структуры дерева в вихревом следе позволяет пред-дожить новые, существенно более
Рис. 6. Трудоемкость операций в методе эффективные алгоритмы выпол- вихревых элементов при использовании
нения остальных операций МВЭ. быстрого метода
В частности, операция реструктуризации вихревого следа, предполагающая поиск близкорасположенных ВЭ, становится значительно менее трудоемкой, поскольку поиск «соседних» ВЭ отныне можно производить только в «близких» ячейках дерева. Данные эффекты вызывают изменение соотношений тру-доемкостей операций на рис. 6 по сравнению с рис. 2, когда все действия выполнялись «прямым» методом.
Описанные выше принципы МР1-распараллеливания вычислений остаются применимыми и в случае использования быстрого метода. При этом разделение задачи по вычислительным ядрам производится не по вихревым элементам, скорости которых предстоит определить, а по вершинам дерева нижнего уровня. Отметим, что в «прямом» методе операции вычисления скоростей каждого ВЭ имеют абсолютно одинаковую вычислительную трудоемкость, что позволяет обеспечить равномерную загрузку вычислительных ядер при проведении расчета в параллельном режиме. В случае же использования быстрого метода,
Операция 3
во-первых, различные ячейки дерева содержат различное число ВЭ (это можно учесть сравнительно просто), а во-вторых, время вычислений существенным образом зависит от расположения соответствующего прямоугольника в вихревом следе. Неучет этого обстоятельства несколько снижает равномерность загрузки вычислительных ядер. Тем не менее реализованные алгоритмы позволяют получить существенное ускорение расчета в параллельном режиме. Ниже приведены полученные в расчетах в ходе решения одной тестовой задачи значения ускорения вычислений при МР1-распараллеливании быстрого метода.
Ускорение алгоритма МВЭ, использующего быстрый метод, при проведении расчета в параллельном режиме:
Число ядер 1 2 4 8 16
Ускорение, раз 1,00 1,58 2,56 4,31 5,97
Суммарное ускорение вычислений за счет использования алгоритма быстрого метода и технологии MPI, таким образом, может составлять несколько сотен и даже тысяч раз, что позволяет существенным образом расширить возможности применения вихревых методов решения задач аэрогидродинамики и аэрогидроупругости.
Заключение
В работе описаны подходы, позволяющие повысить эффективность расчетов при решении задач прямого численного моделирования в гидродинамике и гидроупругости методом вихревых элементов. Указаны способы повышения точности применяемых математических моделей и разработан программный комплекс POLARA-TVEM, позволяющий проводить расчеты в параллельном режиме на вычислительных кластерах.
В программном комплексе реализована возможность применения быстрого метода решения задачи N тел. Построена оценка вычислительной сложности алгоритма, хорошо согласующаяся с результатами практических расчетов.
Работа выполнена в рамках Государственного задания № 2014/104 в сфере научной деятельности в части Государственного задания Ми-нобрнауки России (проект 1712) при финансовой поддержке гранта Президента РФ (проект МК-3705.2014.8).
Библиографический список
1. Андронов П.Р., Гувернюк С.В., Дынникова Г.Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. - М.: Изд-во МГУ, 2006. - 184 с.
2. Cottet G.-H., Koumoutsakos P.D. Vortex Methods: Theory and Practice. - Cambridge University Press, 2008. - 328 p.
3. Accuracy considerations for implementing velocity boundary conditions in vorticity formulations / S.N. Kempka, M.W. Glass, J.S. Peery, J.H. Strickland // SANDIA REPORT. SAND96-0583, UC-700. - 1996. - 52 p.
4. Kuzmina K.S., Marchevsky I.K. On Numerical Schemes in 2D Vortex Element Method for Flow Simulation Around Moving and Deformable Airfoils // Advanced Problems in Mechanics: Proc. of the XLII School-Conference. - Saint Petersburg, 2014. - P. 335-344.
5. Гречкин-Погребняков С.Р. Использование технологии CUDA для моделирования обтекания кругового профиля методом вихревых элементов [Электронный ресурс] // Молодежный научно-технический вестник. - 2014. - № 4. - URL: http://sntbul.bmstu.ru/doc/719651.html (дата обращения: 19.09.2014).
6. Лукин В.В., Марчевский И.К. Учебно-экспериментальный вычислительный кластер. Ч. 1. Инструментарий и возможности // Вестник Моск. гос. техн. ун-та им. Н.Э. Баумана. Сер.: Естественные науки. -2011. - № 4. - С. 28-43.
7. Barnes J., Hut P. A hierarchical O(N log N) force-calculation algorithm // Nature. - 1986. - Vol. 324, № 4. - P. 446-449.
8. Дынникова Г.Я. Использование быстрого метода решения «задачи N тел» при вихревом моделировании течений // Журнал вычислительной математики и математической физики. - 2009. - Т. 49, № 8. -С.1458-1456.
9. Кузьмина К.С., Марчевский И.К. Оценка трудоемкости быстрого метода расчета вихревого влияния в методе вихревых элементов // Наука и образование: электрон. науч.-техн. изд. - 2013. - № 10. -С.399-414.
10. Морева В.С. Способы ускорения вычислений при решении плоских задач аэродинамики методом вихревых элементов // Вестник Моск. гос. техн. ун-та им. Н.Э. Баумана. Сер.: Естественные науки. -2011. - № S. - С. 83-95.
11. Марчевский И.К., Щеглов Г.А. Моделирование динамики вихревых структур высокопроизводительным методом вихревых элементов // Известия вузов. Машиностроение. - 2013. - № 9. - С. 26-36.
References
1. Andronov P.R., Guvernyuk S.V., Dynnikova G.Ya. Vikhrevye me-tody rascheta nestatsionarnykh gidrodinamicheskikh nagruzok [Vortex methods for unsteady hydrodynamic loads computation]. Moskovskiy gosu-darstvennyy universitet, 2006. 184 p.
2. Cottet G.-H., Koumoutsakos P.D. Vortex Methods: Theory and Practice. Cambridge University Press, 2008. 328 p.
3. Kempka S.N., Glass M.W., Peery J.S., Strickland J.H. Accuracy considerations for implementing velocity boundary conditions in vorticity formulations. SANDIA REPORT. SAND96-0583, UC-700, 1996. 52 p.
4. Kuzmina K.S., Marchevsky I.K. On Numerical Schemes in 2D Vortex Element Method for Flow Simulation Around Moving and Deformable Airfoils. Advanced Problems in Mechanics: Proc. of the XLII School-Conference. Saint Petersburg, 2014, pp. 335-344.
5. Grechkin-Pogrebnyakov S.R. Ispolzovanie tekhnologii CUD A dlya modelirovaniya obtekaniya krugovogo profilya metodom vik-hrevykh ehlementov [CUDA technology application for the flow simulation around a circular airfoil by using the vortex element method]. Molodezhnyy nauchno-tekhnicheskiy vestnik, 2014, no. 4, available at: http://sntbul.bmstu.ru/doc/719651.html (accessed 19 September 2014).
6. Lukin V.V., Marchevskiy I.K. Uchebno-ehksperimental'nyy vychis-litel'nyy klaster. Instrumentariy i vozmozhnosti [Computing cluster for training and experiments. Part 1. Instrumentation and opportunities]. Vestnik Moskovskogo gosudarstvennogo tekhnicheskogo universiteta imeni N.E. Baumana. Estestvennye nauki, 2011, no. 4, pp. 28-43.
7. Barnes J., Hut P. A hierarchical O(N log N) force-calculation algorithm. Nature, 1986, vol. 324, no. 4, pp. 446-449.
8. Dynnikova G.Ya. Ispolzovanie bystrogo metoda resheniya "Zadachi N tel" pri vikhrevom modelirovanii techeniy [Fast technique for solving the N-body problem in flow simulation by vortex methods]. Zhurnal vychislitelnoy matematiki i matematicheskoy fiziki, 2009, vol. 49, no. 8, pp. 1458-1456.
9. Kuzmina K.S., Marchevskiy I.K. Otsenka trudoemkosti bystrogo metoda rascheta vikhrevogo vliyaniya v metode vikhrevykh elementov [Es-
timation of computational complexity of the fast numerical algorithm for calculating vortex influence in the vortex element method]. Nauka i obra-zovanie, 2013, no. 10, pp. 399-414.
10. Moreva V.S. Sposoby uskoreniya vychisleniy pri reshenii ploskikh zadach aehrodinamiki metodom vikhrevykh elementov [Approaches to computations acceleration when solving plane problems of aerodynamics by usint the vortex element method]. Vestnik Moskovskogo gosudarstvennogo tekhnicheskogo universiteta imeni N.E. Baumana. Estestvennye nauki, 2011, no. S, pp. 83-95.
11. Marchevskiy I.K., Shheglov G.A. Modelirovanie dinamiki vikhrevykh struktur vysokoproizvoditelnym metodom vikhrevykh ehlementov [Modeling the dynamics of vortex structures by the vortex element method]. Izvestiya vysshikh uchebnykh zavedeniy. Mashinostroenie, 2013, no. 9, pp. 26-36.
Об авторах
Кузьмина Ксения Сергеевна (Москва, Россия) - студентка 5-го курса, лаборант-исследователь кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (105005, г. Москва, ул. 2-я Бауманская, 5, e-mail: [email protected]).
Марчевский Илья Константинович (Москва, Россия) - кандидат физико-математических наук, доцент кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана (105005, г. Москва, ул. 2-я Бауманская, 5, e-mail: [email protected]).
About the authors
Kseniya S. Kuzmina (Moscow, Russian Federation) - 5th year student, Laboratory Assistant, Department of Applied Mathematics, Bauman Moscow State Technical University (5, 2nd Baumanskaya st., Moscow, 105005, Russia, e-mail: [email protected]).
Ilya K. Marchevskiy (Moscow, Russian Federation) - Ph. D. in Physics and Mathematical Sciences, Associate Professor, Department of Applied Mathematics, Bauman Moscow State Technical University (5, 2nd Baumanskaya st., Moscow, 105005, Russia, e-mail: [email protected]).
Получено 1.10.2014
