Вестник нефтегазовой отрасли Казахстана

2707-42262957-806X

KMG Engineering

109025

10.54859/kjogi109025

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

Kudaikulov

Aziz

Leading Researcherazizkudaikulov@gmail.comhttps://orcid.org/0000-0002-4389-4097Amanzholov

Tannur

Senior Researchertannur.amanzholov@gmail.comhttps://orcid.org/0000-0001-7386-9717Alibek

Kuljabekov

PhD, Head of LaboratoryAlibek.kuljabekov@gmail.comhttps://orcid.org/0000-0003-4384-6463Aliuly

Abdurashid

Senior Researcheraliulyabdurashid@gmail.comhttps://orcid.org/0000-0001-6845-6571Seitov

Abzal

seitovabzal9@gmail.comhttps://orcid.org/0000-0002-7317-3047Assilbekov

Bakhytzhan

Chief Researcherassilbekov.b@gmail.comhttps://orcid.org/0000-0002-0368-0131Shabilov

Spartak

Junior Researchersshabilov@gmail.comhttps://orcid.org/0009-0002-3467-5842Syrlybekkyzy

Samal

PhD, Associate Professor, Dean of the Faculty of Engineeringsamal.syrlybekkyzy@yuedu.kzhttps://orcid.org/0000-0002-0260-0611Kaltayev

Aidarkhan

Doctor of Physical and Mathematical Sciences, Professor, Head of Laboratoryaidarkhankaltayev@gmail.comhttps://orcid.org/0000-0003-2180-2785

Al-Farabi Kazakh National Universityl-Farabi Kazakh National UniversitySatbayev Kazakh National Research Technical UniversityInstitute of ionosphereYessenov Caspian University of Technologies and Engineering83

28072026

09092026

Relevance. Kazakhstan is one of the largest oil producers in Central Asia, with major projects concentrated in the Caspian region. Oil spills in marine areas pose a serious ecological threat, making it essential to develop effective forecasting tools for their spread. Materials and Methods. Three-dimensional hydrodynamic models with parameterized turbulent viscosity were applied. Numerical schemes were implemented in explicit and implicit forms, parallelized along x and y coordinates using the OpenCL library. Simulations were performed on CPU and GPU for both uniform and non-uniform grids. Results. GPU-based computations demonstrated significant acceleration: up to 23 times faster on uniform grids and up to 13 times faster on non-uniform grids. The accuracy of simulations with 40 vertical nodes was comparable to that of uniform grids with 500 nodes. Conclusion. Parallel GPU algorithms substantially reduce computation time while maintaining high accuracy in modeling wind-driven currents. This makes them a promising tool for forecasting oil spill dispersion in the Caspian Sea.

GPU, OpenCL, parallel computing, wind velocity field, sea currents

GPUOpenCLпараллельные вычисленияполе скоростей ветраморские течения

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

1. U.S. Energy Information Administration. Oil and Gas – Kazakhstan: Explore oil and gas export opportunities and the regulatory environment in Kazakhstan [Электронный ресурс]. – Режим доступа: https://www.trade.gov/energy-resource-guide-oil-and-gas-kazakhstan?anchor=content-node-t7-field-lp-region-2-1 (дата обращения: 01.06.2026).

2. North Caspian Operating Company. Отчет об устойчивом развитии за 2022 год [Электронный ресурс]. – Режим доступа: https://www.ncoc.kz/ru/sustainability/2022/reports (дата обращения: 01.06.2026).

3. Davies A.M. Application of a sigma coordinate sea model to the calculation of wind-induced currents // Continental Shelf Research. – 1985. – Vol. 4, No. 4. – P. 389–423.

4. Elliott A.J. Shear diffusion and the spread of oil in the surface layers of the North Sea // Deutsche Hydrographische Zeitschrift. – 1986. – Vol. 39. – P. 113–137.

5. Davies A.M., Stephens C.V. Comparison of the finite difference and Galerkin methods as applied to the solution of the hydrodynamic equations // Applied Mathematical Modelling. – 1983. – Vol. 7. – P. 226–240.

6. Zuur E.A.H., Dietrich D.E. The SOMS model and its application to Lake Neuchatel // Aquatic Sciences. – 1990. – Vol. 52, No. 2. – P. 115–129.

7. Korotenko K., Bowman M.J., Dietrich D.E. High-Resolution Numerical Model for Predicting the Transport and Dispersal of Oil Spilled in the Black Sea // Terrestrial, Atmospheric and Oceanic Sciences. – 2010. – Vol. 21, No. 1. – P. 123–136.

8. Korotenko K., Bowman M.J., Dietrich D.E. Modeling of the circulation and transport of oil spills in the Black Sea // Oceanology. – 2003. – Vol. 43. – P. 367–378.

9. Ibrayev R.A. Model of enclosed and semi enclosed sea hydrodynamics // Russian Journal of Numerical Analysis and Mathematical Modelling. – 2001. – Vol. 16, No. 4. – P. 291–304.

10.

10. Ibrayev R.A., Ozsoy E., Schrum C., Sur H.I. Seasonal variability of the Caspian Sea three dimensional circulation, sea level and air sea interaction // Ocean Science. – 2010. – Vol. 6, No. 1. – P. 311–329.

11.

11. Blumberg A.F., Mellor G.L. A Description of a three dimensional coastal ocean circulation model // Coastal and Estuarine Sciences. – 1987. – Vol. 4. – P. 1–16.

12.

12. Mellor G.L. On the two components of wind driven ocean surface stress with extension to scalar fluxes // Ocean Dynamics. – 2019. – Vol. 69. – P. 43–50.

13.

13. Mellor G.L. On theories dealing with the interaction of surface waves and ocean circulation // Journal of Geophysical Research: Oceans. – 2016. – Vol. 121, No. 7. – P. 4474–4486.

14.

14. Korotenko K., Osadchiev A., Melnikov V. Mesoscale Eddies in the Black Sea and Their Impact on River Plumes: Numerical Modeling and Satellite Observations // Remote Sensing. – 2022. – Vol. 14, No. 17. – Art. no. 4149.

15.

15. Korotenko K., Osadchiev A., Melnikov V. Mesoscale Eddy Chain Structures in the Black Sea and Their Interaction with River Plumes: Numerical Modeling and Satellite Observations // Remote Sensing. – 2023. – Vol. 15, No. 6. – P. 1606–1628.

16.

16. OpenCL Homepage [Электронный ресурс]. – Режим доступа: https://www.khronos.org/opencl (дата обращения: 01.06.2026).

17.

17. GitHub repository: oil_trans_in_caspian_sea [Электронный ресурс]. – Режим доступа: https://github.com/azizka85/oil_trans_in_caspian_sea (дата обращения: 01.06.2026).