Дисциплина Computational fluid dynamics(CFD) или, на русском языке, Вычислительная гидродинамика изучает поведение различных потоков, в том числе вихревых. Это и моделирование цунами, и лавовых потоков, и выбрасываемых из жерла вулкана камней вместе с лавой и газами и многое другое. Посмотрим, как можно использовать совместно MantaFlow и ParaView, реализовав на встроенном в MantaFlow языке Python необходимые функции конвертации данных. Как обычно, исходный код смотрите в моем GitHub репозитории: MantaFlow-ParaView.

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

Введение

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

Как и в случае с методами спутниковой интерферометрии, см. предыдущую статью Геология XXI века как наука данных о Земле, методы моделирования потоков позволяют узнать многое о происходивших процессах именно узнать, а не предположить. Действительно, моделируя лавовые потоки на палеорельефе (построенным, к примеру, путем решения обратной задачи геофизики), мы можем сравнить полученную модель с реально существующими [застывшими] лавовыми потоками и полями, поскольку они достаточно прочные и могут хорошо сохраняться многие и многие миллионы лет. Построенная модель позволяет изучить залегание различных слоев там, где у нас нет геологической информации, зачастую, один и то же вулкан извергался многократно из разных жерл, при этом потоки лавы разных извержений [и разного состава] могут перекрываться. Кроме того, эта же модель покажет неувязки между используемой моделью палеорельефа и существующими лавовыми проявлениями и позволит внести уточнения. То есть, вместо сложно формализуемой геологической интуиции, мы можем работать с моделью в том числе, изменять параметры и оценивать, насколько разные геологические предположения вообще разумны. Поскольку я сам не геолог, а физик, для меня путь моделирования представляет вполне понятный интерес. В результате, геолог может точнее оценить возможные области залегания полезных ископаемых и их потенциал разумеется, в геологическом исследовании вовсе не ставится задача обойтись без геолога, но, как и везде, ценность результатов согласуется с известным принципом информатики: "Мусор на входе мусор на выходе", поэтому каждая возможность уточнить, заверить и дополнить имеющиеся данные в буквальном смысле на вес золота (или нефти, или воды,...).

MantaFlow

Как уже сказано выше, MantaFlow это хорошо известный и поддерживаемый программный пакет, а также, на его основе есть несколько интересных производных проектов: с большим акцентом на использование машинного обучения PhiFlow, с возможностью вычисления параметров процесса по набору снимков reconstructScalarFlows и другие.

MantaFlow позволяет моделировать ламинарные и турбулентные потоки, в том числе, пламя и клубы дыма и потоки воды и многое другое. Для примера я построил несколько моделей, чтобы оценить качество моделирования. Например, эта модель потока воды выглядит отлично и даже видны формирование и разрыв воздушных пузырей в потоке:

А это модель заполнения потоком воды заданного рельефа:

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

Добавим рельеф в MantaFlow

Для меня самый удобный вариант это использование ParaView с моим расширением N-Cube ParaView plugin for 3D/4D GIS Data Visualization для построения 3D модели рельефа в ParaView на основе NetCDF или GeoTIFF данных и сохранение нужного участка в формате OBJ для его использования в MantaFlow. Поскольку MantaFlow умеет этот формат загружать и работать с ним в безразмерных координатах, нам потребуется лишь указать нужный размер в безразмерных координатах (скажем, 100% по горизонтальным координатам и 25% по вертикальной чтобы было достаточно места для моделирования столба дыма) и сохранить параметры преобразования для экспорта результатов в физических координатах. Вот скрипт репозитория с реализацией соответствующей функции: mesh2manta.py

Сохраним результаты моделирования в MantaFlow для ParaView

Поскольку мы задаем исходное пространство моделирования в физических координатах (файл OBJ и коэффициенты его масштабирования), у нас есть все необходимое, чтобы и результаты сохранить в физических координатах. По умолчанию, MantaFlow сохраняет безразмерные результаты в формате сжатых массивов Numpy, поэтому мы добавим сохранение в формат с поддержкой физических координат (VTK), см. скрипт репозитория npz2vtk.py. Добавлю, что в скрипте создается массив xarray: N-D labeled arrays and datasets in Python, из которого одной командой можно сохранить данные в формате NetCDF и некоторых других.

Визуализация в ParaView

Как мы уже рассмотрели в предыдущих статьях (с примерами), ParaView поддерживает работу с сериями данных, так что мы можем работать с 4D данными например, в виде 3D анимации. Вот пример анимации вулканического дыма из серии файлов VTK, экспортированных из MantaFlow:

Модели высокой детализации

Увеличение детальности моделей требует и больше ресурсов для их построения. Если модели на сетке 64х64х64 вычисляются за несколько минут на ноутбуке, то при удвоении разрешения по каждой координате время увеличивается в 8 раз (третья степень двойки).

Ниже показана намного более детальная модель на примере турбулентного торнадо с сайта проекта MantaFlow:

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

Заключение

Хотя именно о построении моделей мне и хотелось написать ранее, но пришлось начинать с серии статей про визуализацию данных и потом постепенно переходить уже непосредственно к моделированию. В самом деле, если не рассказать заранее, как хотя бы посмотреть результаты, получится не интересно. Хотя и можно использовать MantaFlow в Blender, но с последним я сам не работаю и не уверен, многие ли из читателей с ним достаточно знакомы. Так что пусть будет как есть знакомство с ParaView ранее и рассказ про MantaFlow и ParaView, а те читатели, кому этого захочется смогут попробовать MantaFlow и Blender.

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

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

О качестве программного обеспечения и инструментарии

В Typeable мы придаем огромное значение качеству программного обеспечения и прикладываем все усилия, чтобы обеспечить это качество. В настоящее время мы искореняем ошибки следующими способами:

1. Анализ и создание спецификаций
2. Устранение простых ошибок с использованием системы типов Haskell
3. Стандартные юнит-тесты и интеграционные тесты
4. Непрерывная интеграция
5. Обязательные ревью кода
6. Тестирование на стендах, проводимое QA инженерами
   (мы используем Octopod для оптимизации процесса разработки и QA)
7. Тестирование в pre-production среде
8. Ведение логов и контроль ошибок на этапе эксплуатации

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

Один из способов сокращения этих затрат заключается в выявлении ошибок на ранней стадии. По грубой оценке, если система типов обнаруживает вашу ошибку, это происходит в течение 30 секунд после сохранения файла. Если ошибка найдена во время CI, получение информации об ошибке займёт до 30 минут. Кроме того, после исправления ошибки вам придётся ждать еще 30 минут, пока CI не отработает снова.

Чем дальше вы продвигаетесь по цепочке, тем длиннее становятся эти перерывы и тем больше ресурсов уходит на исправление ошибок: чтобы достигнуть этапа QA-тестирования могут потребоваться дни, после чего инженер-тестировщик еще должен будет заняться вашей задачей. А если на этом этапе будет обнаружена ошибка, то не только тестировщикам придется снова провести тесты после исправления ошибки, но и разработчики опять должны будут пройти все предыдущие стадии!

Итак, каков самый ранний этап, на котором мы можем выявить ошибки? Удивительно, но мы можем существенно повысить шансы на выявление ошибок ещё до того, как будет написана первая строка кода!

Alloy выходит на сцену

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

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

Давайте приведем пример. Недавно у нас возникла неприятная проблема со следующим куском кода:

newAuthCode  :: (MonadWhatever m)  => DB.Client  -> DB.SessionId  -> m DB.AuthorizationCodenewAuthCode clid sid = do  let codeData = mkAuthCodeFor clid sid  void $ DB.deleteAllCodes clid sid  void $ DB.insertAuthCode codeData  return code

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

Как это получилось?

Моделирование

Проблема, указанная выше, представляет собой хороший пример задачи для Alloy. Давайте попробуем представить ее, построив модель. Обычно мы начинаем моделирование конкретной проблемы с описания нашего представления об операциях newAuthCode для Alloy. Иными словами, необходимо сначала построить модель операций, затем дополнить ее, построив модель базы данных и привязав поведение базы данных к операциям.

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

Процесс, описанный в приведенном фрагменте кода, имеет две интересующие нас части. Он производит удаление в некоторый момент времени, а затем вставляет новый токен в другой момент времени. Вот одна из моделей Alloy, используемая для описания этого поведения:

open util/time  // Импортируем предопределённые объекты Timesig Operation       // У нас есть операции...  { delete : Time   // ...которые удаляют в какой-то момент времени  , insert : Time   // ...и производят вставку в какой-то другой  }  { lt[delete,insert]  // Удаления происходят до вставок    lt[first,delete]   // По техническим причинам в первый                        // момент времени ничего не происходит  }  run {some Operation} for 4 // Показать произвольный пример модели                             // с <= 4 операциями

Приведенная выше модель описывает систему абстрактных объектов и отношений между этими объектами. Выполнение модели далее приведет к образованию произвольного пространства, содержащего некоторые операции Operation, которые распределяются по определенным правилам.

Если вы хотите проследить этот процесс, скачайте alloy и скопируйте в него приведенный выше фрагмент кода. Затем нажмите 'execute' и 'show', чтобы получить модель следующего вида:

Чтобы Alloy показал другие модели, можно нажать 'next'.

Вот один из таких случайных экземпляров, представленный в виде таблицы отношений (нужно несколько раз нажать 'next и выбрать вид 'Table'):

this/OperationdeleteinsertOperation    Time Time   Operation удаляет в момент Time и   вставляет в момент TimeOperation    Time Time   Operation удаляет в момент Time и   и вставляет в момент Time                                                 ОЙ!

Как правило, на данном этапе мы начинаем моделировать таблицы базы данных и семантику операций, но оказалось, что Alloy уже смог показать, почему наши логи содержат нарушение требований!

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

Проблема найдена!

Давайте её исправим!

Как только я нашел причину проблемы, я написал для нее следующее исправление.

code <- run $ do  handleJust constraintViolation    (launchPG $ selectCodeForSession clid scope sid    (launchPG . pgWithTransaction $ newAuthCode clid scope sid)

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

Будет ли это работать сейчас?

Давайте быстро построим модель Alloy для нашего исправления, чтобы проверить его корректность:

open util/time // Импортируем Timesig Token {} // Объекты с названием Tokenone sig DBState // База данных с токенами {userToken : Token lone -> Time}    // В БД не более одного токена в каждый момент врвмени    // (т.к. ограничения БД не позволяют хранить больше одного)sig Operation {   delete : Time , insert : Time , select : Time // Наши операции теперь могут выполнять select}{  lt[first,delete]   // Ничего не происходит в первый момент времени                     // по техническим причинам  lt[delete,insert]  // Первой выполняется операция delete  lte[insert,select] // select выполняется после или во время insert'а  no userToken.(insert.prev) // Если вставка сработала (т.е. таблица  => insert = select         // была пустой во время выполнения),                             // получаем значение в тот же самый                              // момент времени (т.е. у нас запрос                             // 'INSERT RETURNING').                             // В противном случае вызываем обработчик                             // исключения, и select выполняется чуть позже}

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

Затем мы переходим к интересному аспекту модели, который заключается в описании взаимодействия между операциями и состоянием базы данных. К счастью, для нашей модели это довольно просто:

fact Trace {                           // Факт Trace описывает поведение системы all t : Time - first | {              // на всех шагах, кроме первого:   some delete.t => no userToken.t       // Если происходит удаление, таблица пуста   some insert.t => some userToken.t     // Если происходит вставка, таблица не пуста   no delete.t and no insert.t           // Если не происходит ни вставок, ни удалений,    => userToken.t = userToken.(t.prev)  // таблица не меняется  }}

То есть мы описываем, как состояние базы данных изменяется в зависимости от некоторых происходящих событий.

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

Давайте примем это за утверждение и попросим Alloy проверить его.

assert selectIsGood {         // То, что мы хотим проверить all s : Operation.select |   // Всегда, когда выполняется select,  some userToken.s            // в базе присутствуем токен}check selectIsGood for 6 // Проверить, что selectIsGood всегда истинно

К сожалению, запуск этой проверки дает нам следующий контрпример:

DBState userToken   DBStateTokenTime                 Token находится в БД в моменты Time и Time              Time                TokenTime   Token в БД в момент Time.                                    Токены есть в таблице только                 моменты Time, Time и Time                 Заменит, что в момент                 Time токенов нет!Operation     deleteinsertselectOperation     TIME Time TimeOperation     Time Time TIME    Таблица пуста в момент Time и     select не работает для Operation!Operation     Time Time Time                                               Это моменты времени, когда                происходят соответствующие действия

На этот раз контрпример усложняется. К неудаче нашего исправления приводят три операции, происходящие одновременно. Во-первых, две из этих операций производят удаление и очистку базы данных. Далее, одна из этих двух операций вставляет новый токен в базу данных, в то время как другая не может сделать вставку, поскольку в таблице уже есть токен. Неудачная операция начинает выполнять обработку исключений, но еще до ее завершения запускается третья операция и снова очищает таблицу, в результате чего select в обработчике исключений теперь не может ничего выбрать.

Итак, предлагаемое исправление, которое было проверено на согласование типов, протестировано, прошло интеграцию и проверку коллегами, оказалось ошибочным!

Параллельная обработка легко не дается.

Мои выводы

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

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

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

Где взять Alloy?

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

• https://alloytools.org/ < Здесь можно скачать
• https://alloy.readthedocs.io/en/latest/ < Эта документация намного лучше официальной
• https://mitpress.mit.edu/books/software--abstractions-revised-edition < книга об Alloy, написанная понятным языком. В ней описываются большинство основных схем, с которых можно начать.
• https://alloytools.org/citations/case-studies.html < Здесь авторы Alloy перечисляют способы использования своего инструмента. Этот перечень довольно полный и содержит разнообразные примеры использования.

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