muParser

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

Самому писать распознаватель выражений писать не хотелось, и к тому же я усвоил, что существует множество Open Source проектов направленных на решение всевозможных задач. А значит и можно найти готовый и бесплатный распознаватель выражений.

Достаточно быстро нашелся подходящий мне проект под названием muParser. Официальный сайт - http://muparser.beltoforion.de/. Что меня привлекло в этом парсере:

1. Наличие примеров.

2. Наличие документации на сайте.

3. Простота использования.

4. Кроссплатформенность.

Сразу скажу, что muParser я не сравнивал с аналогами. Вопрос высокой эффективности остро не стоял, так как время вычисления выражения muParser'ом существенно меньше временных затрат остальной части программы. 

Здесь я опишу процесс установки в ОС Ubuntu Linux 12.04 64-bit.

1. Скачиваем дистрибутив библиотеки в официального сайта. После скачивания у меня в папке Downloads появился архив muparser_v2_2_3.zip. 

2. Копируем архив во временную папку temp

cp Downloads/muparser_v2_2_3.zip temp/

3. Переходим в папку temp и распаковываем архив 

cd ~/temp

unzip muparser_v2_2_3.zip

4. Переходим в папку извлеченную из архива

cd muparser_v2_2_3

Рекомендую прочитать содержимое файла License.txt. В файле Install.txt описана процедура установки. 

5. Запускаем скрипт configure и указываем в параметрах запуска, что хотим собрать динамически подключаемую библиотеку и примеры.

./configure --enable-shared=yes --enable-samples=yes --enable-debug=no

6 Запускаем утилиту

make

7. Устанавливаем muParser

sudo make install

В папке /usr/local/include появятся заголовочные файлы парсера.

В папке /usr/local/lib появится динамически подключаемая библиотека libmuparser.so

8. Обновляем кэш динамических библиотек

sudo ldconfig

9. Проверяем работоспособность примеров

./sample/example1/example1

example2 написан для Windows

В качестве демонстрации возможностей muParser'а я написал рисовалку графиков для функций одного действительного переменного. Для визуализации использовал freeglut.

Программа требует ввести выражение для функции, далее выводит кривую графика функции. 

Поскольку muParser понимает только основные элементарные функции я зарегистрировал в своей программе первые пять полиномов Лежандра, которые называются P0, P1, P2, P3, P4.

Так как у меня не стояла задача писать универсальный построитель графиков, то программа может вывести только одну кривую на промежутке [-1, 1], подписи осей отсутствуют.

Рисунок 1 - Скриншот терминала

Рисунок 2 - Скриншот построителя графиков. Полином Лежандра 4-й степени

Теперь построим 11 первых слагаемых  разложения функции

в ряд Фурье

Рисунок 3 - Скриншот терминала

Рисунок 4 - 11 первых членов Фурье разложения

Исходники можно скачать по ссылке

https://drive.google.com/file/d/0Byp7G91RaS2aSWdyWURHREhtNmM/edit?usp=sharing. 

Для сборки выполнить в консоли следующее

g++ -O3 main.cpp Legendre.cpp colors.cpp -lmuparser -lGL -lglut -o muParserTest

Текущий релиз 2.2.3 выпущен 22.12.2012. Подозрительно долго не выходило новых версий библиотеки. Тем не менее имеющийся функционал библиотеки может быть полезен при решении определенного круга задач.

О разработке вычислительных программ. Часть 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. Графика помогла мне найти не одну ошибку в алгоритме конвертации сетки, а также расчете некоторых сеточных величин. Одно дело проверить сетку алгоритмически, рассчитать какие-то критерии, но совсем другое увидеть ее своими глазами на экране. Это придает больше уверенности в корректности своего кода.

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

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