О разработке вычислительных программ. Часть 1.

Хочу поделиться своими впечатлениями о разработке своей CFD программы.

Поставили мне интересную задачу – смоделировать движение тела в вязкой жидкости в двумерной и трехмерной постановках. Давайте разберемся какие ограничения можно ввести для уменьшения трудозатрат и какую математику можно будет использовать для описания исследуемого процесса.

Ограничения

1. Тело движется в безграничном объеме жидкости.

2. Тело не деформируется.

3. Тело только одно.

4. Жидкость несжимаемая.

Математическая модель.

1. Введенные ограничения позволяют пользоваться подвижной системой координат жестко связанной с телом, и все уравнения будем записывать в этой же системе координат.

2. Движение тела подчиняется законам Ньютона.

3. Движение жидкости подчиняется уравнениям Навье-Стокса.

Вот тут и начинается самое интересное. Как реализовать данную мат. модель? С чего начать?

Для начала я ограничился двумерной постановкой.

По отдельности уравнения движения тела и уравнения Навье-Стокса, особенно в двумерном случае, решаются без особых проблем. Однако сейчас они связаны друг с другом, что на порядок повышает вычислительную сложность задачи.

Для сокращения трудоемкости работы можно воспользоваться некоторыми открытыми кодами. Неплохим вариантом представляется OpenFOAM. Но этот пакет плохо подходит для параллельных вычисления. В частности, для одной из задач разбивая вычисления на 4 процесса я получил ускорение только в 1,5 раза. Это несерьезно. Любой вычислительный пакет пишется для определенного класса задач, для которых он является достаточно универсальным. Решение нестандартных задач может быть неэффективным. Используя готовый пакет, ты зависишь от чужой реализации, и не всегда можешь что-либо в ней изменить.

Поэтому решил писать все в нуля, чтобы полностью знать что, где и как реализовано, заточить свой вычислительный пакет под свою конкретную задачу, эффективно реализовать распараллеливание.

Разработку веду на языке C++, так как привык к нему, а еще он есть на всех платформах, можно работать с OpenMP, MPI, Cuda. Некоторые утверждают, что этот язык умер, однако разработки вычислительных пакетов ведутся на нем, а не на Java или .NET. Так что утритесь критики.

Стоял выбор между функциональной парадигмой и объектно-ориентированной. Чаша весов упала в сторону ООП. Объектно-ориентированный код для меня получается более понятным, а поскольку я буду основным пользователем своего же кода, то и мне решать каким он должен быть. В дальнейшем я не исключаю перехода к другой парадигме.

Следующим важным элементов решения является выбор типа сетки. Выбор подвижной системы координат и отсутствие деформации тела позволяет использовать фиксированные сетки.

Геометрия тела может быть сложной и мне представляется более удобным использовать единые неструктурированные сетки. Для двумерного случая ячейки будут либо треугольными, либо четырехугольными.

Задача построения сеток является достаточно трудоемкой в плане реализации, поэтому я решил использовать SALOME Mesh и писать конвертер из dat-файла в свой формат.

Далее необходимо было выбрать метод дискретизации пространственных производных. Выбор остановился на разрывных методах Галёркина (DGM) и на методе конечных объемов (FVM).

DGM дает возможность использовать неструктурированные сетки и локально повышать порядок точности. Для ознакомления реализовал схемы 9-го порядка точности для системы жестких ОДУ и уравнения одномерной стационарной теплопроводности с переменными коэффициентами. Однако очень тяжело этот метод программировать с нуля, особенно если не знаком с ним. И временно от него решили отказаться. К DGM вернемся позже.

Остается FVM, с которым более менее все ясно и можно легко получить хотя бы 2-й порядок точности, которого вполне хватает для большинства расчетов.

Подошли к этапу проектирования кода и отладки.

Сетка – ключ к решению задачи и от нее зависит очень много. Поэтому ошибки в проектировании класса сетки могут стоить очень дорого. По-моему мнению лучше хранить избыточные данные, но чтобы они были всегда под рукой.

Некоторые идеи по проектированию класса сетки я нашел при работе с пакетом OpenFOAM и в книге Löhner R. Applied Computational Fluid Dynamics Techniques: An Introduction Based on Finite Element Methods

Я определил четыре основных класса для двумерной сетки: узел сетки Point, ребро Edge, ячейка Cell, граница Boundary. Плюс вспомогательный класс для проверки качества сетки, несколько структур и функций для отбора граничных граней.

Любой узел моей сетки знает свои координаты, исходящие из него грани и окружающие его ячейки.

Пары узлов формируют ребра. Ребро знает номера своих узлов, номера ячеек, соприкасающихся через это ребро, свою длину, нормаль и центр. Ячейка с меньшим номер считается хозяином ребра, ячейка с большим номером считается соседом, нормаль будет внешней по отношению к хозяину. Граничные грани соседей не имеют.

Ячейка знает номера своих узлов, номера граней, номера соседей, площадь и центр масс. Узлы упорядочены против часовой стрелки, это требование существенно и учитывается при расчете площади ячейки.

Граница знает свое имя, количество граней и индекс первой грани в массиве граней.

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

Следующая статья будет посвящена реконструкции сеточных функций.

Продолжение следует...