Высшая математика

Численные методы решения обратных задач теплообмена (ОЗТ)

Поиск причины по ее следствиям

В инженерной теплофизике прямая задача теплопроводности решается просто: если мы знаем начальную температуру тела, его физические свойства (теплопроводность, теплоемкость) и граничные условия (например, мощность нагревателя на поверхности), численная модель (МКЭ или МКР) точно предскажет распределение температуры внутри детали в любой момент времени. Но в аэрокосмической отрасли и металлургии инженеры чаще сталкиваются с обратной ситуацией. Представьте, что спускаемый космический аппарат входит в плотные слои атмосферы. Его внешняя обшивка подвергается экстремальному плазменному нагреву. Установить термодатчики на самой внешней поверхности невозможно — они мгновенно сгорят или нарушат аэродинамику. Датчики устанавливаются глубоко внутри теплозащитного экрана.

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

Итерационная регуляризация и метод Алифанова

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

Задача формулируется как проблема оптимального управления (или минимизации). Мы ищем такой тепловой поток на границе, при котором рассчитанная моделью температура в местах расположения датчиков максимально совпадает с реальными экспериментальными показаниями (минимизируется среднеквадратичная невязка). Для поиска минимума этого функционала используется метод градиентного спуска или метод сопряженных градиентов. Роль регуляризатора (сглаживателя) здесь играет сам номер итерации! На первых итерациях алгоритм восстанавливает плавные, низкочастотные компоненты теплового потока. Если мы не остановим процесс, на поздних итерациях алгоритм начнет «подгонять» решение под высокочастотный шум датчиков, и ответ разрушится. Критерий останова (принцип невязки) приказывает прервать итерации в тот момент, когда отклонение модели от эксперимента становится равно погрешности самих приборов. Так мы получаем устойчивый и физически осмысленный результат.

Применение сопряженных уравнений для поиска градиента

Главная трудность метода Алифанова — это вычисление градиента функционала невязки. В нашей задаче параметром оптимизации выступает не одно число, а целая функция (график теплового потока во времени), содержащая тысячи дискретных точек. Если вычислять градиент по каждой точке конечными разностями, нам придется решать прямую задачу теплопроводности тысячи раз на каждой итерации, что займет вечность.

Эта проблема элегантно решается с помощью метода сопряженных уравнений (Adjoint method). Составляется специальная вспомогательная краевая задача, называемая сопряженной, которая интегрируется «вспять» по времени (от конца эксперимента к его началу). Правой частью для этого сопряженного уравнения служит разность между показаниями датчиков и расчетом. Решив сопряженную задачу всего один раз, мы получаем функцию, которая путем простого умножения мгновенно дает нам точнейший градиент целевого функционала сразу по всем тысячам параметров времени! Этот красивый математический аппарат сегодня применяется не только в теплофизике, но и в обратных задачах сейсморазведки для поиска месторождений нефти по отраженным акустическим волнам.

Подробнее

Методы декомпозиции областей (Domain Decomposition): параллельные вычисления

Преодоление ограничений памяти и распараллеливание

С развитием суперкомпьютеров вычислительная математика столкнулась с новой, сугубо архитектурной проблемой. Когда мы решаем дифференциальные уравнения в частных производных (УЧП) на сетках, состоящих из десятков миллиардов узлов (например, при расчете глобальной климатической модели Земли или аэродинамики космического корабля в сборе), физически невозможно поместить глобальную матрицу системы в оперативную память одного вычислительного узла. Возникает острая необходимость разрезать гигантскую математическую задачу на тысячи более мелких кусков, раздать их независимым процессорам (или узлам кластера) и заставить эти процессоры общаться между собой по сети MPI для получения единого, математически строгого решения.

Эта концепция легла в основу методов декомпозиции областей (Domain Decomposition Methods, DDM). Идея метода невероятно стара: еще в 1870 году математик Герман Шварц придумал альтернирующий метод для аналитического решения уравнения Лапласа в сложных областях. Он показал, что если сложная область состоит из двух перекрывающихся простых фигур (например, круга и прямоугольника), можно поочередно решать задачу в каждой из них, используя на границе перекрытия значения, полученные у соседа на предыдущем шаге. Спустя сто лет эта абстрактная математическая теорема стала фундаментом параллельных суперкомпьютерных вычислений.

Методы с перекрытием (Метод Шварца) и без перекрытия

Современные численные реализации DDM делятся на две большие группы. Первая группа — это методы с перекрытием (Overlapping methods), прямые наследники классического метода Шварца. Расчетная сетка разбивается на субдомены, которые искусственно расширяются так, чтобы они немного наслаивались друг на друга. Итерационный процесс заключается в том, что каждый процессор решает свою локальную СЛАУ внутри своего домена, а затем пересылает значения из зоны перекрытия соседним процессорам, чтобы те использовали их как краевые условия Дирихле на следующей итерации. Этот метод очень надежен, но пересылка больших массивов данных в зонах перекрытия может замедлять сеть кластера.

Вторая группа — это методы без перекрытия (Non-overlapping methods, или методы подконструкций). Здесь граница между доменами представляет собой строгую математическую линию (или поверхность) толщиной в один узел сетки. Эти методы математически более сложны, так как на общей границе необходимо обеспечить не только непрерывность самой функции (температуры), но и непрерывность нормальной производной (теплового потока). Для этого используются специальные алгоритмы на интерфейсах (например, метод Дирихле-Неймана или метод Неймана-Неймана).

Дополнение Шура и предобуславливание FETI

Самым мощным математическим аппаратом для методов без перекрытия является матричное дополнение Шура (Schur Complement). Алгоритм алгебраически исключает все внутренние узлы каждого субдомена, сводя гигантскую глобальную СЛАУ к гораздо меньшей системе уравнений, которая записана исключительно относительно неизвестных значений на самих интерфейсах (границах разреза). После решения интерфейсной задачи каждый процессор абсолютно независимо и быстро восстанавливает решение внутри своего куска.

Для решения сложной интерфейсной системы обычно применяется метод сопряженных градиентов (CG). Чтобы он сошелся быстро, французский ученый Шарль Фархат разработал знаменитый алгоритм FETI (Finite Element Tearing and Interconnecting), который вводит множители Лагранжа (силы взаимодействия) на границах и использует сложные проекционные предобуславливатели. Сегодня алгоритмы типа FETI-DP (с двумя уровнями грубой сетки) обеспечивают идеальную масштабируемость: программа может эффективно задействовать сотни тысяч процессорных ядер, ускоряя решение пропорционально размеру суперкомпьютера.

Подробнее

Метод дифференциальной эволюции: векторная мутация в глобальной оптимизации

Оптимизация в непрерывных недифференцируемых пространствах

Инженерам часто приходится решать задачи глобальной оптимизации, где целевая функция представляет собой «черный ящик» — мы можем вычислить ее значение при заданных параметрах (например, запустив тяжелую аэродинамическую симуляцию), но мы совершенно не знаем ее градиентов (производных). Более того, ландшафт этой функции может быть изрезан оврагами, содержать сотни ложных локальных минимумов и иметь разрывы. Как мы уже выяснили, классические методы Ньютона или градиентного спуска в таких условиях бесполезны: они застревают в первой же яме.

Для таких задач были созданы стохастические генетические алгоритмы (ГА). Классические ГА работают путем кодирования параметров в виде бинарных строк (хромосом) из нулей и единиц. Однако при оптимизации непрерывных вещественных параметров (например, длина волны лазера, концентрация реагента в молях, напряжение в вольтах) бинарное кодирование порождает массу проблем, связанных с потерей точности и эффектом Хэмминга. В 1995 году Кеннет Прайс и Райнер Сторн разработали революционный алгоритм — Дифференциальную эволюцию (Differential Evolution, DE). Этот алгоритм изначально создавался для работы напрямую с векторами чисел с плавающей запятой, устраняя необходимость в громоздких бинарных преобразованиях.

Векторная разность как направляющая сила мутации

Дифференциальная эволюция относится к классу популяционных методов. Программа создает случайное облако (популяцию) из N возможных решений, разбросанных по всему многомерному пространству поиска. Вектор каждого решения хранит конкретные вещественные числа (координаты). Настоящая магия DE кроется в ее операторе мутации, который в корне отличается от слепого случайного шума в классических ГА.

Для того чтобы создать нового «мутанта» из текущей особи (целевого вектора), алгоритм случайным образом выбирает из популяции еще три других, абсолютно независимых вектора: A, B и C. Затем он вычисляет векторную разность между B и C. Эта математическая разность масштабируется (умножается на постоянный коэффициент мутации F) и прибавляется к вектору A. Полученный вектор становится основой для нового потомка. В чем гениальность этого подхода? Векторная разность автоматически адаптирует алгоритм к топологии ландшафта! В начале работы, когда популяция разбросана далеко друг от друга, разности огромны, и алгоритм делает гигантские разведывательные шаги, перепрыгивая через любые локальные минимумы. Когда алгоритм нащупывает глобальную долину, популяция стягивается в плотный клубок, разности между векторами становятся микроскопическими, и алгоритм переходит в режим ювелирной, сверхточной настройки параметров, не требуя ручного изменения радиуса мутации (в отличие от методов имитации отжига).

Скрещивание, селекция и индустриальное применение

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

Всего три простых математических операции — сложение векторов, масштабирование разности и условная замена. Алгоритм дифференциальной эволюции умещается в 20 строк программного кода, требует настройки всего двух интуитивных параметров (F и CR) и феноменально легко распараллеливается на многоядерных процессорах. Благодаря своей беспрецедентной способности находить глобальный минимум сложных невыпуклых функций, метод DE стал основным алгоритмом для настройки гиперпараметров нейронных сетей, автоматического синтеза фрактальных Wi-Fi антенн сложной формы и проектирования молекул новых лекарственных препаратов в хемоинформатике.

Подробнее

Методы Ланцоша и Арнольди: поиск собственных значений гигантских матриц

Охота за экстремальными резонансами

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

К счастью, в таких задачах инженеру почти никогда не нужны абсолютно все миллионы частот. Обычно критически важно знать лишь несколько самых малых собственных значений (они определяют фундаментальную частоту колебаний моста) или несколько самых больших (они определяют алгоритмическую устойчивость системы). Для поиска этих экстремальных собственных значений у сверхбольших разреженных матриц были созданы мощнейшие итерационные алгоритмы: метод Ланцоша (для симметричных матриц) и метод Арнольди (для несимметричных).

Подпространство Крылова и ортогональный базис

Оба метода базируются на элегантной концепции подпространств Крылова, которую мы упоминали в контексте решения СЛАУ. Алгоритм выбирает случайный начальный вектор q1 и начинает многократно умножать его на гигантскую матрицу A. В результате формируется последовательность векторов: q1, A*q1, A^2*q1, A^3*q1. Эти векторы несут в себе колоссальную математическую информацию о доминирующих собственных векторах матрицы (так как при многократном умножении вектор вытягивается вдоль главного собственного направления).

Проблема в том, что эти векторы очень быстро становятся почти коллинеарными (параллельными), что уничтожает их вычислительную полезность (матрица их скалярных произведений становится вырожденной). Гениальность метода Арнольди (Уолтер Арнольди, 1951) заключается в применении процесса ортогонализации Грама-Шмидта прямо на лету! Каждый новый сгенерированный вектор тут же жестко проецируется и очищается от компонентов всех предыдущих векторов. В результате мы получаем идеальный ортонормированный базис подпространства Крылова. В этом новом, маленьком базисе исходная гигантская матрица проецируется (сжимается) в крошечную матрицу Хессенберга размера m x m (где m — число итераций, например, 50 или 100).

Метод Ланцоша и проблема потери ортогональности

Если исходная матрица системы была симметричной (как в задачах сопротивления материалов), процесс Арнольди превращается в Метод Ланцоша (Корнелий Ланцош, 1950). Из-за свойств симметрии матрица Хессенберга магическим образом схлопывается в симметричную трехдиагональную матрицу. Это означает, что для генерации каждого нового ортогонального вектора алгоритму Ланцоша не нужно вспоминать и вычитать все предыдущие 100 векторов — ему достаточно использовать короткую рекуррентную формулу, опирающуюся только на два последних вектора! Это делает метод Ланцоша феноменально быстрым и экономным по памяти.

Получив эту крошечную матрицу проекции, мы просто решаем ее стандартным QR-алгоритмом за миллисекунды. Ее собственные значения (называемые числами Ритца) с фантастической скоростью сходятся к истинным экстремальным собственным значениям гигантской исходной матрицы. Однако у метода Ланцоша есть вычислительный недуг: из-за конечной разрядности компьютеров (ошибок округления float) генерируемые векторы постепенно теряют свою взаимную перпендикулярность. Это порождает появление фиктивных (призрачных) собственных значений, которые являются математической галлюцинацией. Для борьбы с этим современная вычислительная алгебра использует методы полной или выборочной переортогонализации, которые встроены в ядро таких сверхмощных индустриальных пакетов, как ARPACK (используемых функциями `eigs` в MATLAB и SciPy).

Подробнее

Векторизация вычислений (SIMD): аппаратное ускорение численных алгоритмов

Когда закон Мура перестает ускорять код

На протяжении десятилетий программисты и математики привыкли к простому правилу: чтобы программа для численного интегрирования или умножения матриц заработала быстрее, нужно просто подождать пару лет, пока Intel или AMD выпустят процессор с более высокой тактовой частотой. Этот феномен, подпитываемый законом Мура, позволял писать неоптимальный код и все равно получать рост производительности. Однако в начале 2000-х годов рост тактовых частот уперся в непреодолимый физический барьер — тепловой предел кремниевых чипов (частота застыла на уровне 3-5 ГГц). Чтобы продолжать наращивать вычислительную мощь, создателям процессоров пришлось менять архитектуру, внедряя многоядерность и, что еще более важно для математики, векторные инструкции.

Концепция SIMD (Single Instruction, Multiple Data — одна инструкция, множество данных) является фундаментом высокопроизводительных вычислений (HPC). В классической, скалярной модели процессора (SISD), чтобы сложить два массива из четырех чисел (A1+B1, A2+B2, A3+B3, A4+B4), процессору нужно выполнить цикл: загрузить два числа, сложить их, записать результат, и повторить это 4 раза. Архитектура SIMD внедряет в процессор сверхширокие регистры. Загрузив сразу 4 числа из массива A в один 256-битный регистр, и 4 числа из массива B во второй регистр, процессор за один единственный такт времени выполняет сложение сразу всех четырех пар чисел одновременно! Это обеспечивает мгновенный 4-кратный (или даже 8-кратный) прирост чистой математической производительности.

Наборы инструкций AVX и выравнивание памяти

Современные процессоры оснащены мощнейшими блоками векторной математики. Инструкции SSE (128 бит) позволяют складывать 4 числа типа float за такт. Инструкции AVX и AVX2 (256 бит) удваивают этот показатель до 8 чисел. Новейшие серверные процессоры с инструкциями AVX-512 оперируют колоссальными 512-битными регистрами, выполняя по 16 операций сложения или умножения 32-битных чисел с плавающей запятой за один микроскопический такт процессора. Более того, существует специальная инструкция FMA (Fused Multiply-Add), которая за один такт выполняет сразу две операции (A * B + C), удваивая пиковую производительность алгоритма матричного умножения.

Но эта фантастическая аппаратная мощь недоступна программам автоматически. Чтобы векторные блоки заработали, данные в оперативной памяти компьютера должны быть расположены (упакованы) строго последовательно, без разрывов, и желательно быть выровненными (aligned) по границам 32-байтовых или 64-байтовых кеш-линий. Если вы написали алгоритм метода конечных разностей, который скачет по памяти (собирая данные из разных классов или объектов), процессор просто не сможет загрузить их широким вектором. Он откатится к медленным скалярным вычислениям (потеряв до 90% скорости).

Компиляторы, раскрутка циклов и библиотеки BLAS

Чтобы заставить численный код летать, разработчики применяют агрессивную раскрутку циклов (Loop Unrolling) и прагмы для компиляторов (например, #pragma omp simd в OpenMP). Важнейшим требованием для векторизации является отсутствие логических ветвлений (операторов if/else) внутри тяжелого математического цикла, так как ветвление разрушает векторный конвейер.

В мире профессиональной вычислительной математики никто не пишет циклы умножения матриц вручную. Вся базовая алгебра (нахождение скалярных произведений, умножение матрицы на вектор, разложение матриц) делегируется стандартизированным библиотекам BLAS (Basic Linear Algebra Subprograms) и LAPACK. Такие реализации, как Intel MKL (Math Kernel Library) или OpenBLAS, написаны на глубоком ассемблерном уровне. Они досконально учитывают размер L1/L2 кэша процессора и микроархитектуру AVX регистров, нарезая матрицы на микро-блоки. Именно благодаря идеальной SIMD-векторизации внутри библиотек BLAS, скрипты на языке Python (использующие NumPy) умудряются работать со скоростью скомпилированного языка C++ при решении гигантских систем дифференциальных уравнений.

Подробнее

Ускорение граничных элементов: Быстрый мультипольный метод (FMBEM)

Преодоление барьера плотных матриц в акустике и электродинамике

Метод граничных элементов (МГЭ, или BEM) обладает колоссальным преимуществом перед объемными методами: он снижает размерность физической задачи на единицу. Для расчета распределения электромагнитного поля вокруг огромного самолета инженеру нужно покрыть сеткой только поверхность обшивки фюзеляжа, а не кубические километры окружающего пространства. Однако за это математическое изящество приходится платить суровую цену. Переход к поверхностным интегральным уравнениям порождает систему линейных алгебраических уравнений (СЛАУ), матрица которой является плотной (dense). Она заполнена ненулевыми элементами на 100%, потому что каждая точка (заряд) на поверхности самолета непрерывно взаимодействует со всеми остальными точками посредством закона Кулона.

Для сетки всего из 100 000 элементов матрица потребует около 80 гигабайт оперативной памяти для хранения и колоссального количества процессорного времени (порядка O(N^2) операций для умножения на вектор и O(N^3) для прямого решения методом Гаусса). Это означает, что классический МГЭ упирается в «стеклянный потолок» размерности: попытка решить задачу для детализированной сетки современного подводного корабля или антенной решетки из миллионов элементов становится математически невозможной даже для самых мощных суперкомпьютеров министерства обороны.

Интеграция с Быстрым мультипольным методом (FMM)

Спасением для метода граничных элементов стала его интеграция с революционным Быстрым преобразованием мультиполей (Fast Multipole Method, FMM), разработанным Владимиром Рохлиным. Этот синтез получил название Быстрого мультипольного метода граничных элементов (Fast Multipole Boundary Element Method, FMBEM). Идея алгоритма основана на строгом математическом факте: влияние группы далеких источников (точек на поверхности) на приемник можно сгруппировать и заменить одним эквивалентным источником, используя разложение фундаментального решения (функции Грина) в ряды по сферическим гармоникам.

С вычислительной точки зрения алгоритм FMBEM полностью отказывается от сборки и хранения в памяти плотной глобальной матрицы. Алгоритм строит иерархическое дерево кубических ячеек (Octree), заключающее в себе всю поверхность модели. Если два элемента сетки находятся в соседних, близких ячейках, их взаимное влияние вычисляется точно, прямым интегрированием функций формы, и эти матричные коэффициенты образуют сильно разреженную локальную матрицу (ближнее поле). Если же элементы находятся в далеких ячейках (дальнее поле), матричные элементы не вычисляются вообще!

Мультипольные разложения и ускорение до O(N log N)

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

Использование FMM в качестве способа быстрого умножения плотной матрицы на вектор позволяет применять итерационные методы подпространств Крылова (например, GMRES). В результате вычислительная сложность и затраты оперативной памяти радикально падают с O(N^2) до O(N log N) или даже O(N). Это технологический прорыв колоссального масштаба: сегодня коммерческие FMBEM солверы позволяют моделировать рассеяние акустических волн на турбинах подводных лодок, используя сетки из 10-50 миллионов узлов, решая задачи, которые еще 20 лет назад считались абсолютно нерешаемыми в вычислительной акустике.

Подробнее

Разрывный метод Галеркина (DG): точные скачки в вычислительной гидродинамике

Когда непрерывность становится помехой

Классический метод конечных элементов (МКЭ) десятилетиями является непререкаемым стандартом в задачах строительной механики и теории упругости. Его успех базируется на использовании непрерывных базисных функций. Для расчета прогиба балки непрерывность перемещений между конечными элементами — это строгий физический закон. Однако при попытке применить этот же непрерывный базис к задачам газовой динамики, гиперболическим законам сохранения и расчету сверхзвуковых течений классический МКЭ потерпел катастрофическое фиаско. Причина кроется в природе самих физических явлений: сверхзвуковые потоки порождают ударные волны (скачки уплотнения), которые с математической точки зрения представляют собой строгие разрывы плотности и давления. Непрерывные полиномы МКЭ просто не могут адекватно описать разрыв функции; при попытке это сделать они генерируют гигантские нефизичные осцилляции по всей расчетной области, полностью разрушая вычислительный процесс.

Для решения этой фундаментальной проблемы математики Рид и Хилл в 1973 году (а позже Кокберн и Шу в 1980-х) разработали Разрывный метод Галеркина (Discontinuous Galerkin Method, DGM). Этот алгоритм представляет собой элегантный синтез двух миров: он берет теоретическую строгость и высокую точность полиномиальных аппроксимаций из метода конечных элементов и объединяет их с физически корректной обработкой ударных волн, заимствованной из метода конечных объемов (МКО).

Локальная база и численные потоки на границах

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

Но как эти острова обмениваются физической информацией (массой, энергией)? Здесь вступает в игру математический аппарат метода конечных объемов. При интегрировании по частям в слабой формулировке Галеркина возникают контурные интегралы по границам элементов. В эти интегралы подставляются так называемые «численные потоки» (Numerical Fluxes) — например, потоки Роу (Roe) или HLLC, которые получаются путем точного или приближенного решения локальной задачи Римана о распаде разрыва между соседними элементами. Этот поток действует как клей, который обеспечивает строгий закон сохранения массы и энергии во всей глобальной системе, при этом допуская наличие резких физических скачков прямо на границах ячеек.

Матрицы масс и метод Рунге-Кутты (RKDG)

Поскольку базисные функции жестко локализованы внутри каждого элемента и не перекрываются с соседями, глобальная матрица масс системы перестает существовать в привычном виде. Она распадается на множество крошечных блочно-диагональных матриц (размером с полиномиальный базис одного элемента). Эту маленькую матрицу можно обратить аналитически один раз до начала расчетов. В результате мы получаем чисто явную дифференциальную схему по времени, которая вообще не требует решения гигантских СЛАУ на каждом шаге!

Для интегрирования по времени используются специальные TVD (Total Variation Diminishing) методы Рунге-Кутты, которые математически гарантируют невозрастание ложных осцилляций. В сочетании с ограничителями наклона (Slope limiters), которые принудительно «срезают» возникающие локальные пики возле ударных волн, метод RKDG позволяет моделировать аэродинамику гиперзвуковых ракет, детонацию взрывчатых веществ и акустические шумы турбин с невероятной точностью (сходимость 4-го, 5-го и более высоких порядков) на совершенно неструктурированных сетках, обеспечивая идеальную масштабируемость на графических ускорителях (GPU).

Подробнее

Адаптивные методы построения сеток (AMR): динамическое измельчение в расчетах

Борьба за вычислительные ресурсы

В вычислительной гидродинамике (CFD), плазменной физике и астрофизике исследователи сталкиваются с процессами, которые охватывают колоссальные пространственные масштабы, но содержат локальные явления невероятной сложности. Представьте моделирование взрыва сверхновой звезды. Большую часть времени в окружающем пространстве не происходит ничего интересного — там вакуум и плавные градиенты. Но в тончайшей зоне фронта ударной волны происходят катастрофические скачки давления и плотности. Если мы попытаемся накрыть всю область моделирования равномерной сеткой высокого разрешения (чтобы «поймать» узкий фронт волны), нам потребуются триллионы узлов, и памяти не хватит даже у крупнейших суперкомпьютеров.

Идеальное решение — использовать густую сетку только там, где это необходимо (на фронте волны), а в остальном пространстве оставить грубую редкую сетку. Но проблема в том, что ударная волна движется! Зона, которая нуждалась в густой сетке секунду назад, сейчас уже не требует такого разрешения. Сетка должна уметь «дышать» и подстраиваться под физику процесса на лету. Эта технология называется Адаптивным измельчением сеток (Adaptive Mesh Refinement, AMR).

Алгоритмы Бергера-Колеллы и деревья октантов (Octrees)

Фундамент современных блочно-структурированных методов AMR был заложен Маршей Бергер и Филиппом Колеллой в 1980-х годах. Алгоритм начинается с базовой, очень грубой прямоугольной сетки. На каждом временном шаге специальный математический контроллер анализирует текущее решение. Он использует эвристические критерии (например, вычисляет градиенты давления или оценки ошибки усечения). Если в определенной ячейке ошибка превышает допустимый предел, эта ячейка помечается для измельчения (Refinement).

В древовидных структурах (Quadtree в 2D, Octree в 3D) алгоритм делит помеченную ячейку (родителя) ровно на 8 дочерних ячеек меньшего размера. Этот процесс может продолжаться рекурсивно на несколько уровней в глубину (Level of Detail), создавая сверхплотное разрешение точно вокруг фронта горения или завихрения. И наоборот, если градиенты в ячейке исчезают (волна ушла), 8 мелких ячеек сливаются обратно в одну крупную (Coarsening), освобождая драгоценную оперативную память процессора.

Балансировка нагрузки (Load Balancing) на суперкомпьютерах

Алгоритмы AMR — это не просто математика, это шедевр системного программирования. Главный вызов при использовании AMR на кластерах из тысяч процессоров (MPI) — это балансировка нагрузки (Load Balancing). Если ударная волна переместилась из левой части расчетной области в правую, то процессор, обсчитывавший левую часть, внезапно лишится мелких ячеек и начнет простаивать, а процессор, отвечающий за правую часть, захлебнется от миллионов новых сгенерированных узлов.

Чтобы избежать паралича вычислений, библиотеки AMR (такие как открытые пакеты AMReX или p4est) каждые несколько шагов перераспределяют узлы сетки между процессорами. Они используют кривые заполнения пространства (Space-Filling Curves, например кривые Гильберта или Мортона), которые сводят сложную 3D-топологию узлов в один длинный одномерный массив. Этот массив легко нарезать на равные порции и разослать процессорам, гарантируя идеальное распределение работы. Благодаря AMR человечество получило возможность моделировать слияние черных дыр (проект LIGO) и работу сложнейших реакторов термоядерного синтеза.

Подробнее

Численные методы в топологическом анализе данных (TDA): устойчивые гомологии

Геометрия скрытых структур в больших данных

В эпоху Big Data исследователи часто работают с наборами данных (Point Clouds), которые представляют собой миллионы точек в пространствах очень высокой размерности (сотни и тысячи координат). Это могут быть профили экспрессии генов, данные функциональной МРТ мозга или финансовые транзакции. Традиционные методы машинного обучения (кластеризация, PCA) опираются на метрику — точные расстояния между точками. Но в многомерных пространствах метрика теряет свою силу, шумы искажают расстояния, и методы дают сбои. Нужен был инструмент, устойчивый к растяжениям, сжатиям и непрерывным деформациям данных.

Этим инструментом стала алгебраическая топология, адаптированная для компьютеров под названием Топологический Анализ Данных (TDA). Топология изучает форму объектов (связность, наличие дыр, туннелей и пустот). Например, топологически кофейная кружка и бублик — это один и тот же объект, так как у обоих есть ровно одна дырка. Главный численный алгоритм TDA называется Устойчивыми гомологиями (Persistent Homology). Он позволяет компьютеру находить глобальные топологические структуры в хаосе дискретных многомерных точек.

Симплициальные комплексы и фильтрация Виеториса-Рипса

Как из разрозненных точек собрать сплошной топологический объект? Программа создает математические связи (симплексы). Точка — это симплекс 0-мерный. Отрезок между двумя точками — 1-мерный. Треугольник (с заливкой) — 2-мерный, тетраэдр — 3-мерный. Алгоритм начинает «раздувать» вокруг каждой точки данных гиперсферу радиуса R.

Если две сферы пересекаются, компьютер проводит между точками отрезок. Если три сферы попарно пересекаются, алгоритм закрашивает их треугольником. Эта структура называется комплексом Виеториса-Рипса. Фишка метода в том, что радиус R медленно увеличивается от нуля до бесконечности (процесс фильтрации). При малом R у нас просто набор точек. При увеличении R точки сливаются в кластеры, образуются циклы (дырки), а при очень большом R все сливается в один гигантский монолит без дыр. Алгоритм отслеживает «рождение» (появление дыры) и «смерть» (закрашивание дыры треугольниками) каждого топологического свойства.

Вычисление штрих-кодов и матричная редукция

С вычислительной точки зрения процесс отслеживания гомологий сводится к тяжелой линейной алгебре. Для каждого шага фильтрации компьютер строит матрицу границ (Boundary matrix), которая описывает, какие вершины образуют какие грани. Элементы этой матрицы лежат в поле вычетов по модулю 2 (цифры только 0 и 1).

Затем к этой гигантской разреженной матрице применяется специальный алгоритм редукции столбцов (аналог метода Гаусса), который приводит матрицу к нормальной форме. Из этой формы программа извлекает так называемые «числа Бетти» (инварианты) и строит финальный результат TDA — штрих-код (Barcode) или диаграмму устойчивости. Если какая-то дыра «живет» долго (длинный штрих-код), значит, это реальная топологическая макроструктура данных (например, циклическая структура в химической молекуле или вирусной РНК). Если дыра появляется и сразу исчезает — это просто шум. TDA произвел фурор в поиске новых материалов и анализе сложнейших нейросетей, давая инженерам возможность буквально «увидеть» форму многомерного массива данных.

Подробнее

Вейвлет-сжатие матриц в вычислительной математике: алгоритмы Добеши и Хаара

Плотно заполненные матрицы: проклятие интегральных операторов

При численном решении интегральных уравнений (например, уравнений Фредгольма) или при использовании метода граничных элементов (МГЭ) в электродинамике инженеры сталкиваются со страшной проблемой. В отличие от дифференциальных операторов, которые связывают только соседние узлы сетки и порождают удобные разреженные матрицы с нулями, интегральные операторы глобальны. Каждая точка на поверхности антенны электромагнитно взаимодействует со всеми остальными точками антенны. В результате получается СЛАУ, матрица которой заполнена на 100% (dense matrix). Если сетка содержит миллион узлов, матрица будет содержать триллион ненулевых элементов, что потребует 8 терабайт оперативной памяти и сделает матричное умножение катастрофически медленным.

На помощь приходит гениальная концепция из теории обработки изображений — вейвлет-сжатие. В 1991 году Бейлким, Койфман и Рохлин (BCR) предложили алгоритм, который радикально решает эту проблему. Идея основана на математическом свойстве вейвлетов (в частности, ортогональных вейвлетов Добеши и Симлета): они блестяще сжимают гладкие сигналы. Большинство матричных элементов интегрального оператора описывают взаимодействие удаленных частей поверхности друг с другом, которое является плавным и гладким, изменяясь без резких скачков.

От плотной матрицы к разреженной через пороговую фильтрацию

Алгоритм вейвлет-сжатия работает следующим образом. Вместо того чтобы вычислять СЛАУ в стандартном (узловом) базисе, мы применяем к матрице системы двумерное дискретное вейвлет-преобразование (ДВП). Эта операция эквивалентна умножению матрицы слева и справа на ортогональные матрицы вейвлет-фильтров. Преобразование можно выполнить очень быстро, за время O(N^2) (или даже быстрее с помощью специальных алгоритмов).

После преобразования матрица меняет свое физическое лицо. Теперь ее элементы — это не интенсивность взаимодействия узлов, а коэффициенты взаимодействия различных вейвлет-масштабов. Поскольку оператор (например, функция Грина) гладкий вдали от сингулярности, подавляющее большинство (более 95-99%) этих новых вейвлет-коэффициентов оказываются микроскопически малыми числами, очень близкими к нулю. Алгоритм просто обнуляет все коэффициенты, которые по модулю меньше заданного порога отсечения (Thresholding). Плотно заполненная матрица волшебным образом превращается в сильно разреженную (Sparse matrix) без потери физической точности решения! Теперь для ее хранения требуется ничтожный объем памяти, а итерационные решатели (например, GMRES) выполняют умножение матрицы на вектор за время O(N log N) вместо O(N^2), превращая нерешаемые электромагнитные задачи в повседневные инженерные расчеты.

Подробнее