CC BY
26
УДК 519.6
DOI 10.25205/1818-7900-2020-18-3-81-95
Применение криволинейных сеток для моделирования 2D волновых полей в неоднородных областях со сложной топографией
П. А. Титов
Институт вычислительной математики и математической геофизики СО РАН Новосибирск, Россия
Аннотация
Численное моделирование обширно используется при изучении процессов распространения волновых полей в различных сложно построенных средах. Одним из способов является дискретизация области, для которой проводится моделирование, и построение разностной схемы для дальнейшей численной реализации. В данной работе задействуется построение сетки из криволинейных четырехугольников, хорошо согласованной с топографией поверхности. Предложен параллельный алгоритм численного моделирования на основе численного решения линейной 2D-системы теории упругости, записанной в смещениях, с использованием криволинейной сетки и явной разностной схемы. Представлены результаты моделирования. Расчеты проводились при помощи ресурсов ССКЦ СО РАН Ключевые слова
теория упругости, упругие волны, криволинейная поверхность, криволинейная сетка, суперкомпьютер, моделирование Благодарности
Работа выполнена в рамках государственного задания ИВМиМГ СО РАН (проект 0315-2019-0009) (пункты 1, 3.1, 3.2, 8), а также при поддержке грантов РФФИ № 19-07-00085 (пункт 5) и № 20-01-00231 (пункты 2, 4, 6, 7).
Автор выражает благодарность ССКЦ СО РАН за предоставленные вычислительные ресурсы для проведения численных экспериментов. Для цитирования
Титов П. А. Применение криволинейных сеток для моделирования 2D волновых полей в неоднородных областях со сложной топографией // Вестник НГУ. Серия: Информационные технологии. 2020. Т. 18, № 3. С. 8195. DOI 10.25205/1818-7900-2020-18-3-81-95
Utilizing Curvilinear Grid for 2D Wave Field Simulation Heterogeneous Media with Topography
P. A. Titov
Institute of Computational Mathematics and Mathematical Geophysics SB RAS Novosibirsk, Russian Federation
Annotation
The main practical interest of geophysics is the restoration of the structure of the medium. So one can understand whether there are hydrocarbon reservoirs at depth. Studying the structure of the medium by taking numerous samples and drilling wells is a very financially and time-consuming task. Therefore, the use of mathematical modeling seems to be a good alternative. If the calculations show a result similar to the initial data of the medium, then such algorithm for restoring the medium can be applied new regions. Verification of the algorithm should be carried out on real struc-
© П. А. Титов, 2020
tures of the media, but real wave field data are not provided, since they are often the subject of trade secrets. In this case, the result of this work comes to the rescue. Numerical simulation is widely used in studying the processes of wave field propagation in various complex media. One way is to discrete the domain for which the simulation is being performed, and to construct a difference scheme for further numerical implementation. This work involves the construction of a mesh with curved quadrangles, which allows good agreement with the topography of the surface. A parallel algorithm of numerical modeling based on the numerical solution of a linear 2D system of elasticity theory written in displacements using a curvilinear mesh and an explicit difference scheme is proposed. The simulation results are presented. The calculations were carried out using the resources of the SSCC SB RAS. Keywords
elasticity, elastic waves, curvilinear surface, curvilinear mesh, supercomputer, simulation Acknowledgements
This work was supported by the budget project 0315-2019-0009 of ICMMG SB RAS (sections 1, 3.1, 3.2, 8), RFBR grants 19-07-00085 (section 5) and 20-01-00231 (2, 4, 6, 7). For citation
Titov P. A. Utilizing Curvilinear Grid for 2D Wave Field Simulation Heterogeneous Media with Topography. Vestnik NSU. Series: Information Technologies, 2020, vol. 18, no. 3, p. 81-95. (in Russ.) DOI 10.25205/1818-7900-2020-183-81-95
1. Введение
Вкратце опишем суть процесса распространения волн: изначально считаем, что среда находится в состоянии покоя. В среду помещается источник возмущений (на поверхности или заглубленный), который по заранее определенному сценарию создает возмущения. Частицы среды вокруг источника отклоняются от положения равновесия, передавая это возмущение по соседним частицам. В результате создается распространяющееся волновое поле. Свойства этого волнового поля определяются параметрами среды (р - плотность, А, [ - коэффициенты Ламе) и ее строением. Вектор отклонения (он же вектор смещения) каждой частицы можно разложить на две компоненты в двумерном случае. Если восстановить значения вектора смещения в каждой точке области в каждый момент времени, то можно сказать, что мы восстановили процесс распространения волн в рамках математической модели. В нашем случае, восстановление волнового поля происходит по дискретному набору точек в дискретные моменты времени. Для этого задействуются соответствующие математические и численные методы, а также вычислительные мощности суперкомпьютера. Решение задачи основывается на численной реализации явной разностной схемы, созданной с учетом метода баланса. Моделирование будет проводиться на реальных данных, полученных с одного из участков Сибирского региона. Отличительной особенностью используемого подхода является задействование технологии построения криволинейных сеток. Криволинейные сетки активно используются в различных областях науки. Со способами построения сеток, а также их применением для решения некоторых актуальных задач можно ознакомиться в [1-3]. Также есть работы и в области моделирования волновых полей [4; 5], включая работы автора [6-9].
В дальнейшем текущая работа будет применена при создании 3D параллельного алгоритма моделирования волновых полей в неоднородных средах, характерных для Сибирского региона (например, Баженовская свита).
2. Построение криволинейной сетки
Существует много методов построения криволинейной сетки. Часть из них описана в работах [1-3]. Опишем суть технологии: криволинейная сетка в «физической» области (области с топографией, в которой ставится задача) есть результат взаимно-однозначного отображения равномерной прямоугольной сетки «расчетной» области (области, в которой будут производиться численные расчеты). Необходимо установить взаимно-однозначное соответствие между «физической» и «расчетной» областями. «Физическая» область находится в пространстве переменных , а «расчетная» - в пространстве переменных . Для
этого был использован метод трансфинитной интерполяции. Дополнением автора к этому методу является подбор управляющих метрик - функций, позволяющих влиять на параметры сетки.
В данной работе важным отличительным свойством сетки является ее ортогональность к свободной поверхности, что накладывает ограничение на форму свободных поверхностей, с которыми можно работать при помощи данного метода построения сетки, а именно: если поверхность рассматривать как функцию , она должна быть взаимно однозначной
гладкой функцией первого порядка на промежутке, на котором она построена для данной области (рис. 1). В дальнейшем будет показано использование такой локально-ортогональной сетки для реального участка Сибирского региона.
Я
А' В'
q=q(x,y) r=r(x,y)
x=x(q,r)
y=y(q,r)
D' с
Расчетная область
Рис. 1. Пример 2D области с топографией поверхности Fig. 1. Example of 2D domain with free surface topography
G(x}
Физическая область
Таким образом, задачу, поставленную в переменных (х, у), можно переписать в переменных (д, г), используя правила замены переменных. В результате мы получаем видоизмененную задачу, которую требуется решить в области простой геометрической формы (прямоугольник). Для такой области регулярная прямоугольная сетка будет обеспечивать хорошее согласование с границей.
Опишем формулы для установления взаимно-однозначного соответствия между областями X( ц, г) : ф2 ~> X2 в случае, когда в «физической» области 3 грани - прямые вдоль координатных направлений (см. рис. 1).
Сначала опишем отображение граней областей. Общие формулы принимают вид
А'В' -> АВ:х(Ч,г) = х(А) + (х(5) - х(А))(Ч - Ч(А'))/(Ч(В') - ч(А')),у(Ч,г) = С(х(Ч,г)) (1)
В'С' -> ВС\х(ц,г) = х(В),у(ц,г) = у(В) + (у(С) - у(Л))(г - г(5'))/(г(С') - г(5'))
(2)
й'С' -> ЭС: х(Ч, г) = х(£>) + (х(С) - х(0))(д - д(0'))/(д(С') - д(О')), у(д, г) = У(С) (3)
А'Б' -> АБ: х(ц, г) = х(А),у(ц,г) = у (А) + (у (!>) - у (Л)) (г - г(А'))/(гф') -г(А')) (4)
Для удобства положим, что точка А' имеет координаты (0,0), а А '5 'С' является единич-
ным квадратом. Формулы (1)-(4) упрощаются:
А 'В ' - А В : х ( q ,0 ) = х (А ) + (х (В ) -х (А ) ) ( q, 0) = G (х ( q, 0) ) (5)
В ' С' - В С : х ( 1, г) = х (В ) , у ( 1, г) = у (В ) + (у ( С) - у (А ) ) г (6)
D ' С' - D С : х ( q , 1 ) = х (D ) + (х ( С) - х (D ) ) q ,у ( q , 1 ) = у ( С) (7)
А 'D' - А D: х( 0,г) = х (А ) ,у ( 0,г) = у (А ) + (у (D) - у (А ) ) г (8)
Далее, зная отображения границ областей друг на друга, можно эти формулы интерполировать внутрь области, тем самым получив полное отображение самих областей. Для управления свойствами сетки, такими как сгущение и разрежение узлов, контроль над перехлестом ячеек, используются переходные функции:
а1, о ( С) = ( 1 + 2 С?) (q - 1) 2 0(= 1 — < 0(, < ^ = q( 1 - ^ а^=
С с -1 ) с2 *,
а^ 0 (г) = ( 1 + 2 г) (г — 1) 2, а^0(г) = 1 — а^ 0(г) , а2 х(г) = г ( 1 — г) 2 к, а£ 1 (г) = (г — 1)г2^ (9)
Здесь параметр 2£ - четная степень, являются нововведением автора. Стоит отметить, что этот параметры подбираются отдельно для каждой физической области в зависимости от кривизны поверхности и размера области. Контроль над перехлестом ячеек осуществляется при помощи Якобиана / = ХдУг — хгУд. Во всей области должно выполняться условие
/ > 0. Здесь и далее используем обозначения для производных ^ = ^">/у = =
д£ г- дГ г- г- д2г
дд'7г дг'7С1 дС2' дг ддаг
Формулы для отображения единичного квадрата на физическую область
(х(Ч,г),у(Ч,г)у.(22 ^Х2
выглядят следующим образом:
х(ц,г) = Б\ц,г) + а10(г)[х(д, 0) -51(?,0)] + а1г{г)[хг{4,0) - Б^ц, 0)] +
а 20 (г) [х ( ?0) —51 (СО)] + а?, 1 (г) ^ (с0) — 51 (с0)] (10)
у0?,г) = 52((7,г) + а£0(г)[у(<7,0)-Б2(ц, 0)] + а1г{г)\уг{4,0) - Б2(ц, 0)] +
а 2 0 (г) [у (СО) —52 (СО) ] + а^ 1 (г) [у. (?0) —5Г2 (?0) ] (11)
где
5 1(?г) = а110 (С)х (0' г) + а£ 1 (?)х^ ( 0' г) + а£0 ( с) х( 1 г) + а£ 1 (?) х^ ( 1 г) ,
52( сг) = а 1 0( у ( 0' г) + а 1 1 ( Уд ( 0' г) + а1 0 ( У ( 1 г) + а1 1 ( Уд ( 1 г) , #(<7, г) =
а 1 0 ( С) хг ( 0' г) + а ^ 1 ( ?) хдг (0' г) + а£0 (?)х( 1 г) + а^ 1 (с) хдг ( 1 г) , 52((7,Г) =
а1, 0 ( С)Уг ( 0' г) + а1, 1( с)Уд г(0' г)+а2, 0(С)У ( 1 г)+а2, 1 ^Удг ( 1 г).
3. Постановка и перенос задачи 3.1. Постановка задачи в «физической» области
Моделирование процесса распространения волновых полей проводится на основе численного решения полной линейной 2D системы теории упругости, записанной в смещениях. Физическая область предполагается с криволинейной свободной поверхностью, среда - изотропная. Граница области разделяется на две части: дГ - две боковые и нижняя грани и д5 -свободная криволинейная поверхность. В «физической» области (х, у) система принимает вид
д д Р и с С = ^ ((А + 2 /) их + Аиу) + — (/ (их + Му) ) + ^ (12)
д д Р П С = (/ + Му)) + — ((А + 2 /) их + Аиу) + ^2 (13)
Граничные условия на д Г:
и | аг = и I аг = 0 (14)
Начальные условия при t = 0:
и (х,у, 0 ) = и (х,у, 0) = 0, иС = иС = 0 (15)
Условия на свободной поверхности 35:
пх((А + 2 [1)их + Ли/у) + пу(д(и/х + иу)) = 0, + иу)) + пу((А + 2 [1)их +
А иуу=0 (16)
Здесь (и (х, у, £:)(х, у, £:) ) - вектор смещения вдоль координатных осей ОХ и OY, Р (х, у) - плотность, Ц, ( х, у) - скорость продольных волн в среде, ( х, у) - скорость поперечных волн в среде (скорости определяются параметрами среды)
2 К;2 ( х,у)) /р ( х,у), / ( х,у) = К;2( х,у)/р ( х,у) - параметры Ламе, д 5 - свободная поверхность (на рис. 1 изображена как участок функции С( х) , ограниченный точками А и В), д Г - граница области, не включая свободную поверхность, [пх, пу]т - единичный вектор нормали к свободной поверхности д 5. ^ = х, у, £:), F2 = F2 (х, у, £:) - функции, отвечающие за работу источника возмущений в среде.
3.2. Перенос исходной задачи на «расчетную» область
Сформулирована математическая постановка задачи в «физической» области при помощи уравнений (12)-(16). Теперь, используя принципы замены переменных, необходимо эти уравнения переписать в терминах . Учитывая замену , получа-
ем , . Важно сохранить ди-
вергентный вид системы, чтобы в дальнейшем применить метод баланса, описанный в [10]. Аналогично работе [14], система (12), (13) преобразуется в систему вида
/ри С С=^а1+|:а2+/ ^ (17)
dq
dr
(18)
Граничные условия на д Г:
и\дг = w\дг = О
(19)
и начальные условия при t = 0:
u(q, г, 0) = w(q, г, 0) = 0, wt = ut = 0
(20)
остаются неизменными.
На свободной поверхности д S условия примут вид
<г2 = 0 , <г4 = 0 ,
(21)
= уг//, г = — у«?//, qy = — хг//, г = х?// - метрические коэффициенты, / = х «?уг — хгу? - Якобиан преобразования.
= + 2 ß)(qxuq + rxur) + Ä(qywq + rywr)] + ]qy\p.{qxwq + rxwr) +
д ( qyuq + TyUr) ] ,
ö"2 = /г* [(A + 2 ß)(qxuq + rxur) + À(qywq + rywr)] + /гу[д(<7х^ + rxwr) + д ( qyu? + гу^г) ] ,
ö"3 = Jqx\n(.4xWq + w) + K4yUq + ryur)] + /qy[(A + 2 Li)(qxiiq + TxXlr) + A ( qywq + гуИ/у) ] ,
04 = + rxwr) + ß(qyuq + ryur)] + /ry[(A + 2ß)(qxuq + rxur) + Ä(qywt
Для компонент, отвечающих за источник возмущений, также производим замену пере-
менных , .
Далее для системы (17), (18) с граничными и начальными условиям (19)-(21) опишем созданную конечно-разностную схему.
В данном параграфе предлагается явная по пространству и времени конечно-разностная схема для численного решения задачи (17)-(21) внутри «расчетной» области. На (рис. 2) изображена стандартная ячейка с представленной в ней информацией. Численному восстановлению подлежат элементы , расположенные в центре ячейки. Стоит отметить, что для получения большей точности при расчетах все метрические коэффициенты и якобиан рассчитываются на сетке с половинным шагом. Как следствие, для них предусмотрена "удвоенная" индексация. Например, узлу с и^ ^ у' у будет соответствовать /2 ^ (2 у' 2 1 , 2 у. Для нецелого значения узла получаем соответственно .
+ rywr)]
4. Численная реализация с использованием явной разностной схемы
Рис. 2. Расположение элементов конечно-разностной схемы Fig. 2. Elements' positioning in finite-difference scheme
Для удобства считаем пространственный шаг сетки одинаковым по всем направлениям Л = Л^ = Лг (для метрических коэффициентов пространственный шаг берем Л / 2 ), временной шаг обозначим . Вся расчетная область покрывается регулярной квадратной сеткой М XI ячеек ( 2 М X 2 I ) для метрических коэффициентов и якобиана. Для произвольной функции А ( ц, г, £:) обозначим А "у = А ( / Л, у Л, пт) , а для Якобиана и метрических
коэффициентов введем обозначения разностных аналогов дифферен-
циальных операторов:
Ясс [А] "у = т- 2 (А£/ 1 - 2А?у + А?-^ = А „ ( /Л,)Л,пт) , 01 [А]?- 1/2,у = Л- 1 (А£у-А?_ 1,у) = А,(£Л-Л/2,уЛ,пт) , Я [А]?у-1/2 = Л-1 (А"у — А"у _ 1 ) = А г (¿Л,;Л —Л/2,пт),
[А] ?у = (Я [А] ?+ 1/2,у +01 [А] ?- 1/2,у) / 2 , 02° [А ] ?у = (02 [А ] ?,у+1/2+02 [А ] ?,у - 1/2) / 2 ,
01 [0] 2 £, 2 у = Л - 1 (02 I+1 , 2 у — 02 £ - 1,2у0 ,
02 [0 ] 2 £, 2 у = Л - 1 (02 £ , 2 у + 1 — 02 £,2у - 1 ).
В этих обозначениях схема для (17), (18) примет вид
/2£,2уР£,у0«№у = 01°К]£!у +/2£,2У01£,у. I = 2.....М - 1,] = 1.....I -
(22)
J2i.2jPi.jDtt = 01° ВД + 02° К] "у + /2£,2у02 £,у , I = 2.....М~ 1,] =
1,...,! — 1,п > 2 (23)
Граничные условия на 9 Г:
<у = <у = о,<у = W^,y = о,= = О,= wl = 0, i = 1,....M,j =
X..V Д ^ о (24)
Начальные условия:
%°у = = 0, и 1 j = wjj = 0, i = 1,..., M J = 1,..., L (25)
остаются неизменными.
Граница свободной поверхности 95 проходит вдоль граней ячеек с индексом j = 0. На 9 5 условия примут вид
<*2i,l/2 = /2ЦГ,2ц[(Я + г/О^С^дС^М?! + öi №)/2 + rx2U1D2[u\l1/2) +
+öi°Mi!o)/2 +ry2UD2[w]£1/2)] + /2u9y2uKi(9x2u(öiMu +ßiMi!o)/2 +rx2UD2[w][;i/2) +
Mi, i ( qy 2 i, 1 (ö° [и] ¿1 + [U] 1°) / 2 + ry 2 i, ^ [u] ?1/2 ) ] = 0, i = 2.....M , (26)
ßi,i(4y2u(ßiMuy + öl M"o)/2 + ry2UÖ2MÜ1/2)] + /2Ц?У2Ц[(Я + 2M)u(qX2u(öi°Mu + fli №)/2 + Г*2иЯ2[и]£1/2) + A, i(qy 2 i, 1 (Ö1° [w] h + 0°° [w] ) /2 + ry 2 i, 1 02 [w] ¿д /2 ) ] =0, i = 2.....M. (27)
Из (26), (27) можно получить uf° , w™° - значения в «фиктивных» узлах, находящихся вне области. Эти узлы вводятся для вычисления (Г^ _ i / 2 1, Ö31 i _ 1 /2 1 и дальнейшей реализации конечно-разностной схемы вблизи свободной поверхности. Такая реализация граничных условий является нововведением автора.
Для (1 , Г2 , Гз, °4 внутри области разностные аналоги примут вид
°li-l/2,i = J2i-l,2j4x2i-l,2j[(.Ä + 2^)j_1/2y(qx2j_12yD1[u]f_jy2y + rx2j_12y(D2 [u]"y + Ö2°Mf-i,y)/2) + Vl/2,y(qy2i-l,2yßlMr_1/2,y +ГУ2.-1,2У(020№ + ^2° M?_W)/2)] + J2i-12jqy2i-12j\j*i-l/2j(qx2i-12jDi[w]?_1/2j + rx2i-l,2j(D^[w]fj + [w]f_iy)/2) + M i _ 1/2, j ( qy 2 i _ 1, 2jß 1 [u] ?_ 1/2, j + ry 2 i _ 1, 2 j (^ [U¿j + [u]¿_ 1,j) / 2 ) ] , i = 2.....M,
°2i,j-1/2 = J2i,2j-irX2i,2j-l№ + 2M)jJ-l/2(qx2i,2y-l(£,l°Ml!y +ßl°Ml!y-l)/2
+ Гх21,2у-1^2[и]1|у-1/2) + hj-i/2(qy 2i,2y-i(ßi № + Di М^-.о/г + ry 2иу-102М£у_1/2)] + /гиу-^угиу-^О'-^г^хгиу-^^И^у + öiM^-.^/Z + rx2ii2j_1D2[w]fj_1/2) + Mi,j _ 1/2 ( qy 2 i , 2 j _ 1(0°° [u] ¿j + [u] ¿j _ 1) / 2 + ry 2 i, 2 j _ 1^2 [u] ¿j _ 1/2 ) ] , j = 2.....L,
03i-l/2,j = J2i-l,2jqx2i-l,2j[Lli-l/2,j(qx2i-l,2jDl[W]1i-l/2,j + Гх 2i-l,2j (P2. [W]i,j
+ D2°[w]f_1,y)/2 +
21-1,2^1 W-1/2J + ry2i-1<2y(D2°[M]^ + £>2°№-1,у)/2)] + J2i-l,2jqy 2i-l,2j[(.Ä + 2/^)j-i/2,y(i7x2t-l,2yDl[li]f-l/2,y + rx 2i-l,2j(P2. [UVi,j
+ D2°Mf-i,y)/2)) +
A_ 1 / 2, j ( qy 2 i _ 1, 2j£>1 [w] 11 / 2, j + ry 2 i _ 1 , 2 j (D°° [w] ¿j + D2° [w] ¿_ 1,j) /2 ) ] , i = 2.....M,
a4i,j-1/2 - J2i,2j-lrx2i,2j-l\P-i,j-l/2(.4x2i,2j-l(Dl\.WVi,j + Mi!j-l)/2
/21,2 ]—l4y 2i,2 у —1 [(A + ^)i,j-l/2(qX2i,2j-l(Dl[u]fj +£>l°№y-i)/2 + 2i,2y —102 ¿j'—1/2) + 1/2 2 t, 2 у - l(01 M Пу + 0? [w] ! ) /2 + ry 2 i<2/_ x £>2 [w] _ 1/2) ] , j = 2.....L.
Здесь 2 t,2y = 02 [у] 2 t, 2 у//21, 2 у' rx 2 t,2y = - 01 M 2 t , 2 у/72t, 2у'' 2 t, 2у = - 02 [x] 2 1,2у-//2i, 2 у'^ Гу 2 t, 2 у = 01 [x] 2 t, 2 у //2 t, 2 у - разностные аналоги для метрических коэффициентов. /2 t, 2 у = 01 M 2 t, 2 у0 2 [у] 2 t, 2 у - 02 M 2 ¿,2у'01 [у] 2 ¿,2у - Якобиан преобразования.
При этом на свободной поверхности 55 применяются односторонние формулы для вычисления производных со вторым порядком:
2 t, 1 = ( - 3У2 t, 1 + 4У2 t, 2 - У2 t, 3 ) / (У2 t, 1h) ' rx 2 t, 1 = - 01 [У] 2 t, 1 #2 t, 1' <7у 2 t, 1 = ( 3 X2 t, 1 - 4X2 t, 2 + ^2 t, 3) / (/2 t, 1^) , Гу 2 t, 1 = 0 1 M 2 t, 1//2 t, 1, / t, 1 = 01 [x] 2 t, 1 ( - 3 У2 t, 1 + 4У2 t, 2 - У2 t, 3 ) h_ 1 + 0 1 [у] 2 t, 1 ( 3 X2 t, 1 - 4X2 t, 2 + X2 t, 3 ) h _ 1
Для параметров среды:
(Я + 2\i\_xjlij = (Я/2 + fi)Uj + (Я/2 + ^¿_1,у,Я i_1/2J = (Aij + Xi-^j)/ 2,Hi _ 1 / 2, у = (/Чу + /t - 1,у) / 2 ,
(Я + 2д)и_1/2 = (Я/2 + \1)ц + (Я/2 + и_1/2 = + XUj-i)/2,^y_1/2 =
(/,у t, у - 1 ) / 2 .
Представленная схема (22)-(27) обладает вторым порядком точности по пространству и времени во всей области, кроме границ раздела двух сред, в таком случаем схема обладает первым порядком точности по пространству и вторым порядком точности по времени. Критерием устойчивости схемы является видоизмененное условие куранта l^T / h <
/m t п/ V2, где /m t п - минимальное значение якобиана во всей области. Физический смысл условия - за один временной шаг T волна не должна продвинуться по среде дальше, чем диаметр самой маленькой ячейки (максимальное расстояние между двумя точками ячейки).
5. Программная параллельная реализация алгоритма
Учитывая размеры реальной задачи и объем вычислений, расчеты на персональном компьютере займут много времени. Для сокращения времени численных расчетов задействуются средства параллельного программирования, а именно технология MPI (message passing interface) [11].
На рис. 3 изображен способ разбиения области вычислений на подобласти. Каждому процессу (ядру, если не задействована технология hyper threading) определяется его подобласть для расчетов. Все процессы работают одновременно. Ввиду специфики выбранной конечно-разностной схемы, каждому процессу для расчетов на границе подобласти требуются данные с соседней подобласти. Для этого к каждой подобласти добавляются "слои перекрытия" (на рис. 3 обозначены красным цветом), которые используются для обмена данными между соседними процессами. Программно обмен организуется при помощи средств MPI, в данном случае используются блокирующие функции-пересылки MPI_Send, MPI_Recv. Таким образом обеспечивается синхронизация параллельной работы процессов и в слоях перекрытия находятся актуальные данные на каждой временной итерации. В целом, такой способ организации вычислений позволяет существенно сократить время расчетов. Параллельная про-
грамма генерации криволинейной 2D сетки, а также параллельная программа для численного решения исходной 2D задачи созданы автором при помощи языка Fortran.
О
Рис. 3. Одномерная декомпозиция «расчетной» области Fig. 3. ID-decomposition for calculation domain
6. Модель для расчетов
Для численного моделирования распространения волн рассматривается модель неоднородной среды, представленной на (рис. 4), не в масштабе. Данные для модели используются из реальных наблюдений одного из участков Восточной Сибири. Участок имеет ширину 48 910 метров, глубина от поверхности варьируется между 2 400 метров и 2 650 метров. В разных слоях скорость продольных волн варьируется от 800 м/с до 7000 м/с, скорость поперечных волн - от 400 м/с до 3200 м/с, плотность р -от 1500 кг/м 3 до 2600 кг/м3.
Источник волнового поля - типа "центр давления". В качестве функции сигнала используется импульс Рикера Я ( = ( 1 — 2 2ш2( t — ^ ) 2) е- 2 ш ((:-^ с несущей частотой ш = 3 0 Гц и временным сдвигом t 0 = 0 , 1 секунды. Тогда функции, отвечающие за работу источника в (12), (13) и (17), (18), принимают вид
^ = [ъ ( ч 0 ,г0)э * ( «7;' - - * ) + гх я (0 ,
02 0 = [ чу ( ч 0 .Го )^г-о) + Гу ( ч о ,Го)££(£-^] я ( о .
«Физическая» область покрывается криволинейной сеткой, локально-ортогональной возле свободной поверхности (см. рис. 4). Именно такой способ построения сетки позволяет удобно и более точно реализовать граничные условия по нормали (16). Сетка построена при помощи формул (10), (11), для управляющих функций (9) параметр к = 1 0. Для удовлетворения условиям локальной ортогональности сетки возле поверхности накладываются условия в уравнениях (10)-(11) х^ ( ч . 0 ) = уг ( ч. 0 ) . хг ( ч. 0 ) = — у^( ч . 0) . В зависимости кривизны поверхности ячейки будут разниться размерами и кривизной.
Соответствующая «расчетная» область имеет размер 489 10 X 3000 метров и 2445 5 X 1 500 квадратных ячеек (для метрических коэффициентов 489 10 X 3 000 ячеек). Размер ячейки - метра ( км). Временной шаг , время моделирова-
ния - 3 секунды, число временных шагов в схеме - 30000.
Рис. 4. Модель неоднородной среды и соответствующая локально-ортогональная сетка Fig. 4. The Model of heterogeneous medium and corresponding locally-orthogonal mesh
Для области, изображенной на рис. 4, требуется провести серию расчетов с фиксированной системой наблюдения и разным положением источника для каждого расчета. Система наблюдения представлена серией приемников, располагающихся вдоль всей свободной поверхности, в центрах приграничных ячеек и на расстоянии 5 ячеек (около 10 метров) друг от друга в горизонтальном направлении (рис. 5). Результаты расчетов для каждого положения источника записываются в отдельные файлы для каждой компоненты волнового поля (и, и) - трассы. В них помещаются значения смещений и или и, записанные в каждом приемнике с временным шагом 0,002 сек.
Источник для каждого расчета располагается на глубине 5 метров и на горизонтальном смещении 100 метров относительно местоположения источника предыдущего расчета. Таким образом, требуется провести 488 расчетов.
Рис. 5. Система приемников (трассы) и источник Fig. 5. System of receivers (traces) and the source
Для численной реализации необходимо избавиться от отражений волн внутрь области от границы дГ (левой, правой и нижней границ). Для этого можно использовать два способа: создание поглощающих слоев PML [12], которые будут «гасить» отражения от границ, либо расширить расчетную область достаточно, чтобы за время моделирования отражения от границ не успели дойти до системы наблюдения вдоль свободной поверхности. Автором был выбран второй вариант как наименее трудозатратный. Общее число ресурсов, выделенных для проведения численных расчетов, позволяло в общем быстрее провести расчеты именно в такой постановке. «Расчетная» область была расширена на 5 000 ячеек влево, на 5 000 ячеек вправо и на 3 500 ячеек вниз (рис. 6), значения параметров среды в расширении области "продолжаются" те же, что в примыкающей расчетной области слева, справа и снизу, соответственно.
Рис. 6. Слева - исходная «расчетная» область, справа - расширенная «расчетная» область Fig. 6. On the left - original calculation domain, on the right - extended calculation domain
Таким образом, теперь необходимо реализовать разностную схему (22)-(27) на расширенной "расчетной" области размером 34455 х 5000 ячеек (сетка соответственно 68910 х 10000). Система наблюдения и система источников остаются в пределах изначальной «расчетной» области.
7. Результаты численного моделирования
Для расчетов было задействовано 5 узлов кластера НКС-1П ССКЦ СО РАН. Каждый узел включает в себя два 16-ядерных процессора Intel Xeon E5-2697v4 (2.6 ГГц) и 256 Гб оперативной памяти. Каждому узлу было выделено провести 97 или 98 расчетов (488 суммарно на 5 узлах). Каждый расчет проводился на 32 ядрах с использованием технологии avx2,
Рис. 7. Пример трассы для компоненты w, источник смещен на 33,8 км от левой границы области Fig. 7. Example of the trace for w component, source offset from left border of the domain is 33,8 km
hyperthreading отключен (1 поток - 1 вычислительное ядро). Точность для расчетов использовалась одинарная (32 бита на одно число с плавающей точкой). Общее время для всех 488 положений источника заняло чуть меньше 2 суток. Выходные данные - трассы для компонент вектора смещений ( и, w) . На рис. 7 изображен пример такой трассы для компоненты w. По горизонтали слева направо - расстояние от левой границы области в километрах, по вертикали вниз - время в секундах. На цветовом баре - значение отклонения от положения равновесия.
В таком виде данные могут использоваться для верификации алгоритма решения обратной задачи - восстановления строения среды по имеющимся трассам.
8. Заключение
Основной практический интерес геофизики - восстановление строения среды. Так можно понять, есть ли на глубине резервуары углеводородов. Изучение строения среды при помощи взятия многочисленных проб и бурения скважин - очень финансово и время затратное занятие. Поэтому применение средств математического моделирования представляется хорошей альтернативой. Если расчеты покажут результат, схожий с исходными данными среды, то такой алгоритм восстановления среды можно применять для областей, которые не так хорошо исследованы, как представленная на рис. 4, в том числе новые участки. Верификацию алгоритма стоит проводить на реальных характерных данных строения сред, но реальные записи трасс не предоставляются, поскольку зачастую они являются объектом коммерческой тайны. В таком случае на помощь приходит результат данной работы.
В работе был представлен параллельный алгоритм моделирования 2D волновых полей и результаты численного моделирования с использованием этого алгоритма. Новизной работы являются: алгоритм построения локально-ортогональной сетки, новый способ аппроксимации граничных условий с использованием такой сетки. В дальнейшем эта работа послужит основой для решения аналогичной 3D-задачи, результаты которой могут быть применены в реальных исследованиях.
Список литературы
1. Лисейкин В. Д. Метод генерации сетки. 2-е изд. Берлин: Спрингер, 2010. ISBN 978-90481-2912-6
2. Васева И. А., Лисейкин В. Д., Лиханова Ю. В., Мороков Ю. Н. Эллиптический метод построения адаптивных пространственных сеток // Рус. Джей Нумер. Анал. Математика. Моделирование. 2009. Т. 1. C. 65-78. DOI 10.1515/RJNAMM.2009.006
3. Лисейкин В. Д. Разностные сетки. Теория и приложения. Новосибирск: Изд-во СО РАН, 2014. 245 с.
4. Appelo D., Petersson N. A. A stable finite difference method for the elastic wave equation on complex geometries with free surfaces. Communications in computational physics 5, 2007, no. 1, p. 84-107.
5. Jose M. Carcione. The wave equation in generalized coordinates. Geophysics 59, 1994, no. 12, p. 1911-1919. DOI 10.1190/1.1443578
6. Титов П. А. Алгоритм и программа моделирования 2D-волновых полей в областях с криволинейной свободной поверхностью // Материалы конференции «Научный сервис в сети Интернет - 2014». Новороссийск, Абрау-Дюрсо, 2014. С. 446-455.
7. Титов П. А. Моделирование упругих волн в средах со сложной топографией свободной поверхности // Вестник НГУ. Серия: Информационные технологии. 2018. Т. 16, № 4. С.153-166.
8. Титов П. А. Технология моделирования распространения упругих волн в средах со сложной трехмерной геометрией поверхности с помощью высокопроизводительного кластера // Российские суперкомпьютерные дни. 2016. С. 1020-1031.
9. Титов П. моделирование трехмерных волновых полей в сложных топографических средах // Российские суперкомпьютерные дни. CCIS 1129. 2019. С. 451-462. DOI 10.1007/ 978-3-030-36592-9_37
10. Вишневский Д., Лисица В., Чеверда В., Решетова Г. Численное исследование погрешностей сопряжения конечно-разностного моделирования сейсмических волн // Геофизика. 2014. T. 79 (4). C. 219-232. DOI 10.1190/geo2013-0299.1
11. Gropp W., Huss-Lederman S. et al. MPI - The complete reference. Cambridge, MIT Press, 1998, vol. 2, 205 p.
12. Collino F., Tsogka C. Application of the Perfectly Matched Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media. Geophysics, 2001, vol. 66, p. 294-307. DOI 10.1190/1.1444908
References
1. Liseikin V. D. Grid generation method. 2nd ed. Berlin, Springer, 2010. (in Russ.) ISBN 97890-481-2912-6
2. Vaseva I. A., Liseikin V. D., Likhanova Yu. V., Morokov Yu. N. An elliptic method for construction of adaptive spatial grids. Russ. J. Numer. Anal. Math. Modelling, 2009, vol. 1, p. 65-78. (in Russ.) DOI 10.1515/RJNAMM.2009.006
3. Liseikin V. Differential meshes. Theory and applications. Novosibirsk, SB RAS Publ., 2014, 254 p. (in Russ.)
4. Appelo D., Petersson N. A. A stable finite difference method for the elastic wave equation on complex geometries with free surfaces. Communications in computational physics 5, 2007, no. 1, p. 84-107.
5. Jose M. Carcione. The wave equation in generalized coordinates. Geophysics 59, 1994, no. 12, p. 1911-1919. DOI 10.1190/1.1443578
6. Titov P. An algorithm and a program for simulation of 2D-wave fields in areas with a curved free surface. In: Proceedings of "Scientific service on the Internet - 2014" conference. Novorossiysk, Abrau-Dyurso, 2014, p. 446-455. (in Russ.)
7. Titov P. Modeling of elastic waves in media with a complex free surface topography. Vestnik NSU. Series: Information Technologies, 2018, vol. 16, no. 4, p. 153-166. (in Russ.)
8. Titov P. A technology of modeling of elastic waves propagation in media with complex 3D geometry of the surface with the aid of high performance cluster. Russian Supercomputing Days, 2016, p. 1020-1031. (in Russ.)
9. Titov P. The Simulation of 3D Wave Fields in Complex Topography Media. Russian Supercomputing Days, CCIS 1129, 2019, p. 451-462. (in Russ.) DOI 10.1007/978-3-030-36592-9_37
10. Vishnevsky D., Lisitsa V., Tcheverda V., Reshetova G. Numerical study of the interface errors of finite-difference simulations of seismic waves. Geophysics, 2014, vol. 79 (4), p. 219232. (in Russ.) DOI 10.1190/geo2013-0299.1
11. Gropp W., Huss-Lederman S. et al. MPI - The complete reference. Cambridge, MIT Press, 1998, vol. 2, 205 p.
12. Collino F., Tsogka C. Application of the Perfectly Matched Absorbing Layer Model to the Linear Elastodynamic Problem in Anisotropic Heterogeneous Media. Geophysics, 2001, vol. 66, p. 294-307. DOI 10.1190/1.1444908
Материал поступил в редколлегию Received 30.06.2020
Сведения об авторе
Титов Павел Андреевич, младший научный сотрудник лаборатории Суперкомпьютерного моделирования, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Россия)
Information about the Author
Pavel A. Titov, junior researcher, Supercomputing Modeling laboratory, Institute of Computational mathematics and Mathematical Geophysics SB RAS (Novosibirsk, Russian Federation)
