Все большая доступность супервычислителей в настоящее время повышает интерес к методам решения краевых задач, допускающим распараллеливание счета на уровне алгоритма. В связи с этим разработка подобных алгоритмов и их программная реализация на многопроцессорных комплексах находят все большее распространение и применение на практике. В данной работе представлены алгоритм и комплекс программ для решения двухмерных задач теории потенциала. В основу алгоритма положена методика применения метода граничных элементов [1], основанная на точном вычислении интегралов по граничным элементам с помощью полученных авторами аналитических формул и параллельных вычислений на каждом этапе решения. Разработанный подход и созданные на его основе программы показали свою эффективность для решения линейных задач теории упругости [1], а также нестационарных задач диффузии и теплопроводности [2].

 

Краевая задача потенциального течения

Основным дифференциальным уравнением задачи теории потенциала является уравнение Пуассона, в двухмерном случае имеющее вид

Dq(x) = b(x), x(x1, x2) Î Ω,                                     (1)

где θ(x) – значение искомой функции в точке x; ∆=∂2/∂x12+∂2/∂x22 – оператор Лапласа; Ω – рассматриваемая область; b(x) – заданная функция источника, распределенная на области Ω.

Стационарное потенциальное уравнение (1) описывает многие физические процессы, такие как кручение призматического стержня, установившееся течение жидкости, а также стационарное тепловое поле, электростатическое поле и т.д. Задача заключается в определении функции θ(x) (например, электрического или гидравлического потенциала либо температуры), удовлетворяющей уравнению (1) и граничным условиям:

q(x0) = q*(x0), x0(x01, x02) Î S1;

q(x0) = q*(x0), x0(x01, x02) Î S2.                               (2)

Здесь q(x)=∂θ(x)/∂n(x) – поток (соответственно плотность электрического тока, скорость течения жидкости или поток тепла); n(x) – внешняя нормаль к границе S = S1È S2 области Ω; звездочкой отмечены известные величины.

Алгоритм решения

В соответствии с методом граничных элементов [3] осуществляется переход от основных уравнений краевой задачи (1) и (2) к граничному интегральному уравнению для граничной точки x0:

(3)

где u*(x0, x), f*(x0, x) – функции влияния. Для двухмерной задачи теории потенциала

;

;

;                                     (4)

Здесь и далее по повторяющемуся индексу производится суммирование от 1 до 2.

Для решения граничного интегрального уравнения необходима дискретизация границы S, то есть построение модели области Ω с границей, разбитой на граничные элементы. Пусть N + M – общее количество граничных элементов в пред- положении, что часть границы S1 разбита на N граничных элементов, а часть границы S2 – на M элементов. В данной работе используются прямолинейные граничные элементы и принимается, что функции θ(x) и q(x) имеют постоянные значения на каждом элементе

q(xi) = qi = const; q(xi) = qi = const,

i = 1, 2, …, N+M.                                                   (5)

Здесь xi – узел, находящийся в середине элемента [ai–1, ai].

Уравнение (3) в этом случае приводится к виду

            (6)

,  j=1, 2, …, N + M.

Вычисление интеграла по области (последнее слагаемое в правой части уравнения (6)) зависит от вида функции источника b(x). Рассмотрим два случая.

1. b(x) – гармоническая функция в области Ω, то есть имеет место равенство Db(x) = 0. Тогда согласно [3] интеграл по области в правой части (6) преобразуется в интеграл по границе области следующего вида:

(7a)

2. b(x) – произвольная числовая функция. В этом случае область Ω разбивается на конечные элементы ∆Ω1, ∆Ω2, …, ∆ΩK, в каждом из которых функция b(x) полагается постоянной. Тогда интеграл по области в правой части уравнения (6) вычисляется следующим образом:

(7b)

Здесь xk – узел элемента ∆Ωk.

С учетом (7a) и (7b) уравнение (6) принимает соответственно вид

             (8a)

 

или

            (8b)

.                           

Интегралы по граничным элементам от функций влияния вычисляются аналитически по формулам, полученным в [1] с помощью перехода от произвольного граничного элемента к базовому отрезку [0, L], лежащему на оси абсцисс:

;

            (9)

Соотношения (8), записанные для всех граничных узлов, образуют систему линейных алгебраических уравнений относительно узловых значений функций θ(x) и q(x), не заданных граничными условиями:

P*θ + R*q + b* = 0.                                                 (10)

Здесь P* и R* – квадратные матрицы размером (N + M)×(N + M), элементы которых  и  являются интегралами от функций влияния по граничным элементам [ai–1, ai]; θ – (N + M)-мерный вектор с элементами θi – значениями функции θ(x) в узлах, причем первые N значений заданы граничными условиями, а последующие M неизвестны; q – (N + M)-мерный вектор с элементами qi – значениями потока q(x) в узлах границы, при этом первые N значений неизвестны, а последующие M заданы граничными условиями; b* – вектор размерности (N+M), элементы которого   определяются численным интегрированием в соответствии с (8a) или (8b).

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

(11)

Программная реализация

Представленный алгоритм был реализован в виде программного модуля. Языки программирования выбирались с требованием максимальной переносимости комплекта программ под различные операционные системы. Для создания графического интерфейса был использован язык программирования Java. Java-модуль ru.uran.imach.po­tential состоит из одного пакета и включает в себя шесть базовых и несколько вспомогательных классов. Назначение базовых классов:

–      класс FrontFrame отвечает за главное окно программы, описывает вид и стиль программы, панель инструментов и вспомогательные диалоговые окна, отвечает за многооконный интерфейс в целом;

–      класс Expoint2D описывает узловую точку;

–      класс BorderModel отвечает за построение внешней и внутренней границ модели, cодержит методы добавления, перемещения, удаления узловых точек и метод перемещения модели, включает в себя методы сохранения объекта и загрузки объекта из файла;

–      класс MainPanel описывает графическое представление объекта, содержит методы отрисовки отдельных частей модели и работы с масштабом;

–      класс PieceSegment описывает конкретную границу модели и произвольный граничный элемент;

–      класс MGEModel представляет математическую модель объекта, таблицы с результатами расчетов, содержит возможность сохранения результатов в разных форматах на жестком диске.

Для создания программы использовалась среда разработки Eclipse. Модуль, написанный на язы- ке Java, был скомпилирован и экспортирован в jar-архив. Программа компилировалась и выполнялась в среде Oracle Java SE 7 [4]. На рисунке 1 представлен экран ввода данных графического интерфейса.

Расчетный модуль программы написан на языке программирования C++ с использованием библиотеки параллельных вычислений MPI [5, 6] и включает в себя непосредственно алгоритм расчета. C++-модуль mge.potential состоит из четырех базовых классов:

–      класс PieceSegment описывает граничные элементы объекта;

–      класс Triangle содержит методы разбиения внутренней области объекта на конечные элементы;

–      класс MGEModel описывает математическую модель объекта, содержит методы загрузки объекта из файла и методы разбиения границы объекта;

–      класс PotentialSolver решает задачу теории потенциала.

Программный код модуля на C++ после компиляции и проверки на машине разработчика копировался на кластер umt («Уран») ИММ УрО РАН, где компилировался вновь, уже в окружении библиотеки MPI-кластера. При запуске модуля задавалось количество вычислительных ядер на кластере, кратное числу точек разбиения границы области. Выходные данные были сформированы в виде таблицы Excel. При написании модуля использовались также библиотеки BOOST [7] и GSL [8]. Программный модуль компилировался компилятором GCC версии 4.6.2. На рисунке 2 показана общая структура программного модуля.

Примеры

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

Пример 1. Рассмотрим задачу о распространении тепла в пластине, имеющей форму эллипса с полуосями a = 5 м и d = 3 м. Внутри пластины действует источник тепла, заданный в виде гармонической функции b(x)=6x1+6x2+2. На границе задан тепловой поток

.

В результате решения задачи с помощью созданного пакета программ было определено распределение температуры на границе пластины и вычислены значения температуры во внутренних точках пластины. На рисунке 3 показаны результаты сравнения значений температуры вдоль оси абсцисс, рассчитанных для 20 и 100 граничных элементов, с точным решением θ(x)=x13+x23+x22.

Пример 2. Рассмотрим задачу о распространении тепла в квадратной пластине со стороной h=1 м. Пусть на вертикальных участках границы пластины задана температура =100 °С, =200 °С, а на горизонтальных – нуле- вой поток:   соответствующий отсутствию теплообмена с окружающей средой. Внутри пластины действует источник тепла b(x)=–103x12. Задача была решена для 4 и 20 граничных элементов. На рисунке 4 приводится сравнение изменения полученных значений температуры вдоль отрезка AB, концы которого имеют координаты (0.5,0) и (0.5,1), с точным решением θ(x)=100+183.3x1–83.3x14.

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

Литература

1.     Федотов В.П., Спевак Л.Ф. Модифицированный метод граничных элементов в задачах механики, теплопроводности и диффузии. Екатеринбург: Изд-во УрО РАН, 2009. 164 с.

2.     Федотов В.П., Нефедова О.А. Применение модифицированного метода граничных элементов для решения параболических задач // Вестник СамГТУ: Сер. Физ.-мат. науки. Самара, 2011. Т. 25. № 4. С. 93–101.

3.     Бреббия К., Теллес Ж., Вроубел Л. Метод граничных элементов. М.: Мир, 1987. 526 с.

4.     Oracle Technology Network for Java Developers. URL: http://www.oracle.com/technetwork/java/index.html (дата обращения: 11.12.2013).

5.     Pacheco P. An Introduction to Parallel Programming. Morgan Kaufmann, 2011.

6.     The Message Passing Interface (MPI) standard. URL: http://www.mcs.anl.gov/research/projects/mpi/ (дата обращения: 22.04.2013).

7.     Boost C++ Libraries. URL: http://www.boost.org (дата обращения: 09.12.2013).

8.     GSL-GNU Scientific Library. URL: http://www.gnu.org/ software/gsl (дата обращения: 26.11.2013).

References

1.  Fedotov V.P., Spevak L.F.  Modifitsirovanny metod granichnykh elementov v zadachakh mekhaniki,
teploprovodnosti i diffuzii  [A  Modified Boundary Element Method  in  Mechanics, Transcalency  and  Diffusion Prob-lems]. Ekaterinburg, Ural Branch of RAS Publ., 2009, 164 p.

2.  Fedotov V.P., Nefedova O.A.  Vestn. Samarskogo  gos. tekh.  univ.  [The  Bulletin of the Samara State Technical
Univ.]. 2011, vol. 25, no. 4, pp. 93–101 (in Russ.).

3.  Brebbia C.A., Telles J.C.F., Wrobel L.C.  Boundary  Element Techniques. Springer, Berlin, Heidelberg Publ.,
1984, 526 p.

4.  Oracle Technology Network for  Java Developers. Available at:  http://www.oracle.com/technetwork/java/in-dex.html (accessed December 11, 2013).

5.  Pacheco P. An Introduction to Parallel Programming. Morgan Kaufmann Publ., 2011.

6.  The Message Passing Interface (MPI) standard. Available at: http://www.mcs.anl.gov/research/projects/mpi/
(accessed April 22, 2013).

7.  Boost C++ Libraries. Available at: http://www.boost.org (accessed December 9, 2013).

8.  GSL-GNU Scientific Library. Available at: http://www.gnu.org/software/gsl (accessed November 26, 2013).