Содержание к диссертации
Введение
Глава 1. Проблема решения на МВС одномерных краевых задач для уравнении 2-го порядка 23
1. 1. Постановки краевых задач и их разностные аппроксимации 24
1.1 А. Постановки краевых задач 24
1.1.2. Разностные задачи 26
1.1.3. Система алгебраических уравнении 28
1.2. Вазовый алгоритм распараллеливания 29
1.3. Эффективность базового алгоритма при расчетах на МВС 33
Глава 2. Решение с помощью МВС многомерных краевых задач для параболических уравнений 39
2.1. Постановка модельной краевой задачи 39
2.2. Методы численного решения уравнения теплопроводности 42
23. Параллельные алгоритмы численного решения на МВС уравнения теплопроводности 46
2.3.1. Параллельные алгоритмы реализации явной схемы 46
23.2. Параллельные алгоритмы реализации неявной ЛОС 50
233. Параллельные алгоритмы реализации схемы двуцикъического расщепления 52
2.4. Результаты тестовых расчетов 53
Глава 3. Моделирование на МВС течения в недорасширенной струе 64
3.1. Применение МВС для расчета течения в недорасширенной струе с использованием явной схемы 64
3.1.1. Квазигазодинамические уравнения 64
3.1.2. КГДуравнения в r-z геометрии и постановка задачи . 66
3.1.3. Реализация алгоритма на МВС 69
3.2. Неявный метод для КГД-уравнений и его параллельная реализация 72
3.2.1. Линеаризация квазигаюдииамической системы уравнениіі 74
3.2.2. Схема Бима - Уормгтга для решения стационарных задач 78
3.2.3. Граничные условия 87
3.2.4. Способы параллельной реализации 89
3.3. Результаты расчетов 91
3.3.1. Параметры течения иа срезе сопла 91
3.3.2. Сравнение с экспериментом 93
3.3.3. Результаты тестирования многопроцессорных систем на задаче расчета течения в недорасширенной струе 98
Заключение 106
Список литературы 107
- Система алгебраических уравнении
- Параллельные алгоритмы численного решения на МВС уравнения теплопроводности
- КГДуравнения в r-z геометрии и постановка задачи
- Схема Бима - Уормгтга для решения стационарных задач
Введение к работе
Процесс йычислительной техники способствует успешному решению сложных фундаментальных и прикладных задач во многих отраслях науки, технологии и производства. Наиболее сильно эта тенденция проявляется в области высоких технологий, применяемых в аэро-коемической. автомобильной, химической и электронной промышленности, а также в энергетике и медицине. При этом важнейшим фактором эффективного использования современной высокопроизводительной многопроцессорной вычислительной техники являюіся методы математического моделирования, ориентированные на ее использование. В связи с этим одной из задач современной фундаментальной и прикладной пауки стало обеспечение высоких темпов роста в области создания нового и развития традиционного аппарата математического моделирования и информационных технологий в целом, рассматривающих МВС как основной инструмент для вычислений.
В данной работе рассматриваются некоторые аспекты развития методов математического моделирования, связанный с применением массивно-параллельных многопроцессорных вычислительных систем (МВС) с распределенной памятью для решения современных задач механики сплошной среды. В работе изучаются практически все основные этапы решения прикладной задачи (выбор математической модели, построение ее дискретного аналога, разработка численного метода и параллельного алгоритма его реализации на конкретной МВС) с целью получения максимальной отдачи от компьютерного моделирования на МВС. В качестве объекта исследования выбраны методы численного решения с помощью МВС краевых задач для обыкновенного дифференциального уравнения 2-го порядка общего вида, трехмерного уравнения теплопроводности, а также практическая
задача о моделировании струй, истекающих в область низкого давления. Основные усилия в работе сосредоточены на проблемах решения систем сеточных уравнений, возникающих при разностной аппроксимации дифференциальных уравнений. Эти проблемы традиционно считаются ключевыми в вычислительной математике, поскольку построение эффективного численного, а в данном случае еще и параллельного алгоритма решения систем сеточных уравнений в значительной мере определяет возможность успешной реализации всего алгоритма решения сложной задачи на многопроцессорной системе.
Целями диссертации являлись разработка и сравнительный анализ параллельных алгоритмов и эффективно переносимых программ, позволяющих моделировать процессы теплопереноса и динамики вязкого газа на многопроцессорных вычислительных системах с распределенной памятью, а также распространение полученных теоретических и практических результатов на решение систем уравнений механики сплошной среды общего вида.
Выбор в качестве практической задачи расчета струйных течений связан со следующим. Достигнутый в настоящее время прогресс в области вычислительной техники позволяет проводить расчеты на лодробных пространственных сетках, содержащих миллионы и далее миллиарды узлов. Классическими проблемами, для которых такие расчеты всегда играли важную роль, являются аэрокосмические исследования, задачи физики плазмы. Газодинамические расчеты также относятся к данному классу задач, и получили широкое распространение в метеорологии, автомобильной промышленности (создание аэродинамически оптимальных корпусов современных автомобилей и численное моделирование физических процессов внутри двигателей внутреннего сгорания), экологии (исследование распространения загрязнений в атмосфере) и многих других областях современных
наукоемких технологий [45,25]. Практическая необходимость постоянно повышать достоверность результатов численного моделирования в этих областях науки и техники требует проведения расчетов со все более высокой точностью, то есть на более подробных сетках. Использование же таких сеток требует во многих случаях пересмотра традиционных численных подходов и создания новых эффективных численных методов, изначально ориентированных на использования мощных МВС с большим числом процессоров. Задача о струях удачно вписывается в контекст этой проблемы- В ней прослеживается не только вся цепочка известной триады математического моделирования (модель - алгоритм -программа) [1], но и дополнительное звено - МВС. К тому же, для этой задачи имеются многочисленные экспериментальные данные, позволяющие оценить качество модельных расчетов. С другой стороны, среди полученных численных результатов обнаружились и новые детали, не доступные для измерения в натурных экспериментах.
Рассмотрим теперь несколько общих аспектов применения МВС. Многопроцессорная вычислительная техника берет свое начало с создания в 1985 году фирмой INMOS (Великобритания) микроэлектронного прибора, названного транспьютер [35,37], и предназначенного для использования в качестве вычислительного элемента в параллельных архитектурах. Это открыло новые возможности в области повышения производительности самой вычислительной техники, а также в области повышения эффективности методов вычислительной математики. ГТоявлению транспьютера предшествовала разработка фирмой TNMOS в сотрудничестве с Оксфордским университетом языка параллельного программирования высокого уровня Оккам (Occam) [38] для описания модели взаимодействующих последовательных процессов. Транспьютер разрабатывался для реализации этой модели, и в нем имелась аппаратная поддержка основных конструкций языка Оккам, что в свою очередь
позволяло выполнять программы па этом языке с эффективностью Ассемблера [39]. Для реализации численного алгоритма на транспьютерной системе необходимо было выбрать принцип распараллеливания, задать физическую конфигурацию системы, создать описание программ для отдельных процессоров, написать команды обмена информацией для каждого транспьютерного линтса [36,40.41].
В настоящее время параллельные системы создаются на базе мощных процессоров [88] с производительностью порядка 2GFLOPS (IGFLOPS^IO4 floating point instructions per second) [34] и выше. При этом физическая конфигурация системы уже не имеет определяющего значения, поскольку пользователь может писать программу на любом языке программирования, используя для обменов информацией ту или иную библиотеку программ [7,33], установленную на данной МВС.
Адаптация известных численных алгоритмов к вычислениям на МВС остается пока актуальной задачей- Основной причиной такой ситуации была быстрая смена архитектур МВС. В настоящее время системы разделяют по различным признакам: по типам потока команд и потока данных, способам обработки данных, по строению памяти и типу коммуникационной сети, степени однородности компонент системы, степени согласования режимов работы устройств и т.д. [2,32]. Попытки систематизации множества архитектур начались после опубликования М.Флинном [31] первого варианта классификации вычислительных систем в конце 60-х годов- Классификация Флинна по типу потока команд и данных делит все параллельные системы на четыре класса:
SISD (одиночный поток команд - одиночный поток данных),
SJMD (одиночный поток команд - множественный поток данных),
MISD (множественный поток команд - одиночный поток данных),
MIMD (множественный поток команд - множественный поток данных).
R последствии различными авторами уточнялось понятие M1MD-архитектуры (классификации Ванга и Бриггса, Хокни, Джонсона и др. [34J). С этих позиций алгоритмы, представленные в диссертации, можно считать ориентированными на MIMD- и SPMD (Одна программа -множественный поток данных) - архитектуры.
Главный вопрос заключается в том, какие признаки должны лежать в основе классификации МВС. Ответ на него зависит от того, для кого данная классификация создается, и на решение какой задачи она направлена. Основным параметром классификации МВС является наличие обшей (SMP) или распределенной (МРР) памяти. Нечто среднее между SMP и МРР представляют собой NUMA-архитектуры, где память физически распределена, но логически общедоступна. Кластерные системы являются более дешевым вариантом МРР. При поддержке команд обработки векторных данных говорят о векторно-конвейерных процессорах, которые, в свою очередь могут объединяться в PVP-системы с использованием общей или распределенной памяти. Все большую популярность приобретают идеи комбинирования различных архитектур в одной системе и построения неоднородных систем. Приведем особенности перечисленных архитектур:
Массивно-параллельные системы (МРР) - однородные, с локальной памятью, могут содержать до нескольких тысяч вычислительных узлов (например, МВС ASCI Red, Blue Mountain); программирование для них происходит в рамках модели передачи сообщений (MPI);
Симметричные мультипроцессорные системы (SMP) - содержат небольшое число (до 32} однородных процессоров, использующих общую память; программирование осуществляется в рамках модели общей памяти (пример - МВС HP 9000 V-cIass);
Системы с неоднородным доступом к памяти (NUMA) - состоят из однородных базовых модулей (плат), объединенных высокоскоростным
коммутатором (максимальное число процессоров - 256} и работающих в едином адресном пространстве, аппаратно поддерживающим доступ к удаленной памяти (к памяти других модулей); модель проіраммирования аналогична SMP-системам (пример - 1V1BC SGI Origm2G0G);
Векторно-конвейерные компьютеры (PVD) - системы, в процессорах которых предусмотрены команды однотипной обработки векторов независимых данных, выполняющиеся па конвейерных функциональных устройствах; имеют возможность использования общей памяти; программирование предусматривает векторизацию циклов и их распараллеливание (пример - МВС NEC SX-4/SX-5);
Гибридные (кластерные) системы - это объединение векторных или S MP-компьютеров в массивно-параллельные системы через высокоскоростную коммуникационную среду; программирование осуществляется в рамках гибридной модели*
К гибридным системам относятся, например, сети (кластеры) рабочих станций, использующие различную сетевую аппаратуру Ethernet, FastEthcrnct и т.п. При этом возникает неоднородная вычислительная среда (и по производительности процессоров, и по скорости коммуникаций, и по программному обеспечению), позволяющая проводить параллельные вычисления па базе таких стандартов коммуникаций, как PVM, MPI или подобных им. Программирование на таких системах осуществляется в рамках модели передачи сообщений.
Бурно развивающиеся в настоящее время глобальные компьютерные сети (Internet, Internet-2 и д,р.) фактически представляют собой вычислительные системы с распределенной памятью, включающие в себя сотни тысяч взаимодействующих компьютеров, и также могут рассматриваться как МВС.
Учитывая массовый характер использования МРР и гибридных вычислительных систем с распределенной памятью, в настоящей работе исследовались параллельные алгоритмы, использующие модель передачи сообщений. Эти свойства разработанных алгоритмов тем не менее не оітшничивакл их применение на системах с общей памятью, например, с помощью библиотеки коммуникаций MPL
Популярность МРР и гибридных систем в первую очередь определяется малым соотношением их стоимости и суммарной производительности, то есть они выгодны с экономической точки зрения. Кроме того, такая вычислительная техника позволяет получить практически любую производительность путем наращивания числа процессоров (то есть МВС данной архитектуры является масштабируемой в отличие от SMP систем). При этом важным фактором является гибкость настройки аппаратуры на решение конкретных задач, а также общая надежность работы узлов МВС.
Отметим далее, что за время существования МВС накоплен значительный опыт их использования, и они действительно находятся в массовом доступе (МВС имеются в ряде Российских и в большинстве зарубежных научных институтов). Однако, огромные вычислительные возможности, предоставляемые этими системами недостаточно используются. Причина тому - сложность адаптации последовательных алгоритмов и программ к параллельным архитектурам и, как следствие, общее отставание в области разработки прикладного программного обеспечения МВС по сравнению с бурно развивающимся аппаратным обеспечением и вспомогательными программными средствами. Поэтому разработка эффективных параллельных алгоритмов решения фундаментальных и прикладных задач, изначально ориентированных на использование МВС, в том числе с распределенной памятью, является чрезвычайно актуальной проблемой. При этом надо учитывать, что
вместе с возможностями, которые предоставляет программистам
многопроцессорная вычислительная техника, существуют
специфические требования к используемым ею алгоритмам [27]:
обеспечение равномерного распределения вычислительной нагрузки но процессорам;
обеспечение ускорения расчетов при увеличении числа процессоров;
необходимость синхронизации при остановке (завершении) итераций;
снижение затрат времени на синхронизацию;
периодический сбор данных для записи и анализа;
обмен данными между процессорами для продолжения расчета;
снижение объема передаваемой информации при обмене;
снижение числа актов приема-передачи (то есть числа обменов между процессорами);
уменьшение влияние размера задачи на ускорение;
обеспечение масштабируемости и переносимости программы.
Таким образом, разработка алгоритмов, позволяющих решить перечисленные проблемы, и является главной целью параллельного программирования.
Перейдем к проблемам решения с помощью МВС задач механики сплошной среды. Как известно, численное решение краевых задач для обыкновенных дифференциальных уравнений и уравнений в частных производных является одной из центральных проблем математического моделирования в данной области. Важную роль при решении систем линейных уравнений высокого порядка с матрицами специального вида, получающимися при аппроксимации дифференциальных уравнений, играет выбор эффективных методов (прямых и итерационных). При
выборе того или иного численного метода и алгоритма учитывается много обстоятельств, таких, как объем вычислений, требуемый объем оперативной памяти, порядок точности, устойчивость по отношению к погрешностям округления и другие.
При решении задач на МВС в поле зрения сразу попадают явные методы и алгоритмы. Среди них широкое распространение получили явные конечно-разностные схемы для решеттия эволюционных задач. Их основными достоинствами являются наглядность и простата при реализации сложных граничных условий. Кроме того, явные схемы обладают неоспоримыми преимуществами при реализации на МВС. Недостатком явных схем является то, что временной (или итерационный) шаг не может быть произвольным. На значение шага интегрирования по времени из соображений устойчивости накладывается ряд ограничений. Существуют работы, позволяющие для стационарных задач увеличивать шаг по времени [59,86]. Эта меюдика позволяет считать задачу по явной схеме с переменным шагом, значительно его увеличивая в определенных подобластях расчетной области. Однако данный метод пригоден лишь для определенного класса задач и в целом не решает проблему увеличения шага по времени. Расчеты по явным разностным схемам показывают почти идеальную эффективность решения задач на МВС,
Неявные методы и алгоритмы при реализации на МВС требуют больших усилий от исследователей [90]. В сложных многомерных нелинейных задачах их обычно получают при линеаризации конечно-разностных аналогов исходных дифференциальных уравнений относительно приращений искомых функций. В дальнейшем производится приближенная факторизация (расщепление) конечно-разностных операторов, приводящая к цепочкам пространственно одномерных разностных задач. Разработка схем данного типа
представляет собой развитие классических неявных схем метода переменных направлений и метода дробных шагов [60]. Схемы могут строиться па основе расщепления только по пространственным переменным, как это сделано в диссертации и описано в работах [9,10,56] для многомерных задач газовой динамики, либо с применением расщепления как по пространству, так и по физическим процессам [26]. Такое расщепление позволяет без потери абсолютной устойчивости метода получать экономичные разностные схемы, реализуемые при помощи параллельных прогонок и вычислений по явным формулам.
Большое внимание уделяется итерационным методам решения многомерных и нелинейных задач на МВС, Итерационные методы позволяют использовать более крупные шаги по времени и пространству, уменьшая при этом общие вычислительные затраты. Однако при параллельной реализации требования экономичности итераций вступают в противоречия с условиями их сходимости. В результате быстро сходящиеся в скалярном варианте итерационные методы оказываются совершенно непригодными при параллельной реализации. И наоборот, методы, обладающие худшей скоростью сходимости в скалярном случае идеально распараллеливаются.
Разработкой эффективных численных алгоритмов решения различных классов задач на МВС занимаются во всем мире. Большой интерес представляют работы известных авторов Дж, Донгарра и X. Ван дер Ворста [18,63,64,73]. Дж, Донгарра преимущественно развивает прямые методы решения линейных систем уравнений [63]. В это же время X. Ван дер Ворст показывает [64], что для решения одной системы уравнений с большим числом неизвестных, имеющей матрицу с малым числом ненулевых элементов в каждом ряду, при умеренной точности приближения к решению итерационные методы могут быть предпочтительнее, чем прямые, С точки зрения іребуемой оперативной
памяти итерационные методы также выигрывают. Даже в случае плотных матриц итерационные методы могут быть привлекательнее. Основной вопрос, который по мнению Х.Ван дер Ворста в настоящий момент остается открытым - насколько эффективно итерационные методы могут бьнь ^датированы к архитектурам современных суперкомпьютеров.
Подробный анализ существующих алгоритмов и тенденций в решении линейных систем на МВС содержится в обзорной работе X. Ван дер Ворста [18]. Не менее интересна работа [73], в которой сравниваются оценки производительности различных вычислительных систем на решении плотных систем линейных уравнений, В качестве тестов используются программы пакета LTNPACK. [76].
При решении многих задач математической физики определяющими являются ны числительные затраты на решение разностных аналогов одномерных краевых задач для уравнений 2-го порядка. Поэтому вопрос об их эффективной параллельной реализации является особенно актуальным. Первые работы в этой области появились еще в 1978 году [20] и принадлежат отечественным ученым Н.Н, Яненко и А.Н. Коновалову. В основе предложенного ими метода "распараллеливания" прогонки положен известный принцип суперпозиции решений линейных уравнений [20,3,19], В настоящей работе эта методика распространена на задачи о немонотонным оператором, нестационарные и нелинейные задачи, а также задачи большей пространственной размерности (двух- и трехмерные). На базе разработанных алгоритмов создана библиотека программ для решения одномерных краевых задач для уравнений и систем второго порядка.
В 1991 г, Б.Н. Четверушкиным был предложен параллельный вариант (#-/?) итерационного алгоритма, который тоже может быть назван параллельной прогонкой в двумерном случае [24,77], В работах
[24,29,85] рассматриваются многомерные варианты (а - /3) итерационного алгоритма. Этот метод решения систем сеточных уравнений успешно применялся при моделировании задач динамики излучающего газа, ионосферной физики, течений вязкой несжимаемой жидкости и других. Хотя метод не дает оптимальных по скорости сходимости результатов, однако, позволяет решать широкий круг задач.
В работе [67] В.Я, Карповым рассматривается организация параллельных вычислений на МВС на примере решения уравнения Пуассона в трехмерной области в цилиндрических координатах, В направлении оси z и по углу <р выполняются преобразования Фурье, а в радиальном направлении используется метод параллельной прогонки, основанный на разбиении уравнений системы на группы. К сожалению метод не является универсальным. Он не применим в задачах с переменными или нелинейными коэффициентам и.
В 1990 г. Parallel Computing публикует несколько работ о решении трехдиагональпых линейных систем па параллельных компьютерах с распределенной памятью [14-17], Наиболее интересным является алгоритм DAC (divide and conquer) [17], основанный также на принципе суперпозиции решений линейных уравнении. В связи с этим заметим, что алгоритм, представленный в диссертации и имеющий ту же основу, адаптирован к более широкому классу задач. В частности, он может быть использован для решения многомерных задач. Кроме того, с его помощью возможно решать краевые задачи не только для традиционных граничных условий 1-ого, 2-ого и 3-его рода и смешанных, но и задачи с периодическими и нелокальными условиями.
Разработанная методика не является новой. Однако аналогичные алгоритмы практически не использовалась в параллельных вычислениях. Последнее связано с тем, что теоретическая эффективность (не учитывающая реальное время обменов) исходного
(базового) параллельного алгоритма прогонки на системах с распределенной памятью составляет не более 33%. Тем не менее, данная методика имеет свои положительные стороны, а именно;
является экономичным прямым методом;
позволяет использовать неявные схемы для решения нестационарных задач, и тем самым снизить ограничения, налагаемые соотношением Куранта на шаг по времени;
имеет асимптотическую зависимость ускорения и/или эффективности от числа используемых процессоров р (то есть при соблюдении соотношения р<& N (где N - число неизвестных исходной задачи) и быстрых обменах реальная эффективность распараллеливания близка к 33%;
может эффективно использоваться при решении многих больших задач, в том числе нелинейных и многомерных.
Заметим далее, что все параллельные алгоритмы, рассматриваемые в работе основываются на принципе геометрического параллелизма [27] - параллелизма данных или декомпозиции области 162]. В соответствии с тшм при расчетах на МВС производится разбиение расчетной области на подобласти по числу используемых процессоров. При этом каждый процессор производит вычисления Б своей расчетной полобласти и обрабатывает данные только из своей памяти- Расчет во всех подобластях происходит одновременно, и при необходимости процессоры обмениваются информацией между собой по каналам связи. Такой принцип распараллеливания чрезвычайно удобен при решении задач математической физики на МВС с распределенной памятью.
Помимо геометрического параллелизма различают также алгоритмический параллелизм, называемый также функциональным или параллелизмом потока данных [27], когда алгоритм разбивается на
участки, которые могут выполняться параллельно. Каждый процессор выполняет часть общего алгоритма и обменивается с другими процессорами результатами своей работы. Кроме тою, еще выделяют параллелизм типа "Коллективное решение'' [27], при котором одна и та же программа или ее часть выполняется на разных процессорах над разными начальными данными с последующим объединением результатов. В некоторых задачах используется "гибридный" параллелизм - комбинации геометрического и алгоритмического способов распараллеливания.
В настоящее время в связи с приближением к физическим пределам быстродействия вычислительной и коммуникационной аппаратуры в области параллельных вычислений акценты смещаются в сторону исследования алгоритмов и концепций параллельных вычислений. Основы теории параллельных вычислений изложены в книге Дж.Ортега [6], В работе [65], как и в книге Дж.Ортега, основное внимание уделяется параллельным методам решении задач линейной алгебры, а также детально исследуются параллельные методы решения обыкновенных дифференциальных уравнений. Известная монография В,В.Воеводина [32] также посвящена методологии параллельного программирования. В ней содержится систематизированное изложение математических основ совместного изучения параллельных численных методов и параллельных вычислительных систем. Обсуждается влияние алгоритмических языков, как посредников между математиками и техникой. Автор подчеркивает, что для эффективного применения параллельных компьютеров необходимо радикально изменить структуру численных методов.
Встающие перед пользователями параллельной вычислительной техники проблемы, обсуждаются также в книге [66], Здесь
сформулированы критерии оценки параллельных программ, среди которых:
ускорение программы в зависимости от числа процессоров;
величина затрат времени на синхронизацию;
влияние размера задачи на ускорение;
максимальное число занятых процессоров при решении задачи;
детерминизм выполнения программы.
Эти критерии учитывались в диссертации при анализе параллельных алгоритмов и программ.
Кроме методов и алгоритмов созданы и разрабатываются библиотеки параллельных подпрограмм, ь нервую очередь реализующие методы решения систем линейных уравнений. Наиболее известны пакет ScaLAPACK [68,69], включающий в основном прямые методы решения линейных систем, и универсальные библиотеки фирмы NAG Ltd. - NAG Parallel Library [72] - для систем с общей и распределенной памятью и сетей рабочих станций. Особый интерес представляет новый проект при участии NAG, находящийся в стадии разработки и названный P1NFAPL переносимые параллельные библиотеки для трудоемких индустриальных приложений [70,71] - предлагающий широкий спектр методов, в том числе и итерационных. Несмотря на эти работы, количество библиотек параллельных программ невелико, и доступ к ним бывает затруднен по коммерческим соображениям. Как правило, они содержат универсальные к вследствие этого умеренно эффективные реализации. На реальных приложениях их эффективпость оказывается ниже ожидаемой, поскольку они не учитывают ни специфику конкретной задачи, ни особенности используемой для расчетов МВС, В диссертации параллельные библиотеки не использовались.
Приведем несколько замечаний относительно методов решения выбранной газодинамической задачи. Многообразие и сложность
газодинамических явлений (ударные волны, логранслойные, отрывные, осциллирующие течения, турбулентность) стимулируют исследования Б отой области. Растущие запросы аэродинамики, ракетно-космических технологий и других смежных областей требуют совершенствования традиииоппых и разработки новых численных меюдов решения задач газовой и гидродинамики, в том числе с помощью MRC.
Существенный вклад в область математического моделирования течений жидкостей и газов был внесен российскими учеными. Разработка меюдов численного решения уравнений Павье-Стокса ведется весьма интенсивно с 60-ых годов и нашла отражение в работах известных научных школ под руководством академиков К.И.Бабенко [61], А.А.Самарского [4, 28, 91], Н.Н,Яненко [60] и др. Большое число работ посвящено теоретическим и численным исследованиям динамики вязкого газа [43-46]. Для большого числа методов моделирования газодинамических течений объединяющими являются разностные схемы для решения уравнений Эйлера или Навье-Стокса, построенные на базе этих уравнений.
В отличие от такого традиционного подхода, конструкция кинетических разностных схем (использующихся в диссертации) принимает во внимание тот факт, что сами уравнения Эйлера или Навье-Стокса могут быть получены как следствие более сложных кинетических уравнений. В работе [47] была предложена феноменологическая модель, описывающая перенос частин в газе как циклическую последовательность процессов бесстолкновительного разлета молекул с последующей максвеллизацией, В восьмидесятых годах появилась серия публикаций [48-51], в которых эта модель использовалась для построения кинетических численных алгоритмов. В рамках указанною направления и, исходя из рассмотренной феноменологической модели, в 1984 т\ ТТ. Нлизарова и Б.Н,
\9
Четверушкин выписали в [48-51] замкнутую систему уравнений для макропараметров газа, включающую малый параметр г - характерное время бесстолкновительного разлета молекул. Эта система послужила основой для построения класса кинетически - согласованных разностных схем (к.с.р.с.) [45].
Система квазигазодинамических уравнений возникла в рамках этого нового подхода к построению разностных схем для уравнений газовой динамики и основана на специальном принципе формирования искусственной вязкости [52], При численном моделировании задачи об истечении недорасширенной струи конечно—разностный ал гори їм расчета был построен именно на основе квазигазодинамических (КГД) уравнений [11].
Отметим далее, что необходимость изучения струй, истекающих в область низкого давления, называемых недорасширенными [581, возникает при решении целого ряда практических задач, в том числе связанных с проектированием и использованием разнообразных высотных летательных аппаратов. Особенностями рассматриваемых течений является значительный перепад плотности и давлений от среза сопла к периферии струи и соответственно существенное изменение геометрических масштабов струи — от нескольких миллиметров у среза сопла до нескольких сантиметров вдали от него. При численном моделировании педорасширенных струй требуется использовать сетки, позволяющие разрешить как особенности течения у сопла, так и вдали из него, где происходит взаимодействие струй между собой или с препятствиями [21-23]. Моделирование таких процессов представляет собой задачу, трудоемкость которой требует использования высокопроизводительных вычислительных систем.
Приведем далее структуру диссертации. Она состоит из введения, трех глав и заключения. Во введении обосновывается актуальность темы
диссертации, дается обзор литературы но тематике исследования, формулируются цели и структура работы,
В первой главе формулируется базовый параллельный алгоритм прогонки, на основе которого решаются одномерные краевые задачи для уравнений и систем 2-го порядка, обсуждаются различные типы граничные условий, обосновывается применимость метода параллельной проюнки для решения других задач.
Во второй главе рассматриваются параллельные алгоритмы для решения трехмерной модельной задачи о распространении тепла. Сравниваются варианты расчетов с использованием базового алгоритма и расчетов по явной схеме. Рассматриваются различные варианты разбиения расчетной области и влияние такого разбиения на эффективность вычислений. Отдельно исследуется проблема повышения эффективное!и различных параллельных алгоритмов. Ее удаеіся решить за счет оптимальной организации передачи сообщений между процессорами - сокращения количества операций обмена, устранения обменов мелкими порциями информации, синхронизации обменов па различных внутренних границах подобластей.
Третья глава посвящена решению на МВС задачи о струе, истекающей в область низкого давления. Для этого разрабатываются еоответствуЕОщие параллельные алгоритмы и программы. Дается сравнения результатов моделирования с экспериментальными данными. Кроме того, приводятся результаты целого ряда расчетов, которые были использованы при тестировании большой многопроцессорной системы -MVS-1000M.
В заключении сформулированы основные результаты диссертационной работы.
Отметим далее, чш разработанные параллельные алгоритмы решеттия задач газовой динамики и модельной задачи теплопроводности
были адаптированы к МВС с распределенной памятью Parsytec СС, MVS-1000M и 24-процсссорного Linux-кластера на базе процессоров Intel. Соответствующие параллельные программы написаны па языке Fortran 77 с использованием библиотек и средств параллельного программирования PARTX [7] и МРТ [33].
Результаты работы опубликованы в научных журналах, сборниках и трудах конференций [40, 41, 53-55, 87, 89]. Основные результаты диссертации докладывались и обсуждались на
научном семинаре в Центральном Российском Доме Знаний "Вычислительные системы на базе транспьютеров и параллельные вычисления'1. (Москва, 1992 т.);
IV Международной конференции по математическому моделированию, (г-Москва, 27 икшя-1июля 2001 г.);
международной конференции Parallel CFD'2001: "Parallel Computational Fluid Dynamics". (Нидерланды, г. Эгмонд на Зеє, май 2001г.).
международной конференции "Displays and Vacuum Electronics (DVE 2001)", May 2-33 2001, Garmish-Parteankirche (Германия).
В заключение Автор выражает глубокую благодарность своему научному руководителю к.ф.-м.н. Сергею Владимировичу Полякову за постоянную поддержку и помощь в работе. Благодарю также д.ф.-м.н., проф., чл.-корр. РАН Бориса Николаевича Четверушкина за интерес и внимание к работе. Хочу поблагодарить д.ф.-м.н., проф. Татьяну Геннадьевну Елизарову и к.ф.-м.н. Ирину Александровну Граур за тесное научное сотрудничество. Благодарю к.ф.-м.н. Михаила Владимировича Якобовского и всех сотрудников пятого отдела ИММ РАН за содействие в работе и разностороннюю помощь.
Система алгебраических уравнении
Построенные выше разностные схемы приводят к алгебраической задаче, определяющей искомое разностное решение (стаиионарное, либо относящееся к очередному слою по нремени). Простой анализ показывает, 410 в результате построенных аппроксимаций мы получим систему уравнений с трехдиагопалыюи матрицей, или близкую к ІЇСЙ. Этот факт записывают с помощью следующей системы уравнений
В граничных узлах сетки в случае 1-ой, 2-ой, 3-ей или смешанной краевых задач получаются уравнения -С о + воУ\ =-FQ "СыУы + Av Vi = -/V (1 л ]) В случае периодических граничных условий уравнения (1.10) имею і место для всех і = 2,...,/V- 1, а в узлах і = І иі- / справедливы уравнения ЛУн-СіУі+ВіУг -Гг. АыУч CNy, + BKly -F,r (1,12)
В случае нелокальных граничных условий первое и последнее уравнения системы принимают вид -Q.yo + Z = Z -Q =;v 0-13) Относительно свойств коэффициентов системы (1.10)-( 1Л 1 -1.13) предполагается следующее: ДД 0, Ct 0, / = 0,...,ЛГ, (1.14) и либо С Л, + Д, + Ц Д- 0, z=0,...,A , (L15) либо
Последнее условие возникает в случае немонотонного оператора (rc 0) во внутренних узлах сетки. На границе в этом случае должны быть выполнены условия (1Л5). Если краевая задача является нелокальной, то должны выполняться условия ее разрешимости.
Базовый алгоритм распараллеливания
Рассмотрим алгоритм решения полученной алгебраической задачи на примере системы (1Л 0)-(1 ЛI), предполагая использование MB С с распределенной памятью с р процессорами. Для этого используется подход "геометрического параллелизма": всю расчетную область разобьем на число подобластей, равное числу процессоров системы. Внутри каждой подобласти возьмем примерно равное число счетных узлов. Если использовать однородные схемы (алгоритмы, выполняемые по единому правилу во всех узлах), то такое построение автомаїически обеспечивает равномерную загрузку процессоров.
В результате вводится равномерное линейное разбиение множества номеров узлов сетки П = {0,1,...,N} на связные подмножества Ob) ;{Щ) о,= Г І [т -0,...,р-1 )5 соответствующие разбиению вектора неизвестных по процессорам (см, рис.1). х-0 ч -+ р-1 х=1 —I 1 1— te-n . Ю1 І . (0) . (1) 2-І Рис. 1. Разбиение расчетной области по процессорам. &L)=N При таком разбиении процессор с номером т будет обрабатывать (ґ2т—il + 1) точек. Представим далее решение на каждом внутреннем (0 т р-1) процессоре в виде следующей линейной комбинации: v = v n, = v(i№} + v v s,!-mJ I v vilf ll! (1.17) где функции yj ,m (o. = IJIJIl) полностью определены на множестве узлов 1/п и играют роль некоторого базиса (см. рис. 2), а значения искомой функции па границе Пч - _у„„ и _у ,„, - пока не известны. Во внутренних узлах Qfl функция y(tI,fT,) находится из уравнений (1.10), а функции ,(/ "\ уУ/Лт) из уравнений (1.10) с нулевой правой частью.
Параллельные алгоритмы численного решения на МВС уравнения теплопроводности
Для получения равномерной загрузки процессоров количества точек /Vu} = A /J,1 -/ } + 1)Л 3 в диапазонах 1[к] должны быть примерно одинаковы. Расчетные формулы для вычисления решения на верхнем слое выглядят следующим образом: з=Е/С +гК ,X + 4J. (І ОеҐ . (2.11)
Топология межпроцессорных связей в этом случае будет линейной, поскольку в формуле (2,11) для процессора к участвуют значения функции [/,, в гачках (1 - ) = ( 1 l},J3) и (hA \} + 3)-(4 20 h) лежащие в области данных процессоров (Ы) и (А+1). Перед началом очередного шага яычислений по формулам (2.11) на процессоре к необходимо получить требующиеся значения Uh от соседних процессоров- Для этого лучше всею произвести один ірупиовой обмен данными для всех индексов і,. /3. Этот обмен легко организуется в рамках линейной топологии. Количество передаваемых данных одним процессором в общем случае равно М к) =2N ЛГ3 и не зависит от числа процессоров. В случае равного количества узлов сетки по различным направлениям {Nl =N2 = N2-\N ) величина М{к) =2N 2i\
В случае "квадратного" разбиения области, пример которого показан на рис, S, формулы (2.11) остаются в силе, но изменяются диапазон индексов I{ti) и количество точек Nu\ обрабатываемых процессором к: =),/,, ),/, = 0,...,tf„/2 =/ ,..., ,1 = ,..,/1 = /, х/ х/ , = ,(, -, +1)( - +1).
Соответственно изменяегся и топология обменов и межпроцессорных связей, коюрые образуют теперь решетку. Для определения соседей процессора Л, с которыми ему придется обмениваться данными, можно ввести двойную нумерацию процессоров и установить взаимно однозначное соответствие исходных и двойных процессорных номеров. В дашюм случае можно использовать следующие соотношения: к = рък2+к к2 =0,..., -1, 3 =0,...,/ -1, Аа = [ /р3] L=kmodpiri A=0,...,p-ls р р2ръ. Здесь /? - общее число процессоров, образующих решетку р2 х Рз - В этих обозначениях процессор Л - (fe,, къ) перед началом очередного шага вычислений будет обмениваться недостающими данными с четырьмя соседями к ++{к2-икг), r ofe + Uj, кфт (к2 к,-І), к" г (к2,к + 1) в рамках топологии "решетка". Количество передаваемых данных М( ) = 2Nt(N2/р2 + Ni/р2) и уже зависит от числа процессоров. При равном числе узлов сетки по различным направлениям и квадратной решетке процессоров (р2 = рь = ]р ) величина М(Л) =4/V Mp" ;\
В случае "кубического" разбиения области (см. пример на рис. 9) также получим расчетные формулы (2.] 1), в которых = ( - + 1)( - + 1)( -. -1).
Для определения соседей и обмена с ними данными вводится новая тройная нумерация процессоров: к- РгРА+р- +к ,=0,...,/7,-1, А2 =0,...,/-1, =0,..., -1, = 0,..., /7-1, P = P\P,Pj.
В соответствии с ней каждый процессор в кубической решетке имеет в общем случае 6 соседей, с которыми и обменивается необходимыми данными в начале каждого шага вычислений. Количество передаваемых данных Мт =2N[N2f(-plp2) + 2N]N1!(plpi)-r2N2Ni/(p2p3). При равном числе узлов сетки по различным направлениям и равномерной кубической решетке процессоров (/?] = рг=ръ - ifp ) М[к) =6N2np 2П.
Теоретическая эффективность (не учитывающая время обменов) рассмотренных трех алгоритмов одинакова и равна 100%. На практике реальная эффективность алгоритма зависит от схемы обменов и числа процессоров. При прочих равных условиях использование линейной топологии при большом числе процессоров будет менее эффективно, чем использование топологий типа «решетка», поскольку у последних размер передаваемых данных снижаемся с ростом числа процессоров.
Рассмотрим теперь параллельную реализацию неявной ЛОС (2.9). Для этого тоже можно использовать три типа разбиения расчетной области и три типа топологии межпроцессорных связей.
При линейном разбиении области, например, по координате х2, получается следующий параллельный алгоритм. Для вычисления решения на новом временном слое каждый процессор последовательно решает три одномерные задачи; (-rA e/3)t//+e J=t/re",1 3+j/iy",a = l 2,3, (2.12) в своем диапазоне индексов 1{к). При этом оказывается, что первую и третью задачи (первый и третий этапы ЛОС) можно решать независимо от данных па других процессорах и использовать для этого алгоритм обычной (скалярной) прогонки. Во второй задаче (второй этап ЛОС) данные оказываются распределенными между процессорами. Поэтому вторую задачу необходимо решать с помощью какого-либо параллельного алгоритма, В качестве такового предлагается использовать алгоритм параллельной прогонки, разработанный в гл. 1. Он применяется по индексу і2 для всех фиксированных значений
КГДуравнения в r-z геометрии и постановка задачи
Первое из условий представляет собой известное условие прилипания газа к стенке, второе условие обеспечивает отсутствие потока массы через границу, третье условие влечет равенство нулю па стенке нормальной составляющей теплового потока q.
При любых сколь угодно малых числах Кнудсена вблизи поверхности тела существует область, имеющая толщину порядка длины свободного пробега молекулы, так называемый кнудсеновский слой. Для учета этого слоя вблизи стенки в рамках теории Навье- Стокса вводятся так называемые условия скольжения для скорости и скачка для температуры. В литературе существует достаточно много вариантов постановки граничных условий данного типа [80,81]. В настоящей работе использовались условия, полученные в [80]: 1.012 ГГ (дйЛ 2а у-\ т p4lRT Іай J, здесь «t - скорость скольжения на стенке, Т, - температура газа на стенке, Ти. - температура стенки. Условия для скорости скольжения были получены в предположении, что молекулы отражаются от стенки диффузно с максвслловским распределением. В условия температурного скачка входит коэффициент аккомодации а [80]. В расчетах коэффициент аккомодации полагался равным единице, то есть предполагалось, что налетающие на стенку молекулы приходят в термическое равновесие со стенкой. Выпишем условия для скорости и температуры на стенке, использовавшиеся в расчетах: u„ =0, и, =HV, У = 7;, (3,14) где и„ и и. нормальная и тангенциальная компоненты вектора скорости и. Как уже было упомянуто выше, для системы КГД уравнений требуется постановка донолншельного граничного условия, это условие непротекания в виде (3.12) (второе соотношение). Реализация этого условия приводит к дополнительному условию для давления на твердой стенке Б виде: дР = о. (3.15) СП
Расчетная область схематично показана на Рис.19. Сопло находится в точке z=0. Индекс е соответствует параметрам на выходе из сопла. Радиус сопла - гг. Предполагается, что на выходе из сопла профили скорости и температуры соответствуют ламинарному пограничному слою юлщиной 5 = 0.18/;. Стенки сопла рассматриваются как адиабатические. Для задания профилей плотности и температуры на выходе из сопла используются соотношения (3.13): ,п и» Vfr. (3.16) v «, - - У Ш - - Вне пограничного слоя параметры потока полагаются постоянными. С D г=г 2=0 Ма. Е і г=0 2=2 „, ось симметрии твердая стенка внешние границы Рис.19. Расчетная область.
Выпишем условия на границах расчетной области: на левой границе на срезе сопла внутри пограничного слоя (re-5 r rt.) ЩІГ) = н„(г), и, = 0, р{г) = Рс, Т(г) = 7», на срезе сопла вне пограничного слоя (0 r re- S) и,{г) = и„ иг=0, р(г) = Ре! T(r) = Ttt вне сопла (на твердой стенке): условия прилипания и адиабатической стенки 1(.=0, и, = 0, - = 0, — = 0 , дп дп или условия скольжения и температурного скачка в соответствии с (3.13): u:(r) = 0, ur=u„ - = 0, T = Tt. дп на правой границе стоят гак называемые "мягкие" граничные условия: uz & 3z dz на верхней границе стоят условия для температуры и давления окружающего газа: = 0, «,=0, р = р„ Т = ТВ. (3.17) на оси симметрии: - = 0. ur=0, =0, = 0. й- й- й" Индекс ад соответствует параметрам окружающего газа.
С точки зрения разработчика параллельных программ существуют N вычислительных узлов, каждый из которых обладает собственным процессором, оперативной памятью, а также, возможно, дисковой памятью и периферией. Обмен информацией реализуется путем передачи сообщений, при этом любой узел может пересылать сообщения любому другому. Один и тот же выполняемый файл . загружается во вес используемые вычислительные узлы, и создаются однотипные процессы (по принципу "одна программа - разные данные").
Программа расчета течения написана на языке FORTRAN 77. В вычислениях используется библиотека функций, предназначенная для поддержки работы параллельных процессов - MPI. Существует большое число функций MPI для пересылки данных. В данной работе использовались функции MPI_Send, MPI_Recv, MPI_Bcast, MPI_Reducc.
Параллельный алгоритм основывается на принципе геометрического параллелизма. В соответствии с ним расчетная область разбивается на подобласти по направлению z, как это показано на Рис.20.
Схема Бима - Уормгтга для решения стационарных задач
Вычисления проводились для струи азота и параметров па выходе из сопла, приведенных в Таблице 7. Использовались четыре различных значения внешнего давления в соответствии с Таблицей 8. Все использованные в расчетах данные соответствуют условиям эксперимента, опубликованным в [84]. В эксперименте были измерены значения числовой плотности и вращательной температуры вдоль оси струи, а также определено положение диска Маха. Использованная нами вычислительная модель позволяет получать лишь значения средней температуры, в то время как измеряется именно вращательная температура. Однако отличия в этих температурах невелики и заметны лишь R зоне ударной волны, что позволяет использовать для сравнения с экспериментом среднюю температуру.
Лоле течения Основные результаты расчетов представлены на рисунках 22-32. Па рис, 31 представлены линии уровня плотности для варианта А, На этом рисунке хорошо видны все известные из теории особенности истечения сверхзвуковой струи в область низкого давления. Истечение такого типа характеризуется сильными градиентами плотности, температуры и давления в струе. Вытекая из сопла, газ расширяется, образуя висячий скачек уплотнения, имеющий бочкообразную форму. Взаимодействие этого скачка с осью симметрии приводит к образованию центрального скачка (диска Маха), за которым сверхзвуковой поток становится дозвуковым. За диском Маха возникает небольшая рециркуляционная зона, видная на рис. 32. За первым центральным скачком возникает серия вторичных скачков, которая также видна на линиях уровня плотности. Заметим, что в данном эксперименте непосредственная визуализация области возвратного течения невозможна, и имеются лишь косвенные данные, подтверждающие его существование, влияние сетки Для варианта А расчеты были проведены для четырех сгущающихся ссткак (Таблица 9). На рис. 28 изображено распределение числовой плотности, нормированной на свое минимальное значение, вдоль оси струи. Здесь же приведены значения, измеренные в эксперименте [84]. На рис, 2У приведено распределение температуры вдоль оси струи (в градусах Кельвина) для всех трех вариантов сеток вместе с измеренными значениями.
Для всех использованных сеток хорошо разрешается первая, наиболее интенсивная ударная волна. Излом в распределении плотности, полученный в эксперименте, лишь начинает проявляться на самой подробной третьей сетке. Положение всех вторичных скачков хорошо разрешается па всех сетках, однако числовые значения плотности ближе всего к экспериментальным для самой подробной сетки. Для значений температуры измельчение сетки позволяет лучше разрешить температурный минимум до ударной волны и не оказывает практически образом, для правильного разрешения всех особенностей течения необходимо использовать достаточно подробные сетки, с шагом по направлению z не больше, чем 0.25/; , что возможно лишь при использовании никакого влияния на распределение температуры во вторичных скачках. Таким многопроцессорных вычислительных систем, позволяющих проводить расчеты за приемлемое время даже на очень подробных сетках.
Влияние внешнего давления На рис. 22 - 27 изображены распределения плотности и температуры для вариантов В, С и D вместе с соответствующими измеренными значениями числовой плотности и вращательной температуры. Число Кнудсена (Кп) (отношение длины свободного пробега молекулы, посчитанной по параметрам течения в данной точке к диаметру струи) на выходе из сопла достаточно мало и составляет Кп = 3747 Ю-4. В поле течения струи локальное число Кнудсена становится достаточно большим, особенно в области перед диском Маха. Уменьшение внешнего давления от варианта А к варианту D приводит к увеличению локального числа Кнудсена. При этом совпадение расчета с данными эксперимента ухудшается, что особенно заметно на распределениях температуры. Это можно объяснить тем, что континуальные модели, и используемая в расчетах КГД модель, теряет свою применимость при больших числах Кнудсена.
Адекватное описание поля течения в недорасширенной струе требует использования подробных вычислительных сеюк, что связано с большим объемом вычислений и тем самым приводит к необходимости использования высокопроизводительных вычислительных систем.
Проведенные расчеты показали, что предложенный алгоритм (явная по времени и однородная по пространству аппроксимация КГД уравнений) позволяет эффективно использовать кластерные МВС, Сравнение численных и экспериментальных результатов показывает, что предложенный численный метод позволяет описывать основные особенности истечения осесимметричных струй в область низкого давления. При этом результаты расчетов позиолякл получить дополнительную информацию о течении, недоступную в натурном эксперименте.
