Декомпозиция задачи конечно разностного...

8
Декомпозиция задачи конечно-разностного эластического моделирования для расчета синтетических сейсмограмм в гриде Тульчинский В.Г. 1 , Ющенко Р.А. 1 1 Институт кибернетики им. В.М. Глушкова НАНУ, пр. Глушкова 40, Киев, Украина [email protected], [email protected] Аннотация. Ключом к «гридификации» задач, т.е. к приспособлению алгоритмов для распределенного асинхронного выполнения разнородными компьютерами, связанными медленной сетью, является декомпозиция на возможно меньшие независимые подзадачи. Рассматривается задача расчета синтетических данных сейсморазведки. Эта задача естественным образом разделяется на независимые подзадачи моделирования тысяч волновых полей, возбужденных отдельными источниками. Однако для расчета одного трехмерного волнового поля с необходимой детальностью нужен сегмент кластера в десяток узлов, иначе массивы не поместятся в оперативную память. За счет ограничения геометрии среды без снижения параметрической сложности модели удалось уменьшить объем вычислений на два порядка, а за счет разложения в спектр свести трехмерную задачу к набору независимых двухмерных подзадач. Для расчета отдельной двухмерной подзадачи достаточно одного узла кластера или даже обычного однопроцессорного ПК, что обеспечивает эффективность распределенного вычисления этой задачи в гриде. Ключевые слова Распределенные вычисления, грид, моделирование, синтетические трехмерные сейсмограммы. 1 Введение Сейсморазведка является наиболее распространенным поверхностным методом изучения геологического строения Земли, поиска полезных ископаемых и планирования бурения. Сейсмический метод зондирования заключается в возбуждении низкочастотных колебаний земной поверхности с помощью вибраторов или взрыва и записи эхосигнала, отраженного от подземных границ пластов породы с разными свойствами, с помощью геофонов, фиксирующих свои колебания (рис. 1), или гидрофонов, фиксирующих давление воды. Рис. 1. Оборудование полевых работ сейсморазведки: А вибраторы, Б геофон. А Б Міжнародна конференція "Високопродуктивні обчислення" HPC-UA’2011 (Україна, Київ, 12-14 жовтня 2011 року) ____________________________________________________________________________________________________________ -137-

Transcript of Декомпозиция задачи конечно разностного...

Page 1: Декомпозиция задачи конечно разностного ...hpc-ua.org/hpc-ua-11/files/proceedings/1.27(137).pdfДекомпозиция задачи конечно-разностного

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

синтетических сейсмограмм в гриде

Тульчинский В.Г.1, Ющенко Р.А.1

1 Институт кибернетики им. В.М. Глушкова НАНУ, пр. Глушкова 40, Киев, Украина

[email protected], [email protected]

Аннотация. Ключом к «гридификации» задач, т.е. к приспособлению алгоритмов для распределенного асинхронного выполнения разнородными компьютерами, связанными медленной сетью, является декомпозиция на возможно меньшие независимые подзадачи. Рассматривается задача расчета синтетических данных сейсморазведки. Эта задача естественным образом разделяется на независимые подзадачи моделирования тысяч волновых полей, возбужденных отдельными источниками. Однако для расчета одного трехмерного волнового поля с необходимой детальностью нужен сегмент кластера в десяток узлов, иначе массивы не поместятся в оперативную память. За счет ограничения геометрии среды без снижения параметрической сложности модели удалось уменьшить объем вычислений на два порядка, а за счет разложения в спектр – свести трехмерную задачу к набору независимых двухмерных подзадач. Для расчета отдельной двухмерной подзадачи достаточно одного узла кластера или даже обычного однопроцессорного ПК, что обеспечивает эффективность распределенного вычисления этой задачи в гриде.

Ключевые слова

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

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

Рис. 1. Оборудование полевых работ сейсморазведки: А – вибраторы, Б – геофон.

А Б

Міжнародна конференція "Високопродуктивні обчислення" HPC-UA’2011 (Україна, Київ, 12-14 жовтня 2011 року)

____________________________________________________________________________________________________________

-137-

Page 2: Декомпозиция задачи конечно разностного ...hpc-ua.org/hpc-ua-11/files/proceedings/1.27(137).pdfДекомпозиция задачи конечно-разностного

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

В результате рынок сейсморазведки во всем мире растет. Аппаратура, методы наблюдения и обработки постоянно совершенствуются. В настоящее время типичным является выполнение трехмерной сейсморазведки, при которой колебания последовательно возбуждаются в узлах регулярной сетки источников, и эхосигнал записываются регулярной сеткой приемников-геофонов, расположенных вокруг источника. Шаги сеток источников и приемников обычно составляют 20-50 метров, размах сетки приемников – до 5 км от источника (рис. 2). Популярными становятся трехкомпонентные исследования с использованием геофонов, записывающих свои колебания в трех направлениях, и девятикомпонентные исследования, когда в каждом источнике последовательно возбуждаются колебания, направленные вдоль осей X, Y и Z, а отраженная волна записывается трехкомпонентными приемниками. Применение многокомпонентных наблюдений, более густой сетки возбуждения/приема и большего размаха расстановки приемников для отдельного источника повышает стоимость полевых работ. Одновременно растет объем полевых данных (сейсмограмм) и стоимость их обработки.

Полевые работы, как правило, составляют 95% от общих затрат на сейсморазведку, поэтому задача их

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

Рис. 2. Пример системы наблюдений трехмерной сейсморазведки: черные точки – источники, розовые – приемники.

Міжнародна конференція "Високопродуктивні обчислення" HPC-UA’2011 (Україна, Київ, 12-14 жовтня 2011 року)

____________________________________________________________________________________________________________

-138-

Page 3: Декомпозиция задачи конечно разностного ...hpc-ua.org/hpc-ua-11/files/proceedings/1.27(137).pdfДекомпозиция задачи конечно-разностного

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

2 Теория Для синтеза сейсмограмм по заданной модели геологической среды и расстановке применяют методы, основанные на упрощенных физических моделях. Для планирования трехмерных наблюдений широко внедрено лучевое трассирование. Этот метод основан на высокочастотном приближении и использует известные из оптики законы преломления и отражения продольных и поперечных волн на скоростной границе [1]. Метод лучевого трассирования достаточно эффективен для оценки освещенности целевой поверхности в зависимости от ее геометрии, скоростей над ней и расстановки [2]. Его преимуществом является возможность идентификации типа каждой волны, в частности, построения сейсмограмм без кратных и преобразованных волн, или с ограниченной кратностью. Однако синтезированные сейсмограммы достаточно далеки от реальных, т.к. не учитывают изменение формы и расщепление сигнала на неоднородностях, размер которых сопоставим с длиной волны. При типичных частотах сейсмического сигнала 20-50 Гц и скоростях продольных волн 2-6 км/с длина волны изменяется в десятках и даже сотнях метров. А мощность типичных нефтяных коллекторов – до 5 м. Помимо этого лучевое трассирование непригодно для моделирования отражений, возникающих в градиентных средах, и плохо приспособлено к условиям резких границ и точек дифракции (так как использует интерполяцию лучей и амплитуд). C развитием высокоточной широкополосной сейсмологической аппаратуры появились факты регистрации “нелучевых” сейсмических волн, которые не описываются нулевым членом лучевого ряда, и для их вычисления необходимо учесть последующие члены ряда [3].

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

( ) ( ) ( ) ( )ttt ,,, xfxΣxsx +⋅∇=′′ρ , (1)

где ( )t,x – точка в пространстве-времени, ρ – плотность породы, s – вектор смещений, Σ – тензор напряжений, f – внешняя сила (сигнал). Оператор « ⋅∇ » обозначает дивергенцию в пространстве, а штрих – производную по времени. Тензор напряжений Σ – это квадратная симметричная матрица mkkm σσ = , причем

kkσ – напряжение сжатия/растяжения вдоль оси k , а если mk ≠ , то kmσ – это напряжение сдвига в плоскости

km . Внешняя сила f в задачах моделирования сейсморазведки прилагается в точку источника, ее можно выразить через δ-функцию Кронекера: ( ) ( ) ( )tt F, 0xxxf −= δ . Конкретизация системы (1) зависит от характера физической модели. В эластическом изотропном приближении [4]:

( ) ( ) ( )( ) ( ) ( ) ( )( )( )T,,,, tttt xsxsxIxsxxΣ ∇+∇+⋅∇= µλ , (2)

где I – это единичная матрица: ( )mkIkm −= δ ; модель среды представляют коэффициенты Ламе

( ) ( ) ( )xxx 2SVρµ = и ( ) ( ) ( ) ( )xxxx µρλ 22 −= PV ; оператор ∇ обозначает градиент в пространстве, а TM –

транспонированную матрицу M : mkTkm MM = . Таким образом, в трехмерном изотропном эластическом

случае состояние поля описывается девятью массивами параметров: 3 компоненты s и 6 разных компонент Σ . Модель задается тремя массивами параметров (коэффициенты Ламе и ρ ), что в сумме дает 12 массивов.

Для самого простого акустического приближения (идеальной жидкости) ( ) 0=xµ , и тензор напряжений Σ можно сократить до единственной скалярной компоненты σ . В этом случае состояние поля описывается четырьмя массивами параметров: 3 компоненты s и σ , а для модели достаточно двух параметров (всего – 5 массивов). Акустическое приближение используется для ускорения вычислений. Однако в силу отсутствия поперечных, преобразованных волн, поверхностных волн Релея и т.п., полученные синтетические сейсмограммы имеют ограниченную сферу применения.

В более полном эластическом анизотропном приближении [5] тензор напряжений выражается через тензор жесткости Λ и может включать компоненты пары сил (диполей) M :

Міжнародна конференція "Високопродуктивні обчислення" HPC-UA’2011 (Україна, Київ, 12-14 жовтня 2011 року)

____________________________________________________________________________________________________________

-139-

Page 4: Декомпозиция задачи конечно разностного ...hpc-ua.org/hpc-ua-11/files/proceedings/1.27(137).pdfДекомпозиция задачи конечно-разностного

( ) ( ) ( ) ( )tMtxst mn

k q q

kmnkqmn ,,, xxxx +

∂∂

Λ=∑∑σ . (3)

Состояние поля описывается девятью массивами параметров: 3 компоненты s и 6 компонент Σ . В общем анизотропном случае модель включает до 22 коэффициентов ( Λ и ρ ) в каждой точке, что в сумме дает 31 массив. Анизотропия возникает в реальных данных в местах переслаивания тонких пластов, в зонах трещиноватости, аномально высокого порового давления и т.п. Эти явления характерны для большинства месторождений нефти и газа. Недоучет анизотропных свойств среды во время полевых работ, обработки и интерпретации может привести к значительным ошибкам (рис. 3).

Еще более точным является вискоэластическое приближение [6], поскольку для высоких частот и слабых

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

Исторически первый и самый распространенный метод численного решения системы (3) – явная схема метода конечных разностей [8]. Поскольку сейсмические наблюдения проводятся с большими удалениями источник-приемник и на большую глубину зондирования, размеры сетки измеряются в сотнях длин волн. Для сохранения формы сигнала необходимо иметь примерно 10 точек на длину волны. Поэтому сетка конечных разностей при моделировании сейсмических наблюдений в 3D состоит из миллиардов узлов. Учитывая, что в эластическом анизотропном случае необходимо хранить 31 такой массив, и что число с одинарной точностью занимает 4 байта, получаем, что на каждой итерации нужно обрабатывать сотни гигабайт. Это больше, чем объем оперативной памяти большинства современных кластерных узлов и многопроцессорных серверов.

Для сокращения необходимого объема памяти и ускорения вычислений применяется широкий спектр методов: конечно-разностные с аппроксимацией производных схемами высокого порядка, неравномерные сетки, метод конечных элементов, спектральные и псевдо-спектральные подходы. Хороший обзор применяемых вычислительных методов дан в [3, 9]. Однако, в силу того, что размер неоднородностей среды, как правило, гораздо меньше длины волны, а временной шаг решения практически ограничен шагом дискретизации выходных сейсмограмм (обычно 2 мс), эффект от применения сложных методов к детальным моделям редко перекрывает дополнительные вычислительные затраты. Хороший обзор современных конечно-разностных схем трехмерного эластического анизотропного моделирования сделан Лисицей и Вишневским в [5]. Они предложили дискретное решение, оптимизированное с точки зрения объема памяти, и оценили условия его стабильности.

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

Рис. 3. Пример обработки синтетических сейсмограмм с учетом анизотропии (слева) и без ее учета (справа): Скалистые горы, Канада.

Міжнародна конференція "Високопродуктивні обчислення" HPC-UA’2011 (Україна, Київ, 12-14 жовтня 2011 року)

____________________________________________________________________________________________________________

-140-

Page 5: Декомпозиция задачи конечно разностного ...hpc-ua.org/hpc-ua-11/files/proceedings/1.27(137).pdfДекомпозиция задачи конечно-разностного

общего назначения (GPU), появляется возможность моделировать отдельные волновые поля за приемлемое время. Например, в [8] сообщено о достижении 25-50 кратного ускорения анизотропного эластического моделирования на GPU-кластере с картами NVidia. Однако вычислительная стоимость моделирования полного набора синтетических данных трехмерной сейсморазведки (включающего тысячи полей) все еще существенно превосходит практические ограничения.

Трехмерное моделирование может быть упрощено, если зафиксировать свойства среды вдоль одного направления ( 2x ). Модели такого типа (рис. 4) называются «двух с половиной мерными» (2.5D). Ограничение 2.5D не дает возможность синтезировать сейсмограммы, соответствующие реальным объектам. Следовательно, 2.5D не позволяет оптимизировать расстановку с точки зрения освещенности целевых горизонтов (для чего используется лучевое трассирование). Однако 2.5D позволяет исследовать влияние многокомпонентных наблюдений и плотности расстановки на интерпретацию результатов сейсморазведки сложных сред с учетом тонкослоистости, трещиноватости, и анизотропии, что недоступно лучевому трассированию.

Поскольку в 2.5D-среде сейсмограмма не меняется, если источники и приемники одинаковым образом

переместить вдоль оси 2Y x= , в такой модели достаточно моделировать только одну линию источников, а остальные получать копированием данных с заменой координат трасс. Благодаря этому в случае 2.5D общий объем вычислений можно снизить в десятки, иногда – в сотни раз.

3 Метод Сонг и Вильямсон [11] свели задачу акустического моделирования в 2.5D-среде с постоянной плотностью к системе линейных уравнений путем преобразования Фурье во времени и вдоль оси Y. Они решали эту систему методом LU-декомпозиции. Чао и Гринхалг [12] вывели критерий устойчивости и предложили метод реализации поглощающих границ. Нето и Коста [13] предложили метод моделирования 2.5D для эластического изотропной и трансверсально-изотропной среды. Реализации 2.5D для произвольной 3D TTI анизотропии описаны в [14, 15]. Применение моделирования 2.5D описано в [16].

Для линеаризации систему (1)-(3), как правило, переписывают через скорости смещения su ′= . Запишем ее в скалярном виде (по повторяющимся индексам выполняется суммирование):

( ) ( ) ( ) ( )( ) ( ) ( ) ( )

∂∂

+∂

∂Λ=

∂∂

+∂

∂=

∂∂

ttxxxM

xtxxxuxxx

ttxxx

txxxfx

txxxt

txxxuxxx

mn

q

kmnkq

mn

nk

nkn

,,,,,,,,,,,

,,,,,,,,,,,

321321321

321

321321321

321

σ

σρ. (4)

В 2.5D-среде ( ) ( )31321 ,0,,, xxxxx mnkqmnkq Λ=Λ и ( ) ( )31321 ,0,,, xxxxx ρρ = , поэтому систему (4) удобно

записать в спектральном виде, сделав Фурье-преобразование вдоль оси 2x :

( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( ) ( )

∂∂

+Λ+∂

∂Λ=

∂∂

++∂

∂=

∂∂

ttxxMtxxujxx

xtxxuxx

ttxx

txxftxxjx

txxt

txxuxx

mnmnk

q

kmnkq

mn

nnq

nqn

,,,~,,,~,0,,,,~

,0,,,,~

)5(,,,~,,,~,,,~,,,~,0,

31312312

3131

31

313123131

31

ωωωωωσ

ωωσωωσωρ

Рис. 4. Модель 2.5D.

Міжнародна конференція "Високопродуктивні обчислення" HPC-UA’2011 (Україна, Київ, 12-14 жовтня 2011 року)

____________________________________________________________________________________________________________

-141-

Page 6: Декомпозиция задачи конечно разностного ...hpc-ua.org/hpc-ua-11/files/proceedings/1.27(137).pdfДекомпозиция задачи конечно-разностного

где ( ) ( ) dyetxxutxxxu yjnn

ωω∫+∞

∞−

= ,,,~,,, 31321 , ( ) ( ) dyetxxtxxx yjmnmn

ωωσσ ∫+∞

∞−

= ,,,~,,, 31321 , 1−=j , m ,

n и k принимают значения {1, 2, 3}, a q – значения {1, 3}.

Численное решение системы (5) с помощью явной центральной конечно-разностной схемы второго порядка на трех сдвинутых сетках описано в [15]. Две 3D конечно-разностной схемы для (4) проанализированы в [5], а сравнительная оценка сложности вычислительных схем 2.5D и 3D в 2.5D-среде дана в [17].

В целом уравнения (5) аналогичны (4). Шаги сетки в пространстве и во времени ограничиваются теми же физическим предпосылками устойчивости, и при аналогичной схеме решения должны совпадать. Однако, nu~ и

mnσ~ – комплексные числа, в то время, как nu и mnσ – действительные. В результате реализация одного шага

решения (5) требует примерно вдвое больше вычислений. С другой стороны, в силу сопряженности ( )ω−nu~ и

( )ωnu~ , а так же ( )ωσ −mn~ и ( )ωσmn

~ , нет необходимости отдельно моделировать отрицательные частоты. Поэтому моделирование по уравнению (5) пары пространственных частот ω и ω− требует примерно того же числа операций и занимает примерно столько же времени, сколько расходуется при решении (4) на отдельную плоскость сетки: constx =2 . Следовательно, ускорение моделирования 2.5D методом (5) по сравнению с расчетом той же модели методом (4) можно оценить как отношение размера сетки в направлении 2x к количеству неотрицательных пространственных частот ω в частотном представлении волнового поля. Чтобы сохранить форму сигнала, пространственный шаг дискретизации (4): maxmin3 WfVy D =∆ , где 15..5=W – это минимальное число точек сетки на длину волны. При условии симметричной расстановки, максимального удаления приемников в направлении 2x равного maxy и без учета ширины поглощающего края, число

плоскостей сетки в трехмерном расчете minmaxmax2 WVfyN y ≈ . Для спектрального представления шаг

дискретизации определяется частотой Найквиста: ⇒=∆ maxmin5.2 2 fVy D minmaxmax 2 Vfy ππω ≈∆= .

Расстояние вдоль оси 2xy = между настоящим и ближайшим к нему мнимым источником ωπ ∆= 2iy .

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

нужно выполнять условие ( )maxmaxmaxmaxmaxmax 2 VtyVtyyi +<∆⇒>− πω . Отсюда количество частот:

max3

maxmaxmax

min

max

min

maxmaxmaxmaxmax

min

maxmax

21

221 t

yWV

WN

tfVV

VfyVty

VfN

D

yk ∆

+≈+=++

⋅=+∆

πω

ω. (6)

+≈=minmax

maxmax

21

yytVW

NN

Ускорениеk

y . (7)

0

1

2

3

4

5

6

7

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6tmax

Spe

edU

p

Уск

орен

ие

2.5D сложнее 3D!

Рис. 5. Пример соотношения объема вычислений для схем 3D и 2.5D.

Міжнародна конференція "Високопродуктивні обчислення" HPC-UA’2011 (Україна, Київ, 12-14 жовтня 2011 року)

____________________________________________________________________________________________________________

-142-

Page 7: Декомпозиция задачи конечно разностного ...hpc-ua.org/hpc-ua-11/files/proceedings/1.27(137).pdfДекомпозиция задачи конечно-разностного

Ускорение (7) не слишком велико, и при моделировании длительного волнового процесса расчет в спектральной области может потребовать в целом большого объема вычислений, чем 3D (рис. 5). Однако декомпозиция на независимые пространственные частоты полностью решает проблему объемов памяти. 2.5D-схема использует 160 байт на ячейку двумерной сетки для хранения как модели, так и состояния волнового поля. Процесс моделирования для отдельных пространственных частот выполняется полностью независимо, без обменов данными. Это позволяет разделить моделирование на сотни и даже тысячи небольших подзадач, которые можно выполнять одновременно или последовательно в распределенной вычислительной среде, в частности, в гриде. Полученные в результате моделирования 2D-файлы объединяются на хосте в синтетические 3D-сейсмограммы (с помощью обратного преобразования Фурье) и размножаются в соответствии с заданным числом линий источников.

4 Эксперименты Вычислительный эксперимент проводился в рамках виртуальной организации (ВО) GEOPARD [18] Украинского Национального Грида (УНГ). УНГ объединяет 24 грид-сайта, расположенные в 7 регионах Украины [19] и связанные оптоволоконной сетью пропускной способностью 10 Гбит/с. Суммарный объем ресурсов УНГ – 2656 процессорных ядер [20]. Инфраструктура УНГ основана на грид-платформе Advance Resource Connector (ARC) [21].

Тесты подтвердили граничные оценки параметров спектрального разложения (рис.6) и показали принципиальную возможность распределенного решения задачи трехмерного полноволнового сейсмического моделирования в постановке 2.5D в гриде.

B A C

y x

y x

y x

y

t

y

t

y

t

Vp=4000

, 2 . 0 = ε , 1 . 0 − = δ , 03 . 0 = γ , 45 = α

, 30 = ϕ 2500 = P V

2.5D TTI анизотропное эластическое моделирование одного волнового поля от сигнала Риккера 30 Гц в кубе

со стороной 2 км в течение 2 с, потребовало расчета 333 неотрицательных частот. Моделирование было выполнено за 6 ч. с участием грид-узлов ИК НАНУ, НТТУ КПИ, ИКИ НАНУ и КНАУ, ИМБГ НАНУ. Задача аналогичная примеру [5] (куб со стороной 2.5 км и шагом 2.5 м) в оригинале требовала 216 ГБ. Ее 2.5D версия разложена на 1000 подзадач, каждая из которых использовала 432 МБ. Расчет завершен за 13 ч.

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

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

Рис 6. Модель и синтетические сейсмограммы 2.5D (A: недостаточный диапазон частот maxω , B: слишком большой шаг ω∆ , C: правильно

подобранные параметры).

Міжнародна конференція "Високопродуктивні обчислення" HPC-UA’2011 (Україна, Київ, 12-14 жовтня 2011 року)

____________________________________________________________________________________________________________

-143-

Page 8: Декомпозиция задачи конечно разностного ...hpc-ua.org/hpc-ua-11/files/proceedings/1.27(137).pdfДекомпозиция задачи конечно-разностного

Накладные расходы на передачу данных через медленные сети грида в значительной степени компенсируются экономией на синхронизации и обмене данными между устройствами, занятыми моделированием одного поля в 3D-схеме. Более того, схема 2.5D моделирования оказалась эффективной и для вычислений на GPU. Поскольку все данные помещаются в память, отпадает необходимость в синхронизации данных устройства и общей памяти компьютера на каждом шаге. В результате, мы получили примерно вдвое более высокое ускорение СUDA-программы, чем [10].

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

6 Благодарности Авторы выражают благодарность Юрию Вячеславовичу Роганову за разностороннюю помощь и исчерпывающие консультации. Работы по трехмерному сейсмическому моделированию в гриде были поддержаны «Государственной целевой научно-технической программой внедрения грид-технологий на 2009-2013 годы». Мы благодарим руководство Программы и грид-сайты, предоставившие ресурсы ВО GEOPARD.

Литература

[1] Červený V. Seismic Ray Theory. – London: Cambridge University Press, 2001. – 722 p. [2] Gjøystdal H., Iversen E., Laurain L., Lecomte I., Vinje V., Åstebøl, K. Review of ray theory applications in

modelling and imaging of seismic data // Studia geophysica et geodaetica. – 2002. – No 46. – P. 113–164. [3] Активная сейсмология с мощными вибрационными источниками / Отв. ред. Г.М. Цибульчик. –

Новосибирск: ИВМиМГ СО РАН, Филиал “Гео” Издательства СО РАН, 2004. – 387 с. [4] Nilsson S., Petersson N. A., Sjögreen B., Kreiss H.-O. Stable difference approximations for the elastic wave

equation in second order formulation // SIAM J. Numer. Anal. – 2007. – No 45. – P. 1902–1936. [5] Lisitsa, V., Vishnevskiy, D. Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic

elasticity // Geophysical Prospecting. – 2010. – No 58. – P. 619–635. [6] Carcione J. M., Kosloff D., Kosloff R. Wave propagation simulation in a linear viscoelastic medium

//Geophysical Journal. – 1988. – No 95. – P. 597-611. [7] https://www.geophysik.uni-muenchen.de/~kaeser/SS09/Lecture1.pdf [8] Alterman Z., Karal F. C. Propagation of elastic waves in layered media by finite-difference methods // Bull.

Seism. Soc. Am. – 1968. – No 58. – P. 367–398. [9] https://e-reports-ext.llnl.gov/pdf/369608.pdf [10] Komatitsch D., Michéa D., Erlebacher G., Göddeke D. Running 3D Finite-difference or Spectral-element Wave

Propagation Codes 25x to 50x Faster Using a GPU Cluster // Extended Abstracts of 72th EAGE Conference & Exhibition. – 2010.

[11] Song Z., Williamson F.R. Frequency-domain acoustic-wave modeling and inversion of crosshole data: Part 1 – 2.5-D modeling method // Geophysics. – 1995. – Vol. 60, No 3. – P. 784–795.

[12] Cao S., Greenhalgh S. 2.5D modeling of seismic wave propagation: Boundary condition, stability criterion, and efficiency // Geophysics. – 1998. – Vol. 63, No 6. – P. 2082–2090.

[13] Neto F., Costa J. 2.5D anisotropic elastic finite-difference modelling // 76th SEG Annual Meeting, Expanded Abstracts. – 2006. – P. 2275–2279.

[14] Silva Neto, F., Costa, J., Novais, A. 2.5D Elastic Anisotropic Finite-Difference Modeling // 10th International Congress of the Brazilian Geophysical Society, Extended Abstracts. – 2007.

[15] Kostyukevych A., Marmalevskyi N., Roganov Y., Tulchinsky V. Anisotropic 2.5D - 3C finite-difference modelling // 70th EAGE Conference & Exhibition, Extended Abstracts. – 2008. – #P043.

[16] Kostyukevych A., Roganov Y. 2.5D forward modeling: a cost effective solution that runs on small computing systems // ASEG Extended Abstracts. – 2010. – #1–4.

[17] Lavreniuk S., Roganov Y., Tulchinsky V., Kolomiets O. Synergy of 2.5D approach and grid technology for synthesis of realistic 3D/3C seismograms in anisotropic media // 73rd EAGE Conference & Exhibition, Extended Abstracts. – 2011. – #P289.

[18] http://grid.org.ua/voms/ [19] http://grid.kpi.ua/ [20] http://lcg.bitp.kiev.ua/ [21] http://www.nordugrid.org/

Міжнародна конференція "Високопродуктивні обчислення" HPC-UA’2011 (Україна, Київ, 12-14 жовтня 2011 року)

____________________________________________________________________________________________________________

-144-