Введение

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

Численные методы решения уравнения переноса без запаздывания достаточно хорошо изучены, например, в работах [3,2,5]. Вопросы численного решения уравнений переноса с запаздыванием рассматривались, например, в [8,9]. Методы распараллеливания рассматривались ранее, например, в [10,11,12]. В данной работе исследуются вопросы распараллеливания для конкретного алгоритма.

Рассмотрим уравнение переноса с эффектом наследственности:

∂u(x, t) ∂t + a ∂u(x, t) ∂x = f (x, t, u(x, t), u t (x, •)),(1)0 t T, 0 x X, с краевым условием u(0, t) = γ(t), 0 t T,(2)и начальным условием u(x, s) = ϕ(x, s), 0 x X, −τ s 0.(3)

Здесь x, t независимые переменные, u(x, t) искомая функция, u t (x, •) = {u(x, t + s), −τ s < 0}функция-предыстория искомой функции к моменту t, τ > 0 величина запаздывания.

Предполагается, что функционал f , функции γ(t), ϕ(x, s) и константа a > 0 таковы, что задача (1)-(3) имеет единственное решение.

Дискретизация задачи

Пусть шаг h по переменной x такой, что X/h = N , где N -натуральное, тогда обозначим через

x i = ih ∈ [0, X], i = 0 . . . N . Пусть шаг ∆ по переменной t такой, что τ /∆ = m, где m -натуральное, пусть M = T /∆ , обозначим t j = j∆ ∈ [−τ, T ], j = −m . . . M . Сеткой назовем набор пар {x i , t j }. Через u i j обозначим приближенное значение u(x i , t j ).

При всяком фиксированном i = 0, . . . , N введем дискретную предысторию к моменту t j , j = 0, . . . , M : {u i k } j = {u i k , j − m k j}. Оператором интерполяции-экстраполяции назовем оператор, определенный на множестве всех допустимых предысторий и действующий по правилу I:

{u i k } j → v i j (•) ∈ Q[t j − τ, t j + ∆]. Здесь Q[α, β] это набор функций u(s), кусочно-непрерывных на [α, β] с конечным числом точек разрыва первого рода. В точках разрыва будем считать функцию непрерывной справа. Определим норму Q = Q[α, β] как u(•) Q = max s∈[α,β] |u(s)|.

Будем полагать, что, во-первых, оператор интерполяции-экстраполяции липшицев, то есть существует такая константа L I , что для всех предысторий дискретной модели

{u i k } j и {y i k } j выполняется max tj −τ t tj +∆ |v i j (t) − w i j (t)| L I max j−m l j |u i l − y i j |, где v i j (•) = I({u i k } j ), w i j (•) = I({y i k } j ). Во-вторых, положим, что оператор интерполяции-экстраполяции согласован, то есть v i j (t k ) = u i k , k = j − m, . . . , j.

Будем говорить, что оператор интерполяции-экстраполяции имеет порядок p, если существуют константы

C 1 и C 2 такие, что |v i j (t) − u(x i , t)| C 1 max j−m k j |u i k − u(x i , t k )| + C 2 ∆ p для всех i, j и t ∈ [t j − τ, t j+1 ].

Простейший способ интерполяции кусочно-линейная функция. Простейший способ экстраполяции экстраполяция продолжением (см. [7]). Будем использовать кусочно-линейную интерполяцию с экстраполяцией продолжением, которая является липшицевым оператором (L I = 2), согласована и имеет порядок p = 2 (см. [7]).

u i j+1 − u i j ∆ + a u i j+1 − u i−1 j+1 h = f (x i , t j , u i j , u t (x, •))).(4)

Каждая производная в (4) представлена как разность приближенных значений функции в узлах сетки, показанных на рис. 2. Значение u i−1 j+1 может быть выражено из ранее подсчитанных значений, значит, рассматриваемый метод правый уголок является явным для уравнения (1)- (3).

Обозначим ε i j = u(x i , t j ) − u i j , i = 0, . . . , N , j = 0, . . . , M . Будем говорить, что метод сходится, если ε i j → 0 при h → 0 и ∆ → 0 для всех i = 0, . . . , N и j = 0, . . . , M . Будем говорить, что метод сходится с порядком h p + ∆ q , если существует константа C такая, что |ε i j | C(h p + ∆ q ) для всех i = 0, . . . , N и j = 0, . . . , M. В [4] доказано, что рассматриваемый метод сходится с порядком O(∆ + h). Данная схема входит в семейство методов, рассмотренных в [6].

Идея параллельного вычисления

В данной схеме для вычисления очередных неизвестных значений используются значения, стоящие ниже или левее по x и по t, следовательно, мы можем независимо друг от друга вычислять элементы, находящиеся на одной диагонали (далее будем называть эти диагональные элементы фронтом). Это дает возможность для параллельного вычисления результирующей функции. Схематичное расположение фронтов для параллельного вычисления показано на рис. 3.

Рис. 3: Схема Правый уголок

Для реализации данного подхода использовалось многопоточное программирование. В качестве языка было решено использовать C# и платформу .Net, включающую примитивы многопоточного программирования. При реализации использовалась сборка System.Threading.Tasks, в частности метод Parallel.For, позволяющий распараллелить выполнение цикла For между потоками. Все вычисление выполнялось пошагово по диагоналям. На каждом шаге выбиралась следующая диагональ, для узлов которой параллельно вычислялись искомые значения. Так как элементы на диагонали не зависят друг от друга, то данный подход является корректным с точки зрения вычислений. На каждом этапе параллельного вычисления искомых значений в узлах диагонали платформа порождает 8 потоков и распределяет узлы между ними по мере возможности. Из-за того, что одновременно вычисляется не более 8 потоков и количество элементов диагонали в большинстве случаев больше 8, эффект узкого фронта распараллеливания был незначительным, поэтому в большинстве примеров мы рассматривали сетки с одинаковым числом узлов по осям. Так как целью данной работы является исследование возможности ускорения путем распараллеливания вычислений, то мы не прибегали к другим приемам программных и алгоритмических оптимизаций (векторизации, кешированию и др.). В итоге, мы оценивали потенциал ускорения при помощи лишь многопоточного программирования.

Численные эксперименты

Во всех таблицах с результатами экспериментов используем следующие обозначения: Err максимум модуля разности точного и приближенного решений в узлах сетки. T1 время выполнения последовательной версии (в миллисекундах). T2 время выполнения параллельной версии (в миллисекундах).

∂u ∂t + ∂u ∂x = u x, t 2 − 1 2 + 6x + 3t 2 + 6t − 3x 3 − t 2 3 − 5 8 , 0 t 1, 0 x 2, с краевым условием: u(0, t) = 3t 2 + 6t, 0 t 4, с начальным условием: u(x, t) = x 3 + t 3 + 3x 2 + 3t 2 , 0 x 1, − 1 2 t 0.

Точное решение данного уравнения: Отметим, что численные эксперименты для примера 3 проводились без использования рекуррентных соотношений между интегралами по пересекающимся интервалам. Сделано это было для большего числа вычислений на каждом шаге и анализа эффективности использования параллелизма для численного решения уравнений рассматриваемых типов.

u(x, t) = x 3 + t 3 + 3x 2 + 3t

Заключение

Все вычисления выполнялись на ноутбуке с процессором Intel(R) Core i7 4700MQ 2.4GHz x 4, поддерживающем одновременное выполнение 8 потоков. Для чистоты эксперимента все вычисления проводились в максимально стерильных условиях (на компьютере запущены только необходимые приложения), однако операционная система и фоновые процессы в любом случае влияют на время вычислений, поэтому в примерах указано среднее время выполнения среди десяти запусков. Также можно заметить сильное замедление параллельной программы на малом количестве узлов по сравнению с последовательной версией. Это вызвано временем, которое уходит у платформы на подготовку многопоточной экосистемы и переключение между потоками. Последовательная программа лишена данных недостатков. В качестве эффективности в данной работе рассматривается отношение времени выполнения параллельной программы ко времени выполнения последовательной. В дальнейшем планируется более глубокий анализ алгоритмов параллельного вычисления решений описанных выше уравнений: с учетом различных показателей эффективности, таких как слабая и сильная масштабируемость.

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