Численное моделирование сейсмических процессов на высокопроизводительных вычислительных системах

Общая характеристика работы.

Актуальность темы

.

Землетрясения являются одной из наиболее разрушительных сил природы, приносящих человечеству огромные материальные убытки и большое число жертв. Значительная часть поверхности Земли, включая населенные пункты с развитой инфраструктурой, находится в сейсмоактивных зонах. Для обеспечения безопасности при проведении сейсмостойкого строительства используются методы расчета реакции сооружений на сейсмические воздействия. Расчеты должны включать оценки напряженно-деформированного состояния несущих конструкций, связанных с обеспечением их прочности при сейсмических нагрузках.

Наряду с решением проблем сейсмостойкости, большое внимание в последнее время уделяется сейсмической разведке, которая широко используется при поиске месторождений нефти, газа и других полезных ископаемых. Поэтому актуальной задачей является математическое моделирование распространения сейсмического сигнала в земной коре слоистой структуры, с наличием множества неоднородностей.

Работа посвящена созданию программно-моделирующего комплекса для математического моделирования волновых динамических процессов в пространственных конструкциях со сложной структурой и наличием неоднородностей. Подобный комплекс позволяет исследовать сейсмостойкость сложных наземных конструкций при прохождении через них волн, инициированных землетрясением. С другой стороны, становится возможным моделирование распространения упругих волн в слоистой земной коре, возбужденных различными источниками — как естественными, так и искусственными.

При расчете больших пространственных задач, когда рассматривается распространение волны на большие расстояния, требуется высокая точность моделирования. Для этого необходимо использовать численные методы с высокими порядками аппроксимации. Известно множество схем для численного решения уравнений переноса, обладающих большой точностью и отсутствием осцилляций. Чтобы решить указанные проблемы необходимо проанализировать и выбрать схемы, дающие наилучший результат для задач такого типа.

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

Научная новизна.

• Реализован комплекс программ для исследования пространственных динамических волновых процессов в твердых деформируемых неоднородных телах на структурированных сетках, в том числе, в многослойных, перфорированных и сложных наземных и подземных сооружениях. Особенностью комплекса является высокая производительность расчетов и возможность использовать в расчетах сетки с большим количеством узлов (до миллиарда).

• Проведен сравнительный анализ известных лимитеров для ТУБ схем второго порядка точности и схем до 4-го порядка точности. На основе полученных результатов выбрана схема, позволяющая получать достаточно точные численные решения в областях больших градиентов.

• Получено и исследованно поведение поверхностных волн Лява и Рэлея при численном решении динамической пространственной задачи теории упругости.

• Проведено моделирование процесса воздействия фронта упругой волны, инициированной землетрясением, через наземное сооружение, находящееся во вмещающем массиве. Определены области разрушения здания.

• Найдено численное решение задачи распространения волн, инициированных динамическим возмущением, в земной коре через водонасыщен-ные грунты.

• Предложен способ распараллеливания алгоритма для работы в распределенной кластерной среде на основе технологии MPI. Для распараллеливания работы программы в системах с общей памятью использована технология ОрепМР.

• Проведена оптимизация алгоритма решения систем гиперболических уравнений на регулярных сетках для работы на центральных процессорах с эффективным использованием кеша и потоковых инструкций SSE.

• Разработан механизм задания геометрии и входных параметров среды, использующий свободный пакет для преи постпроцессинга Gmsh.

Практическая ценность.

Предложенный в работе комплекс программ позволяет использовать численные эксперименты вместо натурных для ряда задач, связанных с сейсмостойкостью и оценкой последствий различных природных и техногенных катастроф. Исследования динамической стойкости жилых и промышленных наземных сооружений, подземных хранилищ, сложных промышленных сооружений таких как мосты, электростанции, плотины являются одним из приложений данного' программного комплекса. Часто стихийные бедствия сопровождаются большими человеческими жертвами и огромным материальным ущербом. Основными путями снижения этого ущерба является проектирование антисейсмических мероприятий для строительных объектов. Подобное проектирование не мыслимо без математического моделирования воздействий землетрясений на строительные конструкции. Выбор площадок под строительство промышленных стратегических объектов, оценка последствий землетрясений, взрывов и ударов по зданиям — эти и многие другие задачи необходимо решать с учетом предсказания последствий землетрясений на основе компьютерного моделирования.

Реализованный программный комплекс позволяет решать задачи моделирования распространения волн, взрывных и ударных нагрузок на трехмерные объекты со сложной внутренней структурой.

Разработанная параллельная версия комплекса программ позволяет существенно сократить время вычислений, а также увеличить разрешение моделей.

Работа выполнена при поддержке ряда государственных и коммерческих грантов:

• Программа (мероприятие): федеральная целевая программа «Научные и научно-педагогические кадры инновационной России» на 2009;2013 гг., в рамках реализации мероприятия № 1.2.1 «Проведение научных исследований научными группами под руководством докторов наук». Проект: «Разработка вычислительных технологий для моделирования пространственных динамических процессов в проблеме сейсморазведки на высокопроизводительных ЭВМ»;

• Федеральное государственное унитарное предприятие «Российский Федеральный Ядерный Центр — Всероссийский научно-исследовательский институт экспериментальной физики (ФГУП „РФЯЦ-ВНИИЭФ“)». НИР5 «Разработка физико-математических моделей, алгоритмов и эффективных методов решения задач механики сплошных сред на супер-ЭВМ»;

• Грант РФФИ 11−01−12 011;офи-м-2011. Разработка численных методов для решения задач геомеханики и сейсморазведки на многопроцессорных вычислительных системах, 2011;2012 гг.;

• Грант РФФИ 10−01−92 654-ИНДа. Математическое моделирование сложных задач на высокопроизводительных вычислительных системах. 2010.

— 2011 гг;

• Грант РФФИ 11−08−1 240-а. Моделирование ударных, взрывных и сейсмических воздействий на сооружения атомной энергетики. 2011 — 2012 гг;

• Small or medium-scale focused research project (STREP) proposal ICT EU-Russia Coordinated Call. FP7−2011;EU-Russia. 2011;2012 y.

Публикации.

Научные результаты диссертации опубликованы в 23 работах, из которых две [48, 51] - в изданиях, рекомендованных ВАК для публикации основных результатов диссертации.

В работах с соавторами лично соискателем выполнено следующее:

• Создан комплекс программ для численного решения пространственных динамических задач механики деформируемого твердого тела с учетом иеоднородностей среды, наличием подземных и наземных сооружений.

• Проведено моделирование распространения поверхностных сейсмических волн Лява и Рэлея, сравнение полученных результатов с теорией.

• Получено численное решение пространственной задачи о прохождении упругих волн, инициированных землетрясением, через наземные конструкции, определены области возможных разрушений.

• Проведено моделирование динамических процессов в многослойной коре земли, инициированных землетрясением, от гипоцентра землетрясения до земной поверхности.

• Получено численное решение задачи распространения волн, инициированных динамическим возмущением в земной коре, через водонасыщен-ные грунты.

• Проведен сравнительный анализ лимитеров для TVD-схем (второго порядка точности) и схем до четвертого порядка точности применительно к данному типу задач.

• На основе свободного программного обеспечения Gmsh создан пакет программ для задания входных данных (препроцессинг): расчетной геометрии, параметров среды, начальных и граничных условий.

• Реализовано распараллеливание комплекса программ для работы на многоядерных многопроцессорных машинах, в распределенной кластерной среде.

• Алгоритм оптимизирован для эффективной работы на больших сетках (до миллиарда узлов), проведена адаптация расчетного модуля для эффективного использования кеша процессора и ускорения расчетов с использованием потоковых инструкций SSE.

Апробация.

Результаты работы были доложены, обсуждены и получили одобрение специалистов на следующих научных конференциях:

Научные конференции Московского физико-технического институтаВсероссийские молодёжные научные конференции с международным участием «Проблемы фундаментальных и прикладных, естественных и технических наук в современном информационном обществе» (МФТИ, Долгопрудный, 2006 — 2009);

XVI международная научно-практическая конференция «Комплексная безопасность 2011 г.» (МЧС, Москва, 2011);

VI международная конференция по неравновесным процессам в соплах и струях (Санкт-Петербург, 2006);

20-th International Conference On Transport Theory (ICTT-20) (Обнинск, 2007);

26-th International Symphosium on Rarefied Gas Dynamics (Kyoto University, Japan, 2008);

XII международный семинар «Супервычисления и математическое моделирование» (РФЯЦ — ВНИИЭФ, Саров, 2011);

APOS-EU/APOS-RU Project Workshop «Modelling of seismic problems» (ИШ РАН, Москва, 2011);

NPP safety and personal training. XII International Conference. (Обнинск, 2011);

Российско-индийский семинар «Новые достижения математического моделирования» (ИАП РАН, Москва, 2011).

Результаты работы были доложены, обсуждены и получили одобрение специалистов на научных семинарах в следующих организациях:

• Всероссийский научно-исследовательский институт по проблемам гражданской обороны и чрезвычайных ситуаций МЧС России (федеральный центр) (Москва, 2011);

• Институт прикладной математики им. М. В. Келдыша РАН (Москва, 2011);

• Российский федеральный ядерный центр — Всероссийский научно-исследовательский институт экспериментальной физики (Саров, 2011).;

• Московский государственный университет путей сообщения (Москва, 2011).

Структура и объем диссертации

.

Диссертация состоит из введения, пяти глав и заключения. Общий объем диссертации составляет 110 страниц.

Список литературы

содержит 70 публикаций.

Выводы.

Распараллеливание вычислительного алгоритма показало, что наилучшие результаты дает применение технологии MPI, даже в системах с общей памятью, за счет лучшего распределения памяти между ядрами.

При распараллеливании удалось добиться эффективности до 76% при увеличении числа ядер в 64 раза. Применять технологию ОрепМР при расчетах на данном программном комплексе стоит в тех случаях, когда эффективность MPI уже достигла своего минимума. Тогда дальнейшее ускорение может быть достигнуто за счет технологии ОрепМР, при использовании данной технологии внутри одного расчетного узла на гибридных вычислительных машинах.

Применение технологий MPI и ОрепМР позволило существенно уменьшить время расчета и увеличить размеры используемых расчетных сеток.

Программный комплекс Архитектура.

Одним из результатов данной работы является программный комплекс, который позволяет численно решать задачи распространения волн, инициируемых различными динамическими нагрузками на трехмерные объекты со сложной внутренней структурой. Разработанная параллельная версия комплекса программ позволяет существенно сократить время вычислений, а также увеличить разрешение моделей. Программный комплекс применим также для решения задач ЗБ-сейсмики и сейсморазведки.

На рис. 4.9 изображена архитектура программного комплекса, он состоит из:

• вычислительного модуля, основанного на сеточно-характеристическом.

• расчет на блочной геометрии из комбинации сеток, в местах контактов сеток они должны быть согласованы;

• включает в себя монотонные сеточно-характеристические схемы 1−4 порядка точности;

• имеет возможность расчета как в один поток, так и на много ядерных вычислительных системах и на вычислительных системах с распределенной памятью.

Архитектура ядра вычислительного комплекса разрабатывалась с учетом дальнейших расширений возможностей программного комплекса. При написании программ существует эмпирическое правило 80−20, которое говорит о том, что 80% исходного кода выполняется 20% времени и, соответственно, наоборот. Процесс разработки и отладки был значительно ускорен и упрощен, за счет написания частей кода, которые не критичны к скорости выполнения или выполняются один раз за расчет, на скриптовом языке высокого уровня Python.

Вначале все важные модули, критичные к скорости выполнения, были написаны и отлажены на языке С/С-Н-, к ним относится расчетная сетка, граничные и контактные условия, начальные условия, сеточно-характеристические расчетные схемы. Затем на их основе были написаны модуль к языку Python, в результате чего появилась возможность работать с классами из языка С++ прямо из языка Python. Такой подход сильно упрощает отладку приложений, поскольку запуск можно вести в интерактивном режиме, используя оболочку языка Pyton. Помимо этого часть кода, отвечающая за организацию структуры расчетного модуля, частичную подготовку начальных данных, задание геометрии и структуры блоков, частично сохранение вынесены с систему скриптов.

Как оказалось, такой подход себя полностью оправдал. Программный комплекс нисколько не потерял в производительности, однако приобрел большую гибкость и расширяемость, что позволяет его в кратчайшие сроки адаптировать для решения новых задач.

Препроцессинг.

Препроцессинг включает в себя следующие этапы:

• задание расчетной геометрии;

• задание граничных условий;

• задание параметров области начального возмущения;

• задание условий на границе раздела сред и на границах соседних блоков;

• задание механических характеристик среды;

• задание вычислительного алгоритма;

• задание параметров сохранения.

Рассмотрим данные шаги подробнее.

Расчетная геометрия может задаваться несколькими способами. Прежде всего можно прямо в конфигурационном файле указать параметры всех расчетных сеток, а именно число узлов по всем осям, шаги сеток по всем осям и их положение в пространстве. Криволинейные сетки могут быть загружены из файлов, на данный момент поддерживается загрузка таких сеток их текстового варианта формата УТК, однако за счет использования скриптового языка, в программу легко добавить новый функционал по загрузке сеток, и.

Глава 5.

Результаты численного моделирования ряда задач.

Волны Лява Описание.

В слоистых средах возможно возникновение определенных типов волнволн Лява. Вектор смещения у таких воли параллелен границе раздела сред и перпендикулярен направлению распространения, т. е. волны Лява имеют горизонтальную поляризацию. В отличии от волн Рэлея, возникающих в одном полупространстве со свободной границей, волны Лява возникают в структурах типа упругий слой на упругом полупространстве. Теорию этих волн дал Ляв в 1911 г., именем которого они и названы.

Получить выражение для скорости волны в явном виде проблематично [31], определим волну по косвенным признакам и проверим, что соотношение скорости и длинны волны удовлетворяют уравнению волны Лява. Обозначим через cl скорость распространения волны Лява. Тогда для нее выполнено соотношение:

Cil.

На основе всех вышеизложенных результатов, можно сказать что полученная в расчете волна, действительно является волной Лява.

Волны Рэлея Описание.

Волны Рэлея возникают на границе полуплоскости, заполненной однородной изотропной упругой средой. Теоретически эти волны были открыты Рэлеем в 1885 году, и они могут существовать вблизи свободной границы твердого тела, граничащей с вакуумом или достаточно разреженной газовой средой. Фазовая скорость таких волн направленна параллельно границе, а колеблющиеся вблизи нее частицы среды имеют как поперечную, перпендикулярную поверхности, так и продольную составляющие вектора смещения. Фазовая скорость волн Рэлея не зависит от длинны волны, т. е. они являются без дисперсными. Эти волны очень быстро затухают по глубине полуплоскости [32], вследствие наличия экспоненциальных множителей с показателем —где д — волновое число, г — координата, направленная вглубь полуплоскости, а — некий множитель, зависящий от параметров среды и скорости волны Рэлея. Отсюда следует, что чем меньше длина волны (больше волновое число), тем быстрее происходит затухание. Получается, что волны Рэлея поверхностные, т. е. их основная энергия сосредоточена у границы.

Волны Рэлея имеют большое значение в сейсмологии, поскольку наблюдаются вдали от эпицентра землетрясений, и именно они являются основной причиной разрушения наземных объектов.

Обозначим через ср, с^ - продольную и поперечную скорость в среде соответственно, а через сд — скорость волны Рэлея. Соотношение между этими скоростями задается уравнением третей степени [31]: = С3 — 8?2 + 8(2 + х)£ - 8(1 + X) = о, (5.5) где.

Как известно, для всех упругих тел справедливо неравенство 0 < г) < /2, при учете этого анализ (5.5) показывает, что 0,8741 < < 0, 9554. Таким образом, скорость волн Рэлея мало отличается от скорости сдвиговых волн, но всегда меньше ее.

Можно выписать явный вид корней уравнения (5.5) при некоторых частных значениях упругих постоянных среды [32]:

1. х = О, V2 = 2, .

2. Х = 1/3, V2 = 3,? = 2(1 — - среда Коши.

Результаты моделирования.

Данным методом уже осуществлялось моделирование волн Рэлея в двумерном случае [35], однако там была представлена только качественная картина. Основной целью данного моделирования являлось получения волн Рэлея и сравнение полученных численным методом скоростей, со скоростями из уравнения (5.5). В тесте участвовали два частных случая, описанные выше.

Расчетная область для всех тестов представляла собой параллелепипед размерами 1500×500×200 м, соответственно по осям х, у иг. На верхней грани области по оси г было выставлено граничное условие свободной границы, на других гранях устанавливалось граничное условие поглощения. Шаг сетки во всех расчетах брался порядка 10 м, число узлов в сетке — около 160 тыс. Шаг интегрирования по времени был выбран исходя из выполнения условия Куранта, во всех тестах он равен 0,1 615 с. На рис. 5.5 представлена.

Рис. 5.13. Расчетная геометрия скорость и доходит до земной поверхности позднее.

На рис. 5.15 показаны сейсмограммы, взятые на дневной поверхности у данных расчетов. Хорошо видно, что от взрыва приходит сильный фронт продольной волны, а от землетрясений сильный фронт поперечной. Полученный результат имеет хорошее сходство с экспериментальными данными. Также на сейсмограмме хорошо заметны отражения от всех слоев земной породы.

Прохождение волн через водонасыщенные грунты Описание.

В данном разделе рассмотрена задача прохождения сейсмических воли через водонасыщенные грунты. Работа выполнялась совместно со специалистами из Российского Федерального Ядерного Центра — Всероссийского научно-исследовательского института экспериментальной физики (ФГУП «РФЯЦ-ВНИИЭФ»).

Специалистами из Сарова был создан программный комплекс «Нимфа» для численного моделирования задач фильтрации. Этот комплекс ориентирован на решение традиционных гидрогеологических задач, связанных с прогнозом техногенного режима подземных вод при работе водозаборов и ведении в виде фронта упругой продольной волны, направление фронта — к земной поверхности.

На рис. 5.18 показано распределение скорости на поверхности земли, после отражения фронта упругой волны в различные моменты времени.

1.0173 0Я16.

Рис. 5.18. Распределение модуля скорости на поверхности земли в моменты времени i = 0.131,0.138,0.146,0.153,0.160,0.167 с от начала расчета.

На рис. 5.19 показано распределение модуля скорости во вмещающем массиве.

Как видно из расчетов, различная водонасыщенность вносит существенные неоднородности в картинку распространения упругих волн. На дневной поверхности земли возникли сейсмические поверхностные волны, их распространение хорошо видно по рис. 5.18.

Проведена сшивка результатов расчета программного комплекса «Нимфа» для численного моделирования задач фильтрации, с программным комплексом, разработанным в рамках данной работы. На тестовых расчетах показано влияние насыщенности пор водой на изменение в распространении упругих волн.