автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное решение пространственной динамической задачи теории упругости на многопроцессорных вычислительных системах

На правах руь опш и

Кучунова Елена Владимировна

ЧИСЛЕННОЕ РЕШЕНИЕ ПРОСТРАНСТВЕННОЙ ДИНАМИЧЕСКОЙ ЗАДАЧИ ТЕОРИИ УПРУГОСТИ НА МНОГОПРОЦЕССОРНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ

Специальное!ь 05 13 18 - математическое моделирование, численны^ методы и комптексы программ

АВТОРЕФЕРАТ

диссертации па (оис канне ученой с тепени кандидата физико-математических паук

Красноярск - 2008 г

00344Э185

003449185

Работа выполнена в Федеральном государственном образовательном учреждении высшего профессионального образования 'Сибирским федеральный университет'' и в Учреждении Российской академии наук Институте вычислительного моделирования Сибирского отделения РАН

Научный руководитель доктор физико-математических наук,

профессор Садовский Владимир Михайлович

Официальные оппоненты доктор физико-математических наук,

профессор Добронец Борис Станиславович

кандидат физико-математических наук, профессор Распопов Витаний Ешеиьевич

Ведущая организация Учреждение Российской академии наук Институт нефтегазовой геологии и геофизики им А А Трофимука Сибирского отделения РАН

Защита диссертации состоитс я "31" октября 2008 г в 14 часов на заседании диссертационного совета ДМ 212 099 06 при Сибирском федеральном университете по адресу 660074, г Красноярск, ул Киренского, 26, ауд УЛК 115

С диссертацией можно ознакомиться в библиотеке Сибирского федерального университета по адресу г Красноярск, ул Киренского 26, ауд Г 376 Автореферат выставлен па сайте СФУ

01зыв на автореферат в 2-х экземплярах, с подписью составителя и заверенный печатью организации просим направлять в адрес диссертационного с овета 660074 г Красноярск, ул Киренского, 26, ауд УЛК 319

Автореферат разослан '30" сентября 2008 г

Ученый секретарь диссертационного совета, кандидат технических наук, доцент

Л

PK) Царев

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

Диссертация посвящена разработке численною метода моделирования волновых полей для никоторого класса н'офнзическнх сред, созданию комплекса параллельных прикладных программ и исследованию на этой ос нове распространения волн в с ложно построенных средах При этом рассматриваются все основные этапы решения прикладной задачи (выбор магматической модели построение ее дискретного аналога, разработка члененною метода и параллельного алгоритма его реализации на многопроцессорных вычис лительных с ис I ем ах)

Актуальность работы

Масемашческие модели механики деформируемых сред, описывающие процессы распространения упругих волн напряжений и деформаций, всегда пред-с гавляли практический интерес Maie\iaiическое моделирование дает физические представления об образовании волн в различных средах Эти представ к'ния положены в основу сейсмических методов исследования, применяемых при изучении внутреннего строения Земли, поис ко и разведке полезных ископаемых, в инженерно-геологических изысканиях

Начиная с работ Релея и Лзмба, математическое моделирование распространения сеисмичес ких волн проводится на основе с ис гемы уравнений динамической теории упругости На практике важно проанализировать сейсмическое волновое поле с учетом всех типов волн, волн интерференционного характера и возникающих при этом "нелучевых" явлений Совершенствование систем наблюдения и повышение детализации экспериментальной информации приводит к необходимое!и рассмотрения волновых полей с большей точностью и для все более возрастающих времен Иск/ночшелыюе значение для исследования распространения сейсмических волн в сложнопостроепных средах имеет вычислительный эксперимент Однако применение прямых методов- числс иного анализа требует большо! о объема вычис лительных ресурс ов Решение прос трапственпых задач на современном персональном компьютере может ока)атьея невозможным из-за длительности времени счета Кроме того однопроцессорный компьютер зачастую не может обеспечить требуемый обьем оиераппшой памяти для работы с многомиллионными ячейками расчетных сеток использование которых необходимо для достижения требуемой точное гп Эффективным средством решения ресурсоемких простран-

(твенных задач является технология параллельных вычислении на распределенных системах кластерного типа Развитие суперкомпьютеров позволяет успешно решать задачи в трехмерной пос тановке

На пути перехода от персонального компьютера к суперкомпьютеру с параллельной архитектурой имеются определенные трудности, которые во-первых, связаны с необходимостью распараллеливания вычислительного алгоритма, во-вторых, с принципиально более сложным написанием программного кода, недостаточно развитой системой отладки программ на кластере и пр При этом одним из факторов эффективного ис пользования высокопроизводительной многопроцессорной вычислительной техники является применение специальных методов В связи с этим, одной из задач современной науки является развитие технологии математического моделирования рассматривающей многопроцессорные вычислительные системы (МВС) как основной аппарат для вычислений Отметим, что методы моделирования с ис пользованием высокопроизводительных компьютеров часто оказываются едва ли не единственным способом получения информации об исследуемых явлениях, так как возможности теоретического и экспериментального изучения этих явлений ограничены методическими и техническими трудностями

01 мет им далее, чго за время существования МВС накоплен значительный опыг их использования и они находятся в мае с овом доступе (МВС имеются и ряде Российских и в большинстве зарубежных научных иненпугов) Однако, огромные вычислительные возможности предоставляемые этими системами, используются недостаточно Причиной этому является с ложность адаптации последовательных алгоритмов и программ к параллельным архитекторам Поэюму разработка эффективных параллельных алгоритмов решения фундаментальных и прикладных задач, изначально ориентированных на ис пользование МВС (в том числе с распределенной памятью) является чрезвычайно актуальной проблемой

Цели и задачи исследования

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

• разработка параллельного вычислительного алгоритма расчета во пю-вых полей 1! трехмерных кусочно-однородных упру: их средах,

• адаптация разработанного алгоритма к многопроцессорным вычис ли-

тельным системам, которая складывается из выбора методов реализации алгоритма на параллельных ЭВМ выяплепия трудно распараллелпвае-мых блоков их оптимн '.ации и верификации,

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

Объект исследований

Обьекгом исследования диссертации являются трехмерные прямые динамические задачи сеисмики для упругих кусочно-однородных изотропных сред

Методика исследования

В качестве методики исследовании используется вычислительный эксперимент включающий в с ебя следующие ^тапы математичес кая формулировка задачи, построение приближенного (численного) метода решения задачи написание вычислительного алгоритма, программирование на ЭВМ вычислительного алгоритма, проведение расчетов па ЭВМ, анализ полеченных численных результатов и уточнение матемашчес кон модели

Новые научные результаты, выносимые на защиту

1 Параллельный вычислительный алгоритм для решения гиперболической с ш гемы уравнений теории упругости на криволинейных сетках на ос но-ве метода двуциклического расщепления но пространственным переменным с использованием ENO-схемы для решения одномерных систем

2 Универсальный алгоритм распределения многобпочнои области по процессорам основанный па требовании равномерной вычис тигельной нагрузки

3 Комплекс параллельных программ для численного моделирования распространения волн в кусочно-однородной среде в трехмерной геометрии, созданный на основе этих алгоритмов

Научная новизна и практическая ценность

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

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

Апробация работы

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

1 VLIII Международной научной студенческой конференции "Студент и научно-технический прогресс" (НГУ, г Новосибирск, 2005г),

2 V,VI и VII Всероссийских конференциях по математическому моделированию и информационным технологиям (ИВТ СО РАН, г Новосибирск 2004г КемГУ, г Кемерово, 2005г ИВМ СО РАН, г Красноярск, 2006г),

3 IV,V и VI Межрегиональных школах-семинарах "Распределенные и кластерные вычисления" (ИВМ СО РАН г Красноярск 2004г 2005г , 2006i ),

4 Конференциях-конкурсах молодых ученых ИВМ СО РАН (г Красноярск 2005г 2006г, 2007г),

5 XXVIII н XXIX Конференциях молодых ученых механнко-математичес ког факультета МГУ (МГУ г Москва, 2006г, 2007г),

6 Международных научных конференциях "Параллельные вычис литель-ные технологии' ПаВТ'2007 и ПаВТ'2008 (ЮУрГУ г Челябинск, 2007г , СПбГПУ, г Санкт-Пет ep6ypi, 2008г)

Кроме того результаты диссертации докладывались па семинаре "Проблемы математического и численного моделирования" Института вычислительного моделирования СО РАН, на семинаре кафедры информатики Институт математики Сибирского федерального университета, а также на семинаре кафедры волновой и i а юной динамики механико-математического факультета Московскою государственного университета им М В Ломоносова

Публикации

По теме диссертации автором опубликовано 10 работ (из них 2 по списку ВАК) Из работ, выполненных совместно в диссертацию включены результаты, поллченные автором лично Список публикаций помещен в конце автореферата

Рабсил выполялась и]>и финансовой поддержке РФФИ (коды проектов 0401-00267 и 08-01-00184), Комплексной Программы фундаментальных исследований Президиума РАН N 17"Параллельные вычисления на многопроцессорных вычислительных системах" и N 14 'Фундаментальные проблемы информатики и информационных технологий", Красноярского краевого фонда пауки (шифр граи га 17С030)

Структура и объем работы

Диссертационная работа состоит из введения, четырех глав заключения, списка литературы содержащего 123 наименований, и двух приложений Диссертация содержит 66 рисунков Объем диссертации составляет 133 страницы, приложений - 15 страниц

Основное содержание работы

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

В первой главе раскрывается современное с остояние проблем, с вя з.шных с численным решением задач теории упругости на многопроцессорных вы-чис лпгельпых сис темах В силу того, что данная работа включает несколько разделов, относящихся к различным математическим дисциплинам данный обзор разделен на пять час I ей В разделе 1 1 дается обзор работ, посвященных методам моделирования волновых процессов в упругих средах В разделе 12 рассмотрены основные методы дискретизации расчетной области В разделе 1 3 представлены работы, посвященные пос троению разностных схем для решения гиперболических систем уравнении В разделе 1 4 дан обзор методов моделирования неотражающих граничных условии Раздел 1 5 посвящен ос новам организации вычислений на многопроцессорных системах аффективному их исполь зовашпо в решении больших задач

Во второй главе диссертации рассматриваете я алгоритм решения трехмерной упругой задачи с помощью метода двуциклического расщепления по пространственным координатам и применения известной монотонной ЕНО-схемы 'предиктор-корректор"

В и 2 1 дается математическая постановка задачи восстановления волновых полей в кусочно-однородной области $2 £ К3 по заданном силам внешнего воздействия Предполагается что исходная среда имеет сложную структуру - содержит внутренние поверхности раздела материалов с существенно

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

ру = V а, а = Л(У v)/+//(V ® v + v ® V)

(1)

Здесь р - плотность материала, А и ¡1 - параметры Ламе, определяющие ) пругие свойства, V - вектор скорости, а - тензор напряжения, V - градиент по пространственным координатам, I - единичный тензор Точка над символом, как обычно, означает производную по времени Предполагается что в О задано начальное распределение скоростей и напряжений

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

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

В и 2 2 опис ываетс я алгоритм построения расчетной с с1 тки на ос нове алгебраического метода, состоящего в нахождении взаимно-однозначного отображения вычислительной области в виде единичного куба [О I]3 с равномерной сеткой на физическую область Перед построением сетки для аппроксимации граничных поверхностей области и внутренних поверхностей раздела исходная кусочно-однородная область представляется регулярным набором блоков £2 = иГ=1 и"=1 Ш=1 Пи*, 1Де гп\ - количество слоев, гп2 - количество полос и с лое и т3 - количество блоков и полосе Далее расчетные сетки строятс я не завис имо в каждом блоке, таг сетки выбирается 1! ¡апис имос ти от упругих свойств материала и характерной длины возбуждаемых волн Соображение вычислительной области на физическую область блока ищется в виде многомерного кубического сплайна х(£) = ^ г_0си/(^у (£з)' Векторные коэффициенты с,,/ находятся в явной форме и5 условий сопряжения в вершинах блоков заключающихся в том, что положения вершин заданы, а направления выхода координатных линии в вершинах определяются из соображений ортогональности

у(хО)=у°(х), ст(х 0) = ст°(х)

(2)

ст п = я,

(3)

В н 2 3 С1 роится вычислтетьный алгоритм на основе с еточно-характе-ристпчсского метода При построении разностной схемы в задаче (1)~(3) исходная система (1) записывается в мафичнои форме

ди тг-^ ди

1=1

01

' д.гг

(4)

Здесь и = = (уиу,, 03,^11,022,033,023, ет12)Т - вектор-функция

составленная (п компонент вектора скорости и тензора напряжении Аг -матрицы-коэффициенты размерности 9x9, произведение которых на вектор и задается равенством

1

А

А(с, V)/ + //(е, ® V + V ® е,)

где С!, е2, е3 орты декартовой с не темы координат

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

ди у-л ~ ди

(5)

- ,=

С учетом обозначения а' — V получается следующая формула для вы-чис лепия произведения матриц на вектор

А.

1а*

Р

А(а'

г)1 + //(а1 ® V + V ® а1)

Применение -метода двуциклического расщепления по пространственным переменным к системе (5) приводит к серии из шести одномерных задач на интервале t € [¿„ + т] Каждая их полученных одномерных систем является гиперболической обладает полным набором собственных векторов с действи[сльными собственными значениями поэтому каждую из таких систем молено переписать в виде

2н = д-т 4 '

(6)

где <2 - матрица составленная из линейно-независимых собственных векторов Л - диагональная матрица, элементами коюрои являются собственные

значения Для направления £ =

Л = |, —с%|а'|, —с,,|а!|, —с4|а'| 0,0,0},

|а'| = у/{<у\)2 + (»г)2 + (аз)2' (-р = + 2/х)/р - скорость продольных волн, = л/й/р ~~ скорость поперечных волн в среде

К одномерным системам (6) применяется явная монотонная ЕКО-счема предиктор-корректор с предельной реконструкцией решения На шаге "предиктор", исходя из характеристических соображений, по известным значениям на п-и временном с чое вычисляются значения решения на границах ячеек Система (6) после замены неизвестной вектор-функции по формуле ш = фи распадается на девять независимых скалярных уравнений переноса

'дt= ~дS,

для решения которых из узла (£,,/„ + т/2) = 0,1, ЛГ) опускаются характеристики па слой / = Компоненты вектора «г постоянны вдоль характеристик поэтому ш"+1/2(С,) = + 1/2т/2), где Л*:,+1/2 - угловой коэффициент А--ОЙ характеристики в /-ой ячейке Характеристики проведенные из узла + т/2), могут не попасть точно в расчетные узлы 1/2> ^п) в которых определено решение, поэтому применяется процедура предельной реконструкции сеточная функция с узловыми значениями м>"л+1/2 на каждом из отрезков [£л£;+1] заменяется кусочно-линейпым сплайном = + п,к}+1/2' коэффициенты которого подбираются из условия минимума суммы модулей скачков сплайна на границах ячеек В итоге вектор имеет следующий вид

{о",+1/2 т

<,+1/2--— (С+1 -0+2 А^ + 1/2а"^ + 1/2 еСЛИ А^+1/2 > 0

^ А 1—1/2

Ч',-1/2 + -2- ^ ~~ + 2 а^-1/2' РСЛИ < 0

„ 71+1/2 п+1/2

После этого восстанавливается искомое решение

Последнюю формулу нельзя использовать в приграничных ячейках, поскольку н нее входят величины в точках которые находятся за пределами области интегрирования Корректное решение задачи на границе требует постановки граничных условий в количестве, равном числу уходящих характеристик Для левой границы области уходящими являются характеристики отвечающие отрицательным собственным значениям, для правой границы - положительным В качестве недостающих условий на дневной поверхности области И выступают граничные условия (3) На остальной части

поверхности ставятся искусственные граничные условия на основе принципа неотражения воли Следуя этому принципу, поыгаются равными нулю граничные шачения харак1еристических переменных па уходящих характеристиках что эквивалентно отсутствию отраженных волн в одномерной задаче В случае, если 1рапица блока не является внешней, для вычисления

/ н+1/2 П+1/2-,

значении I! граничных ячейках (и0 ' , и^у ) используется алгоритм склеи-ки решения на 1раницах соседних блоков

Для сохранения свойств консервативности схемы па криволинейных сетках, когда матрицы-коэффициенты А, зависят от координат, аппроксимация одномерных систем уравнении па шаге "корректор" производится с помощью иитегро-интерполяционного метода В результате интегрирования системы (5) по криволинейной ячейке в физической области с последующим применением формулы Грина получается следующее равенство

ж=Ь £ ^41++"зМл) (7)

в котором щ - направляющие косинусы внешней нормали, и - объем ячеики, и - среднее интегральное значение вектора-решения, а индекс к относится к граням ячейки в частности 7* - площадь соответствующей грани Далее сумма в правой части (7) разбивается на три пары с лягаемых по противоположным граням каждая из которых отвечает аппроксимации производных по пространственным переменным в одномерных системах (6) Например, слагаемые, относящиеся к координатным поверхностям £2>£з Дают аппроксимацию выражения Ахди/д^х

Заметим что применяемая схема решения одномерных задач является вариантом сс'ючно-характернсшчсской схемы второго порядка точности в областях монотонного поведения исходных функции Метод двуциклического расщепления сохраняет второй порядок Кроме того двойной пересчет одной и той же системы уравнений па третьем и четпертом этапах расщепления гарантирует устойчивость метода при выполнении одномерного условия Кура и I а-Фридрихгд-Л еви

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

В н 3 1 обсуждаются основные проблемы реализации параллельных алгоритмов, связанные с необходимостью осуществления в ходе вычислений межпроцессорных обменов и записи результатов в файлы что влечет за собой замедление, а следовательно снижение эффективноети распараллеливания Формулируются основные критерии качества параллельных реализа-

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

Для сбалансированности вычислений и минимизации межпроцессорных обменов ключевая роль отводится выбору способа распределения данных и вычислений по процессорам Для выбора конкретной технологии распараллеливания алгоритма в п3 2 приводится сравнение двух наиболее распространенных технологий MPI и DVM на примере решения одномерной задачи распространения упругой волны сжатия в неоднородном стержне

Здесь v(f, г) - скорость частиц в направлении оси стержня, a(t, т) - напряжение на площадке, нормальной к оси заданные функшш р(т) > ро > О и с(т) > 0 имеют с мыс л плотности среды и скорости распространения волн Для однозначного определения решения задаются начально-краевые ус ловия В начальный момент времени стержень находится в покое, поэтому v(0,x) = О, <т(0 т) = О На левой границе задано значение скорости «(/,0) = u(f), на правой — напряжение cr(f, 1) = 0 Исходная ели тема уравнений (8) решается численно моноюнной ENO-схемои с предельной реконструкцией

В работе подробно описывается каждая из сравниваемых техноло1 ий Обе эти технологии основываются на принципе геометрического параллелизма Исходную вычислптельщ ю область Ç [0 1] / С [О У)} делится па вертикальные полосы по числу процессоров, каждый hï которых вычисляет неизвестные функции только па своей полосе, обмениваясь при необходимости граничными значениями с соседними процессорами В случае техноло-i ии MPI каждый процесс самостоятельно строит часть сетки в своей полосе, причем сетки соседних процессов пересекаются по двум узлам Эю пересечение сеток позволяет частично продублировать вычис пения на соседних

(8)

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

В случае технологии DVM необходимо написать только один вариант программы и для последовательного и для параллельного выполнения Эта программа содержит правила параллельного выполнения алгоритма, остающиеся "невидимыми ' для стандартных компиляторов для последовательных ЭВМ Все используемые в программе массивы (массив сетки, массивы решения и тп ) определяются как распределенные следующей DVM-директивой DVM(DESTRIBUTE|BLOCKj) double *Х, где X - распределенный массив Распределение определяется однозначным соответствием каждого элемента массива и номером процессора на который элемент распределен На каждом процессоре вычис ляете я только локальная час ть массива при помощи специализированных "распределенных циклов" Для доступа к 'удаленным данным'' (данные используемые на одном вычислительном узле но размещенные на другом) применяются специальные "тс пены с' грани, представляющие собои буфер, непрерывно продолжающий локальную секцию массива на указанное количество ячеек В нашем случае таких теневых граней было по одной с каждой стороны секции массива В процессе расчетов возможны случаи синхронного и асинхронного обновления теневых граней В первом случае обновление теневых граней является непрерывной операцией, и делиться на две запуск обмена и ожидание значении Выполнение обменов осуществляется одновременно всеми вычислительными узлами При асинхронном обновлении вычис ленне и обмен являются независимыми операциями н могут выполняться одновременно Так на фоне ожидания значений теневых граней выполняются вычисления в частности вычисления на внутренней области локальной секции мае с ива

Описанные технологии распараллеливания ис пользование ь при написании двух параллельных программ одна на языке С с применением библиотеки функций передачи сообщений MPI, вторая - на языке CDVM Параллельные программы запускались на кластерах MBC-1000/1G Hucniivia вычислительно моделирования СО РАН (г Красноярск) и МВС-1000/М Межведом-с твенного с уперкомпыотерного центра РАН (г Москва) Параллельные программы тестировались на модельных задачах В качестве модельной задачи выбирается с л\ чаи однородной с реды На -пой задаче сравнивались два варианта D\ М-программ с различными типами обновления тс-невых граней При

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

В диссертации приводятся полученные графики зависимости ускорения алгоритмов и эффективности распараллеливания в зависимости от количества используемых процессоров Были рассмотрены случаи мелкой (104 узлов) и крупной (103 узлов) сетки На основании полученных данных об ускорении и эффективности распараллеливания делается вывод о целесообразности применения MPI-лехнологии при распараллеливании вычислений на достаточно большом числе процессоров Данный результат согласуется с выводами сделанными в ИПМ им Келдыша РАН по сравнению различных подходов по распараллеливанию на широко известных тес тах NAS (NPB 2 3)

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

Для распределения построенной глобальной сетки на заданное количество процессоров поступаем следующим образом Если размерность сетки » одном блоке больше среднего количество ячеек, приходящихся на один процессор, то этот блок обслуживается несколькими процессорами и наоборот, один и тот же процессор обрабатывает несколько соседних блоков, если их суммарная размерность не превышает средней При разбиении блока на несколько процессоров используются 1D-, 2D- или ЗБ-разбиения Конкретный способ разбиения выбирается в зависимости от размерности сетки в блоке по всем координатным направлениям с учетом минимизации количества граничных ячеек определяемого формулой

F(ku к2, h) = {ki - 1) NiN3 + (к2 - l)NiN3 + (h - l)NxN2

где ki, /¡'2 и кз — количества разбиений облака вдоль координатных линий Предполагается что m = к\ X Аг х ^з ~ количество частей, на которое разбивается блок Таким образом, получаем что на каждый процессор распределяется некоторая часть глобальной сетки На этой части производятся все дальнейшие вычислений В п 3 3 диссертации приведено подробное сшис ание

разделения глобальной сетки на локальные для процессора с определенным номером

В п 34 описывается организация межпроцессорных обменов, осуществляемых в ходе численною решения задачи При )том введены îpn группы межпроцессорных обменов обмены впуфи одного блока первого и второго уровней и обмены между блоками Обмены внутри блока первого уровня выполняются на каждом шаге численного интегрирования на этапе предельной реконструкции и осуществляются между процессорами обрабатывающими один блок В п 3 4 1 введены понятия внешней границы локальной секции блока, распределенной на данный процессор, и внуфенней границы При внутриблочиом обмене процессор посылает соседним процессорам значения внутренней границы и принимает значении внешней границы В случае ее ли блок распределен на один процессор внутриблочных обменов первого уровня не выполняется

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

Третий вид обменов, межблочные обмены, осуществляются при осуществлении процедуры склейки решения па граничных поверхностях сое едних блоков При этом один их процессоров, обрабатывающих внутреннюю границу блока все множество граничных значений (собранное при помощи внутри-блочных обменов второго уровня) посылает процессору, обрабатывающему внутреннюю границу соседнего блока и получает от этого процессора множество значении на внешней границе Предполагается, что множество граничных значений соседнего блока также уже собрано на одном из процессоров посредством внутриблочных обменов второго } ровня

В и 3 5 дастся описание принципов создания комплекса параллельных прикладных программ, реализующих разработанный вычислительный алгоритм Программный комплекс на алгоритмическом языке Fort, ian-90 с использованием библиотеки MPI состоит из трех частей программы-препроцессора выполняющей подготовку данных и построение сетки программы-процессора выполняющую непосредственно счет и программы-постпроцессора,

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

Количество вычислительных узлов Количество вычислительных узлов

(а) (б)

Рис.. 1: Эффективность распараллеливания вычислений в случае "крупной сетки "(а) и "мелкой сетки" (б)

В п. 3.6 приводятся результаты оценки эффективности параллельной реализации алгоритма. На рис. 1 представлены графики эффективности распараллеливания от числа используемых вычислительных узлов в случае "крупной сетки" (603 расчетных узлов) и в случае "мелкой сетки" (2003 узлов сетки). Средний уровень эффективности распараллеливания составляет около 80%, что связано с неизбежными затратами времени на организацию межпроцессорных обменов и записи результатов в файл. Приведены обобщенные статистические данные по времени счета одного шага 15 зависимости от объема пересылок и от количества ячеек, в которых значения искомых величин записываются в двоичные файлы. Под объемом пересылок понимается число граничных ячеек сеточной области, через которые осуществляются межпроцессорные обмены. Данные собраны на основании расчетов серии модельных задач па разном количестве процессоров (от 2 до 32) при вычислительной нагрузке на один процессор в 106 ячеек. Результаты показывают, что время счета одного шага колеблется от 1 до 3 минут в зависимости от объема пересылок и количества записей в файл. Увеличение последних, естественным образом ведет к увеличению времени счета.

Четвертая глава посвящена результатам численною исследования разработанного вычш лительпого алгоритма на ряде модельных задач В и 4 1 обсуждаются основные принципы проведения вычислительного эксперимента Во всех предложенных моделях исходная среда ограничена прямоугольной областью для облегчения постановки искусственных граничных условии Исходная нагрузка действует по нормали к верхней грани области Функция источника (син)сообразныи импутье) имеет вид Р(/) = Р<)Ып(и;/), если I < Т, где Р(1 - амплитуда и ш - частота генерируемых волн в среде Верхняя часть области за исключением зоны приложения импульсной нагрузки, считается свободной от напряжения На остальных гранях ставятся неотражающие условия, которые соответствуют бесконечной протяженности массива и формируются для одномерных систем с помощью уравнений на характеристиках

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

В и 4 2 рассматривается задача Дамба о действии сосредоточенной нагрузки на упругое полупространство Рассматриваемая ^ пругая среда ограничена прямоугольной областью АВСОА^С^ (рис 2) Па-

раметры среды /> = 7850 кг/м < р — G000 и cs = 3210 м/с Исходная импульсная нагрузка действует в центре верхней грани области Время действия нагрузки Т = 0 005 с Величина амплитуды колебаний Р0 = 250 МПа Ри( 2 Дрисише сосредоточенной

Размеры расчетной области 200x400 x400 пл упругую

нуго область

м Расчеты проводились на сетке из 200 х 400 X 400 узлов IHai по времени выбирался в соответствии с условием К>рапга-Фридрнхса-Леви и составил 0 0016 с Среднее время расчета одного шага по времени на 32-х вычш лшельных узлах кластера МВС-1000 - 2 мин\ты

Рис. 3: Действие сосредоточенной нагрузки па упругую прямоугольную область

На рис. 3 представлена визуализация щ в фиксированный момент времени 0.016 с. Здесь введены следующие цифровые обозначения: 1 - продольная волна, 2 - поперечная волна, 3 - головная поперечная волна, 4 - поверхностная волна Релея. Результаты соответствуют точному решению задачи для случая плоского деформированного состояния: видны продольная и поперечная волны, менее интенсивные головные волны и волны Релея, которым па рисунках соответствуют движущиеся точки па поверхности х\ = 0. На рис. 4 приведено волновое поле нормального напряжения.

Ш

: i

..............-JLI.1

щыямжли-

"Л

шв . i

а®

щвщ

Рис. 4: Поверхности уровня напряжения и,, в моменты времени 0.02, 0.03 и 0.04 с

<4L

Ч

»X.

у'

/

1х.

л/

(а)

(б) (в) (г)

Рис. 5: Модели с различными границами раздела сред: а горизонтальная, 5 - вертикальная, в - наклонная, г - граница со сбросом

В и. 4.3 представлены задачи о прохождении волн через различные границы раздела сред (рис. 5). Рассматривается горизонтальная (а), вертикальная

(б) и наклонная (в) типы границы раздела сред Также рассматривается среда с границей, имеющей небольшой вертикальный сброс (г) Во всех задачах параметры среды подобраны с учетом условия разрыва акустчегкои жест-кос ти па границе раздела сред, что обеспечивает образование отраженных и преломленных волн при переходе падающих волн через поверхность раздела Результаты расчетов представлены в виде фрагментов визуализаций вертикальной компоненты вектора перемещения и синтетических сейсмограмм Рас четы показывают образование вс ех типов отраженных и преломленных волн, и также головных продольных волн Для этих типов вол приведены уравнения теоретических годографов Рассчшываюгся время и точка выхода отраженных волн и головных волн на дневную поверхность области Полученные численные решения полностью соответствуют теоретическим Кроме того, в задаче с вертикальным сбросом в решении явно видны образовавшиеся дифрагированные волны от отражения падающих волн от угла сброса На рис 6 представлены фрагменты визуализации вертикальной компоненты вектора перемещения вертикальных профилей, секущих сброс под различными углами 0° (а), 11° (б) 26 5" (в), 37° (г) и проходящими через источник импульс а Здесь введены с ледующие обозначения 1 - падающая Р-волна, 2 -отраженная РР-волна, 3 - отраженная РЭ-волна, 4 - преломленная РР-волна 5 - преломленная РБ-волна, 6,7- дифрагированные волны Как видно и з рисунков, с увеличением угла, под которым вертикальная плоскость пересекает линию сброса, сферический фронт дифрагированной волны переходит в эллиптический Это связано с тем, что каждая точка грани сброса порождает конус лучей, не лежащих в одной плоскости Полученные численные результаты согласуются с результатами, полученными в Институте вычислительной математики и геофизики СО РАН

В заключении кратко перечислены основные результаты диссертационной работы, выносимой на защиту

В двух последующих приложениях содержатся иллюстрации к диссертационной работе В Приложении 1 представлены блок-схемы программного модуля 'процессор'' и его фрагмента, относящегося к алгоритму численного решения задачи В Приложении 2 представлены иллюстрации к главе 4 с результатами численных расчетов

Заключение

В диссертационной работе получены следующие результаты

1 Предложен эффективный параллельный вычисли!ельный алгоритм решения динамических задач моделирования сейсмических волновых полей в кусочно-однородных упругих средах, основанный на расщеплении пространственной задачи теории упругости на серию одномерных задач

2 Ал) оритм реализован в виде программного комплекса на я зыке программирования Fortran-90 с использованием функции библиотеки передачи с ообщений MPI Техноло! ия распараллеливания вычислении основана на разделении исходной области по процессорам исходя из требований равномерной загрузки Программный комплекс позволяет получать результаты которые могут быть использованы в прикладных геофизических программах, таких как Ргопшх или SeisView для анализа спнтети-чес ких сейсмограмм

3 Проведены чис ленные расчеты модельных задач на многопроцессорной вычислительной системе МВС-1000 Институт вычислительного моделирования СО РАН Результаты расчетов подтверждаются теоретическими выкладками, основанными на законах геометрической оптики, и согласуются с результатами, полученными в Инг титуле вычислительной математики и математической геофизики СО РАН

Список публикаций по теме диссертации

Статьи в ведущих научных журналах, включенных в перечень ВАК

1 Кучунова, ЕВ Вычислительный алгоритм для расчета волновых полей в блочных средах на многопроцессорных вычис лительных системах /ЕВ Кучунова В М Садовский // Журнал Сибирского федерального университета Математика и физика Journal of Siberian Federal University Matheinetics & Physics - Красноярск, 2008 - N2 -С 210-220

2 Кучунова, E В Численное исследование распространения сейсмических волн в блочных средах на многопроцессорных вычис лительных системах /ЕВ Кучунова В М Садовский // Вычисли 1епьные методы и программирование 2008 - Т9 - N1 С 70-80

Другие публикации

3 Кучунова, ЕВ DVM и MPI-технологии распараллеливания в одномерной динамическом задаче теории упругости /ЕВ Кучунова // Тезисы

докладов V Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям - Новосибирск, 2004 - С 24

4 Кучунови, Е В Комплекс прикладных программ для численного решения пространственных задач динамической теории упругости на многопроцессорных вычислительных с истемах /ЕВ Кучуиова, О В Садовская, ВМ Садовский // Материалы IV межрес иопальной школы-семинара "Распределенные и кластерные вычис пения" - Красноярск 2005г - С 159-172

5 Кучунова, ЕВ Вычислительный алгоритм для расчета ynpyi их волн на многопроцессорных вычислительных системах /ЕВ Кучунона // Материалы VI Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям - Кемерово 2005г - С 22

6 Кучуиова, Е В Вычислительный алгоритм для расчета упругих волн в блочной среде на многопроцессорных вычислительных с истемах /ЕВ Ку чупова // Труды XXVIII Конференции молодых ученых ММФ МГУ им М В Ломоносова - Москва 2006г - С 100-104

7 Кучуиова, Е В Параллельный алгоритм расчета упругих волн в блочной среде /ЕВ Кучуиова // Материалы VI межре1 иональной школы-семинара "Распределенные и кластерные вычисления" - Красноярск, 2006г - С 79-89

8 Кучуиова, Е В Параллельная реализация алгоритма для расчета упругих волн в блочной среде /ЕВ Кучуиова, В М Садовский // Труды Международной научной конференции Параллельные вычислительные технологии (ПаВТ'2007) - Челябинск Изд-во ЮУрГУ, 2007 -Т2 -С 37-43

9 Кучуиова, ЕВ Численное решение прямой динамической задачи динамической теории упругости на многопроцессорных вычислительных системах / ЕВ Кучунова // Сборник трудов XX междупарод копе}) ММТТ-20 в Ют 'Г 1 секция 1 / под ред В С Балакирева Изд Яросл гос техн уп- г, Ярославль 2007 С 161-165

10 Кучуиова, ЕВ Численное исследование распространения сейсмических волн в блочных средах на многопроцессорных вычислительных системах /ЕВ Кучуиова, В М Садовский // Труды Международной научной конференции "Параллельные Вычислительные Тсхнологии"(Санм-Петербург 28 января - 1 февраля 2008г) - Челябинс к Изд ЮУрГУ, 2008 - С 130-141

Подписано и псч<иь Форма! 00 х 841 /! Печать оперативная Ус 71 поч Л //V Тираж 100 ж? Заказ № ЦОЪ"

Ошсчаымо и ПИК СФУ 660041 Красноярск, пр Свободным, 79

Введение

I Современное состояние проблем, связанных с численным решением задач теории упругости на многопроцессорных вычислительных системах

1.1 Обзор основных методов решения прямых задач сейсмики.,

1.2 Обзор численных методов решения задач динамической теории упругости.

1.3 Обзор методов моделирования неотражающих граничных условий

1.4 Обзор основных тенденций в области параллельных вычислений

II Численное решение задачи о распространении упругих волн в блочной среде

2.1 Постановка пространственной задачи.

2.2 Построение сетки в области.

2.2.1 Разбиение области на блоки.

2.2.2 Нахождение направляющих векторов.

2.2.3 Построение сетки в блоке.

2.3 Численный метод решения задачи.

2.3.1 Метод двуциклического расщепления.

2.3.2 Решение одномерных систем

2.3.3 Предельная реконструкция инвариантов решения

2.3.4 Склейка решения на внутренней границе двух блоков

III Технологии реализации параллельного алгоритма на многопроцессорных вычислительных системах

3.1 Проблема эффективной реализации алгоритмов па многопроцессорных вычислительных системах.

3.2 Сравнение MPI- и DVM-технологий распараллеливания на примере решения одномерной задачи теории упругости.

3.3 Технология распределения данных по процессорам

3.4 Обмен граничными значениями.'

3.4.1 Обмен значениями внутри одного блока.

3.4.2 Обмен значениями между блоками.

3.5 Программный комплекс.

3.5.1 Основные идеи программного комплекса

3.5.2 Параллельная реализация алгоритма построения расчетной сетки.

3.5.3 Параллельная реализация алгоритма численного решения задачи

3.5.4 Методы визуализации решения.

3.6 Эффективность параллельной реализации алгоритмов

IV Результаты исследования на модельных задачах

4.1 Вычислительный эксперимент

4.1.1 Модельная задача 1. Действие сосредоточенной нагрузки на упругое полупространство.

4.2 Прохождение волн через различные границы раздела сред

4.2.1 Модельная задача 2. Прохождение волн через наклонную границу раздела двух сред.

4.2.2 Модельная задача 3. Прохождение волн через горизонтальную границу раздела двух сред.

4.2.3 Модельная задача 5. Прохождение волн через границу раздела двух сред со сбросом.

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

Начиная с работ Релея и Лэмба, математическое моделирование распространения сейсмических волн проводится на основе системы уравнений динамической теории упругости. На практике важно проанализировать сейсмическое волновое поле с учетом всех типов волн, волн интерференционного характера и возникающих при этом "нелучевых" явлений. Совершенствование систем наблюдения и повышение детализации экспериментальной информации приводит к необходимости рассмотрения волновых полей с большей точностью и для все более возрастающих времен. Исключительное значение для исследования распространения сейсмических волн в сложнопостроенных средах имеет вычислительный эксперимент. Однако применение прямых методов численного анализа требует большого объема вычислительных ресурсов. Решение пространственных задач на современном персональном компьютере может оказаться невозможным из-за длительности времени счета. Кроме того, однопроцессорный компьютер зачастую не может обеспечить требуемый объем оперативной памяти для работы с многомиллионными ячейками расчетных сеток, использование которых необходимо для достижения требуемой точности. Эффективным средством решения ресурсоемких пространственных задач является технология параллельных вычислений па распределенных системах кластерного типа. Развитие суперкомпьютеров позволяет достаточно успешно решать задачи в трехмерной постановке.

На пути перехода от персонального компьютера к суперкомпьютеру с параллельной архитектурой имеются определенные трудности, которые, во-первых, связаны с необходимостью распараллеливания вычислительного алгоритма, во-вторых, с принципиально более сложным написанием программного кода и недостаточно развитой системой отладки программ на кластере. При этом одним из факторов эффективного использования высокопроизводительной многопроцессорной вычислительной техники является применение специальных методов математического моделирования. В связи с этим, одной из задач современной науки является развитие технологии математического моделирования, рассматривающей многопроцессорные вычислительные системы (МВС) как основной аппарат для вычислений. Отметим, что методы моделирования с использованием высокопроизводительных компьютеров часто оказываются единственным способом получения информации об исследуемых явлениях, так как возможности теоретического и экспериментального изучения этих явлений ограничены методическими и техническими трудностями.

Отметим далее, что за время существования МВС накоплен значительный опыт их использования, и они находятся в массовом доступе. Однако, огромные вычислительные возможности, предоставляемые этими системами, используются недостаточно. Причиной этому является сложность адаптации последовательных алгоритмов и программ к параллельным архитектурам. Поэтому разработка эффективных параллельных алгоритмов решения фундаментальных и прикладных задач, изначально ориентированных на использование МВС является чрезвычайно актуальной проблемой. При этом надо учитывать, что вместе с возможностями, которые предоставляет программистам многопроцессорная вычислительная техника, существует специфические требования к используемым алгоритмам [1]:

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

• ускорение расчетов при увеличении числа процессоров;

• необходимость синхронизации при остановке (завершении) итераций;

• снижение затрат времени на синхронизацию;

• периодический сбор данных для записи и анализа;

• обмен данными между процессорами для продолжения расчета;

• снижение объема передаваемой информации при обмене;

• снижение числа актов приема-передачи (то есть числа обменов между процессорами);

• уменьшение влияния размера задачи на ускорение;

• обеспечение масштабируемости и переносимости программы.

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

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

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

• разработка параллельного вычислительного алгоритма расчета волновых полей в трехмерных кусочно-однородных упругих средах;

• адаптация разработанного алгоритма к многопроцессорным вычислительным системам, которая складывается из выбора методов реализации алгоритма на параллельных ЭВМ, выявления трудно распараллеливаемых блоков, их оптимизации и верификации;

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

Объектом исследования диссертации являются трехмерные прямые динамические задачи сейсмики для упругих кусочно-однородных изотропных сред.

В качестве метода исследований используется вычислительный эксперимент, включающий в себя следующие этапы [2]:

• математическая формулировка задачи,

• построение приближенного (численного) метода решения задачи, написание вычислительного алгоритма;

• программирование на ЭВМ вычислительного алгоритма,

• проведение расчетов на ЭВМ,

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

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

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

4.1 Вычислительный эксперимент

Вычислительный алгоритм, описанный в главе 2, реализован комплексом параллельных программ на алгоритмическом языке Fortran90 с использованием библиотеки MPI (message passing interface). Реализация параллельного алгоритма была основана на основных идеях, представленных в главе 3. В данной главе представлены результаты численных тестов для различных моделей геофизических сред.

Во всех представленных моделях исходная нагрузка действует нормально на верхнюю грань области в некоторой точке Хо. Функция источника (си-пусообразный импульс) имеет вид: P(t) = Pq sin (аЛ), если t < Т, где Pq -амплитуда и ш - частота генерируемых волн в среде. Верхняя грань области, за исключением зоны приложения импульсной нагрузки, считается свободной от напряжений. На остальных гранях поставлены неотражающие условия, которые соответствуют бесконечной протяженности массива и формулируются для одномерных систем с помощью уравнений на характеристиках.

Расчеты производились на суперкомпьютере МВС-1000/16 и МВС-1000 Института Вычислительного Моделирования СО РАН. Суперкомпьютер МВС-1000/16 введен в эксплуатацию в начале 2002 года. Общая производительность МВС-1000/16 составляет 14 млрд. операций в секунду, скорость пересылок - ЮОМЬ в секунду. Суперкомпьютер МВС-1000, введенный в эксплуатацию в 2005 году, имеет пиковую производительность 220.4 млрд. операций в секунду и производительность по тесту LinPack 134 ГФлоп.

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

4.1.1 Модельная задача 1. Действие сосредоточенной нагрузки на упругое полупространство

Изотропная упругая среда занимает прямоугольную область ft = {х G R3\ 0 < xi < Я, -Н < х2 < Н, -Н <х3< Н} о в декартовой системе координат. Параметры среды: р = 7850 кг/м , ср — 6000 и cs = 3210 м/с. В точке хо = (0,0,0) приложена точечная импульсная нагрузка, действующая перпендикулярно к плоскости xi = 0. Размеры расчетной области1: 200 X 400 X 400 м. Расчеты проводились на сетке из 200 X 400 х 400 узлов. Шаг по времени выбирался в соответствии с услови

1 Здесь и далее размеры области перечисляются в следующем порядке: первая цифра - по xi, вторая -но Х2, третья - по х3. ем Куранта-Фридрихса - Л еви и составил 0.0016 с. Время действия нагрузки Т — 0.005 с, то есть около четырех шагов по времени. Величина амплитуды колебаний Pq — 250 МПа. Среднее время расчета одного шага по времени на 32-х процессорах кластера MB С-1000 - 2 минуты. Количество расчетных ячеек, приходящихся на один процессор, равно 106.

Рассмотрим полученную волновую картину в вертикальном сечении области, проходящем через точку приложения нагрузки хц. На рис. 3 приложения Б приведена сейсмограмма вертикальной компоненты вектора перемещения вдоль профиля, проходящего через точку Х(). На рис. 4 приложения Б представлена визуализация вертикальной компоненты вектора перемещения вдоль вертикального сечения области в фиксированный момент времени 0.016 с. Здесь введены следующие цифровые обозначения: 1 - продольная волна, 2 - поперечная волна, 3 - головная поперечная волна, 4 - поверхностная волна Релея.

Результаты соответствуют точному решению задачи для случая плоского деформированного состояния: видны продольная и поперечная волны, менее интенсивные головные волны и волны Релея, которым на рисунках соответствуют движущиеся точки на поверхности х\ = 0. На рис. 4.1 приведено волновое поле нормального напряжения.

Рис. 4.1: Поверхности уровня напряжения в моменты времени 0.02, 0.03 и 0.04 с

4.2 Прохождение волн черёз различные границы раздела сред

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

4.2.1 Модельная задача 2. Прохождение волн через наклонную границу раздела двух сред

Рассмотрим модель двухслойной среды с однородным изотропным верхним слоем и скачком акустической жесткости на границе с подстилающим полупространством. Под однородной покрывающей средой с параметрами р = 500кг/м , ср=500 и cs=300 м/с распложена вторая среда с параметрами ро = 750 кг/м3, Ср—750 и с°=450 м/с (рис. 4.2.а). Горизонт, разделяющий две среды, представляет собой плоскость, наклоненную по линии АС под углом ipi = 5°, а по линии DB - под углом ip2 = 7°. Заметим, что на границе раздела сред выполняются условия рср ф роср и Pcs ф Pocs: гарантирующие образование отраженных волн.

Размеры расчетной области - АВ — AD = 600, АА2 = 400 м. Глубина залегания наклонного горизонта: АА\ — 175, ВВ\ = 87,5, СС\ = 100, DD\ = 187, 5 м. Внешний импульс на среду действует в точке О. Глубина залегания границы раздела в точке О составляет 00\ = Н = 156,25 м.

Рис. 4.2: Модель области с наклонным горизонтом раздела двух сред: а - геометрия среды, б - расположение вертикальных профилей в горизонтальной плоскости

Величина приложенной нагрузки - Pq=250 МП а, время действия - 0,01 с, радиус действия - 3.5 м. Число узлов сетки - 200 х 200 х 200, временной шаг - т — 0.002 с. Среднее время расчета одного шага на 8-ми процессорах составило 2.63 минуты. Количество расчетных ячеек, приходящихся на один процессор - 106.

На рис. 5-7 приложения Б представлены сейсмограммы вертикальной составляющей вектора перемещения вдоль профилей ОС, ОМ, DB, EF и OQ. Здесь введены следующие обозначения: 1 - падающая Р-волна; 2 - падающая S-волна; 3 - отраженная РР-волна; 4 - отраженная PS-волна; 5 - преломленная РР-волна; 6 - преломленная PS-волна; 7 - головная волна, образованная РР-волиой; 8 - отраженная SS-волна; 9 - отраженная SP-волна.

Исходя из законов геометрической оптики, годографы прямой Р-волны вдоль профилей AC, MN и PQ: проходящих через источник импульса, представляют собой прямые линии: ti(г) = г/ср, где г - расстоянии от точки на профиле до источника волн. Годографы Р-волны вдоль профилей DB и EF, не проходящих через источник импульса, описываются уравнением параболы, симметричной относительно центра сейсмограммы: t(r) = \Д2 + г2/ср, где г - расстояние от точки профиля до центра, I - расстояние от центра профиля до источника импульса (I — OL для линии EF, I = О К для линии DB). Картина распространения падающей поперечной волны (S-волпы) полностью эквивалентна картине продольной волны и отличается только скоростью. Годографы S-волны вдоль профилей АС, NM и PQ имеют вид: t(r) = r/cs, вдоль профилей DB и EF: ^г(^) — у/12 + r2/cs.

Годографы отраженной РР-волиы вдоль профилей АС, NM и PQ, характеризуются уравнением: t(r) = — л/г2 - 2Hr sin 2<р + 4Я2 cos2 ip, (4.1)

Ср где ip - угол наклона горизонта (для профиля АС угол ср равен 5°, MN -ср = 8,56°, MN - ip = —1,19°). Годографы (4.1) представляют собой параболу, сдвинутую относительно центра сейсмограммы в сторону подъема горизонта на расстояние г\ = iJ sin время прихода отраженной волны ti = 2Н cos (p/ср. Эта точка отмечена на сейсмограммах как Р\. Для профиля АС: г\ = 27 м, t\ — 0.623 с, для профиля MN: ri=50 м, t\ = 0.618 с, для PQ: п -7, 3 м, ti = 0.625 с.

Вдоль профилей DB и EF уравнения годографов отраженной волны усложняется тем, что угол горизонта зависит от г, расстояния до центра сейсмограммы: ч tan + (г/Л tancpo ip(r) = arctan- , w ; — (4.2) y/1 + r2/l2 K J

Годографы отраженной волны вдоль профилей DB и EF: t(r) = -\fr2 + l2- 2Hy/l2 + г2 sin 2<p(r) + 4Я2 cos2 <p(r).

Ср

На рис. 4.3 представлены годографы отраженной РР-волны для профилей EF и DB.

Рис. 4.3: Годограф отраженной Р Р-волны вдоль профилей EF (а) и DB (б)

Годографы отраженной SS-волны вдоль профилей AC, NM и PQ описываются уравнениями t(r) = — \Jr2 - 2Яг sin 2tp + 4Я2 cos2 (р, cs где (р - угол наклона горизонта. Отраженная волна смещена относительно центра на расстояние r2 = Я sin 2с/?, время прихода отраженной волны t2 = 2Нcos(p/cs. Для профиля АС: г2=27 м, t2 = 1, 04 с; для MN: Г2=50 м, ^2=1,'03 с; для PQ: г2=-7.3 м, t2 = 1,04 с. Годографы отраженной SS-волны вдоль профилей DB и EF: t(r) = -^r2 + l2- 2Ну/12 + г2 sin 2ip(r) + 4Я2 cos2 сp{r), cs где угол <р(г) определяется по формуле (4.2).

При переходе продольной Р-волны через границу раздела двух сред кроме отраженной продольной волны образуется вторая отраженная волна - поперечная PS-волна. Аналогично при прохождении поперечной S-волны кроме отраженной поперечной SS-волны образуется отраженная продольная SP-волна. Эти волны обозначены на сейсмограммах цифрами 4 и 9 соответственно.

Угол отражения поперечной PS-волны (/3) связан с углом падения Р-волны (а) по закону Снеллиуса: sin в = — sin а. скорость поперечных волн в среде (cs) меньше скорости продольных (ср), поэтому (3 < а. Следовательно, при любом угле падения а (0 < а < отраженная волна выйдет на верхнюю границу среды на расстоянии: ч Я (tan а + tan/З) ,, . г(а) = -—, (4.3)

1 + tan <р tan а где ip - угол наклона горизонта. Время прихода отраженной Р5-волны определяется: t(a) -- Н

Cg tan а + Ср tan /3 cpc2s sin а(1 + tan <р tan а)' где ip - угол наклона горизонта, для профилей DB и EF определяется соотношением (4.2). Годограф отраженной Р5-волны представляет собой параболу не симметричную относительно центра сейсмограммы.

Угол отражения продольной 5Р-волпы (5) связан с углом падения S-волны (7) соотношением Снеллиуса: г ср • sm д = — sin 7.

Cs

В этом случае существует критический угол 7* = arcsin — = 36, 87° при котоср ром угол отражения становится максимальным (90°). При 7 < 7* расстояние выхода отраженной 5Р-волны определяется аналогично (4.3). Годограф отраженный 5Р-волны 7 < 7* описывается уравнением: cl tan 7 + tan 5 t(j) = Я£--, (4.4)

CpCs sin 7(1 ~Ь tan (р tan 7)' где (р - угол наклона горизонта, для профилей DB и EF определяется соотношением (4.2).

Критический угол г, при падении под которым луч продольной Р-волны преломляется под углом в 90° к нормали, определяется по закону Снеллиуса: sin г = Ср/Ср. При этом преломленная продольная волна начинает распространяться вдоль границы раздела сред. Фронт проходящей продольной

РР-волны, скользя вдоль границы раздела, возбуждает в верхнем слое колебания, которые вызывают появление головной преломленной волны. Одной стороной фронт головной волны касается фронта отраженной из критической точки волны, другой примыкает к фронту скользящей преломленной волны. На рис. 6 приложения Б представлены фрагменты визуализации вертикальной компоненты вектора перемещения для случаев падения продольной Р-волны под углом, меньшим критического угла i (а) и большим критического (б). При этом на рис. 6 (б) прослеживается образование головной продольной волны. Ее годографы для профилей АС, MN и PQ имеет вид: t(r) = — (г sm(i <р) + 2Н cos <р cos г), (4.5)

Ср где знак "—" берется для поднимающегося горизонта, а знак "+" - для опускающегося. Точка выхода головной волны на поверхность (точка Р2 на сейсмограмме) расположена от источника на расстоянии: гз = 2Н cos <рCos(i-p); время выхода головной волны на поверхность составляет: ts = с^оф-р) • Для профиля АС: гз=261 м, £3=0.776 с; для профиля MN: г3=246 м, £3=0.726 с; для профиля PQ: гз=287 м, ^з=0.866 с. Для профилей DB и EF годографы продольной головной волны аналогичны (4.5) с заменой г на л/12 + г2 и ip -на <р(г) по формуле (4.2): t(r) = — (\//2 + г2 sin(z =f ip) + 2Н cos ip cos г).

Ср

Из этого уравнения можно получить, что головная волна вдоль линии DB (.EF) выходит на поверхность в двух точках: на расстоянии 84 м (225 м) слева от центра сейсмограммы и на расстоянии 45 м (183 м) справа. Расстояние выхода головной продольной волны: 2Н cos (р(ф) sin г

-Г--ГмГ"' cos (г — фуф)} где ф - угол подъема горизонта: р(ф) = arctan tan (f \ + tan ф ■ tan y/l + tan2 ф

Картина распространения головной поперечной SS-волны аналогична. Уравнение годографа поперечной головной волны имеет вид: t = — (г sin(i ip) + 2Н cos <р cos г). cs

Точка выхода (точка на сейсмограмме) на поверхность расположена от источника на расстоянии Г4 = 2iiTcos(^i sin г/соз(г—ipi), время выхода составляет = 2Н cos2 ipi/(cs cos(i — ip)). Для профиля AC: 7-4 = 261 м, = 1.475 с. Головная поперечная волна выходит на поверхность в тех же точках, что и головная продольная волна, отличается только время.

На рис. 4.4 кривой линией представлена зависимость расстояния от точки приложения импульса (т. О) до точки выхода головной продольной волны на поверхность в зависимости от угла подъема горизонта (ф).

X,

Рис. 4.4: График зависимости расстояния до точки выхода головной продольной волны на поверхность в зависимости от угла подъема горизонта ф

4.2.2 Модельная задача 3. Прохождение волн через горизонтальную границу раздела двух сред

Рассмотрим модель двухслойной среды (аналогичную модельной задачи 2) с горизонтальной границей раздела сред (рис. 4.5.а). Параметры верхнего слоя: р=500 кг/м3, ср—500 и cs=300 м/с, параметры нижнего слоя: ро—750 кг/м3, Ср=750 и с®—450 м/с. Поверхность раздела двух сред представляет собой горизонтальную плоскость па некоторой глубине Н от верхней границы. Выполняются условия неравенства акустической жесткости двух сред: рср ф р0Ср и pcs ф р0с°3.

Рис. 4.5: Модель области с горизонтальной поверхностью раздела двух сред: а - геометрия среды, б - расположение вертикальных профилей в горизонтальной плоскости Х2Х3

Размеры расчетной области — АВ = AD = 400 м, АА2 = 200 м. Глубина залегания горизонта: АА\ = Н — 100 м. Внешний импульс действует в точке О. Величина приложенной нагрузки — Pq=250 МПА, время действия — 0.01 с, радиус действия — 2,5 м. Число узлов сетки — 200 х 200 х 200, временной шаг — т=0,002 с. Среднее время расчета одного шага на 8-ми процессорах — 1.26 минуты. Количество расчетных ячеек, приходящихся на один процессор

На рис. 8-9 приложения Б представлены сейсмограммы вертикальной сох. х

106 ставляющей вектора перемещения вдоль профилей ОС, ON, EF, DB и QN. На рис. 10-11 приложения Б представлены визуализации вертикальной компонент вектора перемещения вдоль профилей AC, MN и DB в момент времени 0.4 с. Задача является симметричной относительно оси АС.

Здесь введены следующие числовые обозначения: 1 - падающая Р-волна; 2 - падающая S-волна; 3 - поверхностная волна Релея; 4 - отраженная РР-волна; 5 - отраженная Рб'-волна; 6 - преломленная РР-волна; 7 - преломленная Рб'-волна; 8 - головная волна, образованная РР-волной; 9 - вторичные отраженные волны от верхней границы области; 10 - отраженная б'б'-волна; 11 - отраженная б'Р-волна; 12 - преломленная 5Р-волна; 13 - преломленная 55-волна; 14 - головная волна, образованная SP-волной; 15 - головная волна, образованная 55-волной.

Годографы прямых Р- и S-волн описываются аналогичными уравнениями, представленными в описании задачи 2. Годографы отраженных и головных волн упрощаются за счет того, что поверхность раздела сред является горизонтальной. Годографы отраженной РР-волны вдоль профилей АС и MN описываются уравнением параболы, симметричной относительно положения источника импульса (точка О): t — \/г2 + 4Н2/ср. Вдоль профилей EF, DB и QN годографы отраженной РР-волны характеризуются уравнением: t = у/г2 +12 + 4Я2/ср, где I - расстояние от центра профиля до источника импульса (/ = OL для профиля EF, I — О К для профиля DB, I = OG для профиля QN). Уравнение годографа отраженной 55-волны вдоль профилей ОС, ON: t = у/г2 + 4H2/cs, и вдоль линий EF, DB и QN: t = y/x2 + l2 + 4H2/cs.

Годограф продольной головной волны для профилей AC, MN и PQ описывается уравнением: t = (г sin г + 2Н cos i)/cp, где i - критический угол угол под которым продольная волна преломляется под углом 90° к нормали). Согласно формуле Снеллиуса величина критического угла составит г — arcsin ^ = 42°. Головная волна выйдет на поверхность (точка на сейсмограммах помечена Р) на расстоянии от источника волны, равным хя = 2Н sin г/ cosi= 180 м, вместе с отраженной волной. Однако далее головная волна начинает опережать отраженную волну.

Модельная задача 4. Прохождение волн через вертикальную границу раздела двух сред

Рассмотрим модель среды, состоящей из двух однородных изотропных блоков с вертикальной границей между ними и скачком акустической жесткости на вертикальной границе (рис. 4.6). Параметры одного блока: ^=500 кг/м3, ср=500 и ся=300 м/с, параметры второго блока: ро=750 кг/м3, c^^SO и Сд=450 м/с.

Рис. 4.6: Модель области с вертикальной поверхностью раздела двух сред

Размеры расчетной области - АВ = AD — 300 м, АА\ = 200, АР = BQ = I = 200 м. Внешний импульс действует в точке О. Величина приложенной нагрузки - Pq-250 МПА, время действия нагрузки - 0.01 с, радиус действия нагрузки - 1,5 м. Число узлов сетки - 100 х 300 х 300, временной шаг - т=0.01 с. Среднее время расчета одного шага по времени на 9-ти процессорах -1.13 минуты. Количество расчетных узлов на один процессор - 106.

На рис. 13 приложения Б представлены сейсмограммы вертикальной составляющей вектора перемещения вдоль профилей ОС и О2С2. На рис.12 приложения Б представлена визуализация компонентов вектора перемещения вдоль профиля АС в момент времени 0.195 с. Здесь введены следующие числовые обозначения: 1 - падающая Р-волна; 2 - падающая S-волна; 3 -поверхностная волна Релея; 4 - отраженная продольная РР-волна; 5 - отраженная поперечная PS-волна; 6 - преломленная продольная РР-волна; 7 - преломленная поперечная Р5-волна; 8 - отраженная поперечная волна 55-волна; 9 - отраженная продольная волна 5 Р-волна; 10 - преломленная поперечная б'З'-волна; 11 - преломленная продольная волна 5 Р-волна.

Вдоль профиля АС годографы всех волн являются прямыми линиями (рис. 13.а прил. Б): t = г/ср (Р-волна), t = r/cs (S-волна), t = 1/ср + г/сРр (отраженная РР-волна), t = l/cp + r/c® (отраженная Р5-волна), t = l/cs г/с® (отраженная S Р-волна), t = l/cs + r/c® (отраженная б'б'-волна), где I - расстояние от центра импульса до границы раздела сред (I = OR).

Отраженные РР— и SP—волны распространяются с одинаковой скоростью (с®). Аналогично PS— и SS—волны соответственно со скоростью с®. Поэтому линии годографов для этих пар воли представляют собой параллельные линии. На сейсмограмме хорошо просматриваются преломленные РР- и З^-волпы, в то время как преломленные PS- и 5Р-волны выражены крайне слабо. Большая часть сейсмической энергии при пересечении границы переходи в.волну той же поляризации, какая у падающей волны. Отраженных волн на сейсмограмме практически не видно.

Вдоль профиля О2С2, расположенного на глубине 75м, отраженные волны прослеживаются более явно (6,7,10,11). Кроме того, более четко прослеживается преломленная поперечная Р5-волна.

4.2.3 Модельная задача 5. Прохождение волн через границу раздела двух сред со сбросом

Рассмотрим модель трехмерного сброса (рис. 4.7, а) с амплитудой L = 15 м. Параметры материала верхнего слоя: р=500 кг /м3, ср=500 и cs=200 м/с, подстилающего слоя: pq—750 кг/м3, с^=750 и с^=300 м/с. Размеры расчетной области - АВ = AD = 750, АА2 = 600 м. Глубина залегания поверхности раздела материалов: АА\ — DD\ — 200, ВВ\ = СС\ = 185 м. Сброс проходит по линии ran. Внешний импульс на среду действует в точке О. Величина приложенной нагрузки - Ро—250 МПа, время действия - 0,01 с, радиус действия - 3.5 м. профилей в горизонтальной плоскости Х2Х3

Число узлов сетки - 300 х 300 х 300, временной шаг - т = 0.002 с. Число используемых процессоров - 27, количество расчетных ячеек, приходящихся на один процессор - 106. Число шагов по времени - 1000, среднее время расчета одного шага - 2.1 минуты.

Полученное в ходе расчетов волновое поле представлено на серии вертикальных плоскостей, секущих сброс под различными углами: 0° (EF), 11° (GH), 26.5° (KL), 37° (PQ), 45° (АС) и проходящими через источник импульса.

На рис. 14, 15 приложения Б представлены сейсмограммы вертикальной составляющей вектора перемещения вдоль профилей OF (а), ОН (б), OL (в), OQ (г) и ОС (д). На рис. 16 приложения Б представлены фрагменты визуализации вертикальной компоненты вектора перемещения для тех же профилей. Здесь введены следующие обозначения: 1 - падающая Р-волна; 2 - падающая S-волна; 3 - отраженная РР-волна; 4 - отраженная PS-волна; 5 - преломленная РР-волна; 6 - преломленная PS-волна; 7, 8 - дифрагированные волны.

Как видно из рисунков, с увеличением угла, под которым вертикальная плоскость пересекает линию сброса (nm), сферический фронт дифрагированной волны переходит в эллиптический. Это связано с тем, что каждая точка грани сброса порождает конус лучей, не лежащих в одной плоскости.

Заключение

В диссертационной работе получены следующие результаты:

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

2. Алгоритм реализован в виде программного комплекса на языке программирования Fortran-90 с использованием функций библиотеки передачи сообщений MPI. Технология распараллеливания вычислений основана на разделении исходной области по процессорам, исходя из требований равномерной загрузки. Программный комплекс позволяет получать результаты, которые могут быть использованы в прикладных геофизических программах, таких как Promax или SeisView для анализа синтетических сейсмограмм.

3. Проведены численные расчеты модельных задач на многопроцессорной вычислительной системе МВС-1000 Института вычислительного моделирования СО РАН. Результаты расчетов подтверждаются теоретическими выкладками, основанными на законах геометрической оптики, и согласуются с результатами, полученными в Институте вычислительной математики и математической геофизики СО РАН.

1. Якобский, М.В. Распределенные системы и сети / М.В. Якобский. - М. : МГТУ "Станкин", 2000. - 118с.

2. Самарский, А. А. Теория разностных схем: учебное пособие / А.А. Самарский. М. : Наука : Главная редакция физико-математической литературы, 1983. - 616с.

3. Domain decomposition method for partial differential equations. Proc. of the 1st Int. Symp. (Eds.R.Glowinsti et al.), SIAM, Philadelphia, 1988.

4. Алексеев, А.С. О лучевом методе вычисления интенсивности волновых факторов / А.С. Алексеев, Б.Я. Гельчинский // Вопросы динамической теории распространения сейсмических волн. JL, 1961. - N 5. - С.3-24.

5. Бабич, В.М. Асимптотические методы в задачах дифракции коротких волн. / В.М. Бабич, B.C. Булдырев. М. : Наука, 1982. - 248с.

6. Бабич, В.М. Пространственно-временной лучевой метод /В.М. Бабич, B.C. Булдырев, И.А. Молотков. JI. : Изд. Леи. Универ., 1985.

7. Клем-Mycamoe, К.Д. Теория краевых волн и ее применение в сейсмике / К.Д. Клем-Мусатов. Новосибирск : Наука, 1980.

8. Keller, J.B. Geometrical theory of diffraction /J.B. Keller //J. Opt. Soc. Am. 1981. - Vol.52. - N 4. - P. 175-188.

9. Киселев, А.П. Высшие приближения лучевого метода и "нелучевые явления" в неоднородных вязкоупругих средах /А.П. Киселев // Изв. АН СССР. Физика Земли. 1992. - N 11. - С.35-38.

10. Popov, М.М. A new method of computation of wave fields in the high-frequency approximation /М.М. Popov. Leningrad, 1981. - 20p. -(Preprint AN SSSR)

11. Алексеев, А.С. О лучевом методе вычисления полей волн в случае неоднородных сред с криволинейными границами раздела / А.С. Алексеев, Б.Я. Гельчинский // Вопросы динамической теории распространения сейсмических волн. JL, 1959. - N 3. - С. 16-47.

12. Cerveny, V. Ray method in seismology / V. Cerveny, I.A. Molotkov, I. Psencik. Prague : Varlovar. Univ.,1977.

13. Keller, J.B. Surface waves on waves on water of non-uniform depth /

14. J.B. Keller //J. Appl. Phys. 1980. - Vol.31. - N 6. - PP.1039-1046.

15. Tygel, M. Transient waves in layered media / M. Tygel, P. Hubral //Elsevier, Amsterdam. 1987.

16. Hron, F. Numerical modeling of nongeometrical effects by the Alekseev-Mikhailenko method / F. Hron, B.G. Mikhailenko //Bulletion of the seismological of America. 1981. - Vol.71. - N 4. - P.1011-1029.

17. Алексеев, А.С. "Нелучевые" эффекты в теории распространения сейсмических волн / А.С. Алексеев, Б.Г. Михайленко //Докл. АН СССР. 1982. - Т.267. - Т 5. - С.1079-1084.

18. Numerical methods used un atmospherical models // GARP Publication Series. 1979. - Vol.11. - N17.

19. Fornberg, В. On a Fouriese mehod for the integration of hyperbolic equation / B. Fornberg // Soc. Industr. Appl. Math., J.Numer. Anal. 1975. - N 12. - P.509-528.

20. Kreiss, H. Comparison of accurate methods for the integration of hyperbolic equation / H. Kreiss, J. Oliger // Tellus. 1972. - N 24. - P.199-215.

21. Orszag, S.A. Comparison of pseudospectral and spectral approximation /S.A. Orszag // Stud. Appl. Math. N 51. - P.253-259.

22. Kosloff, D. Forward modeling by a Fourier method / D. Kosloff, E. Baysal // Geophysics. 1982. - N 47. - P.1402-1412.

23. Kosloff, D. Elastic wave calculation by the Fourier method /D.Kosloff, M. Reshef, D. Loewenthal // Bull. Seis. Am. 1984. - N 74. - P.875-891.

24. Михайленко, Б. Г. Численное моделирование сейсмических полей в двумерно-неоднородных средах / Б.Г. Михайленко, В.И. Корнев // Моделирование волновых полей. Новосибирск, 1983 - С.41-70.

25. Mikhailenko, B.G. Synthetic seismograms for complex threedimensional geometries using an analytical-numerical algorithm /B.G. Mikhailenko // Geopys. J. R. Astr. Soc. 1984. - N 79,3. - P. 963-986.

26. Mikhailenko, B.G. Numerical experiments in seismic investigations /B.G. Mikhailenko // J. Geophys. 1985. - N 58. - P.101-124.

27. Алексеев, А. С. О задаче Лэмба для неоднородного полупространства / А.С. Алексеев, Б.Г. Михайленко// Докл. АН СССР. 1974. - Т.214. -N 11. - С.84-86.

28. Алексеев, А. С. Решение задачи Лэмба для вертикально-неоднородного полупространства /А.С. Алексеев, Б.Г. Михайленко // Изв. АН СССР. Сер. Физика Земли. 1976. - N12. - С.11-25.

29. Михайленко, Б.Г. Численное решение задачи Лэмба для неоднородного полупространства / Б.Г. Михайленко // Математические проблемы геофизики. Новосибирск, 1973. - С.273-297.

30. Алексеев, А. С Метод вычисления теоретических сейсмограмм для сложно построенных моделей сред/А.С. Алексеев, Б.Г. Михайленко // ДАН СССР. 1978. - т.240. - N 5. - С.1062-1065.

31. Михайленко, Б.Г. Сейсмические поля в сложнопостроенных средах (Атлас численных снимков и теоретических сейсмограмм) / Б.Г. Михайленко; под ред. А.С. Алексеева. Новосибирск, 1988. - 312с.

32. Петпрашенъ, Г.И. Распространение упругих волн в слоисто-изотропных средах, разделенными параллельными плоскостями / Г.И. Петрашень //Уч. Зал. ЛГУ, 1952. N 162. - С.3-189.

33. Петрашень, Г.И. Общая количественная теория отраженных и головных волн, возбуждающихся в слоистых средах с плоско-параллельными границами раздела /Г.И. Петрашень //Вопросы динамической теории распространения сейсмических волн. JL, 1957. - С. 70-164.

34. Огурцов, К. И. Динамические задачи для упругого полупространства в случае осевой симметрии /К.И. Огурцов, Г.И. Петрашень. Учен. зап. ЛГУ, 1951. - N 149. - С.3-117.

35. Петрашень, Г.И. О распространении волн в сложно-изотропных средах / Г.И. Петрашень, И.Н. Успенский. Учен. зап. ЛГУ, 1956. - N 208.1. C.58-141.

36. Петрашень, Г.И. Элементы динамической теории распространения сейсмических волн / Г.И. Петрашень //Вопросы динамической теории распространения сейсмических волн. Л., 1959. - Сб.З. - С.11-106.

37. Thomson, W. Т. Transmission of elastic waves through a stratified solid medium /W.T. Thomson // J. Appl. Phys. 1959. - N 21. - P.89-93.

38. Haskell, N.A. The dispersion of surface waves a multilayered media / N.A. Haskell // Bull. Seismol. Soc. Amer. 1953. - Vol.43. - N 1. - P. 17-34.

39. Harkrider D.G. Surface waves in multilayered elastic media. 1. Rayleigh and Lowe waves from buried sources in a multilayered elastic half-space. /

40. D.G. Harkrider // Bul.Seimol.Soc.Amer. 1964. Vol.54. - N 2. - P.627-673.

41. Молотков, Л. А. Матричный метод в теории распространения волн в слоистых упругих и жидких средах / J1.A. Молотков. JI. : Наука, 1984. - 366с.

42. Alford, R.M. Accuracy of finite-difference modeling of the acoustic wave equation / R.M. Alford, K.R. Kelly, D.M. Boore // Geophysics. 1974. -Vol. 39. - P.834-842.

43. Alterman, Z. Propagation of elastic waves layered media by finite-difference method / Z. Alterman, F.C. Karal // Bull. Seism.Soc. Am. 1968. - Vol. 58. - P.367-398.

44. Boore, D.M. Finite-difference methods for seismic wave propagation in heterogeneous materials / D.M. Boore // Methods in Computational Physics. 1972. - N 11. - P.l-37.

45. Smith, W.D. The application of finite element analysis to body wave propagation problems / W.D. Smith // Geophys. J.R.Astr.Soc. 1975. - N 42. - P.747-768.

46. Zahradnik, J. Finite-difference solution to certain diffraction problem / J. Zahradnik // Stud, and Gead. 1975. - Vol.19. - P.233-244.

47. Кошур, В.Д. Континуальные и дискретные модели динамического деформирования элементов конструкций /В.Д. Кошур, Ю.В. Немиров-ский. Новосибирск: Наука. Сиб. отд.-ние, 1990. - 198с.

48. Магомедов, К.М. Сеточно-характеристические численные методы /К.М. Магомедов, А.С. Холодов. М.: Наука, 1988.

49. Магомедов, К.М. О построении разностных схем для уравнений гиперболического типа на основе характеристических соотношений /К.М. Магомедов, А.С. Холодов // Журнал вычислительной математики и математической физики. 1969. - Т.9. - N 2. - С. 373-386.

50. Петров, И. Б. Численное решение некоторых динамических задач деформируемого твердого тела сеточно-характеристическим методом / И.Б. Петров, А.С. Холодов // Журнал вычислительной математики и математической физики. 1984. - N.5. - С. 722-739.

51. Рахматулин, Х.А. Применение метода пространственных характеристик к решению задач по распространению упругопластических волн / Х.А. Рахматулин, Т.Д. Каримбаев, Т. Байтелиев // Изв. Кав. СССР. Сер. физ.-мат. 1973. - N. 1. - С. 141-152.

52. Сабодаш, П. О. Применение методы пространственных характеристик к решению осесимметричных задач по распространению упругих волн / П.О. Сабодаш, Р.А. Чередниченко // ПМТФ. 1971. N. 4. - С. 101-109.

53. Годунов, С.К. Численное решение многомерных задач газовой динамики /С.К. Годунов, А.В. Забродин, А.В. Иванов и др. М.: Наука, 1976. -400с.

54. Lax, P.D. Systems of conservation laws / P.D. Lax, Wendroff // Communs Pure and Apppl. Math. 1960. - V.13. - P.217-237.

55. Lax, P.D. Weak solutions of nonlinear hyperbolic equations and their numerical computation /P.D. Lax // Communs Pure and Apppl. Math.- 1954. V.7. - P.159-193.

56. Русанов, В. В. Разностная схема третьего порядка точности для сквозного счета разрывных решений / В.В. Русанов // ДАН СССР. 1968.- N 6. С.1303-1305.

57. Courant, R. On the solution of nonlinear hyperbolic differential equations by finite differences / R. Courant, E. Isaacon, M. Rees // Communs Pure and Apppl. Math. 1952. - vol.5.

58. Roe, P.L. Approximate Riemann solvers, parameters vectors and difference schemes / P.L. Roe // Journal of Computational Physics. 1981. - vol. 43.

59. Roe, P.L. Characteristic based schemes for Euler equations / P.L. Roe // And. Rev. Fluid Mechanics. 1986. - Vol. 18.

60. Roe, P.L. Efficient construction and utilization of approximate Riemann solutions / P.L. Roe, J. Pike // Computing Method in Applied Sciences and Engineering. Amsterdam, North - Holland, 1984. - Vol.6

61. Osher, S. Numerical solution of singular perturbation problems and hyperbolic systems of conservation laws / S. Osher // North Holland Mathematical Studies. 1981. - Vol. 47.

62. Steger, J.L. Flux vector splitting of the inviscrid gas dynamic equations with application to finite difference methods / J.L. Steger, R.F. Warming // Journal of Computational Physics. 1981. - Vol. 40. - No 2.

63. Harten, A. High resolution schemes for hyperbolic conservation laws / A. Harten // Journal of Computational Physics. 1983. - Vol.49. - P.367-393.

64. Swety, P.K. High resolution schemes using flux limiters for hyperbolic conservation laws / P.K. Swety // SIAM J. Numer. Analys. 1984. - Vol.21. - P.995-1011.

65. Yee, H.C. Application of TVD-schemes for the Euler equations of gas dynamics / H.C. Yee, R.F. Warming, A. Harten // Lecture of Applied Mathematics. 1985. - Vol. 22.

66. Harten, A. On a class of high resolution total-variation-stable finite-difference schemes / A. Harten // NYU Report New York, NYU, 1982.

67. Harten, A. The artificial compression method for computation of shocks and contact discontinuities: III. Self-adjusting hybrid schemes / A. Harten // Math. Comput. 1978. - Vol. 32. - P.363-389.

68. Chakravarthy, S.R. Computing with high-resolution upwind schemes for hyperbolic equation / S.R. Chakravarthy, S. Osher // Lectures in Applied Mathematics. 1985. - vol. 22.

69. Годунов, С.К. Разностный метод численного расчета разрывных решений гидродинамики / С.К. Годунов // Математический сборник. 1959. - Т.47. - N 3. - С.271-306.

70. Van Leer В. Toward the ultimate conservation difference scheme. V. A second-order segcul to Godunov's method / Van Leer B. // Journal of Computational Physics. 1979. - V.32. - P. 101-136.

71. Harten, A. Uniformly high-order accurate nonoscillatory schemes / A. Harten, S. Osher // SIAM J. Numer. Analys. 1987. - V.24. - P.279-309.

72. Richtmycr, R.D. A survey of difference method for nonsteady flyid dynamics / R.D. Richtmycr // NCAR Technical Note 63-2-Colorado, Boulder, 1963.

73. Mac Cormack R. W. The effect of viscosity in hypervelocity impact cratering / Mac Cormack R.W. // AIAA Paper N 69 354, 1969.

74. Колган, В. П. Применение принципа максимальных значений производной к построению конечно-разностпых схем для численного анализа разрывных течений гидродинамики / В.П. Колган // Ученые записки ЦАРИ. 1972. - Т 3.

75. Van Leer В. Toward the ultimate conservation difference scheme. I. The quest of monotonicity / Van Leer В. // Lecture Notes in Physics. 1973. -vol.18.

76. Van Leer B. Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second-order scheme / Van Leer B. // Journal of Computational Physics. 1974. - Vol.14.

77. Van Leer B. Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow /

78. Van Leer В. // Journal of Computational Physics. 1977. - Vol. 23. - P. 263-275.

79. Van Leer B. Towards the ultimate conservative difference scheme IV. A new approach to numerical convection / Van Leer B. // Journal of Computational Physics. 1977. - Vol. 23. - P. 276-299.

80. Harten, A. The method of artificial compression / A. Harten // CIMS Report COO 30077-50- New York, Courant Institute, NYU, 1974.

81. Harten, A. Self-adjusting hybrid schemes for shock computation / A. Harten, G. Zwas // Journal of Computation Physics. 1972. - Vol.6.

82. Beam, R. An implicit finite-difference algorithm for hyperbolic systems in conservation-law-form / R. Beam, R.F. Warming // Journal of Computational Physics. 1976. - Vol.22.

83. Boris, J.P. Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works / J.P. Boris, D.L. Book // Journal of Computational Physics. 1973. - Vol. 11.

84. Boris, J.P. Flux-corrected transport. II Generalization of the method / J.P. Boris, D.L. Book, K. Hain // Journal of Computational Physics. -1975. Vol. 18.

85. Boris, J.P. Flux-corrected transport. II Minimal-error FGT algorithms / J.P. Boris, D.L. Book // Journal of Computational Physics. 1976. - Vol. 20.

86. Коновалов, A.H. Численные методы в динамических задачах теории упругости /А.Н. Коновалов. Сиб. матем. журн. - Новосибирск, 1997. - 38:3. - С. 551 - 568.

87. Белоцерковский, О.М. Метод расщепления для исследования течений стратифицированной жидкости со свободной поверхностью / О.М. Белоцерковский, В.А. Гущин, В.Н. Копынин // Журнал вычислительной математики и математической физики. 1987. - Т. 27.

88. Коновалов, А.Н. Вариационная оптимизация итерационных методов расщепления / А. Н. Коновалов. Сиб. матем. эюурн. - Новосибирск, 1997. - 38:2. - С. 312-325.

89. Софронов, И.Л. Искусственные граничные условия адекватные волновому уравнению вне сферы / И.Л. Софронов. Препринт Института Прикладной Математики им. М.В.Келдыша Академии Наук СССР. -Москва, 1992. - No. 42.

90. Софронов, И.Л. Условия полной прозрачности на сфере для трехмерного волнового уравнения / И.Л. Софронов // Доклады Академии Наук.- 1992. Т. 326. -С. 953-957.

91. Keller, J.B. Exact nonreflecting boundary conditions / Joseph B. Keller, Dan Givoli // J.Comput. Phys. 1989. - Vol. 82. - P.172-192.

92. Givoli, D. nonreflecting boundary conditions based on Kirchhoff-type formulae / Dan Givoli, D.Cohen. // J. Comput. Phys. 1995. - Vol. 117.- P.102-113.

93. Grote,M.J. Exact nonreflection boundary conditions for the time-dependent wave equation / M.J. Grote, J.B. Keller. // SIAM J. Appl. Math. 1985.- Vol 55. P. 280-297.

94. Grote M.J. Nonreflecting boundary conditions for time-dependent scattering /M.J. Grote, J.B. Keller. // J. Comput. Phys. 1996. - Vol 127. - P.52-65.

95. Engquist, B. Far field boundary conditions for computation over long time / B. Engquist, L. Halpern. // Appl. Numer. Math. 1988. - Vol. 4• -P.21-45.

96. Engquist, B. Absorbing boundary conditions for the numerical simulation of waves / B. Engquist, A. Majda // Math. Comput. 1977. - Vol.31. -P. 629-651.

97. Higdon, R.L. Absorbing boundary conditions for difference approximations to the multidimensional wave equation / Robert L. Higdon // Math. Сотр. 1986. - Vol. 47(176). - P.437-459.

98. Higdon, R.L. Absorbing boundary conditions for accoustic and elastic waves in stratified media / Robert L. Higdon // J. Comput. Phys. 1992. - Vol. 101(2). - P.386-418.

99. Bayliss, A. Radiation boundary conditions for wave-like equations / Alvin Bayliss, Eli Turkel // Comm. Pure Appl. Math. 1980. - Vol. 33(6). -P.707-725.

100. Givoli, D. Nonreflection boundary conditions for elastic waves / Dan Givoli, Joseph B. Keller // Wave Motin. 1990. - Vol. 12(3), P.261-279.

101. Hedstrom, G.W. Nonreflecting voundary conditions for nonlinear hyperbolic systems / G.W. Hedstrom // J. Comput. Phys. 1979. - Vol. 30(2). -P.222-237.

102. Atkins, H. Nonreflective boundary conditions for high-order methods / H. Atkins, J. Casper // AIAA Journal. 1994. - Vol. 32. - P.512-518.

103. Collino, F. Application of PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media / F. Collino, C. Tsogka // Geophysics. 2001. - Vol.66. - P. 294-307.

104. Крюков, В.А. Разработка параллельных программ для вычислительных кластеров и сетей / В.А. Крюков //Информационные технологии и вычислительные системы. 2003. - N 12.

105. Марчук, Г.И. Методы растепления / Г.И. Марчук. М. : Наука, 1988.

106. Садовская, О. В. Метод сквозного счета для исследования упругопла-стических волн в сыпучей среде / О.В. Садовская // Ж. вычисл. математики и матем. физики. 2004• - N 10. - С. 1909-1920.

107. Куликовский, А.Г. Математические вопросы численного решения гиперболических систем уравнений / А.Г. Куликовский, Н.В. Погорелое, А.Ю. Семенов. М. : Физматлит, 2001.

108. Магометов, К.М. Сеточно-характеристические численные методы / К.М. Магометов, А.С. Холодов. М. : Наука, 1988.

109. Воеводин, В.В. Параллельные вычисления / В.В. Воеводин, Вл.В. Воеводин Спб. : ВХВ-Перетбург, 2002.

110. Ортега, Дою. Введение в параллельные и векторные методы решения линейных систем : пер. с англ.] / Дж. Ортега. М. : Мир, 1991.

111. Кучунова, Е.В. Параллельный алгоритм расчета упругих волн в блочной среде / Е.В. Кучунова // Материалы VI межрегиональной школы-семинара "Распределенные и кластерные вычисления". Красноярск, 2006г. - С. 79-89.

112. Кучунова, Е.В. Численное исследование распространения сейсмических волн в блочных средах на многопроцессорных вычислительных системах / Е.В. Кучунова, В.М. Садовский // Вычислительные методы и программирование. 2008. - Т.9. - N 1. - С.70-80.