GPU-based parallel simulation of wind-driven sea currents for oil spill forecasting

Full Text

Введение

Каспийский регион является зоной активной добычи углеводородов, где оперативный мониторинг течений имеет критическое значение для экологической безопасности. К наиболее важным морским месторождениям относятся Кашаган, Кайран и Актоты, входящие в Северо-Каспийский проект. Разведка, добыча и транспортировка нефти в морской акватории неизбежно сопровождаются риском аварийных выбросов и загрязнения морской среды [1, 2]. По этой причине требуются инструменты, позволяющие заранее прогнозировать направление распространения нефти, скорость её перемещения, площадь загрязнения и потенциальное воздействие на береговую линию и морскую экосистему. Если спутниковое зондирование обеспечивает фиксацию имеющегося нефтяного пятна, то численное моделирование выступает его необходимым прогностическим дополнением, позволяющим непрерывно рассчитывать дальнейший перенос и трансформацию разлива.

Одним из ключевых этапов моделирования распространения разлитой нефти по поверхности моря является расчёт морских течений, вызванных ветровым воздействием. Для решения этой задачи были разработаны различные модели, различающиеся по уровню сложности. В работах [3–5] была предложена гидродинамическая модель однородной жидкости с упрощенным представлением турбулентности, основанным на постоянной турбулентной вязкости. Позже в исследованиях [6–8] этот подход был расширен за счёт включения эффектов стратификации плотности и температуры, а также адвекции, что сделало модель нелинейной, однако процедура решения осталась по существу той же, что и в работе [3].

В работах [9, 10] детально рассматривалось взаимодействие морского льда и атмосферы с морем с целью исследования сезонной циркуляции и изменчивости водных масс Каспийского моря под воздействием реалистичных потоков массы, импульса и тепла. В этих моделях речной сток, осадки, радиационный поток и ветровое напряжение задавались с учётом сезонной изменчивости на основе доступных наборов данных. Кроме того, интенсивность испарения, поток явного и скрытого тепла на поверхности моря рассчитывались в ходе моделирования с использованием подмодели атмосферного пограничного слоя на основе атмосферных данных реанализа ECMWF ERA-15¹ и рассчитанной модели температуры поверхности моря.

Во всех вышеупомянутых работах использовались упрощённые модели турбулентности, в которых турбулентная вязкость задавалась параметрически. Вместе с тем существуют исследования, в которых применяются более совершенные модели турбулентности. Например, в работах [11–13] использовалась модель рейнольдсовых напряжений, требующая решения значительно большего числа уравнений. Однако, поскольку зависимости для определения турбулентной вязкости хорошо апробированы для широкого диапазона условий, на практике зачастую целесообразнее применять простые модели турбулентной вязкости, а не более сложные подходы, такие как модели рейнольдсовых напряжений.

Кроме того, все вышеупомянутые модели являются трёхмерными и поэтому требуют значительных вычислительных затрат и ресурсов [14, 15]. Например, в работе [15] расчёты выполнялись на сетке размером 426 × 238 × 100 узлов, что соответствует более чем 10 млн расчётных точек. Как правило, такие задачи решаются на вычислительных кластерах с большим числом процессоров и использованием параллельных алгоритмов. В настоящем исследовании представлены параллельные алгоритмы решения вышеупомянутых моделей на графических процессорах с применением библиотеки OpenCL [16]. Результаты вычислений показывают, что может быть достигнуто значительное ускорение, поскольку данные модели хорошо подходят для распараллеливания.

Научная новизна и практическая значимость работы заключаются в разработке параллельного OpenCL-алгоритма трёхмерной гидродинамики для графических процессоров (GPU) на вертикально-неравномерной σ-сетке (N_z=40). В отличие от большинства существующих гидродинамических комплексов [14–15], ориентированных на использование дорогостоящих и труднодоступных суперкомпьютерных кластеров (Message Passing Interface / Open Multi-Processing), предложенный подход позволяет выполнять оперативные 3D-расчёты динамики вод на относительно недорогом и доступном оборудовании (включая персональные рабочие станции с GPU).

Цель настоящей работы – разработать и реализовать параллельные алгоритмы для трёхмерного моделирования ветровых морских течений на графических процессорах с использованием OpenCL, обеспечивающие существенное сокращение времени расчётов при сохранении приемлемой точности для задач оперативного прогнозирования нефтяных разливов.

Материалы и методы

Согласно работе [3], линейные трёхмерные уравнения гидродинамики в декартовой системе координат записываются в следующем виде:

$\frac{\partial \zeta }{\partial t}=-\frac{\partial \left(h\overline{u}\right)}{\partial x}-\frac{\partial \left(h\overline{v}\right)}{\partial y}$ (1)

$\frac{\partial \overline{u}}{\partial t}=f\overline{v}-g\frac{\partial \zeta }{\partial x}-\frac{{\tau }_x^h}{\rho h}+\frac{{\tau }_x^s}{\rho h}$ (2)

$\frac{\partial \overline{v}}{\partial t}=-f\overline{u}-g\frac{\partial \zeta }{\partial y}-\frac{{\tau }_y^h}{\rho h}+\frac{{\tau }_y^s}{\rho h}$ (3)

где t – время, сек; x, y, z образуют левостороннюю декартову систему координат, координата z представляет собой глубину ниже невозмущённой поверхности; h – глубина воды, м; ζ – отметка свободной поверхности, метры; ӣ, $\overline{v}$ обозначают усреднённые по глубине компоненты скорости течения в направлениях x и y, соответственно:

$\overline{u}=\underset{0}{\overset{1}{}}ud\sigma , \overline{v}=\underset{0}{\overset{1}{}}vd\sigma$ (4)

где – безразмерная величина (5):

$\sigma =z/h$ (5)

Основные положения, рассмотренные выше, представлены на рис. 1:

 

Рисунок 1. Система координат сигма

Figure 1. Sigma coordinate system

 

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

$u=\overline{u}+u\text{'},v=\overline{v}+v\text{'}$ (6)

$\frac{\partial u\text{'}}{\partial t}=fv\text{'}+\frac{1}{h^2}\frac{\partial }{\partial \sigma }\left(\mu \frac{\partial u\text{'}}{\partial \sigma }\right)+\frac{{\tau }_x^h}{\rho h}-\frac{{\tau }_x^s}{\rho h}$ (7)

$\frac{\partial v\text{'}}{\partial t}=-fu\text{'}+\frac{1}{h^2}\frac{\partial }{\partial \sigma }\left(\mu \frac{\partial v\text{'}}{\partial \sigma }\right)+\frac{{\tau }_y^h}{\rho h}-\frac{{\tau }_y^s}{\rho h}$ (8)

где μ – коэффициент вертикальной вихревой кинематической вязкости, м²/с; g – ускорение свободного падения, м²/с; – параметр Кориолиса, с⁻¹; ρ – плотность, кг/м³; ${\tau }_x^s$ и ${\tau }_y^s$ – компоненты ветровой нагрузки, Н/м² (9):

${\tau }_x^s=-\frac{\rho }{h}{\left(\mu \frac{\partial u}{\partial \sigma }\right)}_0, {\tau }_y^s=-\frac{\rho }{h}{\left(\mu \frac{\partial v}{\partial \sigma }\right)}_0$ (9)

Аналогично на морском дне σ, ${\tau }_x^h$ и ${\tau }_y^h$ обозначают компоненты придонного напряжения:

${\tau }_x^h=-\frac{\rho }{h}{\left(\mu \frac{\partial u}{\partial \sigma }\right)}_1, {\tau }_y^h=-\frac{\rho }{h}{\left(\mu \frac{\partial v}{\partial \sigma }\right)}_1$ (10)

В описанной далее модели компоненты ветрового напряжения ${\tau }_x^s$ и ${\tau }_y^s$ предопределены и был использован линейный закон донного напряжения:

${\tau }_x^h=k\rho u_h,{\tau }_y^h=k\rho v_h$ (11)

где k – коэффициент трения дна; u_h и v_h – компоненты течения на глубине h.

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

$\frac{\rho }{h}{\left(\mu \frac{\partial u^{\text{'}}}{\partial \sigma }\right)}_0=-{\tau }_x^s, \frac{\rho }{h}{\left(\mu \frac{\partial v\text{'}}{\partial \sigma }\right)}_0=-{\tau }_y^s$ (12)

где компоненты ветрового напряжения ${\tau }_x^s$ и ${\tau }_y^s$ задаются заранее.

На морском дне при σ = 0 соответствующие граничные условия имеют вид (13):

${\tau }_x^h=-\frac{\rho }{h}{\left(\mu \frac{\partial u^{\text{'}}}{\partial \sigma }\right)}_1 =k\rho u_h, {\tau }_y^h=-\frac{\rho }{h}{\left(\mu \frac{\partial v\text{'}}{\partial \sigma }\right)}_1=k\rho v_h$ (13)

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

$z\left(\mathrm{0,}y\right)=z\left(L,y\right),z\left(x\mathrm{,0}\right)=z\left(x,W\right)$ (14)

$\overline{u}\left(\mathrm{0,}y\right)=\overline{u}\left(L,y\right),\overline{u}\left(x\mathrm{,0}\right)=\overline{u}\left(x,W\right)$ (15)

$\overline{v}\left(\mathrm{0,}y\right)=\overline{v}\left(L,y\right),\overline{v}\left(x\mathrm{,0}\right)=\overline{v}\left(x,W\right)$ (16)

где L и W – длина и ширина расчётной области соответственно, ориентированной вдоль пространственных координат x и y соответственно.

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

 

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

Figure 2. Staggered grid

 

В данном разделе приводится краткое описание конечно-разностной сетки, используемой в модели.

Уравнения (1)–(3) дискретизируются следующим образом (17)–(19):

${\overline{u}}_{i,j}^{n+1}={\overline{u}}_{i,j}^n+f{\overline{v}}_{i,j}^n△t-g△t\frac{\partial {\zeta }_{i,j}^n}{\partial x}-\frac{{\tau }_x^h}{\rho h}△t+\frac{{\tau }_x^s}{\rho h}△t$ (17)

${\overline{v}}_{i,j}^{n+1}={\overline{v}}_{i,j}^n-f{\overline{u}}_{i,j}^n△t-g△t\frac{\partial {\zeta }_{i,j}^n}{\partial y}-\frac{{\tau }_y^h}{\rho h}△t+\frac{{\tau }_y^s}{\rho h}△t$ (18)

${\zeta }_{i,j}^{n+1}={\zeta }_{i,j}^n-\frac{\partial {\left(h\overline{u}\right)}_{i,j}^{n+1}}{\partial x}△t-\frac{\partial {\left(h\overline{v}\right)}_{i,j}^{n+1}}{\partial y}△t$ (19)

Подставляя значения ${\overline{u}}_{i,j}^{n+1}$ и ${\overline{v}}_{i,j}^{n+1}$, определяемые выражениями (17)–(18), в уравнение (19), получаем:

${\zeta }_{i,j}^{n+1}={\zeta }_{i,j}^n+gh△t^2\left[\frac{{\partial }^2{\zeta }_{i,j}^n}{\partial x^2}+\frac{{\partial }^2{\zeta }_{i,j}^n}{\partial y^2}\right]+S_{i,j}$ (20)

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

${\tau }_g<\frac{△s}{\sqrt{2gh}},$ $△s=min\left(△x,△y\right)$ (21)

При глубине моря h=100 м и шаге сетки ∆s=1000 м, из уравнения (21) следует, что ${\tau }_g~1\mathrm{мин}$.

Подставляя значение ${\zeta }_{i,j}^{n+1}$ определяемое выражением (20), в уравнение (17), получаем (22):

${\overline{u}}_{i,j}^{n+1}={\overline{u}}_{i,j}^n-f^2△t^2{\overline{u}}_{i,j}^n+R_{i,j}$ (22)

Шаг по времени ∆t=τ_f в уравнении (22) определяется из условия устойчивости при $f=1.2\times {10}^{-4} {с}^{-1}$ (23):

${\tau }_f<\frac{\sqrt{2}}{f}\approx 1\mathrm{ч}$ (23)

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

$\begin{array}{l}{u}_{i,j,k}^{\text{'}n+1}=u_{i,j,k}^{\text{'}n}+f{v}_{i,j,k}^{\text{'}n}\Delta t+\\ +\frac{\Delta t}{h^2\Delta {\sigma }_k}\left[\frac{{\mu }_{i,j,k+1}^n\left(u_{i,j,k+1}^{\text{'}n+1}-u_{i,j,k}^{\text{'}n+1}\right)}{\Delta {\sigma }_{k+1}+\Delta {\sigma }_k}-\frac{{\mu }_{i,j,k}^n\left(u_{i,j,k}^{\text{'}n+1}-u_{i,j,k-1}^{\text{'}n+1}\right)}{\Delta {\sigma }_k+\Delta {\sigma }_{k-1}}\right]+\\ +\frac{\Delta t}{h^2\Delta {\sigma }_k}\left[\frac{{\mu }_{i,j,k+1}^n\left(u_{i,j,k+1}^{\text{'}n}-u_{i,j,k}^{\text{'}n}\right)}{\Delta {\sigma }_{k+1}+\Delta {\sigma }_k}-\frac{{\mu }_{i,j,k}^n\left(u_{i,j,k}^{\text{'}n}-u_{i,j,k-1}^{\text{'}n}\right)}{\Delta {\sigma }_k+\Delta {\sigma }_{k-1}}\right]+\\ +\frac{k{u}_{i,j,N_Z-1}^n}{h_{i,j}}\Delta t-\frac{{\tau }_{x,i,j}^{s,n}}{\rho h_{i,j}}\Delta t\end{array}$ (24)

$\begin{array}{l}{v}_{i,j,k}^{\text{'}n+1}=v_{i,j,k}^{\text{'}n}-f{u}_{i,j,k}^{\text{'}n}\Delta t+\\ +\frac{\Delta t}{h^2\Delta {\sigma }_k}\left[\frac{{\mu }_{i,j,k+1}^n\left(v_{i,j,k+1}^{\text{'}n+1}-v_{i,j,k}^{\text{'}n+1}\right)}{\Delta {\sigma }_{k+1}+\Delta {\sigma }_k}-\frac{{\mu }_{i,j,k}^n\left(v_{i,j,k}^{\text{'}n+1}-v_{i,j,k-1}^{\text{'}n+1}\right)}{\Delta {\sigma }_k+\Delta {\sigma }_{k-1}}\right]+\\ +\frac{\Delta t}{h^2\Delta {\sigma }_k}\left[\frac{{\mu }_{i,j,k+1}^n\left(v_{i,j,k+1}^{\text{'}n}-v_{i,j,k}^{\text{'}n}\right)}{\Delta {\sigma }_{k+1}+\Delta {\sigma }_k}-\frac{{\mu }_{i,j,k}^n\left(v_{i,j,k}^{\text{'}n}-v_{i,j,k-1}^{\text{'}n}\right)}{\Delta {\sigma }_k+\Delta {\sigma }_{k-1}}\right]+\\ +\frac{k{v}_{i,j,N_Z-1}^n}{h_{i,j}}\Delta t-\frac{{\tau }_{y,i,j}^{s,n}}{\rho h_{i,j}}\Delta t\end{array}$ (25)

Уравнения (24)–(25) приводятся к трёхдиагональной системе линейных алгебраических уравнений (далее – СЛАУ) следующим образом (26)–(27):

$\begin{array}{l}-\frac{{\mu }_{i,j,k+1}^n}{\Delta {\sigma }_{k+1}+\Delta {\sigma }_k}\frac{\Delta t}{h^2\Delta {\sigma }_k}{u}_{i,j,k+1}^{\text{'}n+1}+\left[1+\frac{\Delta t}{h^2\Delta {\sigma }_k}\left(\frac{{\mu }_{i,j,k+1}^n}{\Delta {\sigma }_{k+1}+\Delta {\sigma }_k}+\frac{{\mu }_{i,j,k}^n}{\Delta {\sigma }_k+\Delta {\sigma }_{k-1}}\right)\right]u_{i,j,k}^{\text{'}n+1}-\frac{{\mu }_{i,j,k}^n}{\Delta {\sigma }_k+\Delta {\sigma }_{k-1}}\frac{\Delta t}{h^2\Delta {\sigma }_k}{u}_{i,j,k-1}^{\text{'}n+1}=\\ =u_{i,j,k}^{\text{'}n}+f{v}_{i,j,k}^{\text{'}n}\Delta t+\frac{\Delta t}{h^2\Delta {\sigma }_k}\left[\frac{{\mu }_{i,j,k+1}^n\left(u_{i,j,k+1}^{\text{'}n}-u_{i,j,k}^{\text{'}n}\right)}{\Delta {\sigma }_{k+1}+\Delta {\sigma }_k}-\frac{{\mu }_{i,j,k}^n\left(u_{i,j,k}^{\text{'}n}-u_{i,j,k-1}^{\text{'}n}\right)}{\Delta {\sigma }_k+\Delta {\sigma }_{k-1}}\right]+\frac{k{u}_{i,j,N_z-1}^n}{h_{i,j}}\Delta t-\frac{{\tau }_{x,i,j}^{s,n}}{\rho h_{i,j}}\Delta t\end{array}$ (26)

$\begin{array}{l}-\frac{{\mu }_{i,j,k+1}^n}{\Delta {\sigma }_{k+1}+\Delta {\sigma }_k}\frac{\Delta t}{h^2\Delta {\sigma }_k}{v}_{i,j,k+1}^{\text{'}n+1}+\left[1+\frac{\Delta t}{h^2\Delta {\sigma }_k}\left(\frac{{\mu }_{i,j,k+1}^n}{\Delta {\sigma }_{k+1}+\Delta {\sigma }_k}+\frac{{\mu }_{i,j,k}^n}{\Delta {\sigma }_k+\Delta {\sigma }_{k-1}}\right)\right]v_{i,j,k}^{\text{'}n+1}-\frac{{\mu }_{i,j,k}^n}{\Delta {\sigma }_k+\Delta {\sigma }_{k-1}}\frac{\Delta t}{h^2\Delta {\sigma }_k}{v}_{i,j,k-1}^{\text{'}n+1}=\\ =v_{i,j,k}^{\text{'}n}-f{u}_{i,j,k}^{\text{'}n}\Delta t+\frac{\Delta t}{h^2\Delta {\sigma }_k}\left[\frac{{\mu }_{i,j,k+1}^n\left(v_{i,j,k+1}^{\text{'}n}-v_{i,j,k}^{\text{'}n}\right)}{\Delta {\sigma }_{k+1}+\Delta {\sigma }_k}-\frac{{\mu }_{i,j,k}^n\left(v_{i,j,k}^{\text{'}n}-v_{i,j,k-1}^{\text{'}n}\right)}{\Delta {\sigma }_k+\Delta {\sigma }_{k-1}}\right]+\frac{k{v}_{i,j,N_z-1}^n}{h_{i,j}}\Delta t-\frac{{\tau }_{y,i,j}^{s,n}}{\rho h_{i,j}}\Delta t\end{array}$ (27)

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

На рис. 3 представлена блок-схема параллельного алгоритма, отражающая последовательность вычислительных этапов и распределение задач между управляющим процессором (CPU) и графическим ускорителем (GPU). Вычислительное пространство OpenCL организуется в виде двумерной глобальной сетки потоков размером N_x × N_y, где каждый рабочий элемент (work-item) с индексами (i, j) непрерывно обрабатывает отдельную вертикальную колонку водных слоёв по координате z. Внутри каждого потока последовательно выполняются расчёт адресации в памяти GPU, формирование коэффициентов трёхдиагональной матрицы для неявных уравнений (24)–(25), а также прямой и обратный ход метода прогонки по глубине. Данная организация вычислений полностью устраняет необходимость межпотоковой синхронизации по вертикальной оси и обеспечивает эффективную прямую запись обновлённых компонент скорости u и v в глобальную память видеокарты.

 

Рисунок 3. Блок-схема параллельного алгоритма расчёта ветровых течений на GPU с декомпозицией по координатам x и y

Figure 3. Flowchart of the parallel GPU algorithm for wind-driven current simulation with decomposition along the x and ycoordinates

 

Измерение времени работы алгоритмов проводилось с использованием высокоточных системных таймеров (std::chrono::high_resolution_clock). Для исключения эффектов инициализации библиотеки OpenCL, первичной компиляции кернелов и передач данных («разогрев» GPU/CPU) первый запуск каждого сценария не учитывался. Для каждой конфигурации выполнялось по 5 независимых запусков, после чего рассчитывалось среднее значение времени выполнения. Относительное стандартное отклонение (разброс) результатов замеров во всех экспериментах не превышало 1,5%, что говорит о высокой воспроизводимости измерений.

Исходный код последовательного и параллельного алгоритмов доступен в работе [17].

Результаты

Численные расчёты проводились при следующих значениях параметров: $f=1.2\times {10}^{-4}{\mathrm{с}}^{-1}$, $g=9.81\mathrm{м}/{\mathrm{с}}^{-2}$, $\rho =1025\mathrm{кг}/{\mathrm{м}}^{-3}$, $k=0.002{\mathrm{мс}}^{-1}$, $\mu =0.4{\mathrm{м}}^2{\mathrm{с}}^{-1}$, H = 260 м, L = 260 м, W = 260 м.

Течения генерировались ветровым напряжением, направленным под углом 45^∘ к осям x и y, при этом значения его компонентов составляли ${\tau }_x^s=1.5\mathrm{Н}{\mathrm{м}}^{-2}$ и ${\tau }_y^s=1.5\mathrm{Н}{\mathrm{м}}^{-2}$ соответственно.

Для оценки точности расчёта течений при различных шагах сетки использовались два варианта распределения расчётных узлов. В первом случае применялась равномерная сетка, содержащая N_z = 500 узлов вдоль координаты z. Во втором случае использовалось N_z = 40 узлов по координате с повышенным разрешением в поверхностном пограничном слое (σ = 0.002). По координатам x и y, применялась равномерная сетка с числом узлов N_x = 20 и N_y = 20 соответственно.

Все расчёты выполнялись на центральном процессоре Intel Core i7 с тактовой частотой 3,6 ГГц и графическом процессоре NVIDIA GeForce GTX 670, содержащем 1344 вычислительных ядра и работающем на частоте 915 МГц.

Сравнение профилей скорости, полученных на указанных двух типах сеток, представлено на рис. 4–5.

Для количественной оценки точности решения на неравномерной сетке (N_z = 40) по сравнению с эталонным расчётом на равномерной сетке (N_z = 500) были рассчитаны среднеквадратичная относительная ошибка (28) и максимальная относительная ошибка (29):

$\text{RMSE}=\frac{\sqrt{\frac{1}{N}{{\displaystyle \sum }}_{k=1}^N{\left(U_k^{\left(40\right)}-U_k^{\left(500\right)}\right)}^2}}{\max \left(U^{\left(500\right)}\right)-\min \left(U^{\left(500\right)}\right)}\times 100\%$ (28)

${\text{E}}_{\text{max}}=\frac{\underset{k}{\max }\left|U_k^{\left(40\right)}-U_k^{\left(500\right)}\right|}{\max \left(U^{\left(500\right)}\right)-\min \left(U^{\left(500\right)}\right)}\times 100\%$ (29)

При расчёте этих погрешностей значения профилей скорости на неравномерной сетке (N_z = 40) были сопоставлены со значениями на равномерной сетке (N_z = 500) с использованием линейной интерполяции в соответствующих узлах по глубине. Результаты сравнения для всех шести рассмотренных профилей скорости представлены в табл. 1.

 

Таблица 1. Количественная оценка расхождения профилей скорости (u, v) на неравномерной (N_z = 40) и эталонной равномерной (N_z = 500) сетках в разные моменты времени

Table 1. Quantitative assessment of the discrepancy between velocity profiles (u, v) on non-uniform uniform (N_z = 40) and reference uniform (N_z = 500) grids at different time points

Рисунок Figure	Момент времени Time	Компонента Component	Среднеквад. относит. ошибка, % Root mean square relative error, %	Макс. относит. ошибка, % Maximum relative error, %
Рис. 4 Fig 4	1 ч 40 мин 1 h 40 min	u	0,31	0,76
Рис. 5 Fig 5	1 ч 40 мин 1 h 40 min	v	0,52	1,23
Рис. 4 Fig 4	6 ч 40 мин 6 h 40 min	u	0,32	0,67
Рис. 5 Fig 5	6 ч 40 мин 6 h 40 min	v	0,3	0,79
Рис. 4 Fig 4	15 ч 15 h	u	0,25	0,56
Рис. 5 Fig 5	15 ч 15 h	v	0,58	1,33

 

Видно, что точность расчётов на неравномерной сетке с N_z = 40 сопоставима с точностью, достигаемой на равномерной сетке с N_z = 500. Таким образом, использование неравномерной сетки с N_z = 40 является более целесообразным по сравнению с равномерной сеткой с N_z = 500.

 

Рисунок 4. Профили компоненты скорости u (м/с) по глубине в точке x = 130 м, y = 130 м

Figure 4. Depth profiles of the velocity component u (m/s) at the point x = 130 m, y = 130 m

 

Рисунок 5. Профили компоненты скорости v (м/с) по глубине в точке x = 130 м, y = 130 м

Figure 5. Depth profiles of the velocity component v (m/s) at the point x = 130 m, y = 130 m

На рис. 4–5 красная линия – неравномерная сетка (N_z = 40, t = 1 ч 40 мин), синие точки – эталонная равномерная сетка (N_z = 500, t = 1 ч 40 мин), фиолетовые штрихи – неравномерная сетка (N_z = 40, t = 6 ч 40 мин), зеленые квадраты – эталонная равномерная сетка (N_z = 500, t = 6 ч 40 мин), черные штрихи с точками – неравномерная сетка (N_z = 40, t = 15 ч), коричневые ромбы – эталонная равномерная сетка (N_z = 500, t = 15 ч).

Figure 4–5. The red line represents the non-uniform grid (N_z = 40, t = 1 h 40 min), the blue dots represent the reference uniform grid (N_z = 500, t = 1 h 40 min), the purple dashes represent the non-uniform grid (N_z = 40, t = 6 h 40 min), the green squares represent the reference uniform grid (N_z = 500, t = 1 h 40 min), the black dash-dotted line represents the non-uniform grid (N_z = 40, t = 15 h), and the brown diamonds represent the reference uniform grid (N_z = 500, t = 15 h).

 

На рис. 6 представлено сравнение времени решения задачи на равномерной сетке с при использовании CPU и GPU. Для данной сетки вычисления на GPU выполнялись в 23 раза быстрее, чем на CPU.

 

Рисунок 6. Зависимость времени вычислений (сек) от количества итераций на равномерной сетке (N_z = 500)

Figure 6. Dependence of computation time (s) on the number of iterations for the uniform grid (N_z = 500)

 

На рис. 7 представлено сравнение времени решения задачи на неравномерной сетке с при использовании CPU и GPU. Для данной сетки вычисления на GPU выполнялись в 13 раз быстрее, чем на CPU.

 

Рисунок 7. Зависимость времени вычислений (сек) от количества итераций на неравномерной сетке (N_z = 40)

Figure 7. Dependence of computation time (s) on the number of iterations for the non-uniform grid (N_z = 40)

На рис. 6–7 сплошная красная линия – CPU (Intel Core i7), штриховая синяя линия – GPU (NVIDIA GeForce GTX 670).

In Figs. 6–7, the solid red line represents the CPU (Intel Core i7), while the dashed blue line represents the GPU (NVIDIA GeForce GTX 670)

 

Обсуждение

Полученное ускорение (23× для равномерной сетки, 13× для неравномерной) согласуется с известными результатами применения GPU в гидродинамике, где типичный выигрыш составляет 10–50× в зависимости от размера сетки и архитектуры. Важно отметить, что для неравномерной сетки с 40 слоями ускорение ниже, но это компенсируется меньшим объёмом вычислений и сохранением точности в поверхностном слое, что критично для моделирования разливов.

Различие в коэффициентах ускорения на сетках N_z = 500 (23 x) и N_z = 40 (13 x) обусловлено балансом между объёмом вычислений и фиксированными накладными расходами параллельной среды. При решении СЛАУ методом прогонки вдоль вертикали количество арифметических операций на один рабочий элемент напрямую пропорционально N_z. Количественный анализ структуры времени выполнения показывает следующий баланс составляющих:

На мелкой сетке (N_z = 500) чистые вычисления на GPU занимают около 97% времени шага, а накладные расходы OpenCL (вызов кернелов, синхронизация и передача граничных условий) – не более 3%. Вычислительные блоки GPU работают с высокой степенью загрузки, эффективно маскируя задержки обращения к глобальной памяти (latency hiding);

На редкой сетке (N_z = 40) объём вычислений в одном потоке сокращается в 12,5 раз, в результате чего доля неизменных накладных расходов на управление и передачу данных возрастает до ≈ 40% от общего времени шага. Кроме того, малая длина цикла прогонки (N_z = 40) не позволяет полностью утилизировать пропускную способность видеопамяти. В результате фиксированные задержки интерфейса OpenCL и шины PCIe оказывают определяющее влияние на итоговое время, снижая относительный выигрыш GPU по сравнению с CPU.

Использование OpenCL обеспечивает переносимость на различные GPU и даже на кластеры, что расширяет область применения. Однако явный шаг по времени (ограниченный условием Куранта) остаётся узким местом. В перспективе целесообразно перейти к полностью неявным схемам, позволяющим увеличить шаг до часов, что ещё больше повысит эффективность.

Следует отметить, что тестовые расчёты проводились в прямоугольной области с плоским дном и постоянным ветровым воздействием. Использование такого идеализированного сценария является общепринятым подходом для верификации численных схем и измерения чистого ускорения (speedup) на GPU. Переход к реальным физико-географическим условиям (батиметрия Каспийского моря, береговая линия и переменное ветровое поле) не снизит эффективности параллельного алгоритма. В рамках используемой -системы координат переменная глубина H(x, y) учитывается через коэффициенты и не меняет структуру трёхдиагональных матриц. Учёт береговой линии посредством маскирования суши (land masking) сокращает количество активных расчётных узлов, что положительно скажется на времени счёта, а обновление ветра требует незначительных вычислительных затрат уровня O(N_xN_y) по сравнению с 3D-блоком O(N_xN_yN_z).

Заключение

В данной работе численное моделирование ветровых течений в море проводилось на равномерной и неравномерной сетках с использованием CPU и GPU. Поскольку глубина моря относительно невелика – около 100 м – по сравнению с горизонтальными размерами расчётной области по осям x и y, число вычислительных узлов по оси z также меньше, чем по направлениям x и y. Как подтверждают полученные результаты, для решения трёхдиагональных СЛАУ более эффективным является применение последовательного алгоритма, а не параллельного. Кроме того, результаты показали, что ускорение вычислений было выше для равномерной сетки с чем для неравномерной сетки с .

При решении более сложных моделей возникает необходимость применения неявных схем, позволяющих использовать более крупные шаги по времени. Полученные неявные разностные уравнения могут решаться как прямыми методами, например, посредством решения трёхдиагональных СЛАУ, использованного в настоящей работе, так и итерационными методами. Преимущество прямых методов заключается в их вычислительной эффективности при любых начальных и граничных условиях.

Оценка времени вычислений (менее 1 ч на 24-часовой прогноз) получена для базового модельного сценария и характеризует высокую вычислительную эффективность разработанного параллельного ядра. Это демонстрирует перспективность использования предложенной схемы в составе оперативных систем гидродинамического прогнозирования после её интеграции с модулями реальной батиметрии и пространственно-временной изменчивости метеополей.

ДОПОЛНИТЕЛЬНО

Источник финансирования. Исследование выполнено при финансовой поддержке Министерства науки и высшего образования Республики Казахстан (грант № BR24992964).

Конфликт интересов. Авторы декларируют отсутствие явных и потенциальных конфликтов интересов, связанных с публикацией настоящей статьи.

Благодарности. Авторы выражают искреннюю благодарность Министерству науки и высшего образования Республики Казахстан за поддержку исследования, выполненного в рамках проекта BR24992964 «Разработка интегрированных энергосберегающих технологий для обеспечения экологической устойчивости и эффективной морской деятельности в казахстанском секторе Каспийского моря».

Вклад авторов. Все авторы подтверждают соответствие своего авторства международным критериям ICMJE (все авторы внесли существенный вклад в разработку концепции, проведение исследования и подготовку статьи, прочли и одобрили финальную версию перед публикацией). Наибольший вклад распределён следующим образом: Кудайкулов А.А. – концепция и дизайн исследования, написание черновика; Аманжолов Т.Е. – разработка численных алгоритмов, программирование на OpenCL; Алиулы А. – проведение расчётов, анализ данных; Сейтов А.Н. – подготовка графических материалов; Асылбеков Б.К. – консультации по гидродинамике; Кульджабеков А.Б. – руководство проектом, редактирование; Шабилов С.М. – верификация результатов; Байымбетов Д.Б. – адаптация к задачам нефтяных разливов; Сырлыбеккызы С. – экологическая интерпретация; Калтаев А.Ж. – литературный обзор.

ADDITIONAL INFORMATION

Funding source. This research was funded by the Ministry of Science and Higher Education of the Republic of Kazakhstan (Grant No. BR24992964).

Competing interests. The authors declare that they have no competing interests.

Acknowledgements. The authors express their sincere gratitude to the Ministry of Science and Higher Education of the Republic of Kazakhstan for supporting this research, which was conducted within the framework of Project BR24992964, “Development of Integrated Energy-Saving Technologies to Ensure Environmental Sustainability and Efficient Marine Activities in the Kazakh Sector of the Caspian Sea.”

Authors’ contribution. All authors made a substantial contribution to the conception of the work, acquisition, analysis, interpretation of data for the work, drafting and revising the work, final approval of the version to be published and agree to be accountable for all aspects of the work. The greatest contribution is distributed as follows: Aziz A. Kudaikulov – conceptualization and study design, writing of the original draft; Tannur Ye. Amanzholov – development of numerical algorithms, OpenCL programming; Abdurashid Aliuly – conducting calculations, data analysis; Abzal N. Seitov – preparation of graphical materials; Bakytzhan K. Assilbekov – consultations on hydrodynamics; Alibek B. Kuldzhabekov – project leadership, editing; Spartak M. Shabilov – validation of results; Dinmukhammed B. Baimbetov – adaptation to oil spill forecasting applications; Samal Syrlybekkyzy – environmental interpretation; Aidarkhan Zh. Kaltayev – literature review.

 

¹ ECMWF ERA-15 – это первый проект глобального реанализа атмосферы, который выполнил ECMWF (англ. European Centre for Medium-Range Weather Forecasts – Европейский центр среднесрочных прогнозов погоды) для сбора и уточнения исторических метеоданных за период с декабря 1978 г. по февраль 1994 г.

About the authors

Aziz A. Kudaikulov

Satbayev University; Al-Farabi Kazakh National University

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

PhD

Kazakhstan, Almaty; Almaty

Tannur Ye. Amanzholov

Satbayev University; Al-Farabi Kazakh National University

Author for correspondence.
Email: tannur.amanzholov@gmail.com
ORCID iD: 0000-0001-7386-9717
Kazakhstan, Almaty; Almaty

Abdurashid Aliuly

Satbayev University

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

Abzal N. Seitov

Satbayev University; Al-Farabi Kazakh National University

Email: seitovabzal9@gmail.com
ORCID iD: 0000-0002-7317-3047
Kazakhstan, Almaty; Almaty

Bakytzhan K. Assilbekov

Ionosphere Institute

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

PhD

Kazakhstan, Almaty

Alibek B. Kuldzhabekov

Tashenev University

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

PhD

Kazakhstan, Shymkent

Spartak M. Shabilov

Satbayev University

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

Dinmukhambet B. Baimbetov

Yessenov University

Email: d2.baimbetov@yu.edu.kz
ORCID iD: 0000-0002-4992-0890

PhD

Kazakhstan, Aktau

Samal Syrlybekkyzy

Yessenov University

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

PhD, Associate Professor

Kazakhstan, Aktau

Aidarkhan Zh. Kaltayev

Satbayev University

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

D. Sc. (Physics & Mathematics), Professor

Kazakhstan, Almaty

References

Supplementary files