ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 8, с. 1320-1328
УДК 519.633
Посвящается светлой памяти А.П. Фаворского
ОБ ОДНОМ АЛГОРИТМЕ РЕШЕНИЯ ПАРАБОЛИЧЕСКИХ И ЭЛЛИПТИЧЕСКИХ УРАВНЕНИЙ1
© 2015 г. Н. Д'Асчензо*, В. И. Савельев**, Б. Н. Четверушкин***
(* Germany, D-22607Hamburg, Notkestaße, 85, Deutsche Elektronen-Synchrotron DESY;
** 236041 Калининград, ул. Александра Невского, 14, Балтийский федеральный ун-т;
*** 125047Москва, Миусская пл., 4, ИПМРАН) e-mail: office@keldysh.ru Поступила в редакцию 19.02.2015 г.
Наблюдаемый в настоящее время быстрый рост производительности вычислительных мощностей, в частности параллельных вычислительных систем сверхвысокой производительности, требует нового подхода к моделям и алгоритмам решения важнейших задач. В работе предложен алгоритм решения параболических и эллиптических уравнений. Возможности метода демонстрируются на примере решения задач астрофизики на высокопроизводительных вычислительных системах с массивным параллелизмом. Библ. 9. Фиг. 5.
Ключевые слова: высокопроизводительные вычисления, численные алгоритмы, параболические и эллиптические уравнения.
DOI: 10.7868/S0044466915080037
1. ВВЕДЕНИЕ
В настоящее время наблюдается быстрый рост производительности вычислительных систем, который достигается за счет параллельной работы огромного числа вычислительных устройств: вычислительных процессоров и ядер. Уже объявлены планы ввода в строй в ближайшие два года вычислительных систем с пиковой производительностью 300 PTFLOPS (1 PTFLOPS = 1015 операций с плавающей запятой в секунду). Следует отметить, что использование уже доступной на сегодняшний день производительности порядка 1 PTFLOPS позволяет применять для дискретизации трехмерных пространственных процессов вычислительные сетки, состоящие из нескольких миллиардов ячеек.
Такая детальная дискретизация важна для моделирования как инновационных индустриальных технологий, так и для многих фундаментальных научных проблем, таких как прямое моделирование турбулентности, задачи физики плазмы, астрофизики и др. Однако несмотря на наличие высокопроизводительных вычислительных систем, из-за потребности в детальной аппроксимации многих задач ее реализация протекает не столь гладко.
Причина подобной ситуации заключается в следующем. Неявные схемы, обладающие хорошей устойчивостью, требуют для своего решения тех или иных итерационных процедур. К сожалению, их адаптация на архитектуру вычислительных систем с экстрамассивным параллелизмом, состоящих из сотен тысяч и более независимых вычислительных устройств, приводит к практически полной потере эффективности параллельной обработки.
Явные схемы идеально адаптируются на архитектуру параллельных многопроцессорных вычислительных систем. Однако они требуют для устойчивости счета жесткого ограничения на шаг по времени. И если для гиперболических уравнений ограничение на шаг курантовского типа имеет вид
А t < h, (1.1)
1) Работа выполнена при финансовой поддержке Российского научного фонда (грант 14-11-00170).
где At — шаг по времени, к — шаг пространственной дискретизации, являющийся жестким, но естественным ограничением (см. [1]), то для параболических уравнений условие
А? < к2 (1.2)
делает невозможным проведение расчетов на детальных пространственных сетках. (Естественность ограничения (1.1) вытекает из того факта, что более дательная пространственная детализация должна быть связана с более детальным описанием времени.) Заметим, что одной из побудительных причин использования высокопроизводительных вычислительных систем стала возможность детального описания физических явлений и технологических процессов.
Таким образом, для многих важнейших физических процессов моделирование, дающее возможность детального пространственно-временного описания, крайне затруднено. Парадоксальность этой ситуации подчеркивается как потребностью точного описания, так и наличием технических средств в виде высокопроизводительных параллельных вычислительных систем.
Целью данной работы является обсуждение явных алгоритмов, которые, с одной стороны, допускали бы простую адаптацию на архитектуру вычислительных систем с экстрамассивным параллелизмом, а с другой стороны — более мягкое условие устойчивости, чем условие (1.2).
Одним из побудительных мотивов разработки таких алгоритмов явилась потребность в решении задач магнитной гидродинамики (МГД), в частности связанных с решениями проблем астрофизики. Некоторые результаты расчета этих задач будут приведены ниже.
2. ВТОРАЯ ПРОИЗВОДНАЯ ПО ВРЕМЕНИ (ГИПЕРБОЛИЗАЦИЯ УРАВНЕНИЙ)
В качестве исходного параболического уравнения рассмотрим линейное уравнение
дии = А кди + Р} (2.1)
д? дх{ дх{
где I — время, XI — пространственная координата в /-направлении, и — искомая функция, Е — источник.
Наряду с исходным параболическим уравнением (2.1) рассмотрим гиперболическое уравнение
дии + х * д2и = д_ к д_и + р (2.2)
д? _?2 дх, дх1
Коэффициент т* выбираем из условия таким, чтобы член с первой производной по времени по порядку величины был существенно больше члена со второй производной по времени:
* д и
т Т . д?2-
<
'и
.д?
(2.3)
Идея использования второй производной по времени вытекает из опыта решения задач газовой динамики с помощью квазигазодинамической системы уравнений (см. [2]), которую для примера выпишем для одномерной пространственной геометрии:
до д2р , дри д д / 2 , ч
+ т—-- + = —т —(ри + р), (2.4)
д? д?2 дх дх дх
_Ри + т __2ш + д( р и + р ) = АД(р и3 + Зри), (2.5)
д? д?2 дх дх дх
дЕ д2 Е , д , гт? . \\ дтд ,дтд
~я'- + тТ + Т(и(Е + Р)) = Ти (Е + 2Р)) + -Г-^-
д? д?2 дх дх 2 дх дх 2 дх
1 ( Е + р) р
(2.6)
Здесь использовались следующие обозначения: р — плотность, и — скорость, р — давления, Е —
2
полная энергия Е = р-и- + 6, 6 — внутренняя энергия, т — характерное время между столкновениями молекул.
Появление членов со второй производной по времени и форма диссипативных членов в правой части системы (2.4)—(2.6) связаны с учетом минимальных масштабов по пространству (длина свободного пробега I) и времени т, меньше которых в газовой динамике дальнейшая детализация решения не имеет смысла (см. [3]).
Анализ системы (2.4)—(2.6) показывает, что она отличается от классической системы уравнений Навье—Стокса на члены второго порядка малости по числу Кнудсена. Кстати, с этой точностью сами уравнения Навье—Стокса получаются из исходного кинетического уравнения Больц-мана. Система (2.4)—(2.6) в течение нескольких десятилетий успешно использовалась для моделирования различных гидро- и газодинамических течений (см. [2]).
Частным случаем уравнения (2.6) при отсутствии газодинамического течения является уравнение теплопроводности. Однако данная теплопроводность будет носить гиперболический характер. Заметим, что еще несколько десятилетий тому назад гиперболическое уравнение теплопроводности начали использовать при изучении процессов взаимодействия мощных коротких лазерных импульсов с веществом.
Вернемся к исходным системам (2.1) и (2.2). Условие (2.3) означает, что по физическим соображениям решения обеих систем должны различаться незначительно. Приведем теоретическую оценку разности е = и2 — и1 (см. [4]), где и2 — решение системы (2.2), а и1 — решение системы (2.1):
/(2 - у2)||У е|| 2йг + е(х, 012йх <(т *)2С- Л
т 2 д и-
дг2
2
йхйг, (2.7)
0 О 0 о
где Т — рассматриваемый временной интервал, О — область интегрирования по пространству,
IV Н' = /V , (2.8>
о
у2 — положительна константа (у2 < 2), С1 — ограниченная константа, зависящая от области О.
Из оценки (2.7) следует, что если константа т* мала, то решения уравнений (2.1), (2.2) различаются незначительно. Это утверждение не относится к высокочастотным компонентам решения. Однако необходимо помнить, что быстрое изменение по времени, как правило, связано с резкими изменениями на коротких расстояниях по пространству. А именно такие изменения исключаются при разностной аппроксимации на сетке с характерным пространственным шагом к как для численного решения уравнения (2.1), так и уравнения (2.2).
3. ЯВНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ РЕШЕНИЯ ГИПЕРБОЛИЗИРОВАННОГО УРАВНЕНИЯ
Для решения уравнения (2.2) выберем следующую трехслойную явную схему второго порядка точности 0(А?2 + к2) по пространству и времени:
и+1 _ и_1 и+1 _ 2 и + и_1
и 2Аи + т*и-Щ±и- = ((ки)^ + Г. (3.1)
-А г А Г ' '
Для аппроксимации второй производной по пространству воспользуемся известными формулами для потоков (см. [1]).
Характерный шаблон для данной схемы в пространственно-одномерном случае приведен на фиг. 1.
Коэффициент т* можно выбрать любым достаточно малым для выполнения условия (2.3). Однако из физических соображений наиболее разумным представляется выбор т*, по порядку величины связанным с характерным шагом пространственной сетки к и характерной для данного процесса скоростью с:
т* - к/с. (3.2)
Схема (3.1) позволяет легко вычислять значения функции и +1 на ] + 1-м слое по времени через уже известные значения этой функции на]-м и] — 1-м слое по времени.
и/
+ 1
т. ,
I - 1
и/
т л
. 1+1
и/ -
Фиг. 1. Шаблон для трехслойной разностной схемы.
Так же, как и явная схема
и(+1 - и{
А t
= ((кЩ-х + ¥
(3.3)
для параболического уравнения (2.1), схема (3.3) хорошо адаптируется на архитектуры высокопроизводительных вычислительных систем с экстрамассивным параллелизмом.
Рассмотрим вопрос об устойчивости трехслойной схемы (3.1). Используя предложенный в [5] алгоритм оценки устойчивости, получаем, что схема (3.1) с т*, определяемая с помощью выражения (3.2), устойчива при следующем ограничении на шаг по времени (см. [6]):
Ат < к
(3.4)
Конечно, условие (3.4) является достаточно жестким. Однако оно представляется более приемлемым, чем условие устойчивости (1.2) для явной разностной схемы (3.3). Особенно наглядно преимущество условия (3.4) наблюдается при использовании детальных пространственных сеток.
В следующем разделе приведена демонстрация результатов расчетов с использованием схем (3.1) и (3.3).
1
4. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ РАСЧЕТОВ
В данном разделе будут представлены результаты чи
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.