=Paper=
{{Paper
|id=Vol-1576/173
|storemode=property
|title=Mathematical modeling of acoustic antenna radiation on multi-processor system
|pdfUrl=https://ceur-ws.org/Vol-1576/173.pdf
|volume=Vol-1576
|authors=Aleksandr Sukhinov,Aleksandr Chistyakov,Oleg Savitski,Alena Semenyakina,Alla Nikitina
}}
==Mathematical modeling of acoustic antenna radiation on multi-processor system==
https://ceur-ws.org/Vol-1576/173.pdf
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
Математическое моделирование излучения акустической
антенной на многопроцессорной системе*
А.И. Сухинов1, А.Е. Чистяков2, О.А. Савицкий 3, А.А. Семенякина2, А.В. Никитина2
Донской государственный технический университет1,
Научно-исследовательский институт многопроцессорных вычислительных систем име-
ни академика А.В. Каляева Южного Федерального Университета2,
Акустический институт имени Н. Н. Андреева РАН3
Разработка оптимальных конструкций акустических антенн в настоящее время явля-
ется актуальной задачей. Для этого необходимо рассчитать характеристики антенных
решеток. В данной работе построена дискретная математическая модель распростра-
нения акустических волн антенной решеткой. Для повышения реальной точности
решений были использованы сетки, учитывающие заполненность расчетных ячеек.
Дискретная модель построена на основе метода баланса. Приведены результаты чис-
ленных экспериментов излучаемых антенн с различными характеристиками направ-
ленности, работающих в заданном диапазоне частот. Для численного решения данно-
го класса задач с точностью, требуемой в задачах проектирования акустических ан-
тенн, с учетом краевых эффектов требуются сетки от 105 узлов для двумерных моде-
лей до 109 узлов в трехмерных моделях. Проведение многовариантных расчетов и
ограниченность вычислительных ресурсов и, в особенности, объема памяти последо-
вательных ЭВМ приводят к необходимости построения эффективных параллельных
вычислительных методов решения данного класса задач на доступных проектным
организациям многопроцессорных системах, содержащих сотни-тысячи вычисли-
тельных ядер.
Ключевые слова: волновая задача, акустические волны, MPI.
1. Введение
Запишем исходную систему уравнений гидродинамики (уравнение Эйлера и неразрывно-
сти) в дифференциальной форме:
vt v v 1p, t v 0 ,
где v – вектор скорости; – плотность; p – давление.
Применим соответствующие допущения [1], тогда система уравнений примет вид:
vt 1p, pt / c 2 v 0 ,
где c – скорость звука.
Применяя к первому уравнению оператор div, а ко второму – / t , вычитая результаты,
получим уравнение звука в неподвижной неоднородной и нестационарной среде:
p / c p .
t
2
t
1
Простейшим случаем является уравнение распространения волны в однородной непо-
движной среде:
ptt c 2 p .
2. Постановка задачи
Требуется найти решение неоднородного волнового уравнения [1–3]:
*
Работа выполнена при частичной поддержке Задания №2014/174 в рамках базовой части
государственного задания Минобрнауки России
699
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
ptt c 2 p f , (1)
удовлетворяющего начальным условиям:
p( x, y,0) 0 ( x, y) , pt( x, y, 0) 1 ( x, y) (2)
и граничным условиям:
p ( x, y , t ) 0 , при ( x, y ) , (3)
где f – функция-источник, – граница расчетной области.
Наиболее эффективными методами для решения подобного вида задач являются сеточные
методы.
3. Решение задачи
Расчетная область вписана в прямоугольник и покрыта равномерной прямоугольной рас-
четной сеткой t x y :
t {t n nht , 0 n N t 1, lt ht N t 1}, x {xi ihx , 0 i N x 1, lx hx N x 1} ,
y { y j jhy , 0 j N y 1, l y hy N y 1} ,
где n, i, j – индексы по временной координате и пространственным координатным направле-
ниям Ox , Oy соответственно; ht , hx , hy – шаги по временной координате и пространственным
координатным направлениям; N t , N x , N y – количество узлов по временной координате и про-
странственным координатным направлениям; lt , lx , l y – длина расчетной области по временной
координате и пространственным координатным направлениям.
Для получения дискретной модели использован интегро-интерполяционный метод [4].
Разностная схема, аппроксимирующая уравнение (1):
pin, j 1 2 pin, j pin, j 1 pi 1, j 2 pi , j pi 1, j pi , j 1 2 pi , j pi , j 1
2
c2 2
c2 2
fi ,nj , (4)
h t h x h y
где pi , j 1 pin, j 1 1 1 2 pin, j 2 pin, j 1 ; 1 , 2 – веса схемы [5].
Для улучшения «гладкости» решения сеточного решения будем предполагать, что ячейки
заполнены не полностью [6]. Областью xy будем называть заполненную часть области
Dxy x xi 1/2 , xi 1/2 , y y j 1/2 , y j 1/2 .
Также введем обозначения для следующих областей:
D1 x xi , xi 1/2 , y y j 1/2 , y j 1/2 , D2 x xi 1/2 , xi , y y j 1/2 , y j 1/2 ,
D x x
3 i 1/2 , xi 1/2 , y y j , y j 1/2 , D x x
4 i 1/2 , xi 1/2 , y y j 1/2 , y .
j
Коэффициенты заполненности q0 , q1 , q2 , q3 , q4 для областей Dtxy , D1 , D2 , D3 , D4 вводятся
следующим образом:
S Dxy S Di
q0 , qi , i 1, 4 ,
S xy S i
где S – площадь соответствующей части области; i – заполненная часть области Di .
Дискретный аналог уравнения (1) в случае граничных условий в форме Дирихле ( p 0 )
примет вид:
pin, j 1 2 pin, j pin, j 1 pi 1, j 2 pi , j pi 1, j
ht2
q 0, i , j
ht2
c 2
min q1, i , j ,q 2, i , j hx2
700
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
pi , j 1 2 pi , j pi , j 1
c 2 min q3,i , j , q4,i , j 2
q0,i , j f i ,nj . (5)
h y
Дискретный аналог уравнения (1) в случае граничных условий в форме Неймона ( pn 0 )
примет вид:
pin, j 1 2 pin, j pin, j 1 pi 1, j pi , j pi , j pi 1, j
q0,i , j c 2
q1, i , j c 2
q 2, i , j
ht2 hx2 hx2
pi , j 1 pi , j pi , j pi , j 1
c 2 q3,i , j 2
c 2 q4,i , j 2
q0,i , j f i ,nj .
h y h y
4. Решение модельной задачи
Схему итерационного двухслойного модифицированного попеременно-треугольного мето-
да (МПТМ) запишем в виде [7–9]:
x n+1 x n n1w n , D R1 D 1 D R2 w n r n , r n Ax n f ,
где xn – вектор решения; w n – вектор поправки; А – оператор сеточного уравнения; D – диаго-
нальная часть оператора A; – итерационный параметр; R1, R2 – верхняя и нижняя треугольные
n
части оператора A; r – вектор невязки; f – правая часть сеточного уравнения.
Для параллельной реализации адаптивного МПТМ были использованы методы декомпо-
зиции области по одному направлению. Наиболее трудоемким расчетом с точки зрения постро-
ения параллельной программной реализации является расчет вектора поправки, который вы-
полняется в два шага:
1) D R1 y r , 2) D R2 w Dy .
n n n n
n
На первом шаге рассчитываются элементы вспомогательного вектора y снизу вверх, а
затем, зная его, на втором шаге находятся элементы вектора поправки w сверху вниз. Схемы n
расчета вспомогательного вектора и вектора поправки представлены на
рис. 1 (стрелками показаны направления счета и передачи, штриховкой отмечены расчетные
узлы).
Рис. 1. Схемы расчета вспомогательного вектора yn (а) и вектора поправки wn (б)
n
В схеме для расчета вектора y только первый процессор не требует дополнительной ин-
формации и может независимо от других процессоров посчитать свою часть области, осталь-
ные процессоры ждут от предыдущего передачи элементов, стоящих в начале строки. Передача
по одному элементу не оптимальна, так как появляются временные затраты, связанные с их ор-
ганизацией. Суммарное время на накладные расходы можно уменьшить путем увеличения объ-
ема передач. Данные рассуждения используются для расчета вектора поправки.
701
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
5. Параллельная реализация
При параллельной реализации использованы методы декомпозиции сеточных областей для
вычислительно трудоемких задач, учитывающие архитектуру и параметры многопроцессорной
вычислительной системы. Максимальная производительность МВС составляет 18,8 терафлопс.
В качестве вычислительных узлов используются 128 однотипных 16-ядерных Blade-серверов
HP ProLiant BL685c, каждый из которых оснащен четырьмя четырехъядерными процессорами
AMD Opteron 8356 2.3GHz и оперативной памятью в объеме 32 ГБ. В таблице 1 приведены
временные затраты выполнения одной итерации МПТМ на различных сетках, а также значения
ускорения и эффективности для различного числа вычислительных ядер. Для расчетов
использована технология распараллеливания MPI.
Таблица 1. Ускорение и эффективность работы параллельного варианта МПТМ
100x100 200x200 500x500 1000x1000 2000x2000 5000x5000
1 Время 0.001183 0.010633 0.026031 0.10584 0.381988 3.700073
Ускорение 1 1 1 1 1 1
Эффективность 1 1 1 1 1 1
2 Время 0.000446 0.003435 0.01932 0.05579 0.264114 1.880677
Ускорение 2.652 3.095 1.347 1.897 1.446 1.967
Эффективность 1.326 1.548 0.674 0.949 0.723 0.984
4 Время 0.000232 0.00106 0.005755 0.026683 0.132585 1.2655
Ускорение 5.099 10.031 4.523 3.967 2.881 2.924
Эффективность 1.275 2.508 1.131 0.992 0.72 0.731
8 Время 0.000179 0.000878 0.004379 0.023322 0.092771 0.489768
Ускорение 6.609 12.11 5.945 4.538 4.118 7.555
Эффективность 0.826 1.514 0.743 0.567 0.515 0.944
16 Время 0.000231 0.000407 0.001869 0.013105 0.085056 0.472151
Ускорение 5.121 26.125 13.928 8.076 4.491 7.837
Эффективность 0.32 1.633 0.87 0.505 0.281 0.49
32 Время 0.000365 0.0005 0.001404 0.008871 0.045913 0.318709
Ускорение 3.241 21.266 18.541 11.931 8.32 11.61
Эффективность 0.101 0.665 0.579 0.373 0.26 0.363
64 Время 0.000642 0.000748 0.001557 0.004189 0.026844 0.182296
Ускорение 1.843 14.215 16.719 25.266 14.23 20.297
Эффективность 0.029 0.222 0.261 0.395 0.222 0.317
128 Время - 0.001612 0.002064 0.003442 0.016437 0.076545
Ускорение - 6.596 12.612 30.75 23.24 48.338
Эффективность - 0.052 0.099 0.24 0.182 0.378
256 Время - - 0.005434 0.005446 0.012521 0.06318
Ускорение - - 4.79 19.434 30.349 58.563
Эффективность - - 0.019 0.076 0.119 0.229
512 Время - - - 0.009793 0.012362 0.058805
Ускорение - - - 10.808 30.9 62.921
Эффективность - - - 0.021 0.06 0.123
Из таблицы также видно, что для каждой из расчетных сеток ускорение принимает
наибольшее значение при определенном значении вычислителей и при дальнейшем увеличении
числа вычислительных ядер ускорение только уменьшается. Это связано с временными затра-
тами на обмен данными между вычислителями. Следует отметить, что адаптивный МПТМ
нашел свое применение при решении задач аэро-гидродинамики [10] и транспорта донных ма-
териалов [11–12].
702
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
6. Результаты численных экспериментов
Программный компонент для расчета картин поля излучения акустических волн включает
в себя следующие компоненты: блок ввода исходных данных; блок расчета геометрии; блок
расчета коэффициентов сеточных уравнений; блок расчета функций правых частей сеточных
уравнений; блок перехода на более грубую сетку; блок расчета сеточных уравнений на основе
модифицированного попеременно-треугольного метода; блок расчета положения расчетного
окна; блок учета граничных условий; блок вывода рассчитываемых функций давления; блок
вывода спектра; блок расчета фазы; блок расчета градиента фазы; блок расчета направленно-
сти; блок вывода направленности.
На рис. 2 представлены результаты математического моделирования излучения акустиче-
ских волн антеннами с различными характеристиками направленности. На рисунках показаны
значения колебаний давления. За единицу принято максимальное значение амплитуды колеба-
ний поля давления.
Рис. 2. Распространение акустических волн от антенн с различными характеристиками направ-
ленности
На рис. 3 представлена влияние наличия препядствия (неоднородностей) на картину
распространения акустической волны.
Рис. 3. Рассеяние акустических волн на препятствии
Для изучения отраженных сигналов удобно моделировать короткоимпульсные сигналы. На
рис. 4 представлены картины распространения одиночного акустического сигнала и его
отраженние от препятствия (неоднородности) в различные моменты времени.
703
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
Рис. 4. Отраженный акустический сигнал
7. Результаты численных экспериментов
Требуется рассчитать акустические поля антенны с рабочими частотами: 1.5 кГц, 2.0 кГц,
2.5 кГц, 3.0 кГц. Приведенное давление: 10 кПа, 20 кПа, 50 кПа, 100кПа. Рабочие поверхности
состоят из пьезокерамики ЦТС-19 и поляризованы радиально, электроды находятся на боковых
поверхностях. На рис. 5 приведена геометрия излучающей акустической антенны.
Рис. 5. Геометрия излучающей антенны
На рис. 6 приведены результаты расчета распространения колебательных процессов от ан-
тенны в случаях: непрерывного сигнала (а) и короткого импульса (б). Рабочая частота 3.0 кГц.
Размер окна моделирования составлял 5×12 м. Временной интервал составлял 3 мс.
704
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
Рис. 6. Результаты расчета распространения колебательных процессов от антенны: непрерывный сигнал
(а) и короткий импульс (б)
На рис. 7 приведены результаты расчета интенсивности и направления распространения
колебательных процессов в ближнем поле антенны и вблизи элементов антенны в случае не-
прерывного сигнала.
Рис. 7. Результаты расчета направления распространения и интенсивности колебательных про-
цессов вблизи элементов антенны (а) и в ближнем поле антенны (б)
Для ускорения работы разработанного программного комплекса при расчете акустических
полей в дальнем поле антенны переход на более грубую сетку.
Рис. 8. Результаты расчета распространения колебательных процессов от антенны: временной интервал 3
мс (а) и 12 мс (б)
На рис. 8(а) представлены результаты расчета акустического сигнала в ближнем поле ан-
тенны. Временной интервал составлял 3 мс. После чего осуществлен переход на сетку с 4 раза
большими шагами по пространственным координатам. На рис. 8(б) приведены результаты рас-
705
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
чета картины распространения звуковой волны через 12 мс после начала работы акустической
антенны. Размер окна моделирования составлял 5×12 м на подробной сетке и 20×48 м на грубой
сетке, что вполне достаточно для расчета интенсивности и направления распространения коле-
бательных процессов в ближнем поле и на элементах конструкции антенны.
Для расчета дальних полей акустической антенны предлагается расчетное окно сделать по-
движным и вычислять его местоположения в пространстве. Данный подход позволяет суще-
ственно сократить время счета распространения звуковых волн на большие расстояния.
Заключение
Работа посвящена изучению волновых колебаний, а также построению комплекса про-
грамм, предназначенного для описания волновых процессов излучения в линейных антенных
решетках с изменяемой геометрией. В основе предложенной математической модели лежит
неоднородное волновое уравнение с соответствующими начальными и граничными условиями.
Для решения поставленной задачи был выбран метод сеток. Дискретная модель была построена
при помощи интегро-интерполяционного метода, при этом осуществлялся учет заполненности
расчетных ячеек, что гарантировало выполнение основных законов сохранения на дискретном
уровне. Полученные сеточные уравнения решены адаптивным модифицированным поперемен-
но-треугольным методом вариационного типа, который имеет наилучшую скорость сходимости
в классе двухслойных итерационных методов.
Литература
1. Марков Г.Т., Чаплин А.Ф. Возбуждение электромагнитных волн. М.: Радио и связь, 1983.
2. Владимиров В.С. Уравнения математической физики. Учебник для физич. и механико-
математ. спец. вузов. 4-е изд., испр. и доп. М.: Наука, 1981.
3. Сухинов А.И., Зуев В.Н., Семенистый В.В. Уравнения математической физики. Таганрог:
ТРТУ, 2005.
4. Самарский А.А. Теория разностных схем. М. Наука, 1989.
5. Сухинов А.И., Чистяков А.Е., Шишеня А.В. Оценка погрешности решения уравнения диф-
фузии на основе схем с весами // Математическое моделирование. 2013. Т. 25. № 11. С. 53–
64.
6. Сухинов А.И., Чистяков А.Е., Фоменко Н.А. Методика построения разностных схем для
задачи диффузии-конвекции-реакции, учитывающих степень заполненности контрольных
ячеек // Известия Южного федерального университета. Технические науки. 2013. № 4. С.
87–98.
7. Самарский А.А. Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978 –588
с.
8. Коновалов А.Н. Метод скорейшего спуска с адаптивным попеременно-треугольным пере-
обуславливателем // Дифференциальные уравнения. 2004. Т. 40. № 7. 953 с.
9. Сухинов А.И., Чистяков А.Е. Адаптивный модифицированный попеременно-треугольный
итерационный метод для решения сеточных уравнений с несамосопряженным оператором//
Математическое моделирование. 2012. Т.24. №1. С. 3–21.
10. Сухинов А.И., Чистяков А.Е., Тимофеева Е.Ф., Шишеня А.В. Математическая модель рас-
чета прибрежных волновых процессов // Математическое моделирование. 2012. Т.24. №8.
С. 32–44.
11. Сухинов А.И., Чистяков А.Е., Проценко Е.А. Математическое моделирование транспорта
наносов в прибрежной зоне мелководных водоемов // Математическое моделирование.
2013. Т. 25. № 12. С. 65–82.
706
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
12. Сухинов А.И., Чистяков А.Е., Проценко Е.А. Математическое моделирование транспорта
наносов в прибрежных водных системах на многопроцессорной вычислительной системе //
Вычислительные методы и программирование: новые вычислительные технологии. 2014. Т.
15. № 4. С. 610–620.
707
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
Mathematical modeling of acoustic antenna radiation on multi-
processor system
A.I. Sukhinov1, AE Chistyakov2, O.A. Savitski3,A.A. Semenyakina2, A.V. Nikitina2
Don state technical University1, Kalyaev Scientific Research Institute of Multiprocessor
Computer Systems at Southern Federal University2, N.N. Andreyev Acoustic Institute of the
Russian Academy of Sciences3
At present the development of optimal constructions of acoustic antennas is a relevant
problem. For this it is necessary to calculate characteristics of antenna arrays. A discrete
mathematical model of distribution of acoustic waves by the antennas array was designed
in this paper. Grids that takes into account the occupancy of computational cells were used
to improve the actual accuracy of solutions. A discrete model is based on the balance meth-
od. Results of numerical experiments for different forms of radiated antennas that operated
in a given frequency range were given. Grids with dimensions from 105 nodes for two-
dimensional models to 109 nodes in three-dimensional models are required for numerical
solution of these problems with the accuracy, required in designing tasks of acoustic anten-
nas in view of edge effects. The need of multiple calculations and the limit of computing
resources and, in particular, the memory of serial computers lead to design effective paral-
lel computing methods for solving these problems on available multiprocessor systems,
contained hundreds or thousands cores and available to design organizations.
Keywords: wave problem, acoustic waves, MPI.
References
1. Markov G.T., Chaplin A.F. Vozbuzhdeniye elektromagnitnykh voln [Excitation of electromagnet-
ic waves]. Moscow, Radio and Communications, 1983.
2. Vladimirov V.S. Uravneniya matematicheskoy fiziki. Uchebnik dlja fizich. i mehaniko-matemat.
spec. vuzov. [Equations of mathematical physics. Textbook for phys. and mech. and math. spec.
universities. 4th ed. and ext.]. Moscow, Science Publishing, 1981.
3. Sukhinov A.I., Zuev V.N., Semenistyy V.V. Uravneniya matematicheskoy fiziki [Equations of
mathematical physics].Taganrog, Publishing of the TRTU, 2005.
4. Samarskiy A.A. Teoriya raznostnykh skhem [The theory of difference schemes]. Moscow, Sci-
ence Publishing, 1989.
5. Sukhinov, A.I. , Chistyakov, A.E., Shishenya, A.V. Error estimate for diffusion equations solved
by schemes with weights. Mathematical Models and Computer Simulations. 1 May 2014. Vol. 6,
No. 3. P. 324–331.
6. Sukhinov A.I. , Chistyakov A.E. , Fomenko N.A. Metodika postroyeniya raznostnykh skhem dlya
zadachi diffuzii-konvektsii-reaktsii, uchityvayushchikh stepen' zapolnennosti kontrol'nykh
yacheyek [The method of constructing difference schemes for problems of diffusion-convection-
reaction, taking into account the degree of filling of the control cells]. Izvestiya Yuzhnogo feder-
al'nogo universiteta . Tekhnicheskiye nauki [Izvestiya of SFEDU. Engineering sciences.]. 2013.
No. 4. P. 87–98 .
7. Samarski A.A., Nikolayev E.S. Metody resheniya setochnykh uravneniy [Methods for solving
grid equations]. Moscow, Science Publishing.1978. 588 P .
8. Konovalov A.N. Metod skoreyshego spuska s adaptivnym poperemenno - treugol'nym pereo-
buslovlivatelem [The method of steepest descent adaptive alternating triangular preconditioner].
Differentsial'nyye uravneniya [Differential Equations]. 2004. Vol. 40, No. 7. 953 P.
708
Параллельные вычислительные технологии (ПаВТ’2016) || Parallel computational technologies (PCT’2016)
agora.guru.ru/pavt
9. Sukhinov, A.I., Chistyakov, A.E. Adaptive modified alternating triangular iterative method for
solving grid equations with a non-self-adjoint operator. Mathematical Models and Computer Sim-
ulations. 1 July 2012.Vol. 4, No. 4 P. 398–409.
10. Sukhinov A.I. , Chistyakov A.E., Timofeeva E.F., Shishenya A.V. Mathematical model for calcu-
lating coastal wave processes. Mathematical Models and Computer Simulations. 1 April 2013.
Vol. 5. No. 2, P. 122–129.
11. Sukhinov A.I., Chistyakov A.E., Protsenko E.A. Mathematical modeling of sediment transport in
the coastal zone of shallow reservoirs. Mathematical Models and Computer Simulations. 1 July
2014. Vol. 6, No. 4. P. 351–363.
12. Sukhinov A.I., Chistyakov A.E., Protsenko E.A. Matematicheskoye modelirovaniye transporta
nanosov v pribrezhnykh vodnykh sistemakh na mnogoprotsessornoy vychislitel'noy sisteme
[Mathematical modeling of sediment transport in coastal water systems on multiprocessor com-
puter system]. Vychislitel’nyye metody i programmirovaniye [Numerical Methods and Program-
ming]. 2014. Vol. 15. P. 610–620.
709