Параллельная реализация полулагранжевого метода для уравнения неразрывности Текст научной статьи по специальности «Математика»

мана - М.: Сов. радио, 1973.

7. Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов / Пер. с англ.; Под ред. Ю.И. Александрова - М.: Мир, 1978.

Seasonal dynamics of Shira lake currents by field observations in 2014-2015 years

Lidiya Kompaniets, PhD, Associate Professor, Senior Researcher, ICM SB RAS Nikolay Golenko, PhD, leading research worker, AB FIO RAS. P. Shirshov Institute Olga Volodko, "Immanuel Kant Baltic Federal University" Olga Lipina, "Siberian Federal University"

In the article the variability of wind currents on the lake Shira was analyzed in period 2014-2015 years based on the data offield observations. The spectral analysis of internal waves, caused by variable wind, was carried out.

Keywords: dynamics of wind current, spectral analysis.

УДК 519.642.2, 519.683

ПАРАЛЛЕЛЬНАЯ РЕАЛИЗАЦИЯ ПОЛУЛАГРАНЖЕВОГО МЕТОДА ДЛЯ УРАВНЕНИЯ НЕРАЗРЫВНОСТИ

Александр Владимирович Вяткин, к.ф.-м.н., научный сотрудник Тел.: 8 391 249 5370, e-mail: vyatkin@icm.krasn.ru Институт вычислительного моделирования СО РАН http://icm.krasn.ru Елена Владимировна Кучунова, к.ф.-м.н., доцент E-mail: HKuchunova@sfu-kras.ru Институт математики и фундаментальной информатики СФУ

http://www.sfu-kras.ru

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

Ключевые слова: уравнение неразрывности, полулагранжевый метод, численное решение, параллельное программирование, OpenMP.

Исследование выполнено федеральным государственным бюджетным учреждением науки "Институт вычислительного моделирования Сибирского отделения Российской академии наук" (ИВМ СО РАН) при поддержке Российского фонда фундаментальных исследований (проект № 14-01-00296).

Введение

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

С 1960-х годов начал активно развиваться другой подход [5; 6], вытекающий из метода характеристик. Его суть состоит в том, что вдоль характеристики уравнение неразрывности можно переписать в виде обыкновенного дифференциального уравнения.

Это позволяет для ряда задач использовать большие шаги по времени и сократить время расчетов. В настоящее время рассмотренный подход называется полулагранжевым методом. Современные версии метода [7-9] основаны на балансовом интегральном соотношении при переходе с одного временного слоя на следующий слой. При этом аппроксимация численного решения на каждом слое по времени раскладывается на три составляющих: аппроксимация интеграла на верхнем слое по времени, на котором решение еще не известно; построение характеристик (траекторий) с верхнего временного слоя на нижний слой; приближенное вычисление интеграла на нижнем слое по времени. Для некоторых из этих алгоритмов получены результаты, позволяющие учитывать краевые условия Дирихле, обосновывать сходимость метода и выполнение дискретного аналога балансового соотношения [9].

Основной недостаток полулагранжевого подхода состоит в больших вычислительных затратах при вычислении интеграла на нижнем слое по времени. Существует большое количество различных версий метода, в которых этот интеграл считается разными способами. При этом, как правило [10; 11], преследуется две цели: сократить вычислительные затраты и выполнить балансовое соотношение. В данной работе мы предлагаем подход, который основан на билинейном преобразовании единичного квадрата на область интегрирования на нижнем слое по времени и использовании метода средних прямоугольников для вычисления интеграла по квадрату. Кроме того, в работе исследуется параллельная реализация описанного метода.

1. Постановка задачи

Пусть О = [0,1] х [0,1] - единичный квадрат с границей дО. На множестве [о, Т]х О рассмотрим двумерное уравнение неразрывности:

ф + д(иР) + д(ур)

Е.В. Кучунова

=/(/,х,у).

(1)

д/ дх ду

Здесь и(/,х,у), ,х,у) и /,х,у) - достаточно гладкие функции, известные в [0, Т]х О. Пусть при х = 0 и х = 1 функция и(/, х, у) не отрицательна

и(,0, у) > 0, и(/,1, у) > 0 V/ е [0, Т], Уу е [0,1], (2)

а при у = 0 и у = 1 выполнены условия прилипания

и(, х,0) = х,0) = и((, х,1) = х,1) = 0 V/ е[0,Т] Ух е [0,1]. (3) Будем считать, что искомая функция р(,х,у) известна в начальный момент времени в

О:

Д0, х у) = д (х у) у(х у) е 0,

и в любой момент времени / при х=0

Д(,0, у) = (/, у) V/ е [0, Т] Уу е [0,1].

(4)

(5)

2. Численная схема

2.1 Базовое соотношение для построения схемы. Построим в О квадратную сетку О^ с шагом к = 1N, N > 1, Ок = {(х., у}.): х. = ¡к, у}. = 'к, ¡,' = 0,..., N}. Исполь-

МЕТОДОЛОГИЧЕСКИЕ ИССЛЕДОВАНИЯ

зуем на отрезке [о, T ] равномерную сетку Тт с шагом т = T/K , K > 1, и узлами tk = кт,

к = 0,.. . Решение задачи (1)-(5) будем искать в виде дискретной функции рк, определенной на сетке ТтхОк. Для удобства введем обозначение рк(/к,х,,уу)= рру . Согласно условию (4) искомая функция р известна во всех узлах сетки О^ в начальный момент времени / = 0 :

РР0 =Рш(х,,У у) V иу = 0,...,N . Для построения схемы будем считать, что при / = ^-1 нам известно численное решение рк во всех узлах сетки О^ и требуется вычислить значение функции ръ на к-том слое по времени.

Полулагранжевый подход к построению численного решения задачи (1)-(5) основан на интегральном равенстве [9]

у у +1 (я^), (6)

о',,у

где

| Ртри! ((, У X ,0, У )Я, ^у П (х = 0)^0

Il|j ) =

Rk j

0 , иначе

Здесь в, у - некоторая замкнутая окрестность узла (х,, у у ) на к -том слое по времени, а

в' у - криволинейный четырехугольник на (к -1) -том слое по времени. Замкнутая область в, у может быть выбрана произвольным образом, а стороны четырехугольника

в' у определяются характеристическими траекториями, выпущенными из границы замкнутой области в, у (рис. 1.а). Если в, у расположена достаточно близко к границе втекания Гшри =( ]х Д)П (х = 0) , то некоторые траектории достигают границы Г1(.. В этом случае

Рис. 1. Основная идея полулагранжевого подхода возникает интеграл по множеству В^у ^ Г^При1, смотрите (рис. 1.Ь). В общем случае,

' к

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

' к полагать, что в, у является четырехугольным множеством, а Я, у - пустое множество,

как это показано на рис. 1.а. Для остальных вариантов рассуждения о построении алгоритма являются аналогичными.

2.2 Интеграл на верхнем слое по времени. В качестве в, у рассмотрим окрестность узла (х,, у у) в следующем виде:

Q, = \ -h/2,x, + h/2]x\yj -h/2,y, + h/2JJПQ.

Для аппроксимации интеграла по Q, j в тождестве (6) мы искомую функцию p(tk, х, y) заменим на билинейную интерполяцию p'h{tk, х, y) сеточной функции рh. Интегрирование по Q, j от функции , х, y) проведем аналитически и получим

к 2 .+1 '+1

ЯР(/к,х,у' ~ ЯР1(/к,х,у' = "67 X X ар,дРр9 . (7)

й,' й,' Р=-19 = '-1

Коэффициенты ард, р = .-1,...,. +1; д = ' -1,...,' +1, одинаковые для всех внутренних узлов (х.., уу), ¡,' = 1,.,N -1, и имеют следующий вид: а. ' = 36,

а/-1,' = а/,' -1 = а/,'+1 = а.+1,' = 6, а.-1,'-1 = а.-1,'+1 = а.+1,'-1 = а.+1,'+1 = 1. Для узлов (х.., у'), лежащих на границе множества О, коэффициенты а рд выводятся аналогично.

2.3 Интеграл на нижнем слое по времени. Обозначим через Лл,'г = (л^,г, ),

w, г = 0,1 вершины 2.' . Каждая точка Ам''г спускается по характеристической траектории на предыдущий слой по времени и образует точку Вм''г. Чтобы вычислить координаты точки Вц,'г, w, г = 0,1, необходимо решить задачу Коши на отрезке / е [/к-1,/к] обратно во времени для системы обыкновенных дифференциальных уравнений

= и(), = V (/) (8)

ш т

с начальными данными

) = Ах^,г, Ж ) = А;,г. (9)

Численное решение системы (8)-(9) найдем методом Эйлера [4], а его значение на (( - 1)-том слое по времени обозначим в виде Вк,^г = (вкм,'г,в^,г), w, г = 0,1. Криво-

I

линейный четырехугольник й.,' с вершинами в четырех точках Bw,г мы аппроксимируем прямолинейным четырехугольником <1' с вершинами в точках

-пк^,г I -пк^,г г>к^,г\ тт

в = \Вх , Ву у Для аппроксимации интеграла по 2.' мы заменим подынтегральную функцию Д/к-1, х, у) на дА(/к-1, х, у):

|р (/к-Ъ х у)2'г,' ~ |р (/к-Ь x, у)Шй\,' ~ к(/к-Ъ x, у)Шй\,' . (10)

2\,} ё'и 2'и

Для вычисления интеграла по замкнутой области 2. ' от функции р/А(/к-1, х, у) построим билинейное преобразование О = (Ох, Оу) квадрата Е = [0,1]х [0,1] на 2.,' (рис. 2).

Рис. 2. Преобразование области интегрирования

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

g fe о)=; ; cpw (feWr {в)в^, Gy fe о)=; ; ^ (fe^ (^)ByA,w,r (11)

где

x

w=0 r=0 w=0 r=0

() (1 ^ = 0 ( (1 -0, г = 0

V™ (^) = 1 д т' Уг (0) = 1 0 ,.

[ д , ™ = 1 [ 0 , г = 1

Таким образом, мы получаем

|р? (к-1, х, у) у = }}р? ((к-1, Ох (д, е), Оу (д, е)) (о (д, е))^е, (12)

й', у 00

где якобиан преобразования О имеет вид (\е\((}(д,0) = д—х -—— - д—х -——. Посколь-

Р Р 4 4 " дд д0 д0 дд

ку преобразование О является билинейным, то мы можем вычислить det (О(д,0)) аналитически. Обозначим через у1 приближенное значение повторного интеграла в правой части соотношения (12). Чтобы вычислить I, у \ мы последовательно применим метод средних прямоугольников, предварительно разбив отрезок интегрирования [0,1] на ^ос равных частей, и получим

^оо-1N1оо -1

Ij =Чг ; £ ph((,Gx(s,Od)Gy(s,Od))((s,Od)), (13)

4,j 2 "k-1>wxxbs^d fi^y^

Nloc s=0 d=0

где s = hloc (s + 0 5), Od = hloc(d + 0 5), hloc = 1 Nloc .

2.3 Система уравнений. Мы аппроксимировали все слагаемые базового соотношения (6). Таким образом, с учетом (7), (10)-(13) из соотношения (6) следует, что

i+1 j+1 h k 64

; ;ap,qphpkq = 64i',j, i=1,-,N; j=0,-,N. (14)

p=-1 q=j-1 h

Система линейных алгебраических уравнений (14) состоит N (n +1) уравнений и имеет N (N +1) неизвестных величин. Полагаем, что система имеет единственное решение. Для поиска решения системы (14) используем итерационный метод Якоби. В качестве начального приближения рассмотрим функцию рh на предыдущем слое по времени.

3. Параллельная реализация алгоритма и результаты численных экспериментов

Основные вычислительные затраты ложатся на вычисление интегралов if- \

i = 1, —, N; j = 0,-, N, и поиск решения системы (14). Поскольку интегралы Ikj1 не

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

#pragma omp parallel for collapse (2).

Для вычисления максимального значения вектора невязки на каждой итерации метода Якоби используется опция reduction(max: norma).

При тестировании метода был осуществлен ряд вычислительных экспериментов. В качестве одного тестов была рассмотрена задача параллельного переноса круга. Начальные данные и значения на границе втекания задавались в виде

р^м-к(х - °.2)2+(ач-°-2>2 - ^ у)-о.

[0, иначе,

Функции скорости были выбраны постоянными и(, х, у)- , х, у)- 0.6 . В этом случае точное решение задачи (1)-(5) имеет следующий вид:

р(, х, у)-]1, (х - Хо )2 +(У - Уо )2 * 0,12, [0, иначе,

где Хо(()-0.2 + и, уо(() = 0 2 + V. За время Т - 1 центр круга перемещается из точки (0.2, 0.2) в точку (0.8, 0.8). Сходимость исследовалась в аналоге нормы пространства Ь1

N

Р

= I mes' i, j=0

(j ),

где теБ^}) - площадь }. При проведении указанного эксперимента мы зафиксировали отношение т к к в виде т - к/2. Для исследования сходимости метода проведены расчеты на последовательности сеток О к с различным параметром N - 40*2п, п - 0,...,5. Обозначим через рр численное решение, полученное на сетке Ок при N - 40* 2п. На рис. 3 представлен график, отражающий отношение

h / h

Р~Рп Ч Р~Рп+1

Эксперименты подтверждают первый порядок сходимости численного решения к точному решению. Расчеты производились в ИВМ СО РАН на вычислительной системе МВС1000/16, состоящей из 28 четырехъядерных процессоров Intel Xeon Core E5345 2.33GHz. Производительность системы по LinPack более 430 Гфлопс, пиковая 1043,84 млрд. оп/сек. Время выполнения последовательной и OpenMP-версий программы приведены в таблице.

Рис. 3. Порядок сходимости схемы

Таблица.

Время выполнения различных версий программы (секунды)

N Последовательная программа Параллельная программа

2 нити 4 нити 6 нитей 8 нитей

160 112 56 29 19 15

320 878 445 225 151 114

640 7149 3580 1799 1200 904

1280 57536 28871 14414 9718 7313

ОрепМР-версия программы запускалась на количестве нитей: 2, 4, 6 и 8. Было вычислено ускорение Б р - Т^Тр , получаемое при использовании параллельного алгоритма

по сравнению с последовательным вариантом и эффективность параллельного алгоритма Ер - Бр /р х 100%, где Т - время работа последовательной программы, Тр - время работы

параллельной программы с использованием р нитей в параллельных секциях. На рис. 4 приведены данные по ускорению (слева) и эффективности (справа) ОрепМР-программы в

зависимости от количества нитей. На больших сетках ОрепМР-программа на 8 нитях дает ускорение порядка 8 раз, т.е. эффективность распараллеливания близка к 100%.

8 ■

0

1

Рис. 4. Ускорение (слева) и эффективность (справа) параллельной программы Заключение

Авторы считают, что в данной работе новым является следующий результат: предложен параллельный алгоритм численного решения двумерного уравнения неразрывности. Исследование алгоритма на тестовых задачах продемонстрировало первый порядок сходимости. Эффективность распараллеливания для систем с общей памятью составляет 95-99 %.

Литература

1. Рябенький В.С. Об устойчивости разностных уравнений / В.С. Рябенький, А.Ф. Филиппов. - М.: Гостехиздат, 1956.

2. Вазов В. Разностные методы решения уравнений в частных производных / В. Вазов, Дж. Форсайт. - М.: Наука, 1963.

3. Годунов С.К. Разностные схемы (введение в теорию) / С.К. Годунов, В.С. Рябенький. -М.: Наука, 1977.

4. Бахвалов Н.С. Численные методы / Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. - 3-е изд., перераб. и доп. - М.: БИНОМ. Лаборатория знаний, 2003.

5. Магомедов К.М. Метод характеристик для численного решения пространственных течений газа // Журнал вычислительной математики и математической физики. 1966. Т. 6. № 2. С. 313-325.

6. Wiin-Nielson On the application of trajectory methods in numerical forecasting // Tellus. -1959.Vol. 11. P. 180-186.

7. Iske A. Conservative semi-Lagrangian advection on adaptive unstructured meshes // Numer. Meth. Part. Diff. Eq. 2004. Vol. 20. P.388-411.

8. Xin Wen, Vyatkin A.V., Shaidurov V.V. Semi-Lagrangian Scheme for soling hyperbolic equation of first order // Молодой учёный. 2013. № 9(56). С.6-13.

9. Вяткин А.В., Ефремов А.А., Карепова Е.Д., Шайдуров В.В. Использование гибридных вычислительных систем для решения уравнения переноса модифицированным методом траекторий // Пятая Международная конференция «Системный анализ и информационные технологии»: Труды конференции. В 2 т. Красноярск: ИВМ СО РАН. 2013. Т. 1. С.45-55.

10.Phillips T.N., Williams A. J. Conservative semi-Lagrangian finite volume schemes // Numer. Meth. Part. Diff. Eq. 2001. Vol. 17. P.403-425.

11.T. Arbogast and W.H. Wang Convergence of a Fully Conservative Volume Corrected Characteristic Method for Transport Problems // SIAM J. Numer. Anal., 2010 Vol. 48(3). P. 797-823

A parallel semi-Lagragian algorithm for advection equation

Alexander Vyatkin, researcher, Institute of Computational Modellng of SB RAS

Helen Kuchunova, associate Professor, School of Mathematics and Computer Science

We construct the algorithm of the family of semi-Lagrangian methods for an advection problem. The algorithm is based on the integral balance equation for the neighborhood of a grid node. This equation involves integrals over two neighboring time levels. We study the effectiveness of the parallel algorithm with OpenMP technology.

Keyword: semi-Lagramgian method, integral balance equation, numerical solution, OpenMP.