Оптимизации вычислений | Computational optimizations

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

One of the most important things in resource intensive computations is to perform calculations not only robust, but also fast. I'd like to show in this post, how effective can be an optimization based on proper physical considerations.

В 3D модели задано множество ячеек (стартовая "песочница" имеет размер 128x256x1), в каждой из которых может условно-произвольным образом меняться химический состав на каждом шаге расчёта (которых сотни и тысячи). Требуется рассчитать фазовый состав системы (сколько расплава и сколько твёрдой фазы) в каждой из ячеек.

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

В силу специфики решаемой задачи и ПО, в двух измерениях из трёх статическое давление практически постоянно. 

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

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

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

Теоретически, ещё большую оптимизацию может дать создание большой look-up таблицы, в которой будут просчитаны все возможные вариации входных параметров с некоторым шагом, чтобы потом интерполировать значения. 

Практически, корректная реализация требует некоторого времени, и, кроме того, полученное решение так или иначе надо будет тестировать, сравнивая с "честной" версией. Другой проблемой является оптимизация памяти. Простая оценка говорит, что только для температур плавления нужно порядка 250 точек по давлению и не менее 10 точек по каждому из 3 независимых химических компонентов (соответственно, при добавлении ещё одного компонента добавляется ещё 2500 значений). Это даёт массив в 150|300 МБ (32|64-битные числа с плавающей точкой), что тоже не очень приятно.

Кроме того, поскольку тогда надо хранить ещё и составы расплава и его количества, это добавляет ещё 150|300 * (3+1) = 600|1200 МБ данных. Конкретный же прирост производительности - это то, что нужно будет тестировать.
Поэтому, если я буду двигаться дальше в сторону оптимизации, я предпочту сглаживать полученные точки простыми эмпирическими полиномами. Однако программирование этого тоже потребует довольно много времени.

Дополнительно учтено (уже касаясь "тонкой" физики) влияния динамической компоненты давления, возникающей при перемещении вещества (простейший пример - закон Бернулли). Простой расчёт показывает, что максимальный градиент температуры по давлению составляет примерно 15 градусов на ГПа, что, при имеющейся оценке максимальной величины стрессового давления в 2 ГПа (как выше, так и ниже статической компоненты), позволяет заведомо учесть это, сравнивая температуру в ячейке не с номинальным солидусом, а беря запас в 50 градусов.

Каков реальный профит от описанной выше процедуры? 

Использование описанного выше подхода позволяет потратить 44 секунды в начале работы программы и от 0.7 секунды на каждый проход цикла (можно заложить 5-10 секунд в случае сложной системы).

Расчёт в каждой точке потребует даже для простейшей модели 1.6 часа на каждом шагу вычислений.

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

We have a 3D model full of cells (e.g. sandbox case has 128x256x1 cells), and each of this cells has it's chemical composition. This composition may change with each computational step (hundreds of them). On each step in each cell there should be calculated amount of liquid and solid phases.

Calculation of a melting diagram requires knowledge of actual pressure and chemical composition, which affects the melting temperatures. In a very simple case, we can take SiO2, FeO, MgO and alkalis.

We can also assume, that at least in two directions static pressure is the same. 

The easiest way to reduce the computational time of such problem is to minimize the number of calculations of a phase diagram. We can compare the actual temperature in each cell with the lowest melting temperature under specific pressure, and then decide whether should we calculate the actual phase diagram, or not.

We can store an array with melting temperatures. Chemical composition of a specific cell may change unpredictably, so we need to store it as well. Finally, we get an array of the lowest melting temperatures (in cells with the most easy-to-melt composition) along the pressure gradient.

During computations we just compare the actual temperature to the stored value and also check whether the actual cell's chemical composition is more easy-to-melt.

Theoretically, we can move even further in out optimizations, and use a look-up precomputed table of values. However, the simplest estimate says that we need around 250 points for the range of pressures, plus at least 10 values at each pressure for each of 3 independent components. 

That means, we need 7500 values, or 150 to 300 Mb of data (32|64 bit floats). And each new chemical component requires more and more RAM.

Data for melt composition and quantity at each point requires 150|300 * (3+1) = 600|1200 Mb. It is really a lot, and performance of that code should be tested independently.
So, if I decide to move further, I rather prefer to select several simple functions to interpolated calculated data and store the coefficients. However, such parametrization also requires time to be implemented.

It is also pretty easy to consider some "subtle" physics. Pressures in a moving material might be lower or higher than hydrostatic values (Bernoulli's law). Straightforward computations give us Clausius-Clapeyron slope at around 15K/GPa, which means, that in the most "optimistic" case with -2 GPa dynamic pressure component it is totally within 50K.

What is the real outcome of all the messy things atop?

This approach requires 44 seconds at program's startup and then 0.7-10 seconds on each timestep.

"Straightforward" way takes 1.6 hours on each computational step.

I think, that shows pretty clear, how efficient are optimizations based on physical models and laws.