Физико-химические гидродинамические процессы часто встречаются в природе и технике. Их экспериментальное изучение является непростым. Во многих случаях проведение экспериментов невозможно. Проектирование новых технических установок с улучшенными экономическими и экологическими характеристиками при помощи экспериментов отнимает много времени, является дорогостоящим, а иногда даже опасным. Поэтому возникает большой интерес к компьютерному моделированию физико-химических процессов, например, к моделированию химических реакций в смеси газов. Но проведение вычислительных экспериментов оказывается сложной задачей Применение прямого численного моделирования сильно ограничено такими компьютерными ресурсами, как память и процессорное время, и будет оставаться таковым в ближайшем будущем. Для того чтобы избежать возникающих ограничений, основные усилия исследователей в настоящее время сосредоточены в двух направлениях: упрощение математических моделей процессов и/или моделирование на параллельных вычислительных комплексах.
Используемые математические модели. Все трехмерные системы уравнений в частных производных, описывающие физико-химические гидродинамические процессы, должны включать в себя законы сохранения полной массы (уравнение непрерывности), частичных масс фракций, импульса и энергии, уравнения состояния, уравнения химических реакций и правила для вычисления коэффициентов реакций. Полные системы уравнений, как правило, чрезвычайно трудны для численного решения. Поэтому обычно используются различные способы их упрощения. Основными направлениями таких упрощений являются: уменьшение пространственной размерности моделей, рассмотрение стационарных моделей, упрощение гидродинамической части моделей и моделей излучения, редукция химических реакций.
Методы дискретизации. Класс применяемых методов дискретизации, использующихся для аппроксимации систем уравнений физико-химической гидродинамики сравнительно невелик. Наиболее популярные способы аппроксимации — это метод конечных объемов на сетке в декартовых координатах [1, 2, 3] и методы конечных элементов [4, 5].
Выбор метода конечных объемов и/или метода конечных элементов позволяет получить высокую степень геометрического параллелизма на последующих стадиях решения задачи. Для численного моделирования применяются и некоторые другие схемы. В [6] пространственная дискретизация осуществляется методом конечных разностей с использованием центральных производных шестого порядка. Применение такого метода снижает схемную диссипацию, но не удаляет ее полностью. В приграничных узлах сетки используются схемы более низкого порядка. В [7] рассматриваются методы введения искусственной сжимаемости. Такой подход приводит к уменьшению жесткости дискретных систем уравнений и делает возможным некоторое увеличение временного шага при решении системы явным методом Рунге-Кутты. Как известно, явные методы могут быть легко и эффективно распараллелены на параллельных вычислительных системах с распределенной памятью.
Численные методы. Системы уравнений, получающиеся после дискретизации, являются жесткими и обладают существенной нелинейностью. Для решения жестких систем уравнений применяются два основных класса численных методов: связанные алгоритмы и алгоритмы операторного расщепления.
Связанные методы. В [6] интегрирование по времени осуществляется с помощью явного метода Рунге-Кутты четвертого порядка. Для устойчивости проверяется выполнение условий Куранта-Фридрихса-Леви и Фурье. Дополнительный контроль точности решения обеспечивается по правилу Рунге. Численное моделирование проводилось на
компьютере CRAY T3E. В качестве основного принципа распараллеливания выбрано геометрическое. В [4, 5] исследовались стационарные состояния. Полученная нелинейная система решалась на параллельном комплексе методом Ньютона. Все уравнения обрабатывались одновременно, взаимосвязи между переменными включались в матрицу Якоби.
Каждый процессор вычислял приращения по методу Ньютона для приписанных к нему узлов. После итерации по нелинейности полученные линейные системы решались итерационным методом с использованием программы из библиотеки AZTEC, реализующей предобусловленный метод Крылова.
Для решения уравнений, связывающих давления и скорости, и уравнений, выражающих законы сохранения, в [1] использован алгоритм SIMPLE. Для решения линеаризованных систем применялись итерационный метод с факторизацией оператора и функция пакета NAG, использующая решение СЛАУ в подпространстве Крылова. Конвективные члены аппроксимировались по схеме с центральными разностями. Для улучшения сходимости применялся многосеточный метод. Многосеточные методы были предложены Р. П. Федоренко [8].
Для декомпозиции расчетной области использовались вложенные сетки. Эта декомпозиция области вычислений была взята за основу для распараллеливания. Для балансировки загрузки процессоров использовалась статическая схема. Все процессоры считались одинаковыми по скорости работы, числа операций на каждом узле сетки равными, а подобласти выбирались так, чтобы включать примерно одинаковое количество узлов. Каждая подобласть приписывалась к одному процессору. Численное моделирование проводилось на компьютере CRAY T3E. Аналогичный подход был рассмотрен в [2], где для решения результирующей системы уравнений использовался метод предиктор-корректор второго порядка по времени и четвертого порядка по пространству. Расчетная сетка разбивалась на подобласти, которые распределялись между процессорами. Для переносимости алгоритма ис- пользовался механизм передачи сообщений MPI. Это позволило использовать полученный код на различных параллельных платформ и оценить на них его базовую производительность. Применение явной схемы допускает простую, но достаточно эффективную балансировку загрузки процессоров. В конце работы как предиктора, так и корректора на каждом шаге процессоры обменивались граничными значениями с соседями. При расчетах требовалась только одна общая коммуникационная операция в начале каждого шага для вычисления его величины. Работа [2] — одна из немногих, где в качестве параллельной платформы использовался кластер из персональных компьютеров. В [9] для решения гидродинамической части задачи реагирующих потоков был использован метод характеристик. Разработка численных методов решения уравнений в частных производных с использованием сеточно-характеристических методов была выполнена М.-К. М. Магомедовым и А . С. Холодовым [10]. Оригинальный путь па- раллельной реализации этого класса методов был разработан М. О. Васильевым [11, 12]. Операторное расщепление. Методы операторного расщепления используются в [3]. Для решения уравнений, описывающих поля давлений и скоростей, предложен SIMPLE-подобный алгоритм. Для возникающих систем алгебраических уравнений использован многосеточный подход, для уравнений сохранения массы фракций применена процедура расщепления на два шага: конвекция-диффузия и система химических реакций. Уравнения конвекции-диффузии решаются раздельно, а на химическом шаге — совместно для всех фракций и объемов с помощью модифицированного метода Ньютона для ЖС ОДУ. За основу для распараллеливания берется декомпозиция расчетной области. В качестве параллельной вычислительной системы с распределенной памятью использовался кластер Beowulf. Другой подход к расщеплению по операторам – это использование локально-одномерных операторов. Основные идеи таких методов были разработаны Н.Н. Яненко [13], Писменом и Рэкфордом [14]. В [15] метод локально-одномерной аппроксимации операторов был объединен с быстрым преобразованием Фурье для решения уравнений Навье–Стокса, что позволило добиться хороших показателей эффективности распараллеливания как на машинах массового параллелизма с распределенной памятью, так и на кластерах из рабочих станций.
Решение жестких систем ОДУ (ЖС ОДУ). Большинство численных методов для решения задачи Коши требуют многократного вычисления вектора правой части системы уравнений. Большая размерность делает это вычисление трудоемким. Часто задачи Коши для систем ОДУ являются жесткими, для устойчивости численного решения тре-
буются только неявные или диагонально-неявные методы. При их реализации приходится многократно вычислять матрицу Якоби и обращать ее. При использовании численных методов для решения задачи Коши необходимо управление шагом по времени. Предпочтительными оказываются методы, которые могут оценивать ошибки и использовать оценки для управления величиной шага.
За последние годы было разработано семейство одноитерационных методов Розенброка для ЖС ОДУ и систем дифференциально-алгебраических уравнений [16–18]. Главной их особенностью является только одно вычисление и обращение матрицы Якоби на каждом шаге.
В [17] приведен алгоритм построения жестко-точного A-устойчивого метода Розенброка порядка 4. Коэффициенты метода подобраны из условия минимизации вычислительной ошибки. Метод является вложенным, оценка погрешности вычислений выполняется без дополнительных затрат. Результаты сравнительных тестов [17] показывают, что построенный метод превосходит классические. В [19] метод типа Розенброка был применен для решения задач атмосферной химии. Неявный метод решения чрезвычайно жестких систем ОДУ предложен в [20] для моделирования многофазных химических реакций в облаке с каплями, разделенными на классы по размерам. Для интегрирования уравнений, описывающих химию фаз, фазовые переходы между газом и жидкостью и перенос масс между различными классами капель, использованы формулы дифференцирования назад (ФДН) высокого порядка. Решение возникающих разреженных систем линейных уравнений проводится модификацией кода LSODE, учитывающей особенности их структуры. В [20] применяется приближенная матричная факторизация с явной генерацией матрицы Якоби. Проведенные расчеты показали, что вычислительная стоимость метода растет линейно по отношению к количеству классов капель, но он не является параллельным. По тому же пути использования приближенной матричной факторизации для уменьшения вычислительной стоимости метода типа Розенброка идут и в [21].
Другой подход к модификации ФДН методов для их параллельной реализации на машинах массового параллелизма рассматривается в [22]. Этот подход развивается для жестких систем ОДУ с низкой степенью связности. Благодаря разреженности матриц вычисление правых частей, матриц Якоби, LU разложение матриц в методе Ньютона и прямая/обратная подстановки при решении могут быть выполнены параллельно.
Использование метода минимальной невязки в приближенном неявном методе для вычисления временного шага [23] на фиксированном числе итераций позволяет свести решение нелинейной задачи к решению набора линейных. Каждая из них может быть решена явно, метод обладает высоким потенциалом для распараллеливания. Существует большое количество параллельных реализаций методов типа Рунге–Кутта: однократно неявный метод [24, 25], диагонально-неявный метод [26], Ньютоновский параллельный итерационный солвер линейных систем для методов Рунге–Кутта [27]. Эти методы являются многостадийными, и вычисления стадий могут проводиться параллельно. В [28] исследуется возможность динамической балансировки загрузки процессоров для таких методов. Для параллельных реализаций методов типа Розенброка [29] основная идея является той же, что и для параллельных методов Рунге–Кутта. Каждая стадия обрабатывается на своем процессоре параллельно с остальными.
Заключение. Проведенный анализ работ в области задач физико-химической гидродинамики показал, что практически все проводимые исследования в настоящее время связаны с использованием численного моделирования на параллельных системах. При этом выделяются два основных подхода: геометрическое распараллеливание и использование многостадийных методов для решения жестких систем ОДУ.
Геометрическое распараллеливание применяется на этапе решения гидродинамических частей задач и легко реализуется. После решения задач в пределах своих подобластей на текущем шаге процессоры производят обмен новыми граничными значениями. Малые обмены данными делают при этом эффективным использование вычислительных комплексов с большим количеством процессоров и распределенной памятью.
Все многостадийные методы обладают достаточным потенциалом для распараллеливания. Для s-стадийного метода вычисления для значений стадий могут быть осуществлены параллельно на s процессорах. Из-за больших обменов информацией между процессорами эти методы являются особенно эффективными при использовании на машинах с
разделяемой памятью с небольшим числом процессоров – по числу стадий.
С точки зрения авторов представляется целесообразным соединить еометрический параллелизм, достигаемый при конечно-разностной аппроксимации, с многостадийным параллелизмом, используемым в методах типа Розенброка, на смешанных параллельных системах – кластерах, состоящих из SMP-узлов. Геометрический параллелизм реализуется на уровне отдельных узлов кластера, в то время как параллельное решение ЖС ОДУ внутри каждой подобласти – по стадиям на процессорах узла. Литература 1. Robbert Verweij L. et al. Parallel Computing for Reacting Flows Using Adaptive Grid Refinement // Contemporary Mathematics, V. 218, 1998. Р. 538546. 2. Stone C., Menon S.. Parallel Simulations of Swirling Turbulent Flames // The Journal of Supercomputing. V. 22. 2002. P. 7–28. 3. Cònsul R., Pérez-Segarra C.D. et al. Detailed numerical simulation of laminar flames by a parallel multiblock algorithm using loosely coupled computers // Combustion theory and modelling. V 7. 2003. P. 525–544. 4. Salinger A.G., Shadid J.N. et al. Parallel Reacting Flow Calculations for Chemical Vapor Deposition Reactor Design // Proc. of the Int. Conf. on Comput. Eng. Sci., San Jose, Costa Rica, May 4–9, 1997. 5. Salinger A.G., Pawlowski R.P. et al. Computational Analysis and Optimization of a Chemical Vapor Deposition Reactor with Large-Scale Computing, February 9, 2004 6. Lange M., Warnatz J. Massively Parallel Direct Numerical Simulation of Turbulent Combustion // NIC Symposium 2001, Proceedings, John von Neumann Institute for Computing, Jülich, NIC Series. V. 9. 2002. P. 419–429. 7. Wang Yi., Trouvé A.. Artificial acoustic stiffness reduction in fully compressible, direct numerical simulation of combustion // Combustion theory and modeling. 8. 2004. P. 633–660. 8. Федоренко Р.П. Релаксационный метод решения разностных эллиптических уравнений. // ЖВМиМФ. 1961. Т. 1. № 5. 9. Pakdee W., Mahalingam S. An accurate method to implement boundary conditions for reacting flows based on characteristic wave analysis // Combustion theory and modeling. V. 7. 2003. P. 705–729. 10. Магомедов М.-К., Холодов А .С. Сеточно-характеристические численные методы – М.: Наука, 1988. – 290 с. 11. Полуосьмак В.В., Васильев М.О. Разработка алгоритмов параллельного счета для решения задач магнитной гидродинамики в применении к задаче о взрыве в верхней ионосфере / Современные проблемы фундаментальных и прикладных наук. Тр. XLVII науч. конф. МФТИ. 2004. Т. 3. С. 199–200. 12. Vasilev M.O., Repin A.Ju. et al. Numerical researches of formation of jet stream of plasma in large-scale geophysical experiment. 13. Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. – Новосибирск, Наука. – 1967. 14. Peaceman D.W., Rachford H.H. The numerical solution of parabolic and elliptic differential equations, J. SIAM, 3(1955). Р. 28–41. 15. Averbuch A., Ioffe L., Israeli M., Vozovoi L.. Two-dimensional parallel solver for the solution of Navier–Stokes equations with constant and variable coefficients using ADI on cells // Parallel Computing. V. 24. 1998. P. 673–699. 16. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. – М.: Мир, 1990. – 512 с. 17. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально – алгебраические задачи. – М.: Мир, 1999. – 685 с. 18. Di Marzo G.A. RODAS5(4), methodes de Rosenbrock d’ordre 5(4) adaptees aux problemes differentiels-algebriques // Memoire de diplome en Mathematiques, Universite de Geneve 1992. 19. Sandu A., Verwer J.G. et al, Benchmarking of stiff ODE solvers for atmospheric chemistry problems II: Rosenbrock solvers. // Atmospheric Environment. V. 31, № 20. 1997. P. 3459–3472. 20. Wolke R., Knoth O. Time-integration of multiphase chemistry in sizeresolved cloud models // Applied Numerical Mathematics. V. 42. 2002. P. 473–487. 21. Botchev M.A., Verwer J.G. A new approximate matrix factorization for implicit time integration in air pollution modeling // J. of Comput. and Appl. Math. V. 157. 2003. P. 309–327. 22. Nordling P., Sјö A. Parallel solution of modular ODEs with application to rolling bearing dynamics // Mathematics and computers in simulation, 44 (1997). P. 495–504. 23. Botchev M.A., van der Vorst H.A. A parallel nearly implicit timestepping scheme Journal of Computational and Applied Mathematics, 137. 2001. P. 229–243. 24. Voss D.A., Muir P.H.. Mono-implicit Runge-Kutta schemes for the parallel solution of initial value ODEs. // J. of Comput. and Appl. Math. V. 102. 1999. P. 235–252. 25. Muir P.H. et al. PMIRKDC: a parallel mono-implicit Runge–Kutta code with defect control for boundary value ODEs // Parallel Computing 29. 2003. P. 711–741. 26. Jackson K.R., Norsett S.P. The Potential for Parallelism in Runge-Kutta Methods, SIAM J. Numer. Anal. 32, No. 1. 1995. P. 49–82. 27. Petcu D.. Experiments with an ODE Solver on a Multiprocessor System // An Int. Journal computers & mathematics with applications, 42. 2001. P. 1189–1199 28. Ruiz J. M.M., Lopera J. O., Carrillo de la Plata J. A. Component-Based Derivation of a Parallel Stiff ODE Solver Implemented in a Cluster of Computers. // International Journal of Parallel Programming. V. 30, No. 2. April 2002. P. 99–148. 29. Voss D.A., Khaliq A.Q.M.. Parallel Rosenbrock methods for chemical systems // Computers and Chemistry, 25. 2001. P. 101–107. |
