ПАРАЛЛЕЛЬНОЕ ВЫЧИСЛЕНИЕ ВЕТРОВЫХ МОРСКИХ ТЕЧЕНИЙ НА GPU ДЛЯ ЗАДАЧ ПРОГНОЗИРОВАНИЯ НЕФТЯНЫХ РАЗЛИВОВ

Полный текст

ВВЕДЕНИЕ
Казахстан является одним из крупнейших производителей нефти в Центральной Азии, а значительная часть крупных нефтегазовых проектов страны связана с Каспийским регионом. К наиболее важным морским месторождениям относятся Кашаган, Кайран и Актоты, входящие в Северо-Каспийский проект. Разведка, добыча и транспортировка нефти в морской акватории неизбежно сопровождаются риском аварийных выбросов и загрязнения морской среды [1, 2]. Поэтому требуется инструменты позволяющие заранее прогнозировать направление распространения нефти, скорость её перемещения, площадь загрязнения и потенциальное воздействие на береговую линию и морскую экосистему. Если спутниковое зондирование обеспечивает фиксацию имеющегося нефтяного пятна, то численное моделирование выступает его необходимым прогностическим дополнением, позволяющим непрерывно рассчитывать дальнейший перенос и трансформацию разлива. 
Одним из ключевых этапов моделирования распространения разлитой нефти по поверхности моря является расчет морских течений, вызванных ветровым воздействием. Для решения этой задачи были разработаны различные модели, различающиеся по уровню сложности. В [3-5] была предложена гидродинамическая модель однородной жидкости с упрощенным представлением турбулентности, основанным на постоянной турбулентной вязкости. Позже в исследовании [6-8] этот подход был расширен за счет включения эффектов стратификации плотности и температуры, а также адвекции, что сделало модель нелинейной; однако процедура решения осталась по существу той же, что и в работе [3].
 В работах [9, 10] детально рассматривались взаимодействия морского льда и атмосферы с морем с целью исследования сезонной циркуляции и изменчивости водных масс Каспийского моря под воздействием реалистичных потоков массы, импульса и тепла. В этих моделях речной сток, осадки, радиационный поток и ветровое напряжение задавались с учётом сезонной изменчивости на основе доступных наборов данных. Кроме того, интенсивность испарения, поток явного и скрытого тепла на поверхности моря рассчитывались в ходе моделирования с использованием подмодели атмосферного пограничного слоя на основе атмосферных данных реанализа ECMWF ERA-15 и рассчитанной моделью температуры поверхности моря.
Во всех вышеупомянутых работах использовались упрощённые модели турбулентности, в которых турбулентная вязкость задавалась параметрически. Вместе с тем существуют исследования, в которых применяются более совершенные модели турбулентности. Например, в работах [11–13] использовалась модель рейнольдсовых напряжений, требующая решения значительно большего числа уравнений. Однако, поскольку зависимости для определения турбулентной вязкости хорошо апробированы для широкого диапазона условий, на практике зачастую целесообразнее применять простые модели турбулентной вязкости, а не более сложные подходы, такие как модели рейнольдсовых напряжений.
Кроме того, все вышеупомянутые модели являются трёхмерными и поэтому требуют значительных вычислительных затрат и ресурсов [14, 15]. Например, в работе [15] расчёты выполнялись на сетке размером 426 × 238 × 100 узлов, что соответствует более чем 10 миллионам расчётных точек. Как правило, такие задачи решаются на вычислительных кластерах с большим числом процессоров с использованием параллельных алгоритмов. В настоящем исследовании представлены параллельные алгоритмы решения вышеупомянутых моделей на графических процессорах с применением библиотеки OpenCL [16]. Результаты вычислений показывают, что может быть достигнуто значительное ускорение, поскольку данные модели хорошо подходят для распараллеливания.
Цель настоящей работы – разработать и реализовать параллельные алгоритмы для трёхмерного моделирования ветровых морских течений на графических процессорах с использованием OpenCL, обеспечивающие существенное сокращение времени расчётов при сохранении приемлемой точности для задач оперативного прогнозирования нефтяных разливов.

МАТЕРИАЛЫ И МЕТОДЫ
Основные уравнения
Согласно работе [3], линейные трёхмерные уравнения гидродинамики в декартовой системе координат записываются в следующем виде:

∂ζ/∂t=-∂(hu ̅ )/∂x-∂(hv ̅ )/∂y                (1)
(∂u ̅)/∂t=fv ̅-g ∂ζ/∂x-(τ_x^h)/ρh+(τ_x^s)/ρh            (2)
(∂v ̅)/∂t=-fu ̅-g ∂ζ/∂y-(τ_y^h)/ρh+(τ_y^s)/ρh            (3)

где t это время, x, y, z образуют левостороннюю декартову систему координат, а координата z представляет собой глубину ниже невозмущённой поверхности, h это глубина воды, ζ это отметка свободной поверхности; u ̅, v ̅ обозначают осреднённые по глубине компоненты скорости течения в направлениях x и y, соответственно:

u ̅=∫_0^1▒u dσ,  v ̅=∫_0^1▒v dσ            (4)
где σ=z/h.

Основные положения, рассмотренные выше, представлены на рис. 1.
 
Рисунок 1 – Система координат сигма

В дополнение к приведенным выше уравнениям, ниже представлены уравнения для отклонений u' and v', которые используются для получения полных компонент скорости u and v:

u=u ̅+u',  v=v ̅+v'            (5)
∂u'/∂t=fv'+1/h^2   ∂/∂σ (μ ∂u'/∂σ)+(τ_x^h)/ρh-(τ_x^s)/ρh        (6)
∂v'/∂t=-fu'+1/h^2   ∂/∂σ (μ ∂v'/∂σ)+(τ_y^h)/ρh-(τ_y^s)/ρh        (7)

где v – коэффициент вертикальной вихревой вязкости, g – ускорение свободного падения, f – параметр Кориолиса, ρ – плотность, τ_x^s и τ_y^s обозначим компоненты ветровой нагрузки:
τ_x^s=-ρ/h (μ ∂u/∂σ)_0,  τ_y^s=-ρ/h (μ ∂v/∂σ)_0    (8)
Аналогично, на морском дне σ=1, τ_x^h и τ_y^h обозначают компоненты придонного напряжения:
τ_x^h=-ρ/h (μ ∂u/∂σ)_1,  τ_y^h=-ρ/h (μ ∂v/∂σ)_1    (9)
В описанной далее модели компоненты ветрового напряжения τ_x^s и τ_y^s предопределены, и был использован линейный закон донного напряжения

τ_x^h=kρu_d,  τ_y^h=kρv_d            (10)

где k – коэффициент трения дна, u_d и v_d – компоненты течения на высоте d над морским дном. 

Граничные условия
Верхние граничные условия задаются на поверхности моря σ=0:

ρ/h (μ ∂u'/∂σ)_0=-τ_x^s,  ρ/h (μ ∂v'/∂σ)_0=-τ_y^s        (11)

где компоненты ветрового напряжения τ_x^s и τ_y^s задаются заранее.
На морском дне, при σ=1, соответствующие граничные условия имеют вид:

τ_x^h=-ρ/h (μ ∂u'/∂σ)_1=kρu_d,  τ_y^h=-ρ/h (μ ∂v'/∂σ)_1=kρv_d    (12)

На боковых границах задаются периодические граничные условия:

z(0,y)=z(L,y),  z(x,0)=z(x,W)    (13)
u ̅(0,y)=u ̅(L,y),  u ̅(x,0)=u ̅(x,W)    (14)
v ̅(0,y)=v ̅(L,y),  v ̅(x,0)=v ̅(x,W)    (15)

где L и W – длина и ширина расчётной области соответственно.

Численная схема
В вертикальном направлении используется разнесённая конечно-разностная сетка: точки μ располагаются в граничных и внутренних узлах стандартной сетки, тогда как точки u и v размещаются посередине между соседними точками μ. Точки ζ расположены на верхней границе расчётной области (рис. 2).

 
Рисунок 2 – Разнесённая сетка

В данном разделе приводится краткое описание конечно-разностной сетки, используемой в модели.
Уравнения (1)–(3) дискретизируются следующим образом:

u ̅_(i,j)^(n+1)=u ̅_(i,j)^n+fv ̅_(i,j)^n△t-g△t (∂ζ_(i,j)^n)/∂x-(τ_x^h)/ρh△t+(τ_x^s)/ρh△t        (16)
v ̅_(i,j)^(n+1)=v ̅_(i,j)^n-fu ̅_(i,j)^n△t-g△t (∂ζ_(i,j)^n)/∂y-(τ_y^h)/ρh△t+(τ_y^s)/ρh△t         (17)
ζ_(i,j)^(n+1)=ζ_(i,j)^n-(∂(hu ̅ )_(i,j)^(n+1))/∂x△t-(∂(hv ̅ )_(i,j)^(n+1))/∂y△t             (18)

Подставляя значения u ̅_(i,j)^( n+1) и v ̅_(i,j)^( n+1) определяемые выражениями (16) и (17), в уравнение (18), получаем:

ζ_(i,j)^(n+1)=ζ_(i,j)^n+gh△t^2 [(∂^2 ζ_(i,j)^n)/(∂x^2 )+(∂^2 ζ_(i,j)^n)/(∂y^2 )]+S_(i,j)     (19)

Шаг по времени ∆t=τ_g в уравнении (19) определяется из условия устойчивости:

τ_g<(△s)/√2gh,  △s=min(△x,△y)         (20)

При глубине моря h=100 m и шаге сетки ∆s=1000 m, из уравнения (20) следует, что τ_g∼1 min.

Подставляя значение ζ_(i,j)^(n+1) определяемое выражением (19), в уравнение (16), получаем:

u ̅_(i,j)^(n+1)=u ̅_(i,j)^n-f^2△t^2 u ̅_(i,j)^n+R_(i,j)         (21)

Шаг по времени ∆t=τ_f в уравнении (21) определяется из условия устойчивости при f=1.2×10^(-4)  s^(-1):

τ_f<√2/f≈1 hours                 (22)

Уравнения (6) и (7) дискретизируются следующим образом:

█(u'_(i,j,k)^( n+1)&=u'_(i,j,k)^( n)+fv'_(i,j,k)^( n) △t@& +(△t)/(h^2△σ_k ) [(μ_(i,j,k+1)^n (u'_(i,j,k+1)^( n+1)-u'_(i,j,k)^( n+1) ))/(△σ_(k+1)+△σ_k )-(μ_(i,j,k)^n (u'_(i,j,k)^( n+1)-u'_(i,j,k-1)^( n+1) ))/(△σ_k+△σ_(k-1) )]@& +(△t)/(h^2△σ_k ) [(μ_(i,j,k+1)^n (u'_(i,j,k+1)^( n)-u'_(i,j,k)^( n) ))/(△σ_(k+1)+△σ_k )-(μ_(i,j,k)^n (u'_(i,j,k)^( n)-u'_(i,j,k-1)^( n) ))/(△σ_k+△σ_(k-1) )]@& +(ku_(i,j,k)^n)/h_(i,j) △t-(τ_(x,i,j)^(s,n))/(ρh_(i,j) )△t) (23)

█(v'_(i,j,k)^( n+1)&=v'_(i,j,k)^( n)-fu'_(i,j,k)^( n) △t@& +(△t)/(h^2△σ_k ) [(μ_(i,j,k+1)^n (v'_(i,j,k+1)^( n+1)-v'_(i,j,k)^( n+1) ))/(△σ_(k+1)+△σ_k )-(μ_(i,j,k)^n (v'_(i,j,k)^( n+1)-v'_(i,j,k-1)^( n+1) ))/(△σ_k+△σ_(k-1) )]@& +(△t)/(h^2△σ_k ) [(μ_(i,j,k+1)^n (v'_(i,j,k+1)^( n)-v'_(i,j,k)^( n))┤)/(△σ_(k+1)+△σ_k )-(μ_(i,j,k)^n (v'_(i,j,k)^( n)-v'_(i,j,k-1)^( n) ))/(△σ_k+△σ_(k-1) )]@& +(kv_(i,j,k)^n)/h_(i,j)  △t-(τ_(y,i,j)^(s,n))/(ρh_(i,j) ) △t)  (24)

Шаг по времени ∆t=τ_f в уравнениях (23) и (24) определяется условием устойчивости:

τ_f<√2/f≈1 hours                     (25)

Параллелизация на GPU
Численные схемы (16)–(18) представлены в явном виде и могут быть легко распараллелены. В то же время численные схемы (23)–(24) сформулированы в неявном виде и могут быть решены с использованием алгоритма решения трёхдиагональных систем линейных алгебраических уравнений (СЛАУ). Несмотря на существование параллельных алгоритмов решения трёхдиагональных СЛАУ, их программная реализация зачастую сопряжена со значительными трудностями. Вместе с тем из численных схем (23)–(24) следует, что вычисления могут быть легко распараллелены по направлениям x и y. В настоящем исследовании анализируется влияние числа точек дискретизации по координате z на время вычислений. При этом схемы (23)–(24) распараллеливаются по направлениям x и y, а вычисления по координате z выполняются последовательно. Исходный код последовательного и параллельного алгоритмов доступен в работе [17].

Результаты
Численные расчёты проводились при следующих значениях параметров: f=1.2×10^(-4) s^(-1), g=9.81 m s^(-2), ρ=1025 kg m^(-3), k=0.002 m s^(-1), μ=0.4 m^2 s^(-1), H=260 m, L=260 m, W=260 m.
Течения генерировались ветровым напряжением, направленным под углом 45^∘ к осям x и y, при этом значения его компонентов составляли τ_x^s=1.5 N m^(-2) и τ_y^s=1.5 N m^(-2), соответственно.
Для оценки точности расчёта течений при различных шагах сетки использовались два варианта распределения расчётных узлов. В первом случае применялась равномерная сетка, содержащая N_z=500 узлов вдоль координаты z. Во втором случае использовалось N_z=40 узлов по координате z с повышенным разрешением в поверхностном пограничном слое (△σ=0.002). По координатам x и y, применялась равномерная сетка с числом узлов N_x=20 и N_y=20, соответственно.
Все расчёты выполнялись на центральном процессоре Intel Core i7 с тактовой частотой 3,6 ГГц и графическом процессоре NVIDIA GeForce GTX 670, содержащем 1344 вычислительных ядра и работающем на частоте 915 МГц.
Сравнение профилей скорости, полученных на указанных двух типах сеток, представлено на рисунках 3–8. Видно, что точность расчётов на неравномерной сетке с (N_z=40) сопоставима с точностью, достигаемой на равномерной сетке с (N_z=500). Таким образом, использование неравномерной сетки с (N_z=40) является более целесообразным по сравнению с равномерной сеткой с (N_z=500).
На рис. 9 представлено сравнение времени решения задачи на равномерной сетке с (N_z=500) при использовании CPU и GPU. Для данной сетки вычисления на GPU выполнялись в 23 раза быстрее, чем на CPU.
На рис. 10 представлено сравнение времени решения задачи на неравномерной сетке с (N_z=40) при использовании CPU и GPU. Для данной сетки вычисления на GPU выполнялись в 13 раз быстрее, чем на CPU.

 
Рисунок 3 – Профили компоненты скорости u при x = 130 м, y = 130 м (t = 1час 40 мин)

 
Рисунок 4 – Профили компоненты скорости v при x = 130 м, y = 130 м (t = 1час 40 мин)

 
Рисунок 5 – Профили компоненты скорости u при x = 130 м, y = 130 м (t = 6 час 40 мин)

 
Рисунок 6 – Профили компоненты скорости v при x = 130 м, y = 130 м (t = 6 час 40 мин)

 
Рисунок 7 – Профили компоненты скорости u при x = 130 м, y = 130 м (t = 15 час)

 
Рисунок 8 – Профили компоненты скорости v при x = 130 м, y = 130 м (t = 15 час)

 
Рисунок 9 – Сравнение времени вычислений для равномерной сетки

 
Рисунок 10 – Сравнение времени вычислений для неравномерной сетки

Обсуждение
Полученное ускорение (23× для равномерной сетки, 13× для неравномерной) согласуется с известными результатами применения GPU в гидродинамике, где типичный выигрыш составляет 10–50× в зависимости от размера сетки и архитектуры. Важно отметить, что для неравномерной сетки с 40 слоями ускорение ниже, но это компенсируется меньшим объёмом вычислений и сохранением точности в поверхностном слое, что критично для моделирования разливов.
Использование OpenCL обеспечивает переносимость на различные GPU и даже на кластеры, что расширяет область применения. Однако явный шаг по времени (ограниченный условием Куранта) остаётся узким местом. В перспективе целесообразно перейти к полностью неявным схемам, позволяющим увеличить шаг до часов, что ещё больше повысит эффективность.
Сравнение с работами [14, 15] показывает, что на сетках аналогичного размера наши результаты сопоставимы, при этом мы использовали более доступное оборудование (один GPU вместо кластера). Предложенная стратегия распараллеливания (по горизонтали) может быть масштабирована на большее число узлов.

ЗАКЛЮЧЕНИЕ
В данной работе численное моделирование ветровых течений в море проводилось на равномерной и неравномерной сетках с использованием CPU и GPU. Поскольку глубина моря относительно невелика — около 100 м — по сравнению с горизонтальными размерами расчётной области по осям x и y, число вычислительных узлов по оси z также меньше, чем по направлениям x и y. Поэтому, как подтверждают полученные результаты, для решения трёхдиагональных систем линейных алгебраических уравнений более эффективным является применение последовательного алгоритма, а не параллельного. Кроме того, результаты показали, что ускорение вычислений было выше для равномерной сетки с N_z=500 чем для неравномерной сетки с N_z=40.
При решении более сложных моделей возникает необходимость применения неявных схем, позволяющих использовать более крупные шаги по времени. Полученные неявные разностные уравнения могут решаться как прямыми методами, например посредством решения трёхдиагональных систем линейных алгебраических уравнений, использованного в настоящей работе, так и итерационными методами. Преимущество прямых методов заключается в их вычислительной эффективности при любых начальных и граничных условиях.
Разработанный параллельный алгоритм может быть интегрирован в систему оперативного прогнозирования нефтяных разливов в казахстанском секторе Каспия, обеспечивая расчёт полей течений за время, приемлемое для принятия решений (менее 1 часа на 24-часовой прогноз).

Об авторах

Азиз Кудайкулов

Казахский национальный университет имени аль-Фараби

Email: azizkudaikulov@gmail.com
ORCID iD: 0000-0002-4389-4097

ВНС

Казахстан, 050040, Республика Казахстан, г. Алматы, пр. аль-Фараби, 71

Таннур Аманжолов

Казахский национальный университет имени аль-Фараби

Автор, ответственный за переписку.
Email: tannur.amanzholov@gmail.com
ORCID iD: 0000-0001-7386-9717

СНС

Казахстан, 050040, Республика Казахстан, г. Алматы, пр. Аль-Фараби, 71

Кульджабеков Алибек

Казахский национальный исследовательский технический университет

Email: Alibek.kuljabekov@gmail.com
ORCID iD: 0000-0003-4384-6463

PhD, Заведующий лабораторией вычислительного моделирования и информационных технологий

Казахстан, 050013, Республика Казахстан, г. Алматы, ул. Сатпаева, 22

Абдурашид Алиулы

Казахский национальный исследовательский технический университет имени К. И. Сатпаева

Email: aliulyabdurashid@gmail.com
ORCID iD: 0000-0001-6845-6571

СНС

Казахстан, 050040, Республика Казахстан, г. Алматы, пр. Аль-Фараби, 71

Абзал Сейтов

Казахский национальный университет имени аль-Фараби

Email: seitovabzal9@gmail.com
ORCID iD: 0000-0002-7317-3047

СНС

Казахстан, 050040, Республика Казахстан, г. Алматы, пр. Аль-Фараби, 71

Бакытжан Асылбеков

ТОО «Институт ионосферы»

Email: assilbekov.b@gmail.com
ORCID iD: 0000-0002-0368-0131

ГНС

Казахстан, 050020, Республика Казахстан, г. Алматы, Каменское плато, Садоводческое товарищество «Ионосфера», д. 117

Спартак Шабилов

Казахский национальный исследовательский технический университет имени К. И. Сатпаева

Email: sshabilov@gmail.com
ORCID iD: 0009-0002-3467-5842

МНС

Казахстан, 050013, Республика Казахстан, г. Алматы, ул. Сатпаева, 22

Самал Сырлыбеккызы

Каспийский университет технологий и инжиниринга имени Ш. Есенова

Email: samal.syrlybekkyzy@yuedu.kz
ORCID iD: 0000-0002-0260-0611

PhD, ассоциированный профессор, декан факультета «Инжиниринг»

Казахстан, 130000, Республика Казахстан, Мангистауская область, г. Актау, 32 микрорайон

Айдархан Калтаев

Казахский национальный исследовательский технический университет имени К. И. Сатпаева

Email: aidarkhankaltayev@gmail.com
ORCID iD: 0000-0003-2180-2785

д.ф.-м.н., профессор, заведующий лабораторей

Казахстан, 050013, Республика Казахстан, г. Алматы, ул. Сатпаева, 22

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

Дополнительные файлы