Разработка алгоритмов решения систем гиперболических уравнений на графических процессорах (1187418), страница 5
Текст из файла (страница 5)
Одна из проблем, которая анализировалась в ходе выполения работы, это различия в вычислительных сетках. В газодинамической задаче, которая использовалась для вычисления характеристик головной волны, в первых расчетах сетка не совпадала с сеткой, задаваемой в сейсмической части. Для устранения этой зависимости, была произведена аппроксимация данных с узлов более мелкой сетки (газодинамической) на более крупную (сейсмическую). Стоит отметить, что были произведены расчеты, доказывающие, что ошибки аппроксимации в данном случае пренебрежимо малы. Сетка в газодинамической задаче была укрупнена для того, чтобы мы могли точно наложить узлы газодинамической сетки на сейсмическую.
Сейсмическая задача состояла в следующем. Входными данными является зависимость давления от времени на площади 40 км x 40 км. Результатом вычисления должны быть данные на расстоянии 1600 км от центра взрыва, где находится Обнинская сейсмостанция, на которой был получен сейсмический сигнал от Челябинского метеорита, упавшего 15 февраля 2013 года. Сигнал был очищен от влияния других событий, таких как землетрясений, произошедхих в этот день.
Рис. 6.1 Сигнал, зафиксированный в Обнинске и очищенный от других сейсмических явлений
На Рис. 6.1 показана зависимость вертикальной скорости от времени, и АЧХ (амплитудно-частотная характеристика). Как видно из Рис. 6.1, основная частота (доминантная частота) приблизительно равна 0.03 Гц. Кроме этого, видно что длительность сигнала составляет около 70 с. При t < 490 с и t > 560 с сигнал определяется уже не головной ударной волной Челябинского метеорита, а слабыми землетрясения. Поэтому мы не будем учитывать эти части сигнала при дальнейшем сравнении полученных сейсмограм с сейсмограммой из Обинской сейсмостанции.
Кроме зависимости давления от времени на области 40x40 км2, в качестве начальных данных необходимо задать продольную и поперечную скорости волн.
Для задания геологических свойств среды были приняты следующие предположения. На основании реальных данных и из оценки средней скорости поверхностной волны Рэлея [10], инициированной головной ударной волной Челябинского метеорита, среднее значение длины волны составило примерно 100 км. Для такой достаточно большой длины волны, учитывая структуру земной коры Русской равнины и Уральских гор, основное влияние оказывают более глубокие, базальтовые и гранитные слои, несмотря на различие в структуре приповерхностного слоя [10]. На основании этого на всей протяженности распространения волны от источника до момента регистрации среда рассматривается как однородная с наиболее приемлемыми средневзвешенными значениями продольной и поперечной скорости распространения волн Vp = 7,18 км/с и Vs = 3,4 км/c [11]. Рассмотрим вопрос передачи данных из газодинамической части в сейсмическую. Передача данных осуществлялась на основании прямого механизма, то есть избыточное давление из распределения на поверхности Земли в газодинамическом расчете передавалось на прямую как функция, определяющая упругие деформации в сейсмической модели, и формировало сейсмический источник [12].
-
Постановка задачи
Физическая расчетная область представляла собой прямоугольный параллелепипед с размерами 2400км x 800км x 150км (длина, ширина, высота). Верхняя граница (дневная поверхность) задавалась как свободная, на остальных гранях устанавливались поглащающие условия. Хотя Обнинск находится на расстоянии в 1600 км от места падения метеорита, для корректной работы необходима значительно большая расчетная область. Это связано со сложностью правильной задачи трехмерных поглащающих граничных условий. Из-за несовершенства граничных условий часть волны отражается и влияет на полезный сигнал. Для решения этого вопроса, было добавлено дополнительное расстояние по краям расчетной области.
Центр области, на которой задается давление, помещался на верхней грани (дневной поверхности) на расстояние 400 км от 3 граней расчетной области - Рис. 6.2.
Рис. 6.2 Расчетная область
Аналогично (5.6) мы решаем трехмерную систему
(
)
Где – вектор, состоящий из 9 компонент, которые необходимо найти.
. Вектор состоит из 3 компонент скорости и 6 состовляющих тензора напряжений в виду его симметричности.
Матрицы принимают следующий вид:
(
)
(
)
(
)
В матрицах -
– постоянные Ламе. Система является гиперболической, а это означает, что систему можно представить как систему из 9 одномерных дифференциальных уравнений. Что значительно упрощает вычисления.
Для решения одномерных уравнений использовалась TVD-схема с ограничителем MC [13]:
(
)
Опишем граничные условия, которые были в постановке задачи.
(
)
Где – матрица, набор строк матрицы
, которые относятся к выходящим характеристикам для границы.
Данные граничные условия хорошо работают только для волн, направляющихся в каком-то одном направлении, соответственно не являются идеальными для работы со сферическими волнами.
(
)
Покомпонентная запись граничных условий:
(57)
(58)
Более подробное описание используемых граничных условий можно найти в [14].
-
Результаты
Перед получением основных результатов, необходимо было провести серию расчетов, которая позволит понять на сколько изменяется результат при варировании параметров задачи. Основные параметры входа были описаны в главе 6.2. Стоит еще раз сказать, что скорость падения метеорита и угол его вхождения в атмосферу известны достаточно точно, поэтому в исследовании не варировались. Однако, при проведении начальных расчетов, было выяснено, что в пределах допустимых погрешностей они слабо влияют на конечный результат. Как конечный анализ, проводилось сравнение полученного сигнала из Обнинской сейсмостанции и амплитудно-частотной характеристики с смоделированным сигналом.
В ходе исследования, было получено, что изменение каждого из факторов в отдельности: увеличение высоты взрыва, уменьшение диаметра и уменьшение скорости приводит к уменьшению амплитуды сигнала, что является интуитивно понятным.
В качестве основных вариантов были выбраны 3, которые более точно описывают сигнал, полученный в Обнинске.
№ | Высота взрыва, км | Плотность метеорита, г/см3 | Характерный диаметр, м |
1 | 30 | 3,6 | 18 |
2 | 24 | 3,6 | 18 |
3 | 30 | 3,0 | 16 |
Таблица 4
При этом во всех расчетах значение скорости вхождения в атмосферу км/ч и под углом 18o . Значения для скорости и угла были взяты из [10, 11]. Так же из [10,11] были взяты оценки высоты взрыва
км, плотность хандрита (материала метеорита) варируется от 3.0 – 3.7 г/см3 , в качестве характерного диаметра можно взять интервал (16; 19) м.
Поясним, почему в качестве основого результата рассматриваются именно эти 3 варианта постановки задачи. 1-ый вариант показывает наилучшее совпадение модельных данных сейсмического отклика с реальными. 2-ой вариант считается наиболее вероятным, опираясь на статьи [10,11]. Вариант 3 является результатом при наборе параметров, который дает минимальную амплитуду - оценка снизу. В ходе предварительных расчетов было выяснено, что оценку сверху проводить не имеет смысла, потому что она сильно завышает амплитуду сигнала.
Рис. 6.3 Сравнение реального сигнала и 3 постановок задачи.
Рис. 6.4 Сравнение реального сигнала с постановками 1 и 2 отдельно.
Как видно из Рис. 6.4 постановка 1 лучше описывает пики 2 и 3, а постановка 2 лучше описывает первый пик. Для более полного описания картины, приведем спектральную характеристику сигналов.
Рис. 6.5 Амплитудно-частотная характеристика сигналов.
Как видно из Рис. 6.5 доминантная частота у всех 3 вариантов моделей практически совпадает с соответствующим значением у реального сигнала. Однако, для постановки 2 и 3, амплитуда оказывается значительно ниже, чем у реального сигнала или у сигнала, полученного в первой постановке.
При сравнении вариантов 1 и 3, а именно, существенном уменьшении плотности при незначительном уменьшении диаметра, можно заметить, что амплитуда сильно уменьшилась. Из данного моделирования можно заключить, что метеорит состоял из однородного вещества, а не из льда и хондрита, как считают некоторые исследователи.
Для более подробного анализа постановки 1, приведем значение поверхностной волны Рэлея, которая была получена из реальных данных. Скорость волны Рэлея состовляла ≈3 км/с. Скорость распространения поверхностной волны легко получить из сейсмограммы для постановки 1 (Рис. 6.6).
Рис. 6.6 Сейсмограмма и сигналы на датчика, находящихся на разном расстоянии от источника сигнала
Кроме этого, волну рэлея можно получить из анализа сейсмограм, полученных с датчиков, находящихся на различном расстоянии от источника сигнала. В случае постановки 1, мы получили скорость км/с, что согласуется с результатами, полученными от реального источника. Это еще одно подтверждение корректности постановки задачи.
Рис. 6.7 Распределение давления на поверхности в моменты t = 13 c, t = 28 c.
-
Выводы
-
Исследованы схемы Лакса, Бима-Уорминга, TVD2 для уравнения переноса на сходимость.
-
Исследовано ускорение на примере уравнения переноса при вычислении на графическом процессоре.
-
Проведено исследование системы уравнений упругости на GPU.
-
Разработана и реализована программа для численного решения системы уравнений упругости на нескольких графических устройствах.
-
Исследовано ускорение при работе программы на нескольких графических устройствах.
-
Разработана модель падения Челябинского метеорита совместно с лабораторией математического моделирования нелинейных процессов в газовых средах, что позволило уточнить его характеристики.
Публикации, конференции:
-
Дашкевич А.Д. “Решение систем уравнений гиперболического типа на графических процессорах с использованием технологии CUDA”, тезисы конференции, конференция “Ломоносов”, МГУ, 2013.
-
Дашкевич А.Д. “Решение систем уравнений гиперболического типа на графических процессорах с использованием технологии CUDA”, конференция МФТИ, МФТИ, 2013.
-
Дашкевич А.Д. , Хохлов Н.И. “Решение систем уравнений гиперболического типа на графических процессорах с использованием технологии CUDA”, сборник трудов МФТИ, 2013.
-
А.Д. Дашкевич, Н.И. Хохлов, В.И. Голубев “Моделирование распространения акустических волн в земной коре при падении метеорита”, сборник МОУ, 2014.
-
“Анализ результатов трехмерного моделирования влияния головной ударной волны Челябинского метеорита на поверхность Земли ” конференция МФТИ, 2014.
-
Совместная статья с лабораторией математического моделирования нелинейных процессов в газовых средах А.В. Астанин, А.Д. Дашкевич, И. Б. Петров, М.Н. Петров, С. В. Утюжников, Н. И. Хохлов “Моделирование влияния головной ударной волны Челябинского метеорита на поверхность Земли” подана в печать.
-
Планы
-
Ведутся работы по переносу алгоритма на 3D задачи
-
Перенос на гибридную архитектуру (MPI+CUDA)
-
Реализация на OpenCL.
-
Поддержка различных граничных условий
-
Поддержка криволинейных сеток
-
Перенос на неструктурные сетки
-
Литература
-
Хохлов Н. И., Петров И. Б. Моделирование сейсмических явлений сеточно-характеристическим методом // ТРУДЫ МФТИ. 2011. Т. 3, ќ3. С. 159167.
-
http://www.paraview.org/
-
RANDALL J. LEVEQUE , Finite Volume Methods for Hyperbolic Problems, ISBN 0-521-00924-3
-
Analysis of Numerical Methods by E. Isaacson, H. B. Keller, ISBN 0-486-68029-0
-
http://www.nvidia.ru/object/what_is_cuda_new_ru.html
-
https://developer.nvidia.com/nvidia-gpu-programming-guide
-
CUDA by example. Jason Sanders, Edward Kandrot.
-
https://developer.nvidia.com/cuda-toolkit-40
-
Kamran Karimi Neil G. Dickson Firas Hamz. A Performance Comparison of CUDA and OpenCL.
-
Brown P. G., Assink J. D., Astiz L., Blaauw R., Boslough M. B. et al. A 500-kiloton airburst over Chelyabinsk and an enhanced hazard from small impactors. Nature, 2013
-
Emel’yanenko V. V. et. al. Astronomical and Physical Aspects of the Chelyabinsk Event (February 15, 2013). Solar System Research, 2013, Vol. 47, No. 4, pp. 240–254. 14
-
Edwards W.N., Eaton D.W., Brown P.G. Seismic observation of meteors: coupling theory and observations. Reviews of Geophysics. 2008. Т. 46. № 4
-
van Leer B. Towards the ultimate conservative difference scheme. III Upstream-centered finite-difference schemes for ideal compressible flow. IV – A new approach to numerical convection // Journal of Computational Physics. — 1977. — mar. — Vol. 23. — Pp. 263–299.
-
Челноков Ф.Б. Численное моделирование деформационных динамических процессов в средах со сложной структурой, Диссертация на соискание ученой степени кандидата физико-математических наук, Москва, 2005.