ДОКЛАДЫ АКАДЕМИИ НАУК, 2014, том 459, № 3, с. 280-284
МАТЕМАТИЧЕСКАЯ ФИЗИКА
УДК 519.6
ВЫСОКОСКОРОСТНЫЕ ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ РЕШЕНИЯ ЗАДАЧ МЕХАНИКИ СПЛОШНОЙ СРЕДЫ
МЕТОДОМ СГЛАЖЕННЫХ ЧАСТИЦ © 2014 г. Н. М. Евстигнеев, Ф. С. Зайцев, О. И. Рябков
Представлено академиком Д.П. Костомаровым 08.04.2014 г. Поступило 16.04.2014 г.
БО1: 10.7868/80869565214290040
1. ВВЕДЕНИЕ
Механика сплошной среды имеет большое число важных приложений в физике, химии, технике и других областях. Глубокий уровень современных исследований требует применения технологий высокоскоростных вычислений. Поэтому разработка эффективных численных алгоритмов решения задач механики сплошной среды является крайне актуальной.
В последние годы широкое распространение получил метод сглаженных частиц — Smoothed Particle Hydrodynamics (SPH). Достаточно полные обзоры метода с развернутым анализом различных аспектов алгоритма представлены, например, в [1, 2]. Метод SPH имеет ряд преимуществ перед традиционными сеточными методами и методами Монте-Карло. К достоинствам SPH относится следующее:
1) возможность сохранения свойств инвариантности относительно преобразования Галилея и обратимости во времени;
2) отсутствие сетки по фазовому пространству. В сеточных методах построение стационарной сетки для объектов сложной геометрии или адаптивной сетки для процессов с высокой динамикой и локализацией — обычно сложная задача;
3) простое и естественное описание предписанных движущихся границ, свободных границ, границ разделов сред, разрывов искомых функций. В сеточных методах задачи с такими услови-
ООО "Нью Инфлоу", Москва Московский государственный университет им. М.В. Ломоносова Институт системного анализа Российской Академии наук, Москва Научно-исследовательский институт системных исследований Российской Академии наук, Москва
ями требуют применения сложной техники для получения адекватной аппроксимации;
4) возможность простого описания множества узко локализованных движущихся взаимодействующих объектов. В сеточных методах для решения подобных задач необходимо использовать динамически перестраиваемую неравномерную сетку, сильно сгущающуюся в области локализации объектов;
5) простое описание взаимодействующих многокомпонентных сред с источниками и стоками, в том числе вызванными ядерными и химическими реакциями, а также систем с фазовыми переходами вещества;
6) сведение задачи к системе обыкновенных дифференциальных уравнений. Отсутствие принципиальных сложностей при решении трехмерных задач и задач большей размерности. Отметим, что в сеточных методах, как правило, аппроксимируют непосредственно уравнения в частных производных в эйлеровом описании среды и сводят задачу к обращению матрицы системы линейных алгебраических уравнений большой размерности, что обычно вызывает значительные вычислительные проблемы;
7) по сравнению с методом Монте-Карло метод 8РИ дает гладкие зависимости искомых функций от фазовых переменных и времени, а также обычно требует значительно меньше элементов дискретизации;
8) относительно простой для программирования численный алгоритм. Эффективное распараллеливание для архитектур ЭВМ с центральными процессорами и для гетерогенных архитектур с графическими ускорителями. Возможность создания для 8РИ специальных аппаратных ускорителей.
Указанные свойства делают 8РИ особенно привлекательным для решения задач динамики сплошных сред.
К недостаткам 8РИ относят несколько меньшую по сравнению с сеточными методами точ-
ность описания ударных волн с резкой границей разрыва [1], хотя для таких случаев в SPH разработаны способы повышения точности. Как и в сеточных методах, в SPH имеется задача обеспечения численной устойчивости и выбора параметров дискретизации исходной задачи, в частности, числа и начального расположения SP-частиц, параметров ядра интегральной аппроксимации. Кроме того, в SPH возможно возникновение специфической численной неустойчивости, которая называется эластичной (tensile).
В реальных системах, например в условиях плазменного вихревого реактора [3] или токамака [4], протекают десятки разнообразных физических и химических процессов. Для каждого процесса из-за его сложности обычно разрабатывается своя математическая модель и свой, как правило достаточно специфический, численный метод. Подход SPH позволяет в рамках единой дискретной модели охватить большое количество явлений. Впервые появляется возможность полноценного самосогласованного описания динамики системы.
Целью исследований является разработка и изучение свойств новых высокоскоростных алгоритмов решения задач механики сплошной среды методом сглаженных частиц (SPH) для двух типов архитектур супер-ЭВМ: с центральными процессорами (CPU) и графическими ускорителями (GPGPU).
CPU-архитектуры обычно используют технологию параллельных вычислений MPI / OpenMP. К таким архитектурам относятся, например, компьютеры "Ломоносов" и "Blue Gene/P" в МГУ им. М.В. Ломоносова, имеющие ~104 CPU-ядер, а также проектируемые эксафлопные вычислительные комплексы с ~106 ядер.
В гибридных архитектурах применяются специализированные графические или другие акселераторы, работающие по технологиям OpenCL или CUDA. Такие архитектуры позволяют за счет усложнения алгоритма добиться высокой производительности на относительно недорогой и компактной аппаратуре.
В настоящей работе также получены оценки погрешности аппроксимации интегралов методом SPH. Отсутствие подробного теоретического изучения погрешностей в общем случае затрудняло распространение метода SPH.
2. ПОГРЕШНОСТЬ АППРОКСИМАЦИИ.
ВЫБОР ПАРАМЕТРОВ SPH
Рассматривается класс задач для уравнений вязкой сплошной среды. Система уравнений приведена, например, в [3]. Там же описан эффективный сеточный метод ее решения, даны тестовые примеры и программное обеспечение. Здесь для решения задач применяется метод SPH [1].
Основная идея метода 8РИ состоит в аппроксимации искомой векторной функции Аг) с помощью интеграла
A(r, t ) = JA( t, r' ) W(r - r', h ) dr' +
(1)
с некоторым заданным в области О ядром Щг — г', к) и сведение операторов от искомой функции к операторам от ядра. Здесь г = (хь х2, х3),
dг' обозначает элемент объема ^г' = dх\ dх2 dх3, у — погрешность интегральной аппроксимации, изученная, например в [1, 2]. Величина к называется параметром сглаживания, в общем случае к может быть функцией времени и пространства к = к(?, г). Способы автоматического выбора к обсуждены в [1]. В методе 8РИ интеграл (1) заменяется суммой [1] А(га(0) =
= £ А(гь«)) ЩГа(0 - Гь((), к)ДГь(0 + у + ,
Ь: Гь(0 е П (2)
где у2 — погрешность аппроксимации квадратурной формулы.
В классическом подходе 8РИ [1] для объема Дгь(0 в (2) используется выражение
р(0, Гь(0))ДГь(0)
Ar»( t) =
(3)
р(Гь(0)
В литературе по 8РИ применение уравнения (3) обосновывается физической интерпретацией — сохранением массы р(1, г(?))Дг(0, объема Дг(?) во времени. Однако интуитивные соображения не дают возможности вычисления погрешности у2. Математически строгим обоснованием правомерности использования выражения (3) является его формальное следование из соотношения для эволюции любого объема сплошной среды [6, с. 123] и уравнения непрерывности в отсутствие источников и стоков.
Обозначив ть = р(0, гь(0))Дгь(0), из (2) получаем
А(I, га(0) = V тьА(?' Гь ( 0 ) х
ь V П ьР('' Гь(<))
ь: Гь(() е П
х Щ(Га(^ - Гь(0, к) + у + у, (4)
Ранее погрешность у2 оценивалась для частных случаев [1]. Обсудим вопрос о величине у2 в общем случае. В начальный момент времени разбиение объема О на элементы Дгь(0) можно выбрать в соответствии с той или иной квадратурной формулой, например повторной формулой трапеций. Величина у2 в начальный момент будет иметь соответствующую этой формуле оценку. Однако с течением времени формы областей Дг(?) могут измениться и разбиение О станет отличным от исходного. Для произвольного разбиения области О в предположении существования интеграла (1) можно гарантировать лишь сходи-
Q
282
ЕВСТИГНЕЕВ и др.
мость частичных сумм в (2) при стремлении диаметра разбиения к нулю, т.е. стремление погрешности vz к нулю при уменьшении диаметра разбиения:
Vz = o (1). (5)
В случае конкретных динамических систем может иметь место сходимость более высокого порядка [1].
На практике, как и в сеточных методах, изучение vz и подбор параметров SPH проводятся на модельных задачах с заранее известными решениями.
Получим общую оценку погрешности yz для другого способа задания элемента объема Arb(t). Рассмотрим сначала одномерный случай.
Лемма. Пусть функции f(x'), W(x — x', h) заданы на отрезке x' е [x — p1h, x + p2h], где h, p1, p2 — не зависящие от x' положительные константы, x е е R — фиксированная точка. При этом произведение fx') W(x — x', h) на концах отрезка обращается в ноль. Рассмотрим на отрезке [x — p1h, x + p2h] точки xb, b = 0, 1, ..., N + 1, расположенные в порядке возрастания с ростом индекса b, так что x'0 = = x — p1h, +1 = x + p2h.
Справедлива квадратурная формула
x + в2 h
I(x) = J f(x') W(x - x', h) dx' =
X - Pih
N
£f(xb) W(x - xb, h) x + + ^,
b = 1
M2 (h) 12 h3
NhM = O
f i2 ^ hM
v h2
M2 (h )= max (2, bMw^ o, bh +
b = 1, 2,..., N + 1
+ 2Mf 1, bMWt 1, bh + M, o, bMw, 2, b),
hM =
ем Нм) в области интегрирования О, связанного, например, с увеличением общего числа 8Р-ча-стиц, погрешность стремится к нулю со вторым порядком по 1/^ (по Нм);
2) стремление параметра сглаживания Н к нулю при фиксированном N (фиксированном числе 8Р-частиц) может приводить к увеличению погрешности в связи с уменьшением числа точек N в области интегрирования О;
3) для обеспечения уменьшения погрешности
при Н —► 0 необходимо одновременно увеличивать общее число 8Р-частиц (уменьшать Нм) так, чтобы возрастало число точек N, попадающих в область интегрирования О.
Повторное применение леммы для трехкратного объемного интеграла дает
А(гй(0) = X А((V Х
Ь: г4(0 е П Х Щга(() - Гь((), А)(АГь(())ь
(Arb (t)) tr = П
ч, bk
/tr + V + Vz,
It) - ^t, bt„ ( t)
Vz
k = 1
= O
(6)
(7)
(8)
max (xb — x'b -1) — максимальный шаг,
b = 1, 2,...,N + 1
константы M, k, b и Mw, k,
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.