Параллельный алгоритм для задачи математического моделирования взаимодействия сильной ударной волны с молекулярным облаком*

Abstract

Results of numerical simulation of (SCI) interaction between strong (SW) shock wave generated after supernova explosion and molecular cloud (MC) have been presented in threedimension gas-dynamical arrangement. Mach number of hypersonic SW is accepted by seven. Statement of problem: the compressed supernova remnants gas with plane frontal shock wave with hypersonic velocity runs into spherical cloud and interacts with more dense gas of the cloud. Developed computational code is based on numerical solution of the Euler equations for compressible flow. Gas flow equations were represented conservatively for velocities components, density and energy. The parallel in-house code using OpenMP for PC hybrid system was developed for modeling. The Intel Vtune Amplifier XE was used to profile the code using GPU. More than a quarter of billion cells mesh (1024x512x512) with effective resolution for cloud radius in 124 nodes was implemented. Peculiarities of molecular cloud forming and density fragmentation in time were analyzed during timeconsuming simulation. Shock wave strikes the cloud, triggering weak reflected shock waves and transmitted shock that advances in to the cloud. Inner-cloud SW compresses and defragment MC media. Post-processing used allow to find out the circumstances of vortex transfer in MC, ablation and erosion phenomena.

-

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

В работе [ 6 ] используется однофазная модель газа, которая достаточно хорошо отражает физическую картину течения. В расчетах предполагалось, что невозмущенная субстанция МЗС состоит из сравнительно теплого вещества (~10000 K) и небольших неоднородно распределенных холодных облаков (~100 K), которые обладают большой плотностью [ 3-5 ]. Изначально облака находятся в динамическом равновесии с фоновым газом [ 6 ].

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

Моделирование динамических процессов в трехмерной постановке на достаточно подробных сетках требует больших вычислительных ресурсов. Для проведения расчетов был разработан параллельный алгоритм решения задач нестационарной газодинамики [ 9 ], детально опробованный при решении задачи взаимодействия набегающей ударной волны с областью меньшей плотности [ 10 ]. Для ускорения расчетов была использована технология OpenMP.

Постановка задачи. Ударная волна, образовавшаяся после взрыва сверхновой звезды, набегает на молекулярное облако, которое находится в межзвездной среде. Используется однофазная модель газа, которая обеспечивает достаточно детальную физическую картину. Изначально облако находятся в динамическом равновесии с фоновым газом. Теплопроводностью и радиационными потерями при взаимодействии будем пренебрегать. В работе используется уравнение состояния идеального газа с γ=5/3.

Рис. 1. Схема постановки задачи.

Слева – ударная волна, порожденная взрывом сверхновой звезды, справа – молекулярное облако. На заднем плане рисунка приводится изображение инфракрасного свечения пыли ароматических углеводородов вокруг газового пузыря в туманности RCW 120 (Photo Ian O’Neill from “Top 10 Spitzer Nebula Stunners”, http://news.discovery.com/space/astronomy/top-10-spitzer-anniversary-10-years-nebulae-pictures130826.htm).

Плотность межзвездной среды принимается равной ρa=2.15*10-25 г/см3, температура Ta=104 K, ua=0.0. Плотность в облаке ρc=1.075*10-22 г/см3, температура Tc=100 K, uc=0.0. Параметры газа за ударной волной определяются из уравнений Ренкина – Гюгонио. В наших расчетах число Маха ударной волны равно М=7, плотность ρsh=8.6*10-22 г/см3, температура Тsh=1.5*105 K, скорость ush=1.04*107 см/с [ 3-5 ]. Радиус облака равен rc=0.05 парсека, время прохождения ударной волной расстояния диаметра облака равно 960 лет. Толщина фронта ударной волны достаточно большая ~ 2-5 парсеков, что намного больше радиуса облака.

При заданных начальных условиях начинается расчет процесса взаимодействия высокоинтенсивной ударной волны с молекулярным облаком. В начальный момент времени ударная волна, примыкает к левой границе облака, и в момент времени t=0 начинает взаимодействовать с облаком. 2. Численный метод решения

Для решения задачи был создан параллельный алгоритм, использующий разностные схемы высокого разрешения, основанный на схемах типа ТВД [ 11 ]. Разностная схема имеет второй порядок точности и позволяет вести расчет с высоким разрешением движения ударных волн и контактных разрывов и предотвращать нефизические осцилляции. Для верификации программного кода было проведено большое количество различных двух и трехмерных тестов [ 9,12 ]. Проведенное тестирование показало хорошее совпадение результатов расчетов с аналитическими решениями и экспериментальными данными. Фронты ударных волн размазывались на величину 3-4 ячейки, контактные разрывы на 3-5 ячеек.

Представленные результаты даются в безразмерном виде. Задача решается на эйлеровой сетке. Размеры ячеек по x, y, z в приведенных расчетах одинаковы. Возникающее трехмерное движение описывается системой уравнений Эйлера, в которой отражены законы сохранения массы, количества движения и энергии. Используется уравнение состояния идеального газа. Более подробно математическая постановка задачи приведена в [ 9,10,12 ].

Для хорошего разрешения процессов, происходящих в МО необходимо, чтобы радиус облака был не менее 64 расчетных ячеек. Размер разностной сетки был выбран из этих соображений, основные расчеты проводились на сетке размерами 1024х512х512. Меньшее количество ячеек отрицательно сказывается на разрешении происходящих процессов. 2.1 Параллельный алгоритм

Математическое моделирование трехмерных задач нестационарной газовой динамики на достаточно больших сетках требуют больших вычислительных ресурсов. Для сокращения времени расчетов был построен параллельный алгоритм и программа, которая позволяет проводить вычисления с помощью OpenMP. Настройка производительности параллельного алгоритма проводилась с помощью Intel VTune Amplifier XE. Он позволяет выполнять профилирование приложений непосредственно на узле кластера. Был использован тип анализа Lightweight Hotspots. С помощью этого средства были определены функции, на выполнение которых тратится больше всего времени. Технология Lightweight Hotspots основана на использовании специальных аппаратных регистров (Performance Monitoring Units) [ 13 ], которые предназначены для учета различных низкоуровневых событий. В процессорах Intel Xeon, которые были использованы в расчетах, на ядро приходится 2 регистра, накапливающих информацию о событиях, специфичных для потока или ядра. Результат анализа приведен на Рис. 2. Основное вычислительное время занимают три подпрограммы, предназначенные для вычисления потоков на границах ячеек. Для оптимизации вычислений было использовано разбиение данных по пространству (Dimensional Splitting). Это позволило оптимизировать распределение числа потоков и памяти. На Рис. 2 изображено «качество» построенного параллельного алгоритма. Более 80% общего времени занимает расчет трех подпрограмм. После оптимизации этих подпрограмм их степень распараллеливания существенно возросла. Рис. 2. Профилирование программы с помощью Intel VTune Amplifier XE.

Зеленым цветом на рисунке показано «идеальное» распараллеливание, желтым – “хорошее”, красным – “плохое”. Плохое распараллеливание означает, что эффективность в этом случае меньше 50%. «Идеальное» распараллеливание означает, что эффективность больше 90%. 2.2 Результаты моделирования

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

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

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

На Рис. 3. приведен численный шлирен в меридиональном сечении расчетной области, показано распределение градиента плотности газа на момент безразмерного времени t=50. В левом нижнем углу приведена изоповерхность модуля ротации скоростного поля в этот момент времени, после прохождения ударной волной облака. Цветная карта на изоповерхности отображает распределение модуля местной скорости течения.

После прохождения ударной волной первоначально сферического облака профиль ударной волны (УВС) становится вогнутым, на передней кромке видны тройные ножки Маха. Падающая ударная волна распадается на проходящую и отраженную ударную волны. На поздней стадии, после прохождения ударной волной всего облака, на оси тыльной поверхности облака, в области отраженной ударной волны, формируется область высокого давления. Эта область высокого давления образует ударную волну, которая начинает распространяться навстречу падающей волне, что приводит, в итоге, к фрагментации облака по плотности. Рис. 3. Цифровой шлирен для градиента плотности при t=50 и изображение одной из изоповерхностей модуля ротации скоростного поля, с приведенной на ней картой распределения локальной скорости потока газа. Рис. 4. Эволюция во времени пространственного распределения плотности облака с картами распределения логарифма энергии на выделенных изоповерхностях логарифма плотности. Рис. 5. Возникновение поверхностных возмущений. Для двух моментов времени представлены изоповерхности логарифма плотности log ρ = 3 (красный цвет), и log ρ = 2 – (зеленый). Показаны отдельные линии тока у граничной поверхности. Возникающая неустойчивость течения в приповерхностной области ведет к возникновению вихревых структур, присущих переходным режимам, к возникновению поперечных и продольных вихревых валов и петлевым образованиям. Это иллюстрируется на Рис. 6, где дается визуализация изоповерхностей для второго инварианта тензора скорости течения – Q критерия [ 14 ].

Рис. 6. Эволюция по времени Q – критерия, приводится изображение изоповерхности Q=10 Цветовая палитра (от бирюзового к синему цвету) для изоповерхностей отображает значение логарифма плотности облака от 2 до 0.2. Для решения задачи моделирования взаимодействия сильной ударной волны с молекулярным облаком был предложен эффективный параллельный алгоритм, который позволил провести детальные газодинамические расчеты в приемлемое время. Исследование работы алгоритма, проведенное с использованием “Intel VTune Amplifier XE”, показало, что порядка 80 процентов времени расчета занимает работа всего трех подпрограмм. Была проведена оптимизация, позволившая сократить время вычислений в 4.8 раза.

Анализ результатов позволил сделать заключение, что в процессе взаимодействия волнового фронта мощной ударной волны с молекулярным облаком образуются интенсивные вихревые структуры, которые способствуют его фрагментации по плотности. Полученные результаты качественно подтверждают оптические наблюдения в астрофизике, полученные с использованием современных космических телескопов. Проведенное моделирование с выявлением вихревых структур фрагментируемых молекулярных облаков расширило и хорошо дополняет результаты других исследователей в этой области. Литература Parallel algorithm for mathematical modeling of interaction of a strong shock wave with a molecular cloud Boris Rybakin, Valery Goryachev, Lyuben Stamov and Elena Mikhalchenko Keywords: parallel computing, shock waves, ablation, fragmentation MC Results of numerical simulation of (SCI) interaction between strong (SW) shock wave generated after supernova explosion and molecular cloud (MC) have been presented in threedimension gas-dynamical arrangement. Mach number of hypersonic SW is accepted by seven. Statement of problem: the compressed supernova remnants gas with plane frontal shock wave with hypersonic velocity runs into spherical cloud and interacts with more dense gas of the cloud. Developed computational code is based on numerical solution of the Euler equations for compressible flow. Gas flow equations were represented conservatively for velocities components, density and energy. The parallel in-house code using OpenMP for PC hybrid system was developed for modeling. The Intel Vtune Amplifier XE was used to profile the code using GPU. More than a quarter of billion cells mesh (1024x512x512) with effective resolution for cloud radius in 124 nodes was implemented. Peculiarities of molecular cloud forming and density fragmentation in time were analyzed during time-consuming simulation. Shock wave strikes the cloud, triggering weak reflected shock waves and transmitted shock that advances in to the cloud. Inner-cloud SW compresses and defragment MC media. Postprocessing used allow to find out the circumstances of vortex transfer in MC, ablation and erosion phenomena.

References

1. Hartquist , T.W. , Dyson , J.E. , Pettini , M. , Smith , L.J. 1986 , MNRAS, 221 , 715

2. Hartquist , T.W. , Dyson , J.E. 1988 , Ap &SS, 144 , 615

3. McKee C.F. Supernova remnant shocks in an inhomogeneous interstellar medium // Supernova Remnants and the Interstellar Medium , 1988 . IAU Coll. 101 . pp. 205 - 222 .

4. Klein R.I. , McKee C.F. , Colella P . On the hydrodynamics interaction of shock waves with interstellar clouds . Nonradiative shocks in small clouds // The Astrophysical Journal , 1994 . Vol. 420 . pp. 213 - 236 .

5. Melioli C. , de Gouveia Dal Pino E., Raga A . Multidimensional hydro dynamical simulations of radiative cooling SNRs-clouds interactions: an application to starburst environments // Astronomy & Astrophysics, 2005 . Vol. 443 . pp. 495 - 508 .

6. Frank A. , Poludnenko A. , Gardiner T. , Lebedev V. , Drake R ( 2003 ), Stellar outflows with new tools: advanced simulations and laboratory experiments . RevMexAA (Serie de Conferencias), 15 , pp. 85 - 91 .

7. Stone J.M. , Norman M.L. The three-dimensional interaction of a supernova remnant with an interstellar cloud / / APJL, 1992 . Vol. 390 . pp. 17 - 19 .

8. J. Kelly Truelove , Richard I. Klein , Christopher F. Mckee , John H. Holliman , Louish H. Howell , Jeffrey A. Greenough , and D. Tod Woods . Self-Gravitational Hydrodynamics with ThreeDimensional Adaptive Mesh Refinement: Methodology and Applications to Molecular Cloud Collapse and Fragmentation, 1997 .

9. Rybakin B . Modeling of III-D Problems of Gas Dynamics on Multiprocessing Computers and GPU // Computers and Fluids, 2013 . Vol. 80 . pp. 403 - 407 .

10. Rybakin B. , Goryachev V. The supersonic shock wave interaction with low-density gas bubble / / Acta Astronautica, 2014 . Vol. 94 , Issue 2. pp. 749 - 753 .

11. Harten A. High resolution schemes for hyperbolic conservation laws // J. Comp. Phys., 1983 . Vol. 49 , pp. 357 - 393 .

12. Rybakin B.P. , Stamov L.I. , Egorova E.V. Accelerated solution of problems of combustion gas dynamics on GPUs // Computers & Fluids, 2014 . Vol. 90 , pp. 164 - 171 .

13. Intel Corporation . Advanced Intel Xeon Phi Coprocessor Workshop , Performance Tuning for Intel Xeon Phi Coprocessors.

14. Hunt J. , Wray A. , Moin P. Eddies , streams, and convergence zones in turbulent flows // Proc. Summer Program Center for Turbulence Research (NASA Ames/Stanford Univ.) , 1988 , pp. 193 - 208 .