English version is below
Вчера вечером после года мучений я нашёл, как убыстрить мой код при экстремальных условиях расчёта примерно в 100 раз. С
одной стороны, я нашёл это, и аз есмь хорошо. С другой стороны, осознание
того, что для этого потребовалось вставить в трёх местах на 6.2
миллионов символов кода цифру "2", меня вымораживает.
Постановка задачи (геология "не на пальцах" в этом посте, физика - в этом) Представьте, что вы сидите в ванне и топите жёлтую пластиковую уточку. Этот гештальт нельзя так просто взять и закрыть: бездушная имитация скотины порывается достичь равновесного уровня - то есть поверхности воды. Скорость процесса зависит от содержимого её внутреннего мира (то есть разности плотностей), вязкости воды и прочего. Теперь мы вносим в студию стробоскоп: именно так, отдельными кадрами, мы видим реальность при численном моделировании. И вот, представим, что в один момент у нас есть кадр, где уточка под водой. Мы посчитаем её скорость, а затем попытаемся угадать, где она окажется на следующем кадре. Если следующий снимок будет через 1 минуту, а ваша ванна не на дне океана, то "в теории" наша уточка улетит куда-то далеко-далеко (потому что мы "не знаем", что её скорость станет 0 при достижении поверхности).
Мы можем научиться опознавать поверхность воды и рассчитывать некоторую "эффективную скорость". Если мы снимаем нашу ванну и забавы с утятоньками в ней с частотой, не достаточной, чтобы за один интервал утка всплыла, то её скорость мы оставляем неизменной: пусть проплывёт столько, сколько сможет, я не жадный. А вот если плёнки у нас мало, и кадры редкие, то мы зачеркнём рассчитанную скорость и напишем сильно меньшую - чтобы как раз в момент снятия следующего кадра уточка достигла поверхности и все были счастливы.
Теперь, чтобы представить более реалистичные условия (с учётом некоторых физических эффектов), представим, что в произвольный момент времени в уточке могут разверзнуться врата ада и накормить её внутреннюю пустоту свинцом. После чего она станет тонуть - то есть направление движения поменяется на противоположное (физически процесс обусловлен изменением химического состава (и плотностей) при плавлении/кристаллизации вещества). И это создаёт большую проблему для расчётов: в "боевых" условиях код за сутки на 4 ядрах продвигался на 1 миллион лет притом, что надо просчитать минимум на 100 миллионов лет. Это очень большая проблема. Три месяца ждать нельзя.
Я потратил много времени на поиск оптимальной схемы (не год, конечно, там было много других проблем, но всё равно). Решение оказалось простым. Используемое значение скорости $v_{pract}$ определялось как:
\[ v_{pract} = \frac{v_{th}}{\max(1,v_{th} \Delta t / Z)}\]
где $v_{th}$ - это теоретическое значение скорости (повторю ссылку на пост с физикой), $\Delta t$ - шаг по времени и Z - расстояние до позиции равновесия. $v \Delta t / Z$ аналогично понятию числа Куранта (Courant number). Так вот, всё решение состояло в том, чтобы написать:
\[ v_{pract} = \frac{v_{th}}{\max(1,2\cdot v_{th}\Delta t / Z)}\]
Поскольку шаг по времени, выбираемый программой, зависит от скорости перемещений субстанций, то предыдущее выражение не позволяло его увеличить всё то время, пока есть плавление. Вторая схема позволяет без потерь точности (при количестве шагов >100, у нас их десятки тысяч) увеличивать шаг по времени.
А я сижу, курю бамбук и вспоминаю лектора по численным методам мехмата МГУ, который сказал золотые слова: "90% успеха численных методов, которые заложены в ноу-хау успешных программных пакетов, - это отнюдь не качественные формулы, а эмпирические численные коэффициенты".
Yesterday night after one year of pain and sorrow I realizes how to speed up my computations by factor of 100. On one hand, I'm glad that it works. On the other hand, it's oppressing to understand that you needed to put three digits "2" into > 6.2 millions of characters to do it.
Problem itself in nuts and bolts (I described geology in another post, and physics is out of there as well). Imagine you are sitting in a bath and you try to sink down a yellow plastic duck. One does not simply perform closure of a gestalt, and there is no exception here: mindless synthetic beast keep trying to reach the equilibrium level at the surface of water. Duck's velocity depends on its inner self (density difference between creature and water), liquid viscosity and so on. Now we take a stroboscope and observe the bath in flashes - that's how do we do it in numerical simulations with discrete time steps. We see duck under the water in one frame. We can calculate its velocity from some formulas, and try to draw its location for the next slide. If next slide is in one minute, and our bath does not stand on the oceanic seafloor, we will "compute" duck's location somewhere very-very in the air above the bath (we don't "know" that velocity becomes zero at the water surface).
We can learn how does water surface work to compute some effective speed value. If we take snapshots of the bath and pleasures with ducks in it so frequently that duck can't reach the surface in one timestep, we can be kind enough not to change its velocity. However, if we take photos rarely, we can just cross out the number and put some value to make duck reaching the water surface exactly at the moment we take the next shot. And all must be happy!
However, situation is a bit more weird. In the real world (in terms of physical effects), at any certain moment of time a portal of the hell can open within the duck and feed its internal emptiness with lead. Dense lead will make the creature sinking down, so direction of velocity vector changes to the opposite one (physically this happens because of changes in chemical composition of melt and crystal phases during solidification or fusing). All these complexities lead to very slow computations: 1 Myr in 24 hours on 4 CPUs. And we need at least 100 Ma, and we can't wait for three months. That's a problem.
I spent a lot of time to figure out a way through. Velocity value used by code $v_{pract}$ is finally determined as:
\[ v_{pract} = \frac{v_{th}}{\max(1,v_{th} \Delta t / Z)}\]
where $v_{th}$ is a theoretical value (here I described how do I calculate it), $\Delta t$ is a timestep, and Z is a distance between current melt location and equilibrium level. $v \Delta t / Z$ is Courant number. The solution was just to put "2" in the this formula:
\[ v_{pract} = \frac{v_{th}}{\max(1,2\cdot v_{th}\Delta t / Z)}\]
Since timestep depends on velocities of substances (determined by code adaptively), the former expression didn't allow to increase it. The latter formula allows to increase it with no changes in precision (is there are mpre than 100 timesteps, and we have thousands of them).
I spent a lot of time to figure out a way through. Velocity value used by code $v_{pract}$ is finally determined as:
\[ v_{pract} = \frac{v_{th}}{\max(1,v_{th} \Delta t / Z)}\]
where $v_{th}$ is a theoretical value (here I described how do I calculate it), $\Delta t$ is a timestep, and Z is a distance between current melt location and equilibrium level. $v \Delta t / Z$ is Courant number. The solution was just to put "2" in the this formula:
\[ v_{pract} = \frac{v_{th}}{\max(1,2\cdot v_{th}\Delta t / Z)}\]
Since timestep depends on velocities of substances (determined by code adaptively), the former expression didn't allow to increase it. The latter formula allows to increase it with no changes in precision (is there are mpre than 100 timesteps, and we have thousands of them).
And now I just try to keep calm and recall some words said by lecturer at math department of Lomonosov MSU: "90% of success of software for numerical simulations are not the formulas itself, but know-hows, which are usually just an empirical coefficients to make them working fast and robust".
Комментариев нет:
Отправить комментарий