CC BY
50
УДК 551.247.1:519.632.4:004.272.26
МОДЕЛИРОВАНИЕ СОЛЯНОГО ДИАПИРИЗМА
РАСЧЕТОМ 3D ПОЛЗУЩИХ ТЕЧЕНИЙ С ИСПОЛЬЗОВАНИЕМ
ТЕХНОЛОГИИ ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ CUDA НА GPU
Борис Валентинович Лунев
Институт нефтегазовой геологии и геофизики им. ак. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Ак. Коптюга, 3, кандидат физико-математических наук, старший научный сотрудник, e-mail: bobvalmail@mail.ru
Тимофей Владимирович Абрамов
Институт нефтегазовой геологии и геофизики им. ак. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Ак. Коптюга, 3, аспирант, младший научный сотрудник, email: abramovtv@ipgg.sbras.ru.
Реализован алгоритм параллельных вычислений на GPU с помощью технологии CUDA для 3D моделирования процесса соляного диапиризма, представляемого ползущим течением однородно-вязкой ньютоновской жидкости с переменной плотностью под действием силы тяжести. Высокая эффективность алгоритма обусловлена использованием функции Грина для задачи о полупространстве со свободной поверхностью. Благодаря этому удается значительно сократить количество операций и эффективно организовать параллельные вычисления аналогию с расчетом поля в задаче взаимодействия N тел.
Ключевые слова: численное моделирование, ползущие течения, соляной диапиризм, параллельные вычисления, GPU, CUDA.
MODELING OF SALT DIAPIRISM BY CALCULATION OF 3D CREEPING FLOW WITH PARALLEL COMPUTING ON GPU WITH CUDA
Boris V. Lunev
A. A. Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, 3 Koptuga, Senior Research Fellow, Candidate of Physics and Mathematics, email: bobvalmail@mail.ru
Timofey V. Abramov
A. A. Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, 3 Koptuga, postgraduate, Junior Researcher, e-mail: abramov@ipgg.sbras.ru
The highly efficient program for three-dimensional numerical modelling of salt diapirism is built due to analytical Green's function for Newtonian fluid in the halfspace with free boundary under the action of gravity. Computing is parallel and uses GPU resources with CUDA. Programming approach for calculating velocity field is the same to calculating acceleration field in N body simulation.
Key words: numerical modeling, creeping flow, salt diapirism, parallel computing, GPU, CUDA.
Процессы соляного диапиризма играют важную структурообразующую роль в большинстве нефтегазовых провинций, в связи с чем давно и
интенсивно изучаются, как в природе, так и методами физического и численного моделирования. Однако, возможности использования численных методов до сих пор были сильно ограничены большим временем, нужным для расчетов ползущих течений, моделирующих этот процесс, несмотря на использование параллельных вычислений. Как правило, расчет таких течений основан на разностных методах. В этом случае, несмотря на действия по приведению системы уравнений к удобному для параллельных вычислений виду, неизбежны частые обмены информацией между вычислительными узлами. В результате ускорение вычислений оказывается ограниченным пропускной способностью памяти, а не пиковой производительностью устройства, как показано на рис. 1 [1, 2].
8 12 16 Количество процессоров
а)
я>
40
30
а
20
10
—0—GMG
—О— AMC.
-О-FFT
У^/ухпх 2Ь0
8<Ю0 GTSЛу
8500 GT Л?
20 40 60 80 100 120 Г.ВЛСС
б)
Рис. 1. Зависимость ускорения параллельных вычислений при решении
гидродинамических задач:
а) от числа процессоров (по [1] - расчет неустойчивости Рэлея-Тейлора с использованием схемы Холесского); б) от пропускной способности памяти (по [2] - расчет неустойчивости в свободном сдвиговом слое с использованием геометрического многосеточного метода -GMG, алгебраического многосеточного метода - AMG и метода редукции - FFT)
Основные вычислительные трудности связаны с решением квазистационарных краевых задач, связанная последовательность которых описывает ползущее течение. Для полупространства однородно-вязкой ньютоновской жидкости со свободной границей решение такой краевой задачи удалось получить аналитически в виде функции Грина [3]. Отыскание квазистационарного течения в этом случае сводится к вычислению интеграла свертки:
где - скорость течения в точке {х}, о^) - возмущение плотности в точке
(f) и Vj^) - функция Грина. В работах [4, 5] обоснована корректность
представления среды однородно-вязкой ньютоновской жидкостью применительно к моделированию соляного тектогенеза и продемонстрирована эффективность такого подхода для расчета 2D течений, даже на одном процессоре.
Моделирование реальных геологических объектов такого типа, как правило, требует трехмерной реализации. Вычислительная сложность при этом возрастает значительно, т.к. в двумерном случае для сетки N X N элементов требуется 0(N4) операций, а для трехмерной сетки N X N X N -0(N6) операций. Необходимая оперативность моделирования в этом случае может быть достигнута организацией параллельных вычислений. В этом отношении использование функции Грина особенно эффективно: при вычислении интеграла свертки каждое из вычислительных устройств может рассчитывать значения поля от источников о^) в некотором выделенном для него множестве точек-«приемников» {я} без дополнительных обменов данными между устройствами. Таким образом ускорение от применения параллельной архитектуры оказывается пропорциональным пиковой производительности системы, а не пропускной способности памяти.
Дополнительное ускорение вычислений дает и вытекающая из (1) возможность расчета только скорости движения границ тел (слоев) с разной плотностью, вместо расчета поля скорости во всей области, обязательного при использовании разностных методов. Благодаря этому, в некотором смысле, вычисляется влияние трехмерного объекта (распределение плотности) на двумерный (поверхность), что снижает размерность задачи с 0(N6) до 0(Ns). При этом явное описание движения границ увеличивает точность общего расчета течения.
Заложенные в методе возможности удалось реализовать в программе параллельных вычислений, ориентированной на использование GPU. В разработанной программе для трехмерного моделирования применяется явное описание границ тел в виде триангулированных поверхностей. Сетка источников при этом на каждом временном шаге восстанавливается из конфигурации границ тел (сетка источников имеет фиксированные размеры и является более грубой аппроксимацией области). Такое описание, понижая сложность задачи, одновременно увеличивает точность расчетов, что позволяет вместо двойной обходиться одинарной точностью (тесты показали, что использование одинарной и двойной точности дает практически идентичные результаты вплоть до глубоких стадий эволюции). Возможность использования одинарной точности позволяет значительно ускорить расчеты при использовании современных GPU, для которых характерна большая разница в производительности вычислений с одинарной и двойной точностью. Для организации параллельных вычислений с производительностью близкой к пиковой для GPU, используется
надлежащим образом модифицированный метод решения задачи о взаимодействии многих тел [2].
Явное определение поверхностей, вместе с эффективной параллельной реализацией, показывает отличные результаты моделирования за приемлемое время. Расчеты производились с помощью устройства Tesla с1060, имеющего 240 вычислительных ядер СЦОА с тактовой частотой 1300 МГц и пиковую производительность 933 GFLOPS для вычислений с одинарной точностью и 78 GFLOPS для двойной.
В качестве примеров, ниже представлены рассчитанные зрелые стадии развития плотностной неустойчивости в среде, представленной несколькими слоями с разной плотностью.
Рис. 2 представляет результат эволюции инициированной одиночным начальным возмущением.
Рис. 3. Развитие одиночного возмущения. Вверху 3Б изображение поверхностей кровли и подошвы всплывающего слоя («соли»), внизу - разрез структуры через ось симметрии.
(Плотности всплывающего, перекрывающего и подстилающего слоев взяты 2.2 г/см3, 2.4 г/см3 и 2.55 г/см3, соответственно)
Результат эволюции более сложной модели, состоящей из 5 слоев с плотностями (снизу вверх) 2.65 г/см3, 2.55 г/см3, 2.2 г/см3 (соль), 2.4 г/см3 и 2.2 г/см3, с хаотичными начальными возмущениями представлен на рис. 4.
Во втором примере для задания источников использовалась сетка 128^128x35 узлов и расстояние между соседними точками, задающими поверхности, составляло величину порядка расстояния между узлами сетки. Расчет в этой модели занял два часа с помощью одного устройства Tesla с1060. Полученные результаты говорят о том, что на более мощных
вычислительных системах разработанная программа пригодна для оперативного расчета гораздо более сложных моделей, близко воспроизводящих реальные геологические объекты.
Рис. 4. Множественные возмущения.
Вверху - начальное состояние модели, внизу - конечная стадия эволюции.
Изображение поверхности, разделяющей верхние два слоя, снято
Работа выполнена в рамках программы VIII.73.2 фундаментальных научных исследований СО РАН.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Исмаил-заде. А. Т., Цепелев И. А., Тэлбот К., Остер П. Трехмерное моделирование соляного диапиризма: численный подход и алгоритм параллельных вычислений // Вычислительная сейсмология, вып 31. 2000, с. 62-76.
2. Боресков А. В., Харламов А. А. Основы работы с технологией CUDA. - М.: ДМК Пресс, 2010. - 232 с.
3. Лунёв Б.В. Изостазия как динамическое равновесие вязкой жидкости. // Доклады АН СССР, 1986, т.290, № 1, с.72-76.
4. Лунёв Б.В., Лапковский В.В. Быстрое численное моделирование соляной тектоники: возможность оперативного использования в геологической практике// Физическая мезомеханика, 2009, т.12, №1, с.63 - 74.
5. Лунёв Б.В., Лапковский В.В. Механизм развития инверсионной складчатости в подсолевом комплексе // Физика земли, 2014, № 1, с. 59 - 65.
6. GPU Gems 3 [Электронный ресурс]. - Режим доступа: https://developer.nvidia. com/content/gpu-gems-3
© Б. В. Лунев, Т. В. Абрамов, 2014
