Моделирование влияния волнения поверхности на спектр собственного излучения океана Текст научной статьи по специальности «Математика»

ISSN 2079-3316 ПРОГРАММНЫЕ СИСТЕМЫ: ТЕОРИЯ И ПРИЛОЖЕНИЯ №2(29), 2016, с. 73-84 УДК 537.874.6,004.942

А. А. Гавриков, Д. Ю. Князьков, А. В. Романова, В. В. Черник,

А. С. Шамаев

Моделирование влияния волнения поверхности на спектр собственного излучения океана

Аннотация. Для нужд радиотомографии поверхности Земли рассматривается задача моделирования излучения морской поверхности. Метод параллельных сечений моделирования дифракции электромагнитного излучения на слое с периодической поверхностью был реализован в виде компьютерной программы на языке программирования python с использованием библиотеки параллельной коммуникации pypar. С помощью этой программы был проведен ряд численных экспериментов. Для расчетов использовался гибридный кластерный вычислитель HybriLIT, находящийся в лаборатории информационных технологий ОИЯИ в г. Дубне.

Ключевые слова и фразы: HPC-вычисления, суперкомпьютерные вычисления, радиометрия океана, дифракция, метод параллельных сечений, дифракция на периодическом слое.

Введение

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

Задача радиотомографии морской поверхности предполагает нахождение формы поверхности, например, коэффициентов |а^}, г = 1,... ,п, в ее Фурье-разложении; по измеренному в заданном направлении набору интенсивностей }, ] = 1,...,то; ее собственного

Работа выполнена при финансовой поддержке НТЦ «Космонит» ОАО «Российские космические системы».

© А. А. Гавриков, Д. Ю. Князьков, А. В. Романова, В. В. Черник, А. С. Шамаев, 2016 © Институт проблем механики им. А.Ю. Ишлинского РАН, 2016 © Программные системы: теория и приложения, 2016

излучения на разных частотах }, ^ = 1,... ,т:

(1) М&)} ^ {*}.

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

(2) {о,и,-,Р ^ ^(Р).

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

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

1. Метод параллельных сечений

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

(3)

(4)

(5)

(6)

ш -

ТоШ = -'1е — Е, с

то^ = -Н,

с

^(еЕ1) = 0, &у(#) = 0,

где е = £о + а — комплексная проницаемость среды, а — проводимость, ш — частота электромагнитной волны, с — скорость света, £о(х) — диэлектрическая проницаемость.

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

(7) К = , ^ \Ъ< = Вр

с краевыми условиями

(8) р + гК Ь = с, при х = хт1П = 0,

(9) р - гКЬ = 0, при х = хтах,

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

При этом, так как диэлектрическая проницаемость £(х) непостоянна (терпит разрыв при переходе через поверхность раздела воды и воздуха), отдельно рассматриваются случаи различной поляризации электромагнитных волн. В случае ^-поляризации векторы Е и Н имеют вид

Е =(ЕХ,ЕУ, 0), Н =(0,0,Нг),

а в случае Деполяризации

Е =(0,0,Ег), Н = (НХ,НУ, 0).

Для решения системы (7) сначала решаются вспомогательные задачи вида

у1'? )' = АЬ1'?,

((р1'!) {(Ъ1'*)

с начальными условиями

Ь1'^ = ¿2 , р1'^ = 0 при X = 0

(11)

с начальными условиями

(р2'!)' = АЬ2'3 (Ъ2'1)' = Вр2'!

Ь2'^ = 0, р2'^ = ¿1 при х = 0.

Здесь

р3'3

8'3 8'3

Ры , $N2-1,

= / }?'3 г.8'3 г.8'3

° = 2 ,"N2-1 ,...,°~ы1-

-т,...,м2

з = -N1, ...,N2

где = 1, если г = 3 и ¿¿^ = 0, если г = 3. (12)

Введем равномерную сетку на отрезке [0,хтаж] с шагом к = ^ы^:

%тах

0,

N

/А Т -¡\хтах

N ''

и будем методом прогонки на этой сетке решать системы (10) и (11).

Найдем решение исходной задачи (7), (8), (9) в виде линейной комбинации решений задач (10) и (11) с набором коэффициентов а. Пользуясь первым граничным условием (8) задачи и подставляя эту линейную комбинацию во второе краевое условие (9), получим систему линейных алгебраических уравнений на коэффициенты линейной комбинации а:

(13)

Са = В.

Таким образом, определены основные этапы алгоритма решения исходной задачи дифракции электромагнитной волны на периодическом цилиндрическом слое:

• решение систем ОДУ (10), (11) методом прогонки на сетке (12) — всего 2(М2 + М1 + 1) систем ОДУ;

• решение СЛАУ (13).

и

3

Описанный выше подход для H-поляризации ранее был реализован в компьютерной программе на языке программирования python [4]. В настоящей работе эта программа была доработана для осуществления расчетов дифракции Е-поляризованных волн, было проведено распараллеливание расчётов для использования суперкомпьютера кластерного типа.

2. Организация параллельных вычислений

Для проведения численных расчётов использовался вычислительный кластер HybriLIT, установленный в Лаборатории информационных технологий ОИЯИ, г. Дубна (http://hybrilit.jinr.ru/). На данный момент он состоит из 7 вычислительных узлов, с двумя двенадцати-ядерными процессорами Intel Xeon E5-2695 v2 2.4 ГГц на каждом узле (всего 168 вычислительных ядер, 3.2 ТФлопс пиковой производительности) и имеет 896 Гб ОЗУ. Также на кластере установлены одна видеокарта NVIDIA Tesla K20X, 12 видеокарт NVIDIA Tesla K40 (всего 37248 CUDA ядер, 18.47 ТФлопс на числах двойной точности), один сопроцессор Intel Xeon Phi 5110P и два сопроцессора Intel Xeon Phi 7120P (всего 182 ядра, 3 ТФлопс на числах двойной точности).

2.1. Время счета однопроцессорного варианта программы

Исследуем зависимость количества времени, требуемого на один расчёт дефекта, от размера расчётной сетки. Эта зависимость приведена в таблице 1. Получается, что для сетки размера 400 время счета на одном вычислительном ядре будет составлять около 2.5 минут, что вполне приемлемо при единичных расчётах, но может превращаться в часы или даже десятки часов при массовых расчётах, когда требуется делать много элементарных расчётов для различных сочетаний параметров, принимающих с некоторым шагом все значения из определённых диапазонов изменения. Более того, оказалось, что для достоверного счета необходимо брать период расчёта в 5-10 раз больший, чем период изменения функции, задающей поверхность. Это требует уже 1-2 часов на проведение единичного численного расчёта на одном вычислительном узле.

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

Таблица 1. Зависимость времени счёта от размера расчётной сетки на 2 вычислительных узлах

Размер расчётной сетки 100 200 300 400 500 600 Среднее время счёта, сек:

1 задачи на 1 ядре 11 43 95 167 260 375

48 задач на 48 ядрах 12 47 105 187 288 415

(формы поверхности, длины волны, угла падения), не имеет смысла осуществлять распараллеливание алгоритма самого проекционного метода (точнее самой затратной по времени его части — решения системы дифференциальных уравнений (7)), а следует организовывать параллельное, одновременное осуществление элементарных расчётов дефекта при фиксированных длине волны, форме поверхности, периоде, угле падения.

2.2. Организация расчетов на кластере HybriLIT ЛИТ ОИЯИ

Программа расчёта по проекционному методу была написана на языке python, для распараллеливания расчётов использовался python-модуль pypar (https://code.google.eom/p/pypar/). В программах, написанных на языке python, этот модуль позволяет пользоваться средствами параллельного интерфейса MPI (Message passing interface — интерфейс передачи сообщений). Такой модуль был установлен сотрудниками ЛИТ ОИЯИ на вычислительном кластере HybriLIT. В программе общая расчётная задача разбивалась на подмножества элементарных расчётов (нахождение дефекта при заданных угле падения, длине волны, периоде, форме поверхности), и каждое из этих подмножеств должно было быть обсчитано на отдельном вычислительном ядре.

Морская поверхность задавалась функцией f (х) следующего вида:

(14) f (х) = I sin(x) + l1 — sin(k1х) + l2k2 sin ( — ) .

«1

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

е = {N, М, Ь, X,a,£,l,l1,k1,l2,k2},

где N, М — количество расчетных точек по осям х и z, Ь — период расчета, Л — длина волны излучения, а — угол падения, е — комплексная диэлектрическая проницаемость морской воды, ki G N, амплитудные коэффициенты I, li > 0.

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

{'defect' :< float >,' passed! :< float >,' reflected' :< float >},

где в поле 'defect' выводится дефект, а в полях 'passed' и 'reflected' доли прошедшего и отраженного излучения соответственно. (Выбиралась такая толщина слоя воды, чтобы доля прошедшего излучения была практически нулевой, например, менее 10-8.) Затем, на локальном компьютере, по этому файлу строились графики с использованием python-библиотеки matplotlib (http://matplotlib.org/).

2.3. Эффективность параллельных вычислений

При описанной выше схеме счёта, ввиду независимости отдельных расчетов, должно получиться линейное уменьшение времени счёта при увеличении количества используемых вычислительных ядер. В ходе работы с параллельной версией программы на кластере HybriLIT было установлено, что эффективность вычислений действительно практически равна 1, хотя немого падает при увеличении количества используемых на узле ядер от 1 до 24, а при одновременном запуске 25 процессов время работы сразу удваивается (по сравнению с одновременным запуском 24 процессов). Среднее время одинарного расчёта практически не изменяется, когда используются более одного вычислительного узла (см. нижнюю строку таблицы 1). Таким образом, следует использовать количество MPI-процессов в точности равное суммарному количеству физических ядер процессоров, установленных на каждом узле (в нашем случае — 24 вычислительных процесса на узел).

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

3. Результаты численных экспериментов

С помощью описанной выше программы был произведён ряд численных расчетов [5]. Исследовалась зависимость интенсивности излучения от угла визирования, влияние на излучение возмущения волны мелкой рябью и длинными волнами, зависимость излучения морской поверхности от периода волнения.

Было установлено, что при угле визирования в 30° (резонансном для заданных периода волнения и длины волны электромагнитного излучения) наблюдается резкий скачок интенсивности излучения, далее до угла визирования 60° интенсивность излучения монотонно зависит от интенсивности волнения при фиксированном угле визирования, однако при дальнейшем увеличении угла эта монотонность теряется. Таким образом, было определено, что пригодными углами визирования являются углы от 30° до 60°.

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

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

• в целом селективность при использовании Е-поляризации оказалась выше;

• если падение интенсивности для Н-поляризации с уменьшением периода от основного значения идет медленно, то в случае Е-поляризации оно идет скачками;

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

4. Заключение

Был реализован кластерный суперкомпьютерный вариант организации расчёта по проекционному методу. Эта параллельная реализация использовалась на кластере HybriLIT ЛИТ ОИЯИ, позволяя получать линейное1 ускорение счёта и существенно снижать время счета при проведении численных экспериментов.

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

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

У авторского коллектива есть положительный опыт переноса расчетных алгоритмов дифракции на архитектуру CUDA, когда было получено почти 5-ти кратное ускорение по сравнению с CPU [6] (сравнивался случай задействования всех ядер процессора Intel Xeon 5570 и всех ядер карты Tesla C1060 на двойной точности). Планируется осуществить адаптацию описанных в настоящей работе численных алгоритмов к архитектуре карт Tesla и сопроцессоров Xeon Phi и написать соответствующие программы. Параллельная эффективность отдельного расчета, включающего в себя решение систем ОДУ и СЛАУ, на SIMD архитектуре не столь очевидна, хотя есть основания предполагать, что решение независимых вспомогательных задач (10) и (11) одновременно на разных ядрах должно дать хорошее ускорение счета. В любом случае, от использования карт CUDA и ускорителей Xeon Phi ожидается прирост производительности из-за возможности одновременного проведения независимых отдельных расчетов при увеличении общей вычислительной мощности. Такая адаптация

1по количеству используемых вычислительных ресурсов

представляется целесообразной, поскольку большинство современных суперкомпьютеров (в том числе, доступные авторам для счета вычислители МСЦ РАН и ЛИТ ОИЯИ) имеют сейчас гибридную архитектуру.

Авторы выражают глубокую признательность И.В. Черному за постановку задачи и полезные обсуждения, а также руководству и сотрудникам ЛИТ ОИЯИ, любезно предоставившим возможность и техническую поддержку расчетов на кластере HybriLIT.

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

[1] С. В. Нестеров, А. С. Шамаев, С. И. Шамаев. Методы, алгоритмы и средства аэрокосмической компьютерной томографии приповерхностного слоя Земли, Научный мир, М., 1996, 272 с. t 73

[2] М. Л. Левин, С. М. Рытов. Теория равновесных тепловых флуктуаций в электродинамике, Наука, М., 1967, 310 с. t 74

[3] А. С. Ильинский. «Метод исследования задач дифракции волн на периодической структуре», Ж. вычисл. матем. и матем. физ., 14:4 (1974), с. 1063-1067. t 75

[4] В. В. Черник. «Применение методов декомпозиции и интегральных преобразований для решения задачи прохождения плоской волны через неоднородную среду», Аэрофизика и космические исследования, Труды 57-й научной конференции МФТИ, МФТИ, М., 2014, с. 23-24. t 77

[5] Д. Ю. Князьков, В. В. Черник. «Моделирование процесса электромагнитного зондирования морской поверхности», Аэрофизика и космические исследования, Труды 58-й научной конференции МФТИ, ред. С. С. Негодяев, МФТИ, М., 2015, с. 37-38. t 80

[6] D. Knyazkov. "Simulation of Holography Using Multiprocessor Systems", Mathematical Modeling and Computational Science, MMCP2011, Lecture Notes in Computer Science, vol. 7125, Springer-Verlag, Berlin-Heidelberg, 2012, pp. 270-275. t 81

Рекомендовал к публикации Программный комитет

Четвёртого национального суперкомпьютерного форума НСКФ-2015

Пример ссылки на эту публикацию:

А. А. Гавриков, Д. Ю. Князьков, А. В. Романова, В. В. Черник, А. С. Шамаев. «Моделирование влияния волнения поверхности на спектр собственного излучения океана», Программные системы: теория и приложения, 2016, 7:2(29), с. 73-84.

URL: http://psta.psiras .ru/read/psta2016_2_73-84.pdf

Об авторах:

Александр Александрович Гавриков К.ф.-м.н., с.н.с. ИПМех РАН. Область научных интересов: теория колебаний, собственные колебания упругих систем, акустика неоднородных сред, теория усреднения. e-mail: topol555@gmail.com

Дмитрий Юрьевич Князьков К.ф.-м.н., н.с. ИПМех РАН. Область научных интересов: прямые и обратные задачи дифракции, математическая теория усреднений, суперкомпьютерные вычисления, алгоритмы быстрой свертки и БПФ, радиометрия океана, сейши. Автор более 20 научных работ и 4 патентов.

e-mail: knyaz@ipmnet.ru

Алла Владимировна Романова

Сотрудник ИПМех РАН, аспирант МГУ им. М. В. Ломоносова. e-mail: avromm1@gmail.com

Виталий Валерьевич Черник

Специалист, выпускник кафедры дифференциальных уравнений механико-математического факультета МГУ им. М. В. Ломоносова, младший научный сотрудник лаборатории механики управляемых систем ИПМех РАН. Область научных интересов: управление в механических системах, математическое программирование, математическое моделирование волновых процессов.

e-mail: gungho424@gmail.com

Алексей Станиславович Ш^амаев Д.ф.-м.н., г.н.с. ИПМех РАН, зам. зав.каф. дифференциальных уравнений механико-математического факультета МГУ им. М. В. Ломоносова. Область научных интересов: управление в механических системах, управление стохастическими динамическими системами, дифференциальные уравнения математической физики, асимптотические методы для дифференциальных уравнений, теория усреднения. e-mail: sham@rambler.ru

Alexandr Gavrikov, Dmitriy Knyazkov, Alla Romanova, Vitaliy Chernik, Alexey Shamaev. Simulation of influence of the surface disturbance on the ocean self radiation spectrum.

Abstract. A problem of sea surface radiation simulation is considered for the needs of Earth surface radiotomography. The method of parallel sections for simulation of diffraction of electromagnetic radiation on a periodical surface layer was implemented as a computer program written on the python programming language with the use of pypar parallel communication library. A number of numerical experiments was performed with the use of this program. A hybrid computational cluster HybriLIT (The Laboratory of Information Technologies at JINR, Dubna, Russia) was used for the computations. (In Russian).

Key words and phrases: HPC-computations, ocean radiometry, diffraction, parallel sections method, diffraction on a periodical layer.

References

[1] S. V. Nesterov, A.S. Shamayev, S.I. Shamayev. Methods, algorithms and tools of aerospace computer tomography of near-surface layer of the Earth, Nauchnyy mir, M., 1996 (in Russian), 272 p.

[2] M. L. Levin, S. M. Rytov. Theory of equilibrium thermal fluctuations in electrodynamics, Nauka, M., 1967 (in Russian), 310 p.

[3] A. S. Il'inskiy. "A method of investigating wave diffraction problems on a periodic structure", USSR Computational Mathematics and Mathematical Physics, 14:4 (1974), pp. 242-246.

[4] V. V. Chernik. "The use of decomposition and integral transformations methods to solve the problem of passing a plane wave through the inhomogeneous medium", Aerophysics and Space Researches, Proc. of the 57th MIPT Scientific Conference, MFTI, M., 2014, pp. 23-24 (in Russian).

[5] D. Yu. Knyaz'kov, V. V. Chernik. "The simulation of electromagnetic sensing of sea surface", Aerophysics and Space Researches, Proc. of the 58th MIPT Scientific Conference, ed. S. S. Negodyayev, MFTI, M., 2015, pp. 37-38 (in Russian).

[6] D. Knyazkov. "Simulation of Holography Using Multiprocessor Systems", Mathematical Modeling and Computational ¡Science, MMCP2011, Lecture Notes in Computer Science, vol. 7125, Springer-Verlag, Berlin-Heidelberg, 2012, pp. 270-275.

Sample citation of this publication:

Alexandr Gavrikov, Dmitriy Knyazkov, Alla Romanova, Vitaliy Chernik, Alexey Shamaev. "Simulation of influence of the surface disturbance on the ocean self radiation spectrum", Program systems: theory and applications, 2016, 7:2(29), pp. 73-84. (In Russian).

URL: http://psta.psiras.ru/read/psta2016_2_73- 84.pdf

© A. A. Gavrikov, D. Y. Knyazkov, A. V. Romanova, V. V. Chernik, A. S. Shamaev, 2016 © A. Ishlinsky Institute for Problems in Mechanics RAS, 2016 © Program systems: Theory and Applications, 2016