Алгоритм динамического распараллеливания решения задачи адаптивного разбиения расчетной сетки для численного решения дифференциальных уравнений Текст научной статьи по специальности «Математика»

УДК 519.682.3

А. А. Кузнецов, В. А. Роганов, Г. А. Матвеев, В. И. Осипов

Алгоритм динамического распараллеливания решения задачи адаптивного разбиения расчетной сетки для численного решения дифференциальных уравнений

Аннотация. При численном решении дифференциальных уравнений непрерывная область решений ДУ заменяется дискретной решеткой, в узлах которой приближенно вычисляется значение искомой функции. В зоне ударных волн, межфазных границ и пограничных слоев при использовании регулярных сеток может быть резкий рост нормы производных и как следствие нормы ошибок аппроксимации, что ведет к потере точности численного решения. Для подавления роста ошибок аппроксимации используются адаптивные алгоритмы сгущения сеток в проблемных областях. В работе кратко описан подход к распараллеливанию такого алгоритма на основе концепции динамического распараллеливания «Т-система».

Ключевые слова и фразы: Т-система, динамическое распараллеливание, ОрепТБ, язык программирования Т+—+, дифференциальные уравнения, сетка, триангуляция области решений.

Введение

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

Работы, положенные в основу данной статьи, были выполнены в рамках НИР «Методы и программные средства разработки параллельных приложений и обеспечения функционирования вычислительных комплексов и сетей нового поколения» (№ 01201354596).

© А. А. Кузнецов, В. А. РогАнов, Г. А. Матвеев, В. И. Осипов, 2015 © Институт программных систем имени А. К. Айламазяна РАН, 2015 © Программные системы: теория и приложения, 2015

Общий вид линейного ДУ второго порядка с двумя независимыми переменными:

Аихх + Виху + Сиуу + Dux + Ей у + Fu = G.

Здесь коэффициенты А, В, С, D, Е, F, G являются функциями только от независимых переменных х и у. Зависимой переменной является и, индексы при и означают частные производные. Частное решение выбирается с помощью краевого и/или начального условия.

Численное решение такого дифференциального уравнения представляет собой тяжелую вычислительную задачу и не всегда выполнимо на одиночном компьютере за разумное время. Поэтому представляют интерес методы распараллеливания алгоритмов численного решения ДУ. В данной работе рассматривается подход к решению этой задачи, основанный на использовании идей «Т-системы» и возможностей системы OpenTS [1] в части динамического распараллеливания.

Поскольку нахождение аналитического решения ДУ в сложной области далеко не всегда возможно, для решения практических задач применяют численные методы. Одним из наиболее часто используемых является метод конечных разностей.

1. Конечно-разностная аппроксимация

Определение производной функции f (х) в точке X) записывается в виде:

w/ ч ,• f (хо + h) - f (хо)

f (хо) = lim---.

h^0 h

При численном дифференцировании на ЭВМ нельзя произвести предельного перехода, но можно придать шагу h некоторое малое положительное значение. Таким образом, производная заменяется разностью. Попробуем применить этот метод к уравнениям в частных производных [2].

Приближенно представим их(хо,уо) в виде конечной разности:

и(хо + h, уо) - и(хо,уо) h '

тогда ихх и иуу можно приближенно представить в виде:

и(хо + h, уо) - 2и(хо, уо) + и(хо - h, уо)

их(хо,уо) мож

с(хо, Уо) =

h2

и(хо, уо + к) - 2и(хо, уо) + и(хо, Уо - к)

)(хо,Уо) = к2

и

где к — величина шага сетки по оси х, к — величина шага сетки по оси у.

Рассмотрим пример — уравнение Лапласа в двумерном пространстве:

д2ф д2ф

дх2 + ду2

О

в некоторой области Д, и ф = /(х, у) на границе этой области.

Найдем численное решение уравнения Лапласа в прямоугольной области решений шириной А и высотой В. Разделим ширину прямоугольника на п интервалов, каждый размером к = А/щ так же разделим высоту В на т частей размером к = В/т. Внутри области получаются при этом (п — !)(т — 1) пересечений сетки (внутренних узлов). Затем для численного решения ДУ составляется разностное соотношение для каждого внутреннего узла и решается система уравнений.

Обозначим и(%к,]к) = , ](%к,]к) = ^^. Тогда граничное условие может быть представлено в виде:

Щ,о = До; г = 0,... ,п;

Щ,т !г,т; ^ 0, . . . , П;

ио,п = /о,1; 3 =0,...,т;

и.

п,3

; 3 =0,...,т.

Пусть теперь точка г,з будет точкой хо,уо из формул разностного представления производных функции и(х,у). Если обозначить А = к/к, то уравнение Лапласа сведется к разностному уравнению вида:

^2Щ+1 л + + «¿,¿+1 + Щ,з-\ — 2(1 + = О

для г = 1, .. . ,п — 1 и ] = 1, .. . ,т — 1.

При А =1 это соотношение означает что значение является средним арифметическим из четырех соседних с ним. Мы получили систему линейных алгебраических уравнений. Всего имеется (т — 1)(п— 1) уравнений относительно (ш+1)(п+1) неизвестных. После того как 2(т + п) неизвестных будут исключены с помощью граничного условия, остается точно (т — 1)(п — 1) уравнений относительно (т — 1)(п — 1) неизвестных. Данная система затем решается на ЭВМ одним из существующих методов.

В случае регулярной сетки параллельный алгоритм решения ДУ сводится к разбиению сетки на равные по площади участки и параллельному решению ДУ на каждом из участков методом конечных разностей. При этом аппроксимация решения будет приемлемой только в случае малого шага к. Кроме того, требуется существенный аппаратный ресурс ЭВМ для численного решения такой задачи. В «особых» зонах решения (в зонах ударных волн, межфазных границ) аппроксимация решения будет достаточно слабой. В таких зонах требуется сгущение сетки.

2. Адаптивная триангуляция сетки

При решении реальных физических задач регулярная сетка [3], состоящая из прямоугольных элементов, не может успешно применяться для аппроксимации решения ДУ, поскольку не учитывает неоднородность распределения частиц вещества в пространстве (особенно на приграничной области). Вместо разбиения на прямоугольные подобласти часто применяется метод адаптивной конечно-элементной триангуляции области решений [4]. Триангуляция области решений О в декартовой плоскости представляет собой разбиение на треугольные подобласти с узлами в О и ребрами, каждое из которых либо принадлежит обоим треугольникам, либо имеет оба узла на границе дО. Чем короче ребра, лежащие у нелинейной границы, тем лучше аппроксимация. Поэтому метод можно с успехом применять не только для приближения прямоугольных областей решений ДУ, но и нерегулярных областей произвольной формы. В частности, треугольники малого размера можно использовать для приближения «проблемной» приграничной зоны области решений, в то время как треугольники большего размера можно использовать для аппроксимации прочих подобластей.

3. Алгоритм адаптивного разбиения сеток

Далее приводится алгоритм [4] адаптивного разбиения сеток для некоторого порогового значения V (скажем, 0.01).

(1) Пусть Т — начальное «грубое» разбиение области решений, полученное путем применения алгоритма триангуляции (множество треугольных подобластей).

(2) Сконструировать матрицу жесткости [4, раздел 12.3] А для области Т и ввести соответствующую функцию ].

(3) Решить систему уравнений Ах = / для вектора х неизвестных

решений ДУ.

(4) Пусть Е — множество ребер из Т.

(5) Пройти по порядку по всем ребрам из Е; для каждого ребра е:

a. пусть г и 3 обозначают вершины ребра е;

b. если е расположено внутри области решений, то существует две треугольных подобласти, которые содержат это ребро; если вершины е расположены на границе области решений, то существует одна треугольная подобласть, содержащая е. Пусть Ь — треугольная подобласть с вершинами г, 3 и к, которая содержит е:

t = А(г,з,к);

c. если \х.1 — Ху | > V, то разбить £ на две треугольных подобласти: = А(г, (г + 3)/2,к), 12 = АЦ, (г + 3)/2,к) где (г + з)/2 есть средняя точка ребра е. Заменить £ во множестве Т на треугольные подобласти и

Так же разбить на две части другую треугольную подобласть, содержащую е, если она существует.

(6) Если разрешение сетки недостаточно высокое, то вернуться на

шаг 2.

(7) х будет численным решением исходной задачи.

4. Динамическое распараллеливание алгоритма адаптивного разбиения сеток

Система ОрепТЯ [1], современная реализация концепции «Т-системы», позволяет с успехом решать вычислительные задачи, в которых заранее неизвестен объем вычислений. Алгоритм адаптивного разбиения сеток представляет собой именно такую задачу. На каждом этапе разбиения в зависимости от результатов решения системы уравнений (системы жесткости) принимается решение о разбиении данной треугольной подобласти. Этот факт диктует необходимость использования систем динамического (а не статического) распараллеливания.

Предлагаемый параллельный динамический алгоритм основан на модифицированном последовательном алгоритме, в котором матрица жесткости составляется не для всей области «грубого» разбиения, а для одной конкретной «неразбитой» (необработанной) треугольной

подобласти. Т-функция (гранула параллелизма) вызывается для каждой такой подобласти, решает для нее систему жесткости и рекурсивно вызывается для каждого из двух треугольников, если разбиение было произведено (выполнен критерий — Xj | > v). Если на этой подобласти нет сильной вариации решений ДУ (критерий не выполнен), то Т-функция завершается. Все Xi, полученные при решении системы уравнений, будут решением задачи численной аппроксимации ДУ.

5. Заключение

В работе кратко описан параллельный алгоритм адаптивного разбиения расчетной сетки с использованием концепции динамического распараллеливания «Т-система».

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

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

[1] С. М. Абрамов, А. А. Кузнецов, В. А. Роганов. «Кроссплатформенная версия T-системы с открытой архитектурой», Вычислительные методы и программирование, 8:1(2) (2007), с. 175-180, URL http://num-meth.srcc.msu.su/ t 54, 57.

[2] Н. Г. Бураго. Вычислительная механика, М., 2012, 274 с. t 54.

[3] Д. Мак-Кракен, У. Дорн, Численные методы и программирование на ФОРТРАНе, пер. с англ. Б. Н. Казака, 2-е изд., ред. Б. М. Наймарк, Мир, М., 1977, 584 с. t 56.

[4] Ya. Shapira. Solving PDEs in C++: numerical methods in a unified objectoriented approach, SIAM, 2006 t 56.

Рекомендовал к публикации

д.ф.-м.н. С.В.Знаменский

Об авторах:

Антон Александрович Кузнецов

Научный сотрудник ИПС им. А.К. Айламазяна РАН, разработчик системного и прикладного ПО. Один из разработчиков системы параллельного программирования «OpenTS». Область научных интересов: параллельное программирование, компиляторы, распределенные вычисления в гетерогенных средах, геоинформационные системы. e-mail: tonic@pereslavl.ru

Владимир Александрович Роганов

Научный сотрудник ИПС им. А.К. Айламазяна РАН. Разработчик современных версий Т-системы, ведущий разработчик системы «OpenTS». Принимал активное участие в суперкомпьютерных проектах Союзного государства России и Беларуси, в том числе в проектах «СКИФ» и «СКИФ-ГРИД».

e-mail: var@pereslavl.ru

Герман Анатольевич Матвеев

Ведущий инженер-исследователь ИЦМС ИПС им. А.К. Айламазяна РАН. Один из разработчиков системы «OpenTS». Принимал участие в суперкомпьютерных проектах Союзного государства России и Беларуси.

e-mail: gera@prime.botik.ru

Валерий Иванович Осипов К.ф.-м.н., научный сотрудник ИПС им. А.К. Айламазяна РАН. Один из разработчиков системы «OpenTS». Принимал участие в суперкомпьютерных проектах Союзного государства России и Беларуси.

e-mail: val@pereslavl.ru

Пример ссылки на эту публикацию:

А. А. Кузнецов, В. А. Роганов, Г. А. Матвеев, В. И. Осипов. «Алгоритм динамического распараллеливания решения задачи адаптивного разбиения расчетной сетки для численного решения дифференциальных уравнений», Программные системы: теория и приложения, 2015, 6:3(26), с. 53-60.

http://psta.psn.ras .ru/read/psta2015_3_53-60.pdf

Anton Kuznetsov, Vladimir Roganov, German Matveev, Valeriy Osipov. Dynamic parallel algorithm of mesh adaptive refinement for numerical solution of differential equations.

Abstract. The article describes the dynamic methods to refining meshes used to obtain numerical solution of differential equations. A new parallel algorithm of mesh refining using the "T-system" dynamic approach is described. (In Russian).

Key Words and Phrases: T-system, dynamic parallelism, OpenTS, T++ programming language, differential equations, mesh, triangulation.

References

[1] S. M. Abramov, A. A. Kuznetsov, V.A. Roganov. "Cross-platform version of the T-system with an open architecture", Numerical Methods and Programming, 8:1(2) (2007), pp. 175-180 (in Russian), URL http://num-meth.srcc.msu.su/.

[2] N. G. Burago. Computational Mechanics, M., 2012 (in Russian), 274 p.

[3] D. McCracken, W. Dorn. Numerical methods and FORTRAN programming, with applications in engineering and science, Wiley, New York, 1964.

[4] Ya. Shapira. ¡Solving PDEs in C++: numerical methods in a unified objectoriented approach, SIAM, 2006.

Sample citation of this publication:

Anton Kuznetsov, Vladimir Roganov, German Matveev, Valeriy Osipov. "Dynamic parallel algorithm of mesh adaptive refinement for numerical solution of differential equations", Program systems: theory and applications, 2015, 6:3(26), pp. 53-60. (In Russian.)

URL http://psta.psiras.ru/read/psta2015_3_53-60.pdf

© A. A. Kuznetsov, V. A. Roganov, G. A. Matveev, V. I. Osipov, 2015 © Ailamazyan Program System Institute of RAS, 2015 © Program systems: Theory and Applications, 2015