GPU-BASED PARALLEL COMPUTATION OF WIND-INDUCED SEA CURRENTS FOR OIL SPILL PROPAGATION MODELING IN THE CASPIAN SEA

Толық мәтін

ВВЕДЕНИЕ
Казахстан является одним из крупнейших производителей нефти в Центральной Азии, а значительная часть крупных нефтегазовых проектов страны связана с Каспийским регионом. К наиболее важным морским месторождениям относятся Кашаган, Кайран и Актоты, входящие в Северо-Каспийский проект. Разведка, добыча и транспортировка нефти в морской акватории неизбежно сопровождаются риском аварийных выбросов и загрязнения морской среды [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-часовой прогноз).

Авторлар туралы

Aziz Kudaikulov

әл-Фараби атындағы Қазақ ұлттық университеті

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

Жетекші ғылыми қызметкер

Қазақстан, 050040, Қазақстан Республикасы, Алматы қ., Әл-Фараби даңғ., 71

Tannur Amanzholov

әл-Фараби атындағы Қазақ ұлттық университеті

Хат алмасуға жауапты Автор.
Email: tannur.amanzholov@gmail.com
ORCID iD: 0000-0001-7386-9717

Аға ғылыми қызметкер

Қазақстан, 050040, Қазақстан Республикасы, Алматы қ., Әл-Фараби даңғ., 71

Kuljabekov Alibek

Қ. И. Сәтбаев атындағы Қазақ ұлттық зерттеу техникалық университеті

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

PhD, зертхана меңгерушісі

Қазақстан, 050013, Қазақстан Республикасы, Алматы қ., Сәтбаев көш., 22

Abdurashid Aliuly

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

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

Аға ғылыми қызметкер

Қазақстан, 050040, Қазақстан Республикасы, Алматы қ., Әл-Фараби даңғ., 71

Abzal Seitov

Al-Farabi Kazakh National University

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

Аға ғылыми қызметкер

Қазақстан, 050040, Қазақстан Республикасы, Алматы қ., Әл-Фараби даңғ., 71

Bakhytzhan Assilbekov

Institute of ionosphere

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

Chief Researcher

Қазақстан, 050020, Republic of Kazakhstan, Almaty, Kamenskoye Plateau, Gardening Partnership “Ionosphere”, 117

Spartak Shabilov

Satbayev Kazakh National Research Technical University

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

Junior Researcher

Қазақстан, 050013, Republic of Kazakhstan, Almaty, 22 Satbaev str.

Samal Syrlybekkyzy

Yessenov Caspian University of Technologies and Engineering

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

PhD, Associate Professor, Dean of the Faculty of Engineering

Қазақстан, 130000, Republic of Kazakhstan, Mangystau Region, Aktau, District #32

Aidarkhan Kaltayev

Satbayev Kazakh National Research Technical University

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

Doctor of Physical and Mathematical Sciences, Professor, Head of Laboratory

Қазақстан, 050013, Republic of Kazakhstan, Almaty, 22 Satbaev str.

Әдебиет тізімі

Қосымша файлдар