Содержание к диссертации
Введение
Глава 1. Универсальная технология параллельных вычислений для численного моделирования задач, описываемых системами уравнений гиперболического типа 25
1.1 Постановка задачи 25
1.2 Формализация задачи и унификация расчетных методик 27
1.3 Конечно-разностные схемы 34
1.4 Обобщение монотонной гибридной схемы второго порядка для неравномерных сеток 40
1.5 Применение метода крупных частиц 46
Глава 2. Реализация параллельного программного пакета для численного моделирования задач, описываемых системами уравнений гиперболического типа 51
2.1 Описание программного пакета 51
2.2 Объектная модель 55
2.3 Декомпозиция задач численного моделирования, описываемых уравнениями гиперболического типа, на параллельных компьютерах 61
2.4 Оценки эффективности при пространственной декомпозиции задач на параллельных компьютерах 65
2.5 Тестовые расчеты 68
Глава 3. Численное моделирование пространственных течений в атмосфере, вызванных крупномасштабными возмущениями 72
3.1 Введение 72
3.2 Математическая модель 73
3.3 Тестовые расчеты 76
3.4 Моделирование пожаров 78
3.5 Течения с вихревой структурой 81
Глава 4. Численное моделирование прямой задачи геофизики 85
4.1 Обзор прикладной области 85
4.2 Математическая модель 86
4.3 Тестовые расчеты 88
4.4 Моделирование прямой задачи геофизики в неоднородной среде 89
Заключение 93
Таблицы 94
Рисунки 95
Список литературы 128
- Обобщение монотонной гибридной схемы второго порядка для неравномерных сеток
- Оценки эффективности при пространственной декомпозиции задач на параллельных компьютерах
- Моделирование пожаров
- Моделирование прямой задачи геофизики в неоднородной среде
Введение к работе
Математика всегда черпала свои задачи из физики и других естественных наук. Но и методы исследования этих наук сильно обогатились за последнее время. Там где сложно, дорого или вообще невозможно производить физический эксперимент, прибегают к математическому моделированию. И, если раньше моделирование ограничивалось только аналитическими методами, то с появлением вычислительной техники появился новый способ исследования -вычислительный эксперимент [1]. То, что невозможно получить аналитически, зачастую легко моделируется на ЭВМ. Кроме того вычислительный эксперимент гораздо дешевле физического, а если требуется серия исследований, то применение ЭВМ становится еще более привлекательным. Вычислительный эксперимент необходим там, где есть угроза жизни человека (например, летчика-испытателя или физика-ядерщика). С помощью вычислительного эксперимента можно смоделировать поведение экологически опасного объекта в различных условиях, дать практические рекомендации обеспечения для условий безопасной работы [2].
Большинство задач механики сплошных сред невозможно решить явно в аналитическом виде. Связано это с нелинейностью рассматриваемых уравнений. Так, например, для большинства задач газовой динамики вообще не доказано существование и единственность решения соответствующих систем уравнений [3]. Поэтому применяют также эмпирические и полуэмпирические методы. Сложности также встречаются на пути формализации задачи и построения математической модели, поскольку физическая постановка задачи не всегда формализована и может включать множество дополнительных факторов. Адекватность построенной модели и использованных методов моделирования должны проверяться путем сравнения с результатами экспериментов, других расчетов, аналитических выкладок.
Модели в основном строятся как численное решение дифференциальных уравнений в частных производных. Для решения таких задач используется численное моделирование [1,4]. При численном моделировании осуществляется переход от непрерывной модели, выражающейся в записи исходной системы дифференциальных уравнений и граничных условий, к дискретной модели, представляющей собой систему алгебраических уравнений большой размерности, получаемую на основе различных разностных схем исследуемых уравнений. В дискретных моделях дифференцирование заменяется разностными операторами, интегрирование - суммированием и т.д.
Для широкого круга задач система исследуемых уравнений в частных производных оказывается гиперболической [3] (т.е. в каждой точке может быть расщеплена на N независимых одномерных уравнений переноса в характеристических переменных) или содержит гиперболическую часть (при использовании принципа расщепления по физическим процессам один из этапов численного решения описывается гиперболической системой уравнений). Это связано с тем, что одним из основных механизмов в механике сплошных сред является конвекция - перенос основных физических компонент: массы, импульса, энергии.
Основными представителями таких уравнений являются: уравнения гидро- и аэродинамики, акустики, упругости, магнитной гидродинамики, уравнения "мелкой воды" и другие. Каждая из этих областей включает несколько разделов, каждый из которых содержит множество задач, как решенных, так и открытых для исследования. Таким образом, рассматриваемый в диссертации класс гиперболических уравнений охватывает большое число важных задач, как с теоретической, так и с практической точек зрения.
В последнее время бурное развитие вычислительной техники позволило производить ресурсоемкие расчеты [5] и получать физически оправданные результаты для 3-мерных задач [6,7]. Развитие аппаратных средств в свою очередь стимулирует появление новых программных технологий и ориентацию на языки высокого уровня. Используя эти новые технологии можно обеспечивать большую гибкость и надежность программного обеспечения, разрабатывать сложные программные системы с меньшими затратами, как вследствие высокой степени абстрактности и простоты используемых средств разработки, так и за счет повторного использования программ.
Кроме того, использование суперкомпьютеров и параллельных программ позволяет повышать производительность вычислений в порядки раз [5]. Высокая производительность вычислений выводит численное моделирование на качественно новый этап развития, на котором возможно решение качественно новых задач: либо вследствие достижения необходимой точности вычислений, либо вследствие возможности расчета на больших сетках (или, там, где не используются сетки, с большим числом расчетных параметров).
К сожалению, развитие программного обеспечения для суперкомпьютеров отстает от общего развития технологий и в области, где доминируют персональные компьютеры. Это связано в первую очередь с разницей в алгоритмах решения одних и тех же задач на многопроцессорных ЭВМ и однопроцессорных. Многопроцессорные вычислительные системы позволяют проводить вычисления одновременно на нескольких процессорах, что обеспечивает большую производительность вычислений, но только в том случае, если алгоритм реализует параллельную модель вычислений. Поэтому приходится проводить дополнительную работу по распараллеливанию последовательного алгоритма, либо придумывать новый [8]. Кроме того, используемые средства для разработки параллельных программ сами по себе довольно сложны [8,9]. Высокая цена многопроцессорных вычислительных систем и их небольшой объем использования по сравнению с другими ЭВМ также порождают отставание технологий в области параллельных вычислений от общего развития технологий разработки программного обеспечения [10,11]. Поэтому разработка специальных технологий для параллельных вычислений и эффективных параллельных алгоритмов для решения обширного круга прикладных задач в различных областях численного моделирования является чрезвычайно актуальной задачей и способствует развитию технологий параллельной обработки информации в целом.
Данная диссертационная работа и программный пакет, созданный по ее результатам, призваны заполнить некоторый пробел в области параллельных вычислений и технологий программного обеспечения для многопроцессорных вычислительных систем. Целью настоящей диссертационной работы является создание универсальной технологии параллельных вычислений для численного моделирования задач, описываемых системами уравнений гиперболического типа, создание на ее основе прикладного программного пакета, и проведение с помощью этого пакета исследования нескольких задач численного моделирования в двух различных прикладных областях. Первый класс задач включает численное моделирование пространственных атмосферных течений, вызванных крупномасштабными возмущениями (термики, пожары, взрывы, огненные смерчи, торнадо) [194,182]. Вторая задача - численное моделирование прямой задачи геофизики [12], которая заключается в рассмотрении прохождения звуковых волн (продольных и сдвиговых) в неоднородной земной коре.
Численное моделирование - сравнительно молодая наука. Однако первые работы по этой тематике появились еще до возникновения компьютеров. В 1910 г. Л. Ричардсон [13] опубликовал статью, посвященную итерационным методам решения дифференциальных уравнений в частных производных. Кроме того, что Ричардсон разработал сами етоды, получил оценки сходимости, изучил задание граничных условий, он также произвел по ним первый расчет практической задачи - расчет напряжений в каменной дамбе. Эта статья послужила отправной точкой исследований по численному моделированию. В 1923 г. Филлипс и Винер [14] получили первое строгое доказательство сходимости метода Либмана, в 1918 г. усовершенствовавшего метод Ричардсона для решения эллиптических уравнений [15]. В 1928 году Курант, Фридрихе и Леви [16] опубликовали теоретическую работу по конечно-разностным методам, ставшую основополагающей и определившую развитие этих методов. В этой работе исследовались такие вопросы, как сходимость решений разностных уравнений к решениям дифференциальных, устойчивость решения, доказывались теоремы о существовании и единственности решения как для эллиптических и параболических, так и для гиперболических систем уравнений.
Первое численное моделирование в гидродинамике было проведено Томом в 1933 г. [17]. Был произведен расчет течения вязкой несжимаемой жидкости в следе при обтекании кругового цилиндра. В 1955 г. Аллен и Саусвелл [18] применили метод релаксации Саусвелла [19] для расчета обтекания циллиндра вязкой жидкостью при различных числах Рейнольдса. Результаты вычислений показали наличие неустойчивости течения в определенном диапазоне чисел Рейнольдса. Следом за Саусвеллом различные методы релаксации были разработаны Фоксом в 1948 [20], Франкелом в 1950 [15] и Янгом в 1954 [21]. В
1950 г. были впервые проведены расчеты крупномасштабных метеорологических задач [51]. В той же статье были опубликованы результаты работ фон Неймана, проведенные в Лос-Аламосской лаборатории во время второй мировой войны, в том числе критерий устойчивости для параболических задач. С развитием ЭВМ и внедрением их в научные лаборатории, стали появляться методы, ориентированные в первую очередь на решение задач на ЭВМ, например, методы типа Монте-Карло. Были разработаны новые методы для решения задач, которые ранее не могли быть рассчитаны, например, для многомерных нестационарных задач гидро- и аэродинамики. В книге Рихтмайера 1957 г. [22] было приведено более десяти различных численных методов для решения задач нестационарной гидродинамики. После опубликования статьи Харлоу и Фромма в 1965 г. [38] за вычислительной гидродинамикой закрепляется статус становится самостоятельной науки. Термин вычислительная гидродинамика (computational fluid dynamics, CFD) устоялся с выходом в 1976 г. одноименной книги Роуча [24].
Для гиперболических систем уравнений создано большое количество разнообразных методов решения [1,25,26]. Наряду с прямыми аналитическими методами нахождения решения используются точные методы, использующие разложение в ряд Фурье и другие функциональные ряды, метод разделения переменных в многомерных задачах, также используется метод потенциалов, представляющий решение в виде интеграла, зависящего от параметра. Однако, они дают решение в удовлетворительном виде только для ограниченного количества задач, не включающего наиболее практически ценные.
Среди численных методов решения можно выделить несколько отдельных групп: характеристические методы [27], в которых аппроксимация системы осуществляется на характеристической сетке, строящейся в процессе счета; наиболее многочисленный класс конечно-разностных методов [28,3], включающий как эйлеровы, так и лагранжевы методы [29]; сеточно-характеристические методы [30], в которых расчет осуществляется в характеристических переменных на обычной нехарактеристической сетке (эйлеровой или лагранжевой). метод конечных элементов [31], в котором ищется приближение точного решения в определенном функциональном пространстве; полуаналитические методы, такие как: метод прямых [32], в котором сетка строится только для пространственных переменных, а временная переменная остается непрерывной, таким образом исходные уравнения аппроксимируются большой системой обыкновенных дифференциальных уравнений, обобщающий его метод интегральных соотношений [33], методы динамики вихрей [34], в которых вихревое поле моделируется большим количеством дискретных вихревых частиц, а исходные уравнения аппроксимируются ОДУ, описывающими положения дискретных вихрей; методы квазичастиц, в которых определенные объемы сплошной среды моделируются конечным числом квазичастиц, такие как: метод частиц в ячейках [35], метод крупных частиц [36], метод сеток и частиц [37], метод маркеров в ячейке [38], во всех данных методах имеется сетка ячеек, через которую осуществляют движение различного рода квазичастицы; статистические методы, такие как, методы типа Монте-Карло [39] или статистический метод частиц в ячейках [40].
Конечно-разностные методы и схемы можно классифицировать по шаблону: явные, неявные, полунеявные; по количеству временных слоев в шаблоне: двухслойные, трехслойные итд; по размерности: одномерные, двумерные, трехмерные; по типу шаблона: с регулярным шаблоном, с нерегулярным; по однородности: схемы сквозного счета [3,43] и схемы с явным выделением особенностей [41]); по алгебраическому виду: линейные, квазилинейные, гибридные, нелинейные [42]; и по способу получения: с помощью разностных аппроксимаций [28], с помощью интегро-интерполяционных методов [43], например, консервативный метод потоков [44], и с использованием метода неопределенных коэффициентов [30].
Одной из базовых конечно-разностных схем является схема Куранта-Изаксона-Риса [45], являющаяся также схемой бегущего счета, вследствие того, что при расчете по ней вычисления бегут от одной границы к другой, что удобно для решения смешанной задачи Коши. Фон Нейман и Рихтмайер [46] в схеме типа "крест" впервые предложили метод использования псевдовязкости для обеспечения устойчивости решения путем сглаживания контактных разрывов. Лаке в 1954 [47] предложил центрально-разностную схему, получившую имя Лакса-Фридрихса. В этой работе Лаке также впервые использовал консервативную форму для записи уравнений гидродинамики. В широко известной схеме Лакса-Вендрофа [48] вязкость аппроксимации не приводит к размазыванию контактного разрыва и обеспечивает второй порядок аппроксимации. Рихтмайер [49] сформулировал более экономичную двухшаговую схему, приводящую к той же схеме Лакса-Вендрофа, на первом шаге которой используется схема Лакса-Фридрихса, а на втором - схема "чехарда". Мак-Кормак [50] предложил двухшаговый метод для квазилинейного уравнения, на разных шагах которого вычисляются разности вперед и назад. Метод Мак-Кормака совпадает с методом Лакса-Вендрофа на линейных уравнениях, а для нелинейных - дает лучшие результаты по сравнению с методом Лакса-Вендрофа на движущихся разрывах.
Одновременно развивался и математический аппарат исследования численных методов. Фон Нейманом, Рихтмайером, а также Лаксом в 50-х годах [46,51, 52] были исследованы вопросы устойчивости разностных схем, в частности, был получен фундаментальный критерий устойчивости фон Неймана [53]. Фундаментальная работа Годунова [54] ответила на главный вопрос о несуществовании линейной монотонной схемы для уравнения переноса выше первого порядка аппроксимации и дала направление исследований для построения монотонных схем. В той же работе была сформулирована монотонная схема Годунова, имеющая первый порядок, на каждом шагу которой в каждых соседних ячейках точно решалась задача Римана о распаде разрыва. Годунов также изучал класс схем, обладающих свойством сохранения энтропии (введенным им же) [55,56]. Ошер описал общий критерий сохранения энтропии для конечно-разностных схем для уравнений с постоянными коэффициентами [57].
Вслед за работой Годунова, первым введшего аналитическое решение задачи Римана в кончено-разностный метод, было разработано множество методов, базирующихся на точном или приближенном решении задачи Римана. Мак-
Намара [58] модифицировал метод Годунова, введя подвижную сетку, перестраивавшуюся под положение контактного разрыва. Схожий с методом Годунова метод для решения уравнения Бюргерса был предложен Глиммом [59], который также решал задачу с кусочно-постоянным распределением на соседних ячейках точно, а затем использовал процедуру реконструкции решения. Еще до работ ван Лира Колган [60] предложил метод, использующий процедуру ограничения потоков, полученных с помощью метода Годунова. В середине 70-х годов ван Лир публикует серию статей [61-65], в которых описывается техника создания монотонных схем путем построения или модификации существующих схем с помощью процедуры линейной реконструкции решения и ограничения потоков. В последней из этих статей [65] рассматривается известная схема второго порядка по пространству MUSCL, основанная на приближенном решении задачи Римана в Лагранжевых координатах. В этой схеме был также впервые реализован способ вычисления частичных разностей для системы уравнений. Методы, реализующие данный подход к вычислению разностей, получили в дальнейшем название методов MUSCL-типа, например, метод 4-го порядка, предложенный Ямамото и Дайгуджи [66].
В работе [67] Магомедов и Холодов изложили общий сеточно-характеристический подход для методов решения гиперболических систем уравнений. В этой работе была также предложена векторная схема для решения векторной задачи Римана, обобщающая схему Куранта-Изааксона-Рисса. Аналогичная запись была использована Роу [68,69], однако здесь уже речь шла о методе второго порядка. Роу также предложил осреднение матриц Якоби, при котором решение задачи Римана для линеаризованного уравнения с осредненной матрицей совпадает с решением нелинейного уравнения. Роу и Пайком было доказано [70], что такое осреднение единственно. Кроме того, Роу предложил ограничитель SuperBee [71] и ввел понятие функции-ограничителя, формализующее уже использовавшиеся ранее методы ограничения потоков в других схемах: схемы ван Лира [65], Колгана [60] и Хартена [72] с широко известным ограничителем MinMod. Различные авторы использовали большое количество различных ограничителей, среди которых следует отметить ограничители ван Албады [73], ограничитель, использованный Хартеном и
Ошером [74] в схеме UNO второго порядка, "острый" ограничитель, предложенный Тишкиным [92,93], дающий схему третьего порядка. В такой же форме с ограничителем можно представить и кусочно-параболический метод РРМ Колелла и Вудворда [75].
Альтернативный подход к методу Роу также на основе приближенного решения задачи Римана был предложен Ошером [76,77]. Так в методе Энгквиста-Ошера [77] для вычисления потока использовалось интегральное осреднение матрицы Якоби. Методы Роу и Ошера базируются на представлении решения в виде суперпозиции дискретных волн, движущихся в прямом и обратном направлениях, поэтому они получили название методов с разностным расщеплением потоков. Другой тип методов с векторным расщеплением потоков основывается на представлении решения в виде взаимодействующих псевдочастиц. Так, например, расщепление ван Лира [78] основано на Больцмановском распределении скоростей частиц, а расщепление Стигера-Уорминга [79] исходит из предположения, что половина частиц движется с одной характеристической скоростью, а другая половина - с другой скоростью.
Имеется несколько схожих подходов к построению монотонных схем. Так, Ии, Уорминг и Хартен [80] выделяют 4 принципа построения таких схем: гибридные схемы, подход ван Лира с использованием реконструкции решения для решения задачи Римана, подход Роу на основе метода расщепления потоков, подход Хартена для TVD методов с модификацией потоков. Первыми из них были получены и стали использоваться монотонные гибридные методы, переключающиеся между линейными методами высокого порядка в области гладкости решения и методами первого порядка на разрывах. Среди них следует отметить гибридные методы Хартена [81], Хартена и Цваса [82], ван Лира [61] и Бима-Уорминга [83]. Хотя метод коррекции потоков FCT Бориса и Бука [84,85,86] первоначально не был отнесен к этому классу методов, было показано, что он может трактоваться как гибридный метод. Различные гибридные схемы строились и изучались Магомедовым и Холодовым с помощью подхода с использованием неопределенных коэффициентов и пространств схем [30]. В работе Белоцерковского, Гущина и Конынина [87] была предложена гибридная монотонная разностная схема второго порядка по времени и по пространству.
Развитие TVD схем связано с именем Хартена. TVD-свойство, т.е. свойство схем уменьшать полную вариацию функций, является более слабым требованием, чем монотонность. Однако за счет ослабления требований к схеме удается построить схемы более высокого порядка, тем не менее, не осциллирующие на разрывах решения. В работах [72,88,80] Хартен ввел класс TVD-схем и изучил их свойства, а также предложил известную TVD-схему Хартена второго порядка. Широко известен также метод энтропийной коррекции Хартена [89], применяемый во многих TVD-схемах. Было опубликовано много TVD-схем, в том числе и более высоких порядков. Среди них стоит упомянуть схему Чакраварти-Ошера [90], модифицированную впоследствии Ямамото и Дайгуджи [66]. Дальнейшее развитие этот подход получил в схемах ENO [91], UNO [74], отметим UNO-схему третьего порядка Янга [94]. Схемы ENO и UNO также не являются в общем случае TVD-схемами, как TVD-схемы не являются монотонными, но имеют более слабые свойства, не допускающие нефизичных осцилляции.
До появления компьютеров вычисления производились вручную и с применением арифмометров. Скорость таких вычислений была крайне невелика. С появлением первых ЭВМ в начале 1940 гг. численное моделирование стало стремительно развиваться. Развитие этой области стимулировалось необходимостью проведения сложных расчетов для военных задач, таких как наведение зенитных орудий по быстродвижущимся самолетам и т.п. Первая ЭВМ общего назначения ENIAC, появившаяся в 1946 г., могла производить около 5000 сложений или 14 умножений в секунду и обладала оперативной памятью в 20 слов. К тому времени, как она была готова к работе, задачи наведения, для решения которых ENIAC разрабатывался, уже потеряли актуальность. Однако уже в 1945 г. в ходе испытаний этой машины С. Франкелем и Н. Метрополисом из Лос-Аламосской лаборатории была проведена серия расчетов термоядерных реакций при взрыве водородной бомбы. В СССР в 1953 г. появилась машина "Стрела" со скоростью 2000 операций в секунду и памятью 1024 слов. Было выпущено 7 таких машин, на которых в ВЦ АН СССР, ИПМ АН СССР, МГУ им. М. В. Ломоносова и в вычислительных центрах некоторых министерств было проведено множество научных расчетов. Крупным событием стало появление в СССР в 1967 г. машины БЭСМ-6 с быстродействием до 1 млн оп./сек и памятью
32 Кб. В этой ЭВМ было реализовано множество революционных архитектурных особенностей, в том числе, конвейер команд, параллельная и асинхронная работа процессора и модулей оперативной памяти, ассоциативная кеш-память. Механизмы прерывания, защиты памяти, преобразования виртуальных адресов в физические и привилегированный режим работы для ОС позволили использовать БЭСМ-6 в мультипрограммном режиме и режиме разделения времени.
Концепция параллелизма появляется в архитектуре компьютеров сначала на уровне разрядно-параллельной памяти и арифметики в компьютерах IBM 701 в 1953 г., затем вводятся независимые процессоры ввода/вывода (IBM 709, 1958 г.) и расслоение памяти на несколько банков (IBM Stretch, 1961 г.), которые могут обслуживаться одновременно. Конвейерная обработка команд была впервые использована в 1963 г. в машине ATLAS, разработанной в Манчестерском университете. Данный компьютер оказал огромное влияние на архитектуру ЭВМ и ОС: в нем впервые использована мультипрограммная ОС, основанная на использовании виртуальной памяти и системы прерываний. В 1964 г. фирма CDC выпускает CDC-6600 - первый компьютер, в котором использовалось несколько независимых функциональных устройств. За ним в 1969 г. последовал CDC-7600, в котором 8 независимых функциональных устройств могли осуществлять конвейерную обработку. В 1972 г. Сеймур Крэй покидает CDC и основывает компанию Cray Research, которая в 1976 г. выпускает первый векторно-конвейерный компьютер CRAY-1. Для данного типа параллельных архитектур характерно наличие векторных команд, работающих с целыми массивами независимых данных и позволяющих эффективно использовать конвейерные функциональные устройства. В то же время (1967 - 1974 гг.) разрабатывается компьютер ILLIAC-IV, в котором имелось 64 вычислительных элемента и центральный управляющий процессор, и была реализована концепция SIMD (Single Instruction Multiple Data). Вычислительные элементы (процессор и оперативная память) были соединены сетью в виде матрицы и работали в синхронном режиме, одновременно выполняя одни и те же инструкции, получаемые от центрального процессора. Похожая идея организации вычислительных элементов в двух- или трехмерную матрицу была положена в основу транспьютерной техники. В 1985 г. фирма INMOS начала выпуск транспьютеров, объединявших в одной микросхеме микропроцессор, коммуникационную систему и оперативную память. Из транспьютеров собирались относительно дешевые, но мощные, надежные и масштабируемые суперкомпьютеры. В отличие от матричных процессоров, в транспьютерных системах была реализована модель MIMD. В 1987 г. фирмой Sun была разработана рабочая станция Sun-4 на основе процессоров с масштабируемой архитектурой SPARC. Данная архитектура позволяет строить симметричные мультипроцессоры с общей памятью (SMP). В настоящее время уже многие персональные компьютеры имеют несколько процессоров, соединенных общей шиной данных, и работающих с общей памятью.
С развитием возможностей коммуникационного оборудования и сетей передачи данных большее применение находят массивно-параллельные компьютеры (МРР) и компьютеры с кластерной архитектурой. Уже упоминавшиеся транспьютерные системы были родоначальниками массивно-параллельных архитектур с использованием матричной топологии сети связи. Однако, наиболее оптимальной топологией является n-мерный куб, в котором максимальная длина пути доставки сообщения от одного процессора другому и соответственно время доставки пропорциональны логарифму количества процессоров. Впервые такая система связей была использована в компьютере Caltech Hypercube, созданном в 1984 г. в Калифорнийском технологическом университете. По всей видимости, первой анонсировала концепцию кластерной системы компания DEC в 1983 г. на основе VAX-компьютеров. Однако, еще в начале 80-х годов Билл Джой (один из основателей компании Sun) сказал фразу, впоследствии ставшую лозунгом Sun: "Сеть - это компьютер". Подтверждением этому становится доминирование в настоящее время кластерных гибридных архитектур с использованием высокоскоростных коммуникационных сетей, таких как Fast Ethernet, Myrinet, Paramnet итп. Так, наиболее мощный на сегодняшний день суперкомпьютер NEC Earth-Simulator (по данным списка Тор 500 [95] на ноябрь 2002 г.) представляет собой кластер, составленный из 640 узлов, каждый из которых - 8-процессорный мультипроцессор с общей памятью на основе векторных процессоров. Все основные типы параллельных архитектур: векторно-конвейерные процессоры, SMP и МРР/кластеры находят достойное применение в данном суперкомпьютере-победителе.
Для каждого типа архитектур параллельных компьютеров характерна своя модель программирования. Так, для SIMD-архитектур используется модель data parallel, для векторно-конвейерных процессоров характерно наличие векторных инструкций, осуществляющих операции одновременно над целым массивом регистров. В мультипроцессорах с общей памятью (SMP) вычисления обычно производятся в нескольких нитях или процессах под управлением одной операционной системы (ОС) с явной (или неявной) синхронизацией доступа к разделяемой памяти. Кроме стандарта POSIX программного интерфейса ОС здесь имеется также специализированный для SMP-компьютеров стандарт ОрепМР. Для МРР и компьютеров с кластерной архитектурой стандартной моделью является программирование с использованием сообщений. Впервые сообщения как средство межпроцессного взаимодействия ввел Хоар в 1978 г. в работе [96]. Эта модель является настолько удобной и универсальной, что может с успехом употребляться на всех типах MIMD-архитектур. На основе сообщений было создано много систем параллельного программирования, например, Parallel Virtual Machine [97], High Performance Fortran [98], P4 [99], Express [100] и, конечно же, стандарт MPI для работы с сообщениями [157,101,102]. В настоящее время большинство параллельных программ реализуется с использованием тех или иных реализаций МРІ. Более высокоуровневые параллельные языки и программные среды, такие как трС [103], DVM [104], также часто реализуются с использованием МРІ, добавляя свой слой программного обеспечения, реализующего параллельные конструкции языка. Сборник [105] освещает основные парадигмы и модели программирования для различных классов параллельных ЭВМ.
Параллельные вычисления, как это ни странно, были открыты еще в докомпьютерную эпоху. Пионером в параллельной обработке потоков данных был академик Самарский, выполнявший в начале 50-х годов расчеты, необходимые для моделирования ядерных взрывов. Непосредственно расчеты совершали несколько десятков человек с арифмометрами, которые передавали данные друг другу просто на словах. Таким образом, в частности, была рассчитана эволюция взрывной волны. Это, можно сказать, и была первая параллельная система, в которой к тому же была реализована модель передачи сообщений.
С параллельным программированием связано много важных как теоретических, так и чисто практических вопросов. Как практическим, так и некоторым теоретическим аспектам программирования для параллельных компьютеров уделяется внимание в обзоре [106]. В центре же этой области науки находится построение и анализ алгоритмов распараллеливания. Как уже было сказано, не любой последовательный алгоритм можно распараллелить, часто необходимо строить новый параллельный алгоритм, эквивалентный существующему. Анализу алгоритмов на предмет их возможного распараллеливания посвящены работы Воеводина [8,107,108]. Параллельные алгоритмы для решения задач линейной алгебры наиболее широко освещены в книге Ортега [109]. В книге Фримана и Филипса [ПО] также основное внимание уделяется параллельным методам решения задач линейной алгебры, а также рассматриваются параллельные методы решения обыкновенных дифференциальных уравнений и задач оптимизации. Параллельным методам решения алгебраических задач посвящены также книги Саада [111], Донгарры [112]. Данные методы также используются и при решении уравнений в частных производных, поскольку при их решении задача часто сводится к линейным разностным уравнениям. Параллельные алгоритмы в приложении к задачам численной гидродинамики рассматривались Ван дер Ворстом в сборнике [113]. Основной метод распараллеливания, используемый для задач математической физики и, в частности, численной гидродинамики - декомпозиция расчетной области [114,115] на подобласти, обрабатываемые различными процессорами. Такой метод удобен при распараллеливании программ для компьютеров с распределенной памятью, поскольку каждый процессор может производить вычисления в своей подобласти, находящейся в своей локальной памяти. Параллельные алгоритмы в применении к задачам численной гидродинамики продолжают рассматриваться в литературе и в последние годы [116,117,118].
Для решения задач математической физики на параллельных компьютерах разработано большое число библиотек и программных комплексов. Часть из них, такие, как пакеты UNPACK [121], ScaLAPACK [122], NAG Parallel Library [123], PINEAPL [124], нацелены на решение систем линейных уравнений различного типа, в том числе, с ленточными и разреженными матрицами, получающимися при дискретизации дифференциальных уравнений в разностные. Другая часть пакетов направлена сугубо на численное решение уравнений в частных производных. Наиболее популярными математическими программами являются Mathematica [125] и Maple [126], однако в них можно решать ограниченное число типов уравнений, так, например, в Maple стандартный пакет DETools умеет решать только квазилинейные уравнения в частных производных первого порядка с двумя переменными. Более универсальными и гибкими являются: пакет PDELab [127] для численного исследования систем нелинейных дифференциальных уравнений в частных производных, а также, FlexPDE [128] - программная система для получения численных решений систем дифференциальных уравнений, встречающихся при решении задач во многих областях фундаментальной и прикладной науки. Крупнейший мировой архив Netlib [129] свободно распространяемых программных библиотек (многие из которых предназначены для параллельных компьютеров) для разнообразных математических исследований, в том числе, для численного решения дифференциальных уравнений в частных производных, поддерживается AT&T Bell Laboratories.
Следующий класс программ и пакетов предназначен для численного моделирования в определенных прикладных областях. Это довольно обширный класс, для примера можно привести такие программы, как FLOW-3D [130] -программа для численного моделирования в гидро- и газодинамике; CARDINAL [131] - программа для моделирования динамики различных водных объектов, оценка загрязнений воды; SME [132] - моделирование в гидрологии; Argus [133] -программная среда для численного моделирования в геологии и геофизике; пакет программ KIPARIS [134], разрабатываемый в ИММ РАН, для численного моделирования сверхзвуковых течений вязкого газа в 2-D и 3-D геометриях на многопроцессорных транспьютерных системах; GDT [135] - пакет для моделирования газодинамических задач, включающий различные модели горения, в сложных конфигурациях и геометриях, задаваемых пользователем. Программы этого класса имеют ограниченное число реализованных в них моделей и ограниченное число параметров, доступных для задания пользователю.
Промежуточное положение между программами решения дифференциальных уравнений в частных производных общего типа и пакетами для моделирования в определенной прикладной области занимают так называемые каркасы приложений [136]. Каркасы приложений предоставляют некоторую поддержку для решения задач по определенному шаблону, но реализация основного ядра программы остается за пользователем. В качестве каркасов приложений можно рассматривать следующие пакеты: CLAWPACK [137] для расчета законов сохранения и гиперболических систем, разработанный Ле Веком; GEMS, Goddard Earth Modeling System [139] для моделирования климата Земли, распространения загрязнений; Flux Coupler [138] для вычисления потоков в различных компонентных моделях климатической системы Земли; FMS, Flexible Modeling System [140] для моделирования атмосферных явлений, например, ураганов, а также химических взаимодействий в атмосфере.
Особенно интересны пакеты для параллельных компьютеров. Пакеты DDB, Distributed Data Broker [141] и PAWS, Parallel Application Workspace [142] предоставляют средства для организации параллельных приложений в компонентно-ориентированную инфраструктуру, работающую в распределенной среде. Пакеты DAGH, Distributed Adaptive Grid Hierarchy [143] и SAMRAI, Structured Adaptive Mesh Refinement Applications Infrastructure [144] являются каркасами для реализации параллельных вычислений на адаптивных вложенных сетках. Но более всего указанной категории программных продуктов соответствуют три развитых каркаса параллельных приложений для численного моделирования: CACTUS [145] - каркас приложений для численного решения уравнений в частных производных, он более ориентирован на приложения в области астрофизики и релятивистской физики; РООМА, Parallel Object-Oriented Methods and Applications [146] - большая библиотека классов, реализующих различные абстракции уравнений в частных производных и параллельных вычислений вообще; PETSC, Portable Extensible Toolkit for Scientific Computing [147] - расширяемый каркас для решения уравнений в частных производных с использованием параллельных вычислений. Отличный обзор имеющихся каркасов и программ для численного моделирования на параллельных компьютерах был проведен при предпроектном обследовании группой Earth System Modeling Framework TASC в рамках NASA Computational Technology Project [148]. Сравнительный обзор каркасов и инструментальных средств для разработки параллельных приложений в области численного моделирования был также представлен ДеЛюка и Деннисом на конференции NCAR по информационным технологиям [149].
В данной диссертации рассматриваются вопросы, связанные с решением задач численного моделирования с системами уравнений гиперболического типа. Процесс построения программной модели может быть символически описан как нахождение точки в трехмерном пространстве, координатами которого являются три оси, которые мы сейчас кратко опишем. Как уже указывалось, гиперболические уравнения встречаются в различных прикладных областях. Поэтому, первое измерение, в котором мы перемещаемся - это измерение физических моделей и конкретных задач моделирования. Здесь изменяются: основные расчетные компоненты, конкретный вид системы уравнений, массовые силы, матрицы Якоби итп.
Многие методы численного решения этих уравнений могут быть использованы независимо от конкретного вида системы, причем различные методы обычно комбинируются для совместного использования. Некоторые из этих методов ортогональны в том смысле, что они являются независимыми и могут быть заменены другими методами данного типа независимо от других частей результирующего метода. Это многообразие численных методов составляет второе измерение, в котором мы находим результирующий численный метод для решения нашей задачи.
Последнее, что мы должны выбрать на третьей оси - это алгоритм и программная реализация выбранного метода. Для того чтобы получать наиболее точные результаты, нужно использовать параллельные вычисления на многопроцессорных компьютерах. Следует заметить, что параллельные и последовательные алгоритмы для одного и того же метода могут сильно отличаться. Существует разнообразные программные пакеты и библиотеки методов, но многие исследователи пишут свои собственные программы - в основном, потому что большинство этих пакетов не настолько гибки, как того бы хотелось исследователям, пакеты не являются переносимыми, они поддерживают ограниченное количество методов, их очень тяжело интегрировать с другими программами, и многие из них не имеют параллельных версий.
Для облегчения разработки программ для численного моделирования в различных прикладных областях, с использованием различных численных методов и существующего программного обеспечения предлагается использовать унифицированный каркас приложений, который можно представлять себе в виде общего ядра, встречающегося во всех программах для численного моделирования, лежащих в описанном пространстве моделей. Программный продукт должен представлять собой универсальный, гибкий, расширяемый, открытый, ориентированный на стандарты каркас для построения солверов для параллельных компьютеров. Каждый альтернативный алгоритм (описывающий метод, модель и т.п.) должен быть представлен определенным строительным блоком, встраиваемым в каркас. Наиболее полезные и часто используемые блоки должны предоставляться наряду с каркасом в том же программном продукте, но каркас должен иметь открытый интерфейс для пользователей, которые могли бы подключать через данный интерфейс другие физические модели, методы, солверы, визуализаторы и другие компоненты. Для построения данного каркаса нужно базироваться на универсальной методике, которая также должна быть гибкой и расширяемой. В диссертации описывается такая методика, представляется ее реализация в программном пакете и рассматриваются результаты приложения данного пакета к численному моделированию нескольких актуальных задач газодинамики атмосферных потоков и сейсмического моделирования.
Основными целями диссертации являются: разработка универсальной технологии параллельных вычислений для численного моделирования задач, описываемых системами уравнений гиперболического типа, разработка программного пакета, позволяющего упростить создание параллельных программ для численного моделирования, на основе разработанной технологии, проведение численного моделирования серии задач, связанных с конвективными течениями в стратифицированной атмосфере, вызванными пожарами или другими высокоэнергетическими источниками, проведение серии численных экспериментов, заключающихся в решении прямой задачи геофизики.
Диссертация состоит из введения, четырех глав, заключения и списка литературы.
Во введении обосновывается актуальность темы диссертации, дается обзор работ по численному моделированию, численным методам, обзор существующих программных пакетов для численного моделирования. Затем кратко описывается постановка задачи, формулируются цели диссертации и приводится краткое содержание глав диссертации.
В первой главе, состоящей из пяти разделов, излагается универсальная технология численного моделирования. В разделе 1.1 излагается постановка задачи. Далее, в разделе 1.2 строится формализация задачи и описываются расчетные методики. Использованные конечно-разностные схемы рассматриваются в разделе 1.3. В разделе 1.4 строится обобщение монотонной гибридной схемы второго порядка [87] для неравномерных сеток. В разделе 1.5 описывается применение метода крупных частиц и приводится трехмерный аналог неявной схемы бегущего счета [156] для эйлерова этапа метода крупных частиц для уравнений Эйлера с силой тяжести.
Во второй главе, состоящей из пяти разделов, описывается программный пакет, построенный на основании методики, описанной в первой главе. Раздел 2.1 посвящен описанию программной реализации описываемой методики на параллельных компьютерах. В разделе 2.2 приводится объектная модель программного пакета. В разделе 2.3 описывается метод пространственной декомпозиции задач на параллельных компьютерах. Далее, в разделе 2.4 даются оценки эффективности распараллеливания с применением пространственной декомпозиции. В разделе 2.5 приводятся тестовые расчеты с использованием указанного пакета.
В третьей главе, состоящей из пяти разделов, проводится численное исследование пространственных течений в атмосфере, вызванных крупномасштабными пожарами или другими высокоэнергетическими источниками. Во введении к третьей главе (раздел 3.1) проводится обзор литературы и исследований других авторов на данную тему. Математическая модель, использованная во всех численных экспериментах третьей главы, описывается в разделе 3.2. В разделе 3.3 приводятся результаты тестовых расчетов моделирования атмосферных течений. Раздел 3.4 посвящен моделированию пожаров. В разделе 3.5 проводится моделирование течений с вихревой структурой: огненного смерча и торнадо.
В четвертой главе, состоящей из четырех разделов, проводится численное исследование прямой задачи геофизики, заключающейся в изучении характера распространения акустических волн в неоднородной земной коре. В разделе 4.1 проводится обзор прикладной области (сейсмического моделирования). Далее, в разделе 4.2 описывается математическая модель на основе уравнений упругости. Тестовые расчеты: распространение сигнала в однородной среде и акустическое приближение приводятся в разделе 4.3. Основные результаты моделирования прямой задачи геофизики в неоднородной среде представлены в разделе 4.4.
В заключении сформулированы основные результаты работы и выводы.
Результаты диссертации опубликованы в научных журналах, сборниках и трудах конференций [194-203].
Результаты работы докладывались и обсуждались на конференциях и семинарах: VII Российско-Японский симпозиум "Вычислительная Гидродинамика" (г. Москва, МГУ, 31 июля - 6 августа 2000 г.);
IV международная конференция "Лесные и степные пожары: возникновение, распространение, тушение и экологические последствия" (г. Иркутск, 25-29 сентября 2001 г.);
Российско-Японский симпозиум "Актуальные Проблемы Вычислительной Механики" (г. Санкт-Петербург, 5-10 августа 2002 г.); V международный конгресс "Математическое Моделирование" (г. Дубна, ОИЯИ, 30 сентября - 6 октября 2002 г.); XXXXV научная конференция МФТИ "Современные проблемы фундаментальных и прикладных наук" (г. Долгопрудный, МФТИ, 29-30 ноября 2002 г.); международная школа-семинар "GO++ winter school on numerical methods for HJ/HJB problems" (Франция, г. Роканкур, INRIA, 9-12 декабря 2002 г.); научный семинар ИАП РАН "Математическое моделирование в механике сплошных сред и физике плазмы".
Опыт автора, приобретенный в ходе работы над диссертацией, нашел отражение в проводимых автором курсах лекций "Теоретические основы конвейерных и параллельных операционных вычислений в математической физике" для студентов кафедры "Математические и информационные технологии" Московского Физико-Технического Института и курсе "Программирование распределенных систем на основе Java-технологий", проводимом Московским Центром Непрерывного Математического Образования по программе дополнительного профессионального образования.
Автор выражает глубокую благодарность за постоянное внимание, поддержку и совместную научную работу своим научным руководителям, предложившим также и тему диссертации: академику Олегу Михайловичу Белоцерковскому и к.ф.-м.н. Алексею Михайловичу Опарину, к.ф.-м.н. Андрею Викторовичу Конюхову за ценные обсуждения и полезные советы, н.с. Максиму Николаевичу Антоненко за тесное научное сотрудничество, н.с. Марине Сергеевне Белоцерковской за помощь в подготовке статей, а также всему коллективу ИАП РАН за доброжелательное отношение и содействие в работе.
Обобщение монотонной гибридной схемы второго порядка для неравномерных сеток
С параллельным программированием связано много важных как теоретических, так и чисто практических вопросов. Как практическим, так и некоторым теоретическим аспектам программирования для параллельных компьютеров уделяется внимание в обзоре [106]. В центре же этой области науки находится построение и анализ алгоритмов распараллеливания. Как уже было сказано, не любой последовательный алгоритм можно распараллелить, часто необходимо строить новый параллельный алгоритм, эквивалентный существующему. Анализу алгоритмов на предмет их возможного распараллеливания посвящены работы Воеводина [8,107,108]. Параллельные алгоритмы для решения задач линейной алгебры наиболее широко освещены в книге Ортега [109]. В книге Фримана и Филипса [ПО] также основное внимание уделяется параллельным методам решения задач линейной алгебры, а также рассматриваются параллельные методы решения обыкновенных дифференциальных уравнений и задач оптимизации. Параллельным методам решения алгебраических задач посвящены также книги Саада [111], Донгарры [112]. Данные методы также используются и при решении уравнений в частных производных, поскольку при их решении задача часто сводится к линейным разностным уравнениям. Параллельные алгоритмы в приложении к задачам численной гидродинамики рассматривались Ван дер Ворстом в сборнике [113]. Основной метод распараллеливания, используемый для задач математической физики и, в частности, численной гидродинамики - декомпозиция расчетной области [114,115] на подобласти, обрабатываемые различными процессорами. Такой метод удобен при распараллеливании программ для компьютеров с распределенной памятью, поскольку каждый процессор может производить вычисления в своей подобласти, находящейся в своей локальной памяти. Параллельные алгоритмы в применении к задачам численной гидродинамики продолжают рассматриваться в литературе и в последние годы [116,117,118].
Для решения задач математической физики на параллельных компьютерах разработано большое число библиотек и программных комплексов. Часть из них, такие, как пакеты UNPACK [121], ScaLAPACK [122], NAG Parallel Library [123], PINEAPL [124], нацелены на решение систем линейных уравнений различного типа, в том числе, с ленточными и разреженными матрицами, получающимися при дискретизации дифференциальных уравнений в разностные. Другая часть пакетов направлена сугубо на численное решение уравнений в частных производных. Наиболее популярными математическими программами являются Mathematica [125] и Maple [126], однако в них можно решать ограниченное число типов уравнений, так, например, в Maple стандартный пакет DETools умеет решать только квазилинейные уравнения в частных производных первого порядка с двумя переменными. Более универсальными и гибкими являются: пакет PDELab [127] для численного исследования систем нелинейных дифференциальных уравнений в частных производных, а также, FlexPDE [128] - программная система для получения численных решений систем дифференциальных уравнений, встречающихся при решении задач во многих областях фундаментальной и прикладной науки. Крупнейший мировой архив Netlib [129] свободно распространяемых программных библиотек (многие из которых предназначены для параллельных компьютеров) для разнообразных математических исследований, в том числе, для численного решения дифференциальных уравнений в частных производных, поддерживается AT&T Bell Laboratories.
Следующий класс программ и пакетов предназначен для численного моделирования в определенных прикладных областях. Это довольно обширный класс, для примера можно привести такие программы, как FLOW-3D [130] -программа для численного моделирования в гидро- и газодинамике; CARDINAL [131] - программа для моделирования динамики различных водных объектов, оценка загрязнений воды; SME [132] - моделирование в гидрологии; Argus [133] -программная среда для численного моделирования в геологии и геофизике; пакет программ KIPARIS [134], разрабатываемый в ИММ РАН, для численного моделирования сверхзвуковых течений вязкого газа в 2-D и 3-D геометриях на многопроцессорных транспьютерных системах; GDT [135] - пакет для моделирования газодинамических задач, включающий различные модели горения, в сложных конфигурациях и геометриях, задаваемых пользователем. Программы этого класса имеют ограниченное число реализованных в них моделей и ограниченное число параметров, доступных для задания пользователю.
Промежуточное положение между программами решения дифференциальных уравнений в частных производных общего типа и пакетами для моделирования в определенной прикладной области занимают так называемые каркасы приложений [136]. Каркасы приложений предоставляют некоторую поддержку для решения задач по определенному шаблону, но реализация основного ядра программы остается за пользователем. В качестве каркасов приложений можно рассматривать следующие пакеты: CLAWPACK [137] для расчета законов сохранения и гиперболических систем, разработанный Ле Веком; GEMS, Goddard Earth Modeling System [139] для моделирования климата Земли, распространения загрязнений; Flux Coupler [138] для вычисления потоков в различных компонентных моделях климатической системы Земли; FMS, Flexible Modeling System [140] для моделирования атмосферных явлений, например, ураганов, а также химических взаимодействий в атмосфере.
Особенно интересны пакеты для параллельных компьютеров. Пакеты DDB, Distributed Data Broker [141] и PAWS, Parallel Application Workspace [142] предоставляют средства для организации параллельных приложений в компонентно-ориентированную инфраструктуру, работающую в распределенной среде. Пакеты DAGH, Distributed Adaptive Grid Hierarchy [143] и SAMRAI, Structured Adaptive Mesh Refinement Applications Infrastructure [144] являются каркасами для реализации параллельных вычислений на адаптивных вложенных сетках. Но более всего указанной категории программных продуктов соответствуют три развитых каркаса параллельных приложений для численного моделирования: CACTUS [145] - каркас приложений для численного решения уравнений в частных производных, он более ориентирован на приложения в области астрофизики и релятивистской физики; РООМА, Parallel Object-Oriented Methods and Applications [146] - большая библиотека классов, реализующих различные абстракции уравнений в частных производных и параллельных вычислений вообще; PETSC, Portable Extensible Toolkit for Scientific Computing [147] - расширяемый каркас для решения уравнений в частных производных с использованием параллельных вычислений. Отличный обзор имеющихся каркасов и программ для численного моделирования на параллельных компьютерах был проведен при предпроектном обследовании группой Earth System Modeling Framework TASC в рамках NASA Computational Technology Project [148]. Сравнительный обзор каркасов и инструментальных средств для разработки параллельных приложений в области численного моделирования был также представлен ДеЛюка и Деннисом на конференции NCAR по информационным технологиям [149].
Оценки эффективности при пространственной декомпозиции задач на параллельных компьютерах
С другой стороны, некоторые солверы имеют параметры, отсутствующие в других солверах. Так, например, функция-ограничитель применяется только для приближенных Римановских методов. При применении объектного подхода это обстоятельство легко учитывается и не нарушает логической структуры системы. Ключевым моментом является идея рассматривать солверы как абстрактные типы данных [159] и, соответственно этому, построить дерево наследования классов, реализующих всевозможные солверы. Используя абстрактные типы данных, можно одинаково рассматривать и объединить в рамках одной программы даже такие разные методы как метод крупных частиц и какой-либо из сеточно-характеристических методов, например, описанный в разделе 1.4. Несмотря на то, что данные методы имеют разные алгоритмы и совершенно разные наборы параметров, они являются реализациями одного абстрактного типа данных представленного классом HMethod, и для других частей программы имеют один и тот же интерфейс.
На рисунке 1 представлено дерево наследования солверов, реализующих расчетные методики. Наследование можно понимать как в концептуальном смысле, так и в чисто техническом - как наследование [159] классов, представляющих солверы на объектном языке программирования. Большое число показанных методов составляют приближенные Римановские солверы с ограничителями различного вида, сильно влияющими на свойства метода. Кроме этих схем часто используются гибридные схемы, обеспечивающие высокий порядок аппроксимации и по пространству, и по времени. Отдельные ветви составляют такие методы, как метод крупных частиц, метод частиц в ячейках, статистические методы и другие методы, направленные на решение уравнений определенного вида.
Структура всего пакета показана на рисунке 2 Представлена диаграмма классов на визуальном языке моделирования UML [160]. Главные классы пакета -HApplication, HProblem, HMethod и HState. HApplication отвечает за инициализацию/деинициализацию пакета и предоставляет метод calcLoopO -входную точку для выполнения вычислений. HProblem описывает физическую постановку задачи. Класс HMethod представляет численный метод для решения задачи моделирования. Сами данные хранятся в классе HState, агрегирующим поля компонент (объекты класса HField), представляющих решение в настоящий и предыдущий моменты времени, а также параметры, определяющие текущий момент времени (класс HMoment).
Класс HProblem, описывающий физическую постановку задачи, включает следующие компоненты: координаты, размер сетки, количество расчетных компонент, параметры декомпозиции, начальные условия, граничные условия, вид уравнений, правая часть уравнений. Координаты описываются базовым классом HCoord и бывают двух типов: генерируемые - наследники класса HGeneratedCoord, и вычисляемые - наследники класса HCalculatedCoord. Начальные условия чаще всего задаются пользователем, который для этого наследует абстрактный класс HInitialValues. Правую часть пользователь либо задает таким же образом сам, наследуя класс HRightHandTerms, либо использует предоставленные пакетом классы: HGravity, определяющий силу тяжести Земли, Н Viscosity, учитывающий вязкие члены в уравнении Навье-Стокса, и HHeatTransfer, задающий теплопроводность. Кроме того, можно комбинировать правые части, используя класс HCombinedRightHand, задавая любое количество слагаемых. Вид уравнений задается через матрицы Якоби частных производных потоков по расчетным переменным. Для этого предоставляются три объекта класса HJacobiMatrix, вычисляющие матрицы Якоби потоков по трем направлениям, их собственные числа и собственные вектора. Для постановки граничных условий используются две иерархии классов. Классы, наследованные от HBoundaryValues, являются главными для задания граничных условий. В этих классах для всей грани или для каждой ее точки отдельно (что регулируется с помощью методов selectXXXBndType()) задается тип граничного условия, который должен быть представлен объектом, наследованным от HBoundaryType. Пакетом предоставляются стандартные типы граничных условий: непротекание (класс HWallBoundary), прилипание (класс HAdhesionBoundary), снос (класс HShiftBoundary), свободная граница (класс HFreeBoimdary), но пользователь может и сам задавать новые типы граничных условий, наследуя абстрактный класс HBoundaryType.
Диаграмма классов, использующихся для реализации численного метода, представлена отдельно от других классов на рисунке 3. Базовым классом для всех численных методов является абстрактный класс HMethod, имеющий виртуальный метод calcStep(), вызываемый на каждой временной итерации программы. Повышение порядка схемы, представленной объектом originalMethod, с помощью метода Рунге-Кутта осуществляется замещением объекта originalMethod в главном объекте приложения HApplication новым объектом класса HRungeKuttaMethod, построенным на основе originalMethod. В методе calcStepO классов, наследующих от HRungeKuttaMethod, должен быть закодирован лишь определенный метод Рунге-Кутта с вызовом originalMethod.calcStepO с нужными параметрами (так как схема в данном случае получается многошаговая, в данном классе требуется введение дополнительных полей данных, память под которые выделяется статически или в момент первого обращения на все время работы программы). В пакете реализовано несколько методов использующих расщепление по пространственным переменным, для всех них используется единственная реализация класса HMethod - класс HSpatialVariablesSplittingMethod. Отдельно от этих методов стоит метод крупных частиц, использующий совершенно другой тип расщепления: на эйлеров и лагранжев шаг. Данный метод представлен классом HLargeParticlesMethod, агрегирующим в себе два указанных метода для расчета на каждом этапе, описываемых классами HEulerStep и HLagrangianStep. Для решения уравнений переноса на лагранжевом этапе может быть применен основной метод с расщеплением по пространственным переменным со скалярными матрицами Якоби. Для этого используется адаптер HHyperbolicLagrangianStep, реализующий интерфейс HLagrangianStep и делегирующий выполнение расчета к HSpatialVariablesSplittingMethod.
Несмотря на обилие методов, которые можно реализовать с помощью данного пакета, во всех них используется один и тот же класс HSpatialVariablesSplittingMethod. Новые численные методы реализуются путем задания необходимых новых алгоритмических параметров для класса HSpatialVariablesSplittingMethod. Среди них: тип расщепления по пространственным переменным, осреднение, аппроксимация, расщепление потоков, энтропийная коррекция, конечно-разностный метод для решения уравнения переноса. Каждый из этих параметров определяется как объект определенного абстрактного класса. Так, осреднение, аппроксимация и энтропийная коррекция определяются абстрактными классами HAveraging, HSpline и HEntropyFix, у которых есть стандартные реализации, например, для осреднения расчетных компонент чаще всего используются: HMeanAverage (для выбора среднего арифметического компонент на границах ячеек) и HRoeAverage (при использовании осреднения по Роу).
Моделирование пожаров
В задаче о крупномасштабном пожаре [173,194] сам пожар моделируется источником энергии и примеси. Источник горения - поверхностный, представляет собой круг с радиусом порядка нескольких километров, с мощностью выделения энергии Q = 0.05 МВт/м2. В модели вместо данного источника использовался объемный источник, с высотой, равной одной ячейке, соответственно, объемная мощность составляла Q/h0. Облако горения визуализируется концентрацией легковесной примеси, извергаемой в источнике. Источник работает в течение часа, затем выключается.
Вначале приведем результаты двумерного моделирования задачи с источником радиуса R = 10 км. В трехмерном случае это соответствует горению протяженной полосы шириной 2R, что может являться модельной постановкой для горения протяженного лесного массива. В двумерном же случае можно взять подробную сетку и получить качественные результаты за приемлемое время. Так, в представленных ниже двумерных расчетах сетка составляет 400 х 150 точек, и на 16 процессорах многопроцессорного комплекса МВС-1000/17 IAP расчет данной задачи можно проводить в реальном масштабе времени. На данной задаче также проводилось экспериментальное исследование сходимости решения. Один и тот же расчет выполнялся на нескольких сетках разного размера. Было получено, что характерные крупномасштабные структуры не зависят от сетки, а при увеличении подробности уточняются более высокочастотные характеристики этих структур. На рисунке 20 показано распределение облака примеси в моменты t = 5, 10, 20, 40 и 60 минут. Ввиду того, что задача - двумерная, облако поднимается выше и распространяется дальше, чем в трехмерном случае. Для двумерного расчета характерная высота растекания - 10 км.
Также следует отметить проявляющуюся неустойчивость Релея-Тейлора [6], возникающую на границе двух сред, когда градиент плотности и градиент давления направлены в разные стороны. В данном случае, сила тяжести, обеспечивающая градиент давления направлена вниз, а более теплый и, следовательно, более легкий газ находится на поверхности Земли. Характерные струи горячего воздуха сперва поднимаются вверх, затем охлаждаются и стекают на некоторое расстояние вниз, образуя грибообразные пузыри. Заметно явление укрупнения масштаба пузырей с течением времени. В конце концов остается один крупный пузырь, который, во-первых, начинает растекаться на определенной высоте, а во-вторых, из него выдвигается на определенное расстояние вверх конвективная колонка. На дальних моментах времени наблюдается колебание облака относительно средней высоты растекания с частотой Вяйсяля-Брента минуты.
В трехмерном расчете размер сетки составляет 128 х 128 х 64. Сетка увеличивается по мере удаления от центра и от поверхности Земли. Размер одной ячейки в центре области составляет порядка 300 м. Боковые границы отнесены на расстояние в 30 км, верхняя граница - 24 км. При радиусе источника 10 км общая мощность энерговыделения источника q = 1.57ХІ07 МВт. На рисунках 21-28 представлено положение облака примеси на моменты времени t = 5, 10, 15, 20, 25, 40, 60, 90 минут. В начальные моменты времени на краю источника в вертикальной плоскости образуется вихрь, приводящий к началу подъема облака по периметру источника. Вихрь хорошо виден на рисунке 29, где в вертикальном сечении, проходящем через центр источника, показано распределение примеси и поле скоростей на момент времени 5 минут. Затем, в силу того, что выделяется большая энергия, течение быстро становится неустойчивым, а в области над очагом пожара наблюдается интенсивное перемешивание горячего воздуха, приходящего из источника, и холодного атмосферного воздуха. К 15 минутам в центре облака образуется конвективная колонка, в которой течение имеет направленное движение вверх. Облако достигает стратосферы, однако основная масса облака, сосредоточенная у Земли, продолжает подниматься и перемешиваться в процессе подъема.
К моменту времени 25-30 минут облако выходит на характерные высоты растекания, сосредоточенные в диапазоне от 3.5 до 9 км. При дальнейшем распространении облако развивается неравномерно по всем направлениям в горизонтальной плоскости, что обусловлено определенным типом неустойчивости, нарушающей осевую симметрию. Наличие выделенных направлений декартовой сетки определяет эти направления развития языков облака - вдоль лучей, параллельных координатным направлениям, и лучей, идущих под 45 к ним. При этом, языки чередуют несколько выделенных высот растекания: 4.5, 6.5 и 8.5 км в зависимости от направления. Наличие нескольких языков облака было замечено также и в осесимметричной постановке в работе [173]. Асимметрия также проявляется и в ножке пожара, имеющий на порядок меньший радиус, чем радиус очага. Ножка устойчиво существует вплоть до выключения источника в момент времени 1 час.
На рисунке 30 приведены графики зависимости высоты подъема конвективной колонки, высоты растекания и радиуса облака от времени. Высота растекания определяется как высота с максимальным интегральным содержанием примеси. При определении высоты подъема конвективной колонки и радиуса облака также определяются интегральные характеристики плотности примеси (на данной высоте или с данным радиусом). Поскольку примесь непрерывно распределена в пространстве, то граница облака не является четкой. В качестве характерного уровня берется точка, в которой вторая производная максимальна. Высота подъема конвективной колонки, определяемая таким образом, составляет на момент времени 1.5 часа порядка 13.5 км, что соответствует результатам статьи [173]. Данный расчет подтверждает результаты [173], в которых утверждается, что для источников с радиусом R 8 км высота подъема отличается (в меньшую сторону) от эмпирической формулы [186] для точечного источника z = AqVA, где q - общая мощность энерговыделения, измеряемая в МВт, а г - высота подъема в км, приА=0.31 [187].
Моделирование прямой задачи геофизики в неоднородной среде
Результаты расчетов обоих вариантов программы (с уравнениями упругости и уравнениями акустики) совпадают. Во всех задачах размер сетки - 800 х 800 точек, размеры расчетной области - 1км х 1км, и добавлены дополнительные точки слева, справа и внизу, чтобы избежать отражения волны от этих границ. Параметры основной среды: плотность р - 2200 кг/м2, скорость звука Vp = 2000 м/сек. Точечный источник имеет координаты (300м, 200м). Сигнал в источнике представляет собой вторую производную функции Гаусса, локализованной во времени т = 1/30 сек, т.е. с частотой v = 30 Гц:
Результаты тестового расчета распространения акустических волн с отражением от поверхности Земли представлены на рисунке 48 в моменты времени 0.1 и 0.2 сек. Цветом визуализируется распределение давления. Белый цвет соответствует положительным значениям, черный - отрицательным. Поскольку в сигнале на одном периоде есть и положительные, и отрицательные значения, положение волны легко заметить. В принципе, можно визуализировать и компоненты скорости, но за счет того, что скорость направлена в разные стороны в разных точках среды, профиль волны может опрокидываться, и определить фронт волны визуально сложнее.
Двумерный расчет модельной акустической задачи в неоднородной среде представлен на рисунке 49 в моменты: 0.1, 0.2, 0.3 и 0.4 сек. Неоднородность представлена углом размера 500м х 600м, находящегося в точке (500м, 400м) с более быстрой средой, Vp = 3000 м/сек. Видно, как часть звуковых волн отражается от границы раздела сред, а сама точка (500м, 400м) начинает играть роль самостоятельного источника. Длина волны равна Vp, что хорошо видно на рисунке - длина волны в более быстрой среде в 1.5 раза больше чем в основной.
В двумерном расчете модельной задачи для уравнений упругости в неоднородной среде берутся такие же условия, что и в предыдущих задачах. Единственная разница заключается в том, что в данном численном эксперименте коэффициент Пуассона А в обоих средах равен 0.25, т.е. к = ju. Результаты расчета в те же моменты времени, что и в акустическом случае, представлены на рисунке 50. Визуализируется компонента azz тензора напряжений. В данной постановке появляются сдвиговая s-волна, которые двигаются в ,/ А = v3 раз медленнее. S-волна двигается сама по себе, отражается от границ и интерферирует с другими s- и р-волнами. Отличие от акустического случая состоит также и в том, что при отражении и преломлении р-волны образуется еще и новая s-волна. Точно так же, угол в точке (500м, 400м) порождает теперь две волны: р и s. Порожденную s-волну хорошо видно в более быстрой среде при преломлении первой пришедшей р-волны - ее легко распознать, поскольку компонента czz на этой волне меняет знак, происходит опрокидывание профиля волны.
В трехмерной постановке размер основной области - 1км х 1км х 1км, при этом сетка состоит из 50 х 50 х 50 точек на основную область и по 10 точек с каждой стороны на дополнительные ячейки. Параметры среды - те же, что и раньше. Более быстрая среда сосредоточена в углу размером 700м х 700м х 600м, координаты угла - (300м, 300м, 400м), координаты источника - (300м, 300м, 200м). Изолинии (Ты, соответствующие первой акустической волне, показаны на рисунке 51 в моменты времени 0.1, 0.15, 0.2, 0.25, 0.3 и 0.35 сек. (На этом рисунке земная поверхность находится внизу, а угол, моделирующий более быструю среду, нарисован в виде параллелепипеда, обращенного к наблюдателю.) В трехмерных расчетах частота сигнала 3 Гц (что намного меньше, чем в двумерных расчетах) выбрана специально, для того чтобы на относительной редкой сетке разрешить профиль сигнала, и чтобы сигнал меньше затухал. Однако, частота источника не играет роли и не влияет на геометрическое расположение фронтов волн, по крайней мере, пока длина волны несоизмерима по порядку величины с характерными размерами тел и неоднородностей.
Расчет с несколько видоизмененными начальными данными: координаты источника - (250м, 250м, 200м), координаты угла - (350м, 300м, 300м), показан на рисунке 52 в трех сечениях х = 250м, у = 200м, z - 400м. Размер сетки составляет 280 х 280 х 240 точек, из них 200 х 200 х 200 точек приходится на основную область. Здесь поверхность Земли находится сверху, горизонтальная плоскость пересекает угол с другой средой, вертикальные плоскости не пересекают его, одна из вертикальных плоскостей содержит источник. На рисунке цветом показано распределение компоненты а в этих плоскостях. Представлены моменты времени 0.1, 0.25, 0.4, 0.55, 0.7 и 0.85 сек. Расстояние от источника до угла с более быстрой средой составляет 150 м и проходится продольной акустической волной за 0.2 сек. На четвертом кадре (t = 0.55 сек) кроме отраженной от поверхности Земли и преломленной в более быстрой среде волн, можно заметить еще и отраженную от границы раздела сред волну. Даже на далеких временах около самого источника наблюдается нефизическое поведение, связанное с точечным характером источника, жестким заданием сигнала в этой точке (отраженные волны не могут пройти через данную ячейку напрямую, а огибают ее в соседних ячейках), и наличием выделенных направлений сетки.
В последнем из представленных здесь расчетов моделируется прохождение сигнала в реальной земной коре. Распределение земных слоев по глубине и их параметры (Vp иУ5) приведены на рисунке 53. Источник находится на глубине 90 м, приемник перемещается в скважине на расстоянии 80 метров от источника и замеряет сигнал по всей глубине скважины от 0 до 100 м. Реальная сейсмограмма, полученная на основании опытных данных [193], приведена на рисунке 54.
Сейсмограмма, полученная в результате численного моделирования на сетке размером 200 х 40 х 200 точек (40 точек берется в направлении, ортогональном плоскости, проходящей через источник и ось приемника), представлена на рисунке 55. Здесь, как и на предыдущем рисунке, изображена вертикальная компонента скорости (в цвете). Видно характерное преломление волны на глубине 80 м. На глубине 60 м образуются преломленная и отраженная волны, также заметные на реальной сейсмограмме. Отраженная от земной поверхности волна в реальном случае намного слабее, чем в модельном. Хотя она тоже видна, но очень быстро затухает. Видимо, это связано с искусственным характером граничных условий в модельной задаче - производится полное отражение волны. В реальности же, часть энергии поглощается в поверхностных слоях Земли, часть рассеивается из-за неровностей земной поверхности.
