=Paper=
{{Paper
|id=Vol-1662/mod10
|storemode=property
|title=Параллельный вариант численного метода решения многомерного уравнения переноса с запаздыванием(Parallel variant of numerical algorithm for solving a multidimensional advection equation with time delay)
|pdfUrl=https://ceur-ws.org/Vol-1662/mod10.pdf
|volume=Vol-1662
|authors=Svyatoslav I. Solodushkin,Arsen A. Sagoyan,Irina F. Iumanova
}}
==Параллельный вариант численного метода решения многомерного уравнения переноса с запаздыванием(Parallel variant of numerical algorithm for solving a multidimensional advection equation with time delay)==
https://ceur-ws.org/Vol-1662/mod10.pdf
Параллельный вариант численного метода решения
многомерного уравнения переноса с запаздыванием
С.И. Солодушкин1,2 А.А. Сагоян1 И.Ф. Юманова1
solodushkin_s@mail.ru arsensag@mail.ru yuirfa@mail.ru
1 – УрФУ (Екатеринбург) 2 – ИММ УрО РАН (Екатеринбург)
Аннотация
Построено семейство разностных схем для численного решения
многомерного уравнения первого порядка с запаздывающим ар-
гументом. Изучены порядок аппроксимации, устойчивость и по-
рядок сходимости. Приведен алгоритм решения, допускающий па-
раллельную реализацию. На тестовых примерах показана хорошая
слабая и сильная масштабируемость параллельной версии. Экспе-
риментальные оценки порядков сходимости соответствуют теоре-
тическим оценкам.
1 Введение и формальная постановка
Уравнения в частных производных первого порядка (известные также как уравнения переноса) с запаз-
дыванием возникают при моделировании в различных областях естествознания [1, 2, 3, 4]. Качественная
теория функционально-дифференциальных уравнений в частных производных (ФДУ в ЧП) развита до-
статочно полно (см., например, [4] и ссылки в ней). Существует обширная библиография, посвященная
вопросам существования, единственности, устойчивости и асимптотическому поведению решений ФДУ в
ЧП. В то же время аналитически ФДУ в ЧП решаются лишь в исключительных ситуациях, а потому
разработка, обоснование и программная реализация численных методов для этих уравнений вызывают
значительный интерес.
Метод прямых [5, 6, 7] сводит ФДУ в ЧП к системе ФДУ в обыкновенных производных, которая может
быть решена специальными методами, однако, системы уравнений, возникающие после дискретизации
по пространству, являются жесткими. Неявные разностные схемы для ФДУ в ЧП [8, 9, 10] позволяют
избежать проблемы жесткости, но для того, чтобы найти решения в каждом следующем временном слое,
необходимо решать систему нелинейных уравнений большой размерности.
Эффективные разностные схемы для ФДУ в ЧП, в том числе, уравнений переноса с запаздыванием
были разработаны в [11, 12]. Основная идея этих работ — принцип разделения конечномерной и бесконеч-
номерной (функциональной) составляющих в структуре ФДУ. Чтобы учесть запаздывание, используется
интерполяция и экстраполяция дискретной предыстории. Экстраполяция, в частности, позволяет авторам
строить неявные схемы и позволяет избежать необходимости решать нелинейные уравнения на каждом
временном слое. Подход, используемый в настоящей работе, близок к [11, 12] и основан на применении
методов проверки устойчивости двуслойных разностных схем [13] и принципе разделения конечномерной
и бесконечномерной составляющих.
Copyright c by the paper’s authors. Copying permitted for private and academic purposes.
In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the 47th International Youth School-conference “Modern Problems in
Mathematics and its Applications”, Yekaterinburg, Russia, 02-Feb-2016, published at http://ceur-ws.org
315
Увеличение производительности суперкомпьютеров приводит к повышению интереса к разработке па-
раллельных вариантов численных методов. Ниже мы приведем обзор некоторых подходов. Первый класс
методов основан на декомпозиции области решения, классическим является альтернирующий метод Швар-
ца для эллиптических уравнений, последовательный по своей природе. Чтобы распараллелить метод Швар-
ца, применяются разнообразные подходы: аддитивный метод Шварца, параллельные многосеточные ме-
тоды с предобуславливанием, параллельный метод Шварца с весами и т. д. [14]. Параллельные методы
декомпозиции области зарекомендовали себя как очень эффективные для решения эллиптических урав-
нений [15] и были обобщены и адаптированы для решения уравнений диффузии-конвекции [16, 17].
Параллельная реализация метода полу-Лагранжевого типа для уравнения переноса рассматривалась
в [18]. Разностная схема с переменным шаблоном построена на основе интегрального неравенства между
соседними временными слоями. Данный подход позволяет авторам избежать ограничений типа Куранта,
накладываемых на соотношение между шагами по времени и по пространству.
Настоящая работа обобщает результат [19, 20, 21] на многомерный по пространству случай и представ-
ляет параллельную версию соответствующего численного метода.
Рассмотрим многомерное уравнение переноса с запаздыванием:
p
∂u X ∂u
+a = f (x, t, u(x, t), ut (x, ·)), (1)
∂t α=1
∂xα
здесь x = (x1 , . . . , xp ), xα ∈ [0, Xα ], и t ∈ [t0 , θ] — пространственная и временная переменные; u(x, t) —
искомая функция; ut (x, ·) = {u(x, t + ξ), −τ 6 ξ < 0} — предыстория искомой функции к моменту t, τ > 0 –
величина запаздывания,
Qp коэффициент a > 0.
Пусть G = α=1 [0, Xα ] — p-мерный прямоугольник с границей Γ. Заданы начальные и граничные
условия
u(x, t) = ϕ(x, t), x ∈ G, t ∈ [t0 − τ, t0 ], (2)
u(x, t) = g(x, t), x ∈ Γ, t ∈ [t0 , θ]. (3)
Пусть f, ϕ, g таковы, что (1)–(3) имеет единственное решение, обладающее той степенью гладкости,
которая требуется для построения численного метода [4].
Обозначим через Q = Q[−τ, 0) множество функций u(ξ), кусочно-непрерывных на [−τ, 0) с конечным
числом точек разрыва первого рода, непрерывных справа в точках разрыва.
Введем в Q норму kukQ = supξ∈[−τ,0) |u(ξ)|. Для построения численного метода и обоснования его схо-
димости дополнительно предположим, что функционал f (x, t, u, v(·)) определен на G × [t0 , θ] × R × Q и
Липшицев по двум последним аргументам:
∃ Lf ∈ R ∀ x ∈ G, t ∈ [t0 ; θ], u1 ∈ R, u2 ∈ R, v 1 ∈ Q, v 2 ∈ Q :
|f (x, t, u1 , v 1 (·)) − f (x, t, u2 , v 2 (·))| 6 Lf |u1 − u2 | + kv 1 (·) − v 2 (·)kQ .
2 Построение разностной схемы
Рассмотрим разбиение отрезка [0, Xα ] на части с шагом hα = Xα /Nα , где Nα — число отрезков разбиения
по соответствующей координате. На G введем сетку ω h = {(i1 h1 , . . . , ip hp ), iα = 0, . . . , Nα }.
Мы используем обозначение x(i) = (xi1 , . . . , xip ) = (i1 h1 , . . . , ip hp ), iα = 0, 1, . . . , Nα , здесь индекс i — это
p-мерный вектор. Обратим внимание, что xα ∈ R — это α-я координата вектора x, в то время как x(i) ∈ Rp
— узел в сетке ω h .
Введем ω ∆ = {tj = t0 + ∆, j = 0, . . . , M }, где ∆ = (θ − t0 )/M. Для простоты изложения предполагаем,
что τ /∆ = m целое.
Наконец, определим сетку ω h∆ = ω h × ω ∆ = {(x(i) , tj ), x(i) ∈ ω h , tj ∈ ω ∆ }. В узлах сетки ω h∆ будем
искать приближенное решение задачи (1)–(3).
Обозначим через uij приближение функции u в узле: uij ≈ u(x(i) , tj ), где (x(i) , tj ) ∈ ω h∆ .
Функционал f (x(i) , tj , u(x(i) , tj ), utj (x(i) , ·)) может зависеть от значений u между узами, а потому для его
вычисления необходима интерполяция сеточной функции uij . Для каждого момента времени tj и величины
запаздывания ξ ∈ [−τ, 0) либо tj +ξ 6 t0 , в этом случае полагаем u(x(i) , tj +ξ) = ϕ(x(i) , tj +ξ), либо tj +ξ > t0 ,
в этом случае используем интерполяцию.
316
Для каждого узла (x(i) , tj ) определим его дискретную предысторию:
i
{ul }j = uil | max{0, j − m} 6 l 6 j .
Определение 1 Отображение I, определенное на множестве Aij всех возможных дискретных предыс-
торий и действующее по правилу
i
I : Aij → Q[− min{τ, tj }, 0] : {ul }j 7→ v i,j (·) = v i,j (tj + ξ),
называется оператором интерполяции дискретной предыстории.
Для построения численного метода, обладающего нужными свойствами, достаточно применять кусочно-
i
линейную интерполяцию. Именно, для дискретной предыстории {ul }j определим
1
v i,j (tj + ξ) = (tl+1 − tj − ξ) uil + (tj + ξ − tl ) uil+1 , tl 6 tj + ξ 6 tl+1 . (4)
∆
Оператор кусочно-линейной интерполяции имеет 2-й порядок на точном решении, т. е. существуют кон-
станты C1 и C2 такие, что для всех i = (i1 , . . . , ip ), iα = 0, 1, . . . , Nα , j = 1, . . . , M, и t ∈ [max{0, tj − τ }, tj ]
выполняется:
v i,j (t) − u(x(i) , t) 6 C1 max uil − u(x(i) , tl ) + C2 ∆2 .
max{0,j−m}6l6j
Рассмотрим разностные операторы, которые аппроксимируют первые производные по пространствен-
ным переменным,
!
p
i[−1α ] 2h i[−1α ] ∂g(x(i)[−1α ] ,tj ) ∂g(x(i)[−1α ] ,tj )
+ 4uij
P
−4uj − a − −
α
fj ∂xβ ∂t
β=1,β6 =α
, для iα = 1,
Ωα uij = 2hα
i[−2α ] i[−1 ]
u − 4uj α + 3uij
j
, для iα > 2.
2hα
Здесь x(i)[−1α ] — сосед узла x(i) , сдвинутый влево на один шаг hα по соответствующей координате, т. е.
i[−1 ]
узел (i1 h1 , . . . , iα−1 hα−1 , (iα − 1)hα , iα+1 hα+1 , . . . , ip hp ), uj α — сеточное приближение функции u(x, t),
вычисленное в узле (x(i)[−1α ] , tj ).
Рассмотрим семейство разностных схем (0 6 s 6 1) для j = 0, . . . , M − 1 и i = (i1 , . . . , ip ), iα = 1, . . . , Nα :
p
uij+1 − uij X
s Ωα uij+1 + (1 − s) Ωα uij = fji
+a (5)
∆ α=1
с начальными и граничными условиями
ui0 = ϕ(x(i) , t0 ), v i,0 (t) = ϕ(x(i) , t), t < t0 ,
uij = g(x(i) , tj ), x(i) ∈ Γ.
Здесь fji = f x(i) , tj , uij , v i,j (·) — значение функционала f, вычисленное на приближенном решении,
v i,j (·) — результат интерполяции. Для построения численного метода мы дополнительно предполагаем,
что g(x, t) дифференцируема.
3 Исследование сходимости
Определение 2 Обозначим εij = u(xi , tj )−uij , i = (i1 , ..., ip ), iα = 1, ..., Nα , j = 0, . . . , M. Будем говорить,
что метод (5) сходится, если εij → 0 при h → 0 и ∆ → 0 для всех i и j. Будем говорить, что метод
p
p
hqα1 + ∆q2 , если существует константа C такая, что kεij k 6 C
P q1
hα + ∆q2
P
сходится с пордяком
α=1 α=1
для всех i и j.
317
Определение 3 Невязкой метода (5) назовем сеточную функцию
p
u(x(i) , tj+1 ) − u(x(i) , tj ) X
Ψij = +a (s Ωα u(x(i), tj+1 ) + (1 − s) Ωα u(x(i), tj )) − f¯ji , (6)
∆ α=1
где f¯ji = f (x(i) , tj , u(x(i) , tj ), utj (x(i) , ·)) – значение функционала f , вычисленное на точном решении.
Теорема 1 Пусть точное решение u(x, t) задачи (1) − (3) трижды непрерывно дифференцируемо по xα ,
дважды непрерывно дифференцируемо по времени t и все первые производные решения по xα непрерывно
p
h2α + ∆.
P
дифференцируемы по t, тогда невязка метода (5) имеет порядок
α=1
Доказательство теоремы сводится к разложению функции u в ряд Тейлора в окрестности точки (x(i) , tj );
для одномерного случая см. [20].
В силу нелинейного характера зависимости функционала f от решения и его предыстории, мы
применим аппарат общих схем с последействием, разработанный в [22] для случая функционально-
дифференциальных уравнений с обыкновенными производными и примененый в дальнейшем для уравне-
ний в частных производных [11, 12].
Без ограничения общности мы рассмотрим случай однородных граничных условий u(x, t) = 0, x ∈ Γ,
t ∈ [t0 , θ]. Действительно, замена вида ũ(x, t) = u(x, t) − ĝ(x, t), где ĝ(x, t) — произвольная функция,
обладающая нужной степенью гладкости, совпадающая с g(x, t) при x ∈ Γ, t ∈ [t0 , θ], сводит изначально
поставленную задачу к задаче с однородными граничными условиями.
Доказательство сходимости будем проводить методом вложения семейства (5) в общую разностную
схему с последействем [22]. Идея близка к [19, 20] и основана на методе повышения размерности.
Пусть p фиксировано. Для каждого tj ∈ ω ∆ определим значения дискретной предыстории вектором ~yj ,
для этого упорядочим компоненты вектора uij в лексикографическом порядке:
>
(1,...,1) (2,...,1) (N ,...,1) (N ,...,Np )
~yj = uj , uj , . . . , uj 1 , . . . , uj 1 ∈ Yp ,
здесь j = 0, . . . , M, знак ·> означает транспонирование, Yp — линейное пространство, dim Yp = N1 × N2 ×
>
× . . . × Np . Например, при p = 1, имеем ~yj = u1j , u2j , . . . , uN
j
1
.
p
P
Построим разностный оператор A : Yp → Yp , который соответствует Ωα . Для этого рассмотрим
α=1
последовательность матриц Dα , α = 1, . . . , p, где каждая следующая (начиная со второй) рекуррентно
определяется через предыдущую.
Рассмотрим N1 × N1 матрицу D1 , которая соответствует разностному оператору Ω1 :
4 0 0 0 ··· ··· ··· 0
−4 3 0 0 ··· ··· ··· 0
1 −4 3 0 ··· ··· ··· 0
a 0 1 −4 3 ··· ··· ··· 0
D1 = .. .
.. .. .. ..
2h1 . . . . .
0 ··· ··· 1 −4 3 0 0
0 ··· ··· 0 1 −4 3 0
0 ··· ··· 0 0 1 −4 3
Определим матрицу E1 = 2ha2 I1 , где I1 — единичная N1 × N1 -матрица. Далее, определим матрицу D2 ,
которая имеет блочно-трехдиагональную форму, каждый блок имеет размер N1 × N1 , и в каждой строке
N2 блоков:
318
D1 + 4E1 0 0 0 ··· 0 0 0
−4E1 D1 + 3E1 0 0 ··· 0 0 0
E1 −4E1 D1 + 3E1 0 ··· 0 0 0
D2 =
0 E1 −4E1 D1 + 3E1 ··· 0 0 0 .
.. .. .. .. ..
. . . . .
0 0 0 0 ··· −4E1 D1 + 3E1 0
0 0 0 0 ··· E1 −4E1 D1 + 3E1
Далее аналогично. Таким образом определим линейный оператор A через его квадратную матрицу размера
N1 × N2 × · · · × Np :
Dp−1 + 4Ep−1 0 0 ··· 0 0 0
−4Ep−1 Dp−1 + 3Ep−1 0 ··· 0 0 0
E p−1 −4E p−1 D p−1 + 3E p−1 · · · 0 0 0
A= ,
.. . . . . .
.
. . . .
0 0 0 · · · −4Ep−1 Dp−1 + 3Ep−1 0
0 0 0 ··· Ep−1 −4Ep−1 Dp−1 + 3Ep−1
где Ep−1 = 2hap Ip−1 , а Ip−1 — единичная матрица размерности N1 × · · · × Np−1 .
Лемма. 1 Для всякого p, матрица Dp положительно определена.
Доказательство. Матрица D1 имеет положительные собственые числа λ1 (D1 ) = 2a/h, λ2 (D1 ) = . . . =
= λN1 (D1 ) = 3a/2h, следовательно, положительно определена. Для всякого p матрица Dp нижнетреуголь-
ная с положительными диагональными элементами. Все ее главные миноры равны произведению соответ-
ствующих диагональных элементов, следовательно, положительны. Матрица Dp положительно определена
по критерию Сильвестра. Лемма доказана.
Перепишем систему (5) в форме
~yj+1 − ~yj
+ sA ~yj+1 + (1 − s)A ~yj = F~j . (7)
∆
Определим линейный оператор B = I + s∆A (I — единичный оператор соответствующей размерности)
и запишем (7) в виде двуслойной разностной схемы в каноническом виде [13]:
~yj+1 − ~yj
B + A ~yj = F~j . (8)
∆
В силу леммы 1, оператор A положительно определен, следовательно, B положительно определен. Таким
образом, B обратим, и можно записать (8) в виде
~yj+1 = S ~yj + ∆B−1 F~j ,
где S = (I − ∆B −1 A) — так называемый оператор продвижения на шаг.
В пространстве Yp определим скалярное произведение и энергетическую норму
N1 ×···×Np
X p
(y, u) = h1 h2 · · · hp y i ui , kykYp = (Ay, y),
i=1
далее определим соответствующую индуцированую операторную норму.
Определение 4 Разностную схему (8) будем называть устойчивой, если kSkYp < 1.
Отметим, что в [13, стр. 324–330] дана эквивалентная формализация устойчивости двуслойных разностных
схем.
Теорема 2 Если s > 1/2, разностная схема (8) устойчива.
319
Доказательство. Рассмотрим (8) с позиции операторно-разностных уравнений и применим методы
проверки устойчивости двуслойной разностной схемы [13]. Проведем симметризацию в уравнении (8), до-
множив обе части уравнения на A−1 , который обратим в силу леммы 1, получим
~yj+1 − ~yj
(A−1 + s∆E) + E ~yj = A−1 F~j .
∆
ˆ
Обозначим B̂ = A−1 + s∆E, Â = E и F~j = A−1 F~j , получим
~yj+1 − ~yj ˆ
B̂ + Â ~yj = F~j . (9)
∆
Метод (9) усточив в энергетической норме тогда и только тогда, когда 2B̂ > ∆ [13, стр. 333, теорема 1].
Это требование эквивалентно тому, что A−1 +∆E(s−0.5) > 0. Поскольку A−1 положительно определенный,
последнее неравенство выполняется для любых ∆, как только s > 0.5. Теорема доказана.
Важной особенностью предлагаемого метода является то, что условие устойчивости s > 1/2 не накла-
дывает ограничений на шаги разбиения, как это имеет место быть в условиях Куранта–Леви.
Методом вложения в общую схему [22] доказывается
Теорема 3 Пусть выполняются условия теоремы 1 и используется оператор кусочно-линейной интер-
p
h2α + ∆.
P
поляции (4). Тогда если s > 0.5, то метод (5) сходится с порядком
α=1
4 Параллельная реализация и численные эксперименты
Для простоты изложения все конструкции и обозначения в этом разделе вводятся для случая p = 2; по
этой причине результаты численных экспериментов также приводятся для случая p = 2.
Для нахождения uij1 ,i2 при i1 > 2, i2 > 2, требуется найти решение в 9 узлах сетки (рис. 1а). Такой
тип рекуррентных информационных зависимостей позволяет последовательно находить решения. Соответ-
ственно, метод (5) может быть запрограммирован в виде трех вложенных циклов: внешний — по времени,
два внутренних — по пространству.
Рис. 1: а) сетка и информационные зависимости; б) расположение блоков в пространстве и информационная
зависимость между блоками.
Реорганизуем вычисления таким образом, чтобы решения в некоторых узлах сетки можно было
находить параллельно. Разрежем сетку на K1 × K2 пространственных колодцев. Назовем (k1 , k2 )-м,
k1 = 1, ..., K1 , k2 = 1, ..., K2 , пространственным колодцем множество узлов сетки
(i1 , i2 , j) i1 = (k1 − 1)N1 /K1 , ..., k1 N1 /K1 − 1, i2 = (k2 − 1)N2 /K2 , ..., k2 N2 /K2 − 1, j = 1, ..., M .
320
Задачу нахождения сеточного решения разделим между процессами; именно, пусть (k1 , k2 )-й процесс вы-
числяет решение в (k1 , k2 )-м пространственном колодце, k1 = 1, ..., K1 , k2 = 1, ..., K2 , т. е. требуется орга-
низовать K1 × K2 потоков.
Поток (1, 1) не имеет информационных зависимостей от других потоков. Для нахождения решения на
(j +1)-м временном слое, процесс (k1 , k2 ), k1 > 2, k2 > 2, должен получить от процесса (k1 −1, k2 ) значения
uij1 ,i2 , где i1 = (k1 − 1)N1 /K1 − 2, (k1 − 1)N1 /K1 − 1, i2 = (k2 − 1)N2 /K2 , ..., k2 N2 /K2 − 1, и от процесса
(k1 , k2 −1) — значения uij1 ,i2 , где i1 = (k1 −1)N1 /K1 , ..., k1 N1 /K1 −1, i2 = (k2 −1)N2 /K2 −2, (k2 −1)N2 /K2 −1,
в этом состоит информационная зависимость между потоками.
Для каждого фиксированного набора k1 , k2 , j назовем блоком Bjk1 ,k2 cледующее множество узлов:
(i1 , i2 , j) i1 = (k1 − 1)N1 /K1 , ..., k1 N1 /K1 − 1, i2 = (k2 − 1)N2 /K2 , ..., k2 N2 /K2 − 1, j = 1, ..., M .
Поскольку блоки являются объединением узлов сетки, они наследуют информационные зависимости, имев-
шиеся между узлами; именно, для нахождения решения в блоке необходимо, чтобы был завершен счет в
блоках, от которых он имеет информационные зависимости: Bjk1 −1,k2 и Bjk1 ,k2 −1 . Кроме того, блоки есте-
ственным образом упорядочиваются в пространстве (x1 , x2 , t) (рис. 2б).
Для обеспечения возможности параллельного счета необходимо построить такие сепарабельные гипер-
плоскости [21, 23], которые проходят через блоки таким образом, что блоки, лежащие на данной гипер-
плоскости, имеют информационные зависимости только от блоков, лежащих по одну сторону от данной
гиперплоскости. Блоки, в которых параллельно находятся решения, лежат на одной сепарабельной гипер-
k0 ,k0
плоскости. Не сложно показать, что блоки Bjk1 ,k2 и Bj 01 2 лежат на одной сепарабельной гиперплоскости
тогда и только тогда, когда k1 + k2 + j = k10 + k20 + j 0 . Так, например, в блоках B21,1 , B12,1 и B11,2 реше-
ния можно искать параллельно. Поскольку задача нахождения решения в блоке трудоемкая, а передача
данных между соседними потоками требует мало времени, эффективность распараллеливания высокая.
Численные эксперименты проводились на кластере ИММ УрО РАН «Уран». Использовалось следующее
оборудование и програмное обеспечение: процесоры INTEL XEON E5450, 2 × 4 ядра, 3 ГГц; кэш память
второго уровня 2 × 6 МБ (5400 Sequence); оперативная память 16 ГБ DDR2; операционная система Linux
2.6.32. Код написан на языке C, использован компилятор Intel C++ compiler (ICC) v14.0.0, использована
библиотека MPI MVAPICH2 Intel 13.0. Используемую вычислительную архитектуру можно охарактеризо-
вать как гибридную, так как на каждом узле используется как распределенная, так и общая память.
Выше было использовано понятие «процесс» для обозначения абстрактной единицы, выполняющей вы-
числения, конкретизируем данное понятие. Согласно принципам MPI, число процессов равно числу вы-
числительных узлов.
Чтобы проиллюстрировать качество распараллеливания, рассмотрим тестовое уравнение
∂u ∂u ∂u
+ + = 2u(x, y, t) − u(x, y, t − 1) − x − y − t + 2,
∂x ∂y ∂t
где 0 6 x 6 2, 0 6 y 6 2, 0 6 t 6 4, с начальными и граничными условиями
u(x, y, t) = x + y + t, 0 6 x 6 2, 0 6 y 6 2, −1 6 t 6 0,
u(0, y, t) = y + t, u(x, 0, t) = x + t, 0 6 x 6 2, 0 6 y 6 2, 0 6 t 6 4.
Легко проверить, что поставленная начально-краевая задача имеет точное решение u(x, y, t) = x + y + t.
В силу того, что решение линейно зависит от всех аргументов, а метод (5) имеет порядок второй по
пространству и первый по времени, погрешность численного решения во всех экспериментах будет мала,
порядка 10−14 .
Таблица 1 содержит результаты численных экспериментов. Число отрезков разбиения по пространствен-
ным переменным выбрано одинаковым и изменяется от 64 до 1024, число отрезков разбиения по времени
постоянно и равно 100. Число процессов указано в первой строке таблицы и меняется от 1 до 64. Ячейки
таблицы содержат медианное время счета в минутах. Ускорение рассчитывалось по отношению к после-
довательному, т. е. однопроцессорному, варианту.
Для оценки качества распараллеливания сделана оценка сильной масштабируемости. Сильная масшта-
бируемость характеризует изменение времени решения задачи с увеличением количества процессоров (или
вычислительных узлов) при неизменном общем объеме задачи [23]. Сильная масштабируемость определя-
ется и может быть численно найдена как отношение T (N )/(N0 T (N0 )), где T (N ) — время решения задачи
321
Таблица 1: Время решения тестового уравнения для сеток различного размера на разном числе процессоров
и оценка ускорения
Число процессоров
Размеры сетки 1 4 9 16 25 36 49 64
64 × 64 × 100 0.039 0.022 0.02 0.014 0.015 0.015 0.018 0.017
128 × 128 × 100 0.152 0.043 0.023 0.018 0.018 0.019 0.02 0.019
256 × 256 × 100 0.478 0.16 0.079 0.049 0.038 0.031 0.031 0.028
512 × 512 × 100 2.103 0.503 0.274 0.174 0.12 0.088 0.074 0.068
1024 × 1024 × 100 6.934 2.267 0.879 0.595 0.399 0.299 0.238 0.201
Ускорение 3.058 7.889 11.647 17.399 23.165 29.179 34.471
Рис. 2: Сильная масштабируемость на примере тестового уравнения: сплошная — экспериментальная оцен-
ка, пунктирная — идеальная сильная масштабируемость
с использованием N узлов, N0 — начальное число узлов. При N0 = 1 производится сравнение параллель-
ной версии и последовательной, что может приводить к некорректным результатам сравнения, а потому
в экспериментах положили N0 = 4. В случае идеальной масштабируемости время вычислений обратно
пропорционально числу вычислительных узлов. В эксперименте сильная масштабируемость близка к иде-
альной (см. рис. 2).
Чтобы проиллюстрировать качество численного метода, рассмотрим уравнение
Z 0
∂u ∂u ∂u x+y π π
+ + = −2(u − sin t) + u(x, y, t − π) − u(x, y, t + ξ) dξ + + 1 sin t,
∂x ∂y ∂t 1 + x2 + y 2 2 −π/2 2
где −2 < x < 2, −5 < y < 5, 0 < t 6 4π, с начальными и граничными условиями
1
u(x, y, t) = + sin t, −2 6 x 6 2, −5 6 y 6 5, −π 6 t 6 0,
1 + x2 + y 2
1 1
u(−2, y, t) = + sin t, u(x, −5, t) = + sin t, −2 6 x 6 2, −5 6 y 6 5, 0 6 t 6 4π.
5 + y2 26 + x2
Эта начально-краевая задача имеет точное решение u(x, y, t) = 1+x12 +y2 + sin t. В таблице 2 представ-
лена зависимость погрешности численного решения diff = max |uij − u(xi , tj )| методом (5) при s = 0.8 в
i,j
322
зависимости от выбора шагов h и ∆. Для простоты положим h1 = h2 и обозначим этот шаг h.
Таблица 2: Погрешность численного решения в зависимости от шагов
Вариант 1 2 3 4 5 6 7 8
h 1/5 1/10 1/10 1/20 1/5 1/10 1/20 1/40
∆ π/20 π/20 π/40 π/40 π/200 π/200 π/200 π/200
diff 0.0902 0.0919 0.0481 0.0495 0.0530 0.0135 0.0038 0.0031
В экспериментах № 5—7 погрешность, обусловленная дискретизацией по времени, мала в сравнении с
погрешностью, обусловленной дискретизацией по координате; анализ поведения погрешности показывает
квадратичную сходимость по пространственным переменным: уменьшение шага в два раза приводит к
уменьшению погрешности примерно в два раза. Анализ таблицы показывает, что только согласованное
уменьшение шагов приводит к уменьшению погрешности. Так, в экспериментах № 7—8 уменьшение h в
два раза не приводит к соответствующему уменьшению погрешности, так как общая погрешность уже в
большей степени обусловлена дискретизацией по времени.
Согласно теореме 2, при s = 0.8 разностная схема (5) является устойчивой при любом соотношении
шагов, однако в силу некорректности операции численного дифференцирования при уменьшении h ап-
проксимация производных по пространству в (5) становится более чувствительной к ошибкам машинного
округления, что приводит к росту погрешности. Согласованное с h уменьшение ∆ является своего рода
регуляризатором, который не позволяет вычислительным ошибкам расти и накапливаться. Эксперименты
№ 1—4 служат тому иллюстрацией.
Отметим, что порядок численного интегрирования для вычисления операторов распределенного запаз-
дывания должен быть согласован с порядком метода (5).
Благодарности
Работа поддержана грантом РФФИ 14-01-00065, РНФ 14-35-00005 и программой 02.A03.21.0006 от
27.08.2013.
Список литературы
[1] B. A. van Tiggelen, S. Skipetrov (Eds.). Wave Scattering in Complex Media: From Theory to Applications.
Proceedings of the NATO Advanced Study Institute, Corsica, France, 10-22 June 2002.
[2] S.A. Gourley, R. Liu, J. H. Wu. Spatiotemporal Distributions of Migratory Birds: Patchy Models with
Delay. SIAM J. Appl. Dyn. Syst., 9(2): 589–610, 2010.
[3] T. Luzyanina, J. Cupovic, B. Ludewig and G. Bocharov. Mathematical models for CFSE labelled lymphocyte
dynamics: asymmetry and time-lag in division. J. Math. Biol., 69(6-7): 1547–83, 2013.
[4] J. Wu. Theory and applications of partial functional differential equations. New York: Springer–Verlag,
1996.
[5] L. Tavernini. Finite difference approximation for a class of semilinear Volterra evolution problems. SIAM
J. Numer. Anal., 14(5): 931–949, 1977.
[6] B. Zubik-Kowal. Stability in the Numerical Solution of Linear Parabolic Equations with a Delay Term. BIT
Numerical Mathematics, 41(1): 191–206, 2001.
[7] Z. Kamont and M. Netka. Numerical method of lines for evolution functional differential equations. J.
Numer. Math., 19(1): 63–89, 2011.
[8] Z. Kamont and K. Przadka. Difference methods for first order partial differential-functional equations with
initial-boundary conditions. Computational Mathematics and Mathematical Physics, 31(10): 1476–1488,
1991.
[9] Z. Kamont and K. Kropielnicka. Implicit difference methods for evolution functional differential equations.
Siberian Journal of Numerical Mathematics, 14(4): 361-379, 2011.
323
[10] W. Czernous and Z. Kamont. Comparison of explicit and implicit difference methods for quasilinear
functional differential equations. Appl. Math., 38(3): 315-340, 2011.
[11] V. G. Pimenov, A. B. Lozhnikov. Difference schemes for the numerical solution of the heat conduction
equation with aftereffect. Proceedings of the Steklov Institute of Mathematics, 275: 137–148, 2011.
[12] V. G. Pimenov, S. Sviridov. Numerical methods for advection equations with delay. American Institute of
Physics. Conference Proceeding. Proceedings of 40th International Conference Applications of Mathematics
in Engineering and Economics, 1631(114): 114–121, 2014.
[13] A.A. Samarskii. Theory of differential schemes. Moscow: Nauka, 1989 (in Russian). = А.А. Самарский.
Теория разностных схем. М.: Наука, 1989.
[14] J.H. Bramble, J.E. Pasciak, J. Xu. Parallel multilevel preconditioners. J. Math.Comput., 55: 1–22, 1990.
[15] B. F. Smith, P. E. Bjorstad and W. Gropp. Domain Decomposition: Parallel Multilevel Methods for Elliptic
Partial Differential Equations. Cambridge University Press, 1996.
[16] P.N. Vabishchevich. Parallel domain decomposition algorithms for parabolic problems. Mat. Model., 9(5):
77-86, 1997.
[17] D.P. Yang. Parallel domain decomposition procedures of improved D-D type for parabolic problems. J.
Comput. Appl. Math., 233: 2779–2794, 2010.
[18] Efremov A. et al. A Computational Realization of a Semi-Lagrangian Method for Solving the Advection
Equation. Journal of Applied Mathematics 2014, Article ID 610398, 2014.
[19] S. I. Solodushkin. A difference scheme for the numerical solution of an advection equation with aftereffect.
Russian Mathematics, Allerton Press, 57: 65–70, 2013.
[20] S. I. Solodushkin, I. F. Yumanova and R. H. De Staelen. First order partial differential equations with time
delay and retardation of a state variable, Journal of Computational and Applied mathematics, 289: 322–330,
2015.
[21] S. I. Solodushkin, A. A. Sagoyan, I. F. Iumanova. A Parallel Algorithm for Solving the Advection Equation
with a Retarded Argument. CEUR Workshop Proceedings 1513:57–66, 2015.
[22] V. G. Pimenov, General Linear Methods for the Numerical Solution of Functional-Differential Equations,
Differential Equations, 37(1): 116–127, 2001.
[23] M. McCool, A. D. Robison, J. Reinders. Structured Parallel Programming. Waltham, USA: Elsevier, 2012.
324
Parallel variant of numerical algorithm for solving a multidimensional
advection equation with time delay
Svyatoslav I. Solodushkin1,2 , Arsen A. Sagoyan1 , Irina F. Iumanova1
1 – Ural Federal University (Yekaterinburg, Russia)
2 – Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia)
Keywords: parallel numerical methods, difference scheme, advection equation, time delay, parallel
programming.
We describe a family of difference schemes for the numerical solving of first-order multi-dimensional equations
with partial derivatives and time delay. The order of approximation, stability and order of convergence are
studied. The proposed numerical algorithm allows parallel implementation. High weak and strong scalability of
parallel version are illustrated in test examples. Experimental evaluation of convergence orders are in consistent
with theoretical estimates.
325