Очередные сложности со столкновениями тел. Краткое содержание предыдущих серий:

Сначала мы столкновения реализовывали упругими силами, расталкивающими тела. Недостатки -- силы должны быть большими, чтобы не допускать глубокого проникновения тел друг в друга, а это приводит к жесткости задачи и довольно большим погрешностям интегрирования. Как правило, даже в не очень сложных задачах энергия флуктуировала на уровне 10^-5 что довольно много.

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

Схема отлично работала пока мы не начали ее тестировать в гравитационных полях. Оказывается, если "положить" неупругий шарик на "пол" --плоскость (в поле тяжести), то он начинает (медленно но неуклонно) в нее просачиваться. Причина простая -- на каждом шаге интегрирования шарик приобретает небольшую скорость вертикально вниз, и чуть-чуть смещается тоже вниз. Т.е. к концу шага он оказывается чуть внутри плоскости. Движок collision detection фиксирует взаимопроникновение, движок столкновений отрабатывает физику столкновения, и "отражает" скорость шарика -- теперь он летит вверх. Если столкновение абсолютно упругое, этой скорости точь-в-точь хватает, чтобы к концу следующего шага оказаться там же, где он был в начале предыдущего -- у поверхности пола -- и все начнется с начала. Т.е. шарик будет просто колебаться с периодом 2*dt. Но если столкновение неупругое, то после отражения у шарика будет меньшая скорость, и он потеряет ее в гравитационном поле раньше, чем "выскочит" из взаимопроникновения с плоскостью. И следующий шаг (микропадение-отражение-микроподскок) начнет уже с точки чуть-чуть под поверхностью. Это чуть-чуть будет каждый раз увеличиваться, постепенно шарик будет просачиваться в пол.

Решение было довольно очевидным -- нужно делать коррекцию положения (position correction -- штука, хорошо известная в игровых движках, хоть там и несколько другая физика). Т.е. вместе с модификацией импульсов/моментов тел в ходе столкновения нужно "откатить время" назад, ровно чтобы оказаться в моменте, где объекты еще только-только соприкоснулись. Строго говоря, откатить, конечно, нужно всю систему, но для упрощения жизни мы начали откатывать только пару сталкивающихся тел. Тела перестали проваливаться.

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

И вот, очередной раунд. Создаем стакан (стенки+пол), поле тяжести, +20 шариков, столкновения неупругие. Через несколько секунд шарики скапливаются на дне стакана. И начинают довольно быстро просачиваться друг в друга и в пол.

Почему? Кратные столкновения. Наше приближение -- при коррекции положения откатывать надо всю систему, но мы откатываем только текущую пару тел -- работает удовлетворительно только если в системе достаточно "свободно", если столкновения только парные, если отматывая время для одной пары тел и убирая их взаимопроникновение мы, тем самым, не создаем взаимопроникновение другой пары тел. В случае стакана с шариками на дне это не так. Каждый шарик в контакте с 2-3-4 другими, а то еще и с полом или стенками. Получается, что для данного шарика коррекция положения будет суперпозицией коррекций для всех его столкновений с соседями. По факту, более-менее правильно пройдет только последняя коррекция, остальные превратятся черт знает во что. Вот это "черт знает" и видно на экране :)

Что будет следующей итерацией пока не знаю. Надо думать. Честно рассчитывать кратные столкновения нереально -- во-первых, только из законов сохранения более чем парное столкновение не рассчитать, во-вторых под это дело придется сильно корежить движок collision detection, в третьих сам алгоритм position correction для этого случая неочевиден.

В раздумьях

Решил перед уходом в отпуск протестировать "потолок" скорости расчета, на который можно рассчитывать в джаве. Для примера взял простейшую систему из N материальных точек (без столкновений, без вращений), в однородном поле тяжести. Все координаты развернул в линейные массивы double[], и реализовал стандартный Рунге-Кутты 4-го порядка. Итог: real-time (т.е. когда симуляционное время не отстает от реального) считается примерно 4500 объектов (jre1.6.0_14, -server -XX:+DoEscapeAnalysis, P4 3.2Ггц).

Для сравнения: когда у нас еще не было вращения и столкновений, в real-time обсчитывалось примерно 500 объектов. Т.е. я оцениваю так, что потенциал чисто низкоуровневых оптимизаций (развернуть структуры данных, улучшить cache-locality, b т.п.) у нас еще где-то на 300-500% прироста.

В качестве ориентировки переписал (матерясь -- забыл уже все) код на С++, запустил из под Visual Studio. Включил все оптимизации, которые нашел (__inline, глобальную оптимизацию, оптимизацию скорости, SSE2). Результат (примерный) -- процентов на 5-9% быстрее, чем в джаве

Вчера, наконец-то, вчерне закончил работу над вращением. Теперь у нас есть объекты с честными 3-мя вращательными степенями свободы, тензор момента инерции, и уравнения Эйлера для вращательной динамики. Но это еще только треть дела. У нас есть идеальные абсолютно упругие столкновения таких тел, в результате которых они могут приобретать момент импульса. Причем с точным (в пределах машинной точности) сохранением всех задействованных инвариантов -- энергии, импульса и момента импульса.

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

В планах неупругие столкновения и трение. А пока ухожу в отпуск

Наконец-то раскочегарились работы по физконструктору. Изучаю физику заново -- ощущение, что забыл все. Точнее, саму физику на качественном уровне я еще помню, но вот выписывать уравнения я напрочь разучился. Сижу туплю по полчаса над школьными уравнениями движения. Решение задачи о столкновении шаров заняло часа три и вызвало бурю радости :)

Float vs Double

Тестировал производительность плавающей запятой в java. Обычно я везде использую double, но пока изучал jbullet заметил, что там вся арифметика на float-ах, и мне стало интересно.

1. import junit.framework.TestCase;
2.  
3. public class FloatVsDoublePerfomance extends TestCase {
4.     private static final int[] SIZES = new int[]{ 1000, 10000, 100000, 300000 };
5.     private static final int ITER = 10000;
6.  
7.     public void test() {
8.         for ( int i = 0; i < SIZES.length; i++ ) {
9.             final int size = SIZES[i];
10.  
11.             _testFloat( size );
12.             System.gc();
13.             System.gc();
14.  
15.             _testDouble( size );
16.             System.gc();
17.             System.gc();
18.         }
19.     }
20.  
21.     private static void _testFloat( final int size ) {
22.         final float[] arr = createFloatArray( size );
23.  
24.         System.out.printf( "float[%d * %d]\n", size, ITER );
25.         //warm-up -- cache load, etc
26.         testFloatCompare( arr );
27.  
28.  
29.         long startedAt = System.currentTimeMillis();
30.         //compare
31.         for ( int i = 0; i < ITER; i++ ) {
32.             testFloatCompare( arr );
33.         }
34.         System.out.printf( "\tcompare: %d ms\n", ( System.currentTimeMillis() - startedAt ) );
35.  
36.         //summ
37.         startedAt = System.currentTimeMillis();
38.         for ( int i = 0; i < ITER; i++ ) {
39.             testFloatSumm( arr );
40.         }
41.         System.out.printf( "\tsum: %d ms\n", ( System.currentTimeMillis() - startedAt ) );
42.  
43.         //mul
44.         startedAt = System.currentTimeMillis();
45.         for ( int i = 0; i < ITER; i++ ) {
46.             testFloatMul( arr );
47.         }
48.         System.out.printf( "\tmul: %d ms\n", ( System.currentTimeMillis() - startedAt ) );
49.  
50.         //div
51.         startedAt = System.currentTimeMillis();
52.         for ( int i = 0; i < ITER; i++ ) {
53.             testFloatDiv( arr );
54.         }
55.         System.out.printf( "\tdiv: %d ms\n", ( System.currentTimeMillis() - startedAt ) );
56.     }
57.  
58.     private static void _testDouble( final int size ) {
59.         final double[] arr = createDoubleArray( size );
60.  
61.         System.out.printf( "double[%d * %d]\n", size, ITER );
62.         //warm-up -- cache load, etc
63.         testDoubleCompare( arr );
64.  
65.  
66.         long startedAt = System.currentTimeMillis();
67.         //compare
68.         for ( int i = 0; i < ITER; i++ ) {
69.             testDoubleCompare( arr );
70.         }
71.         System.out.printf( "\tcompare: %d ms\n", ( System.currentTimeMillis() - startedAt ) );
72.  
73.         //summ
74.         startedAt = System.currentTimeMillis();
75.         for ( int i = 0; i < ITER; i++ ) {
76.             testDoubleSumm( arr );
77.         }
78.         System.out.printf( "\tsum: %d ms\n", ( System.currentTimeMillis() - startedAt ) );
79.  
80.         //mul
81.         startedAt = System.currentTimeMillis();
82.         for ( int i = 0; i < ITER; i++ ) {
83.             testDoubleMul( arr );
84.         }
85.         System.out.printf( "\tmul: %d ms\n", ( System.currentTimeMillis() - startedAt ) );
86.  
87.         //div
88.         startedAt = System.currentTimeMillis();
89.         for ( int i = 0; i < ITER; i++ ) {
90.             testDoubleDiv( arr );
91.         }
92.         System.out.printf( "\tdiv: %d ms\n", ( System.currentTimeMillis() - startedAt ) );
93.     }
94.  
95.     private static float[] createFloatArray( final int size ) {
96.         final float[] arr = new float[size];
97.         for ( int i = 0; i < arr.length; i++ ) {
98.             arr[i] = ( float )( Math.random() * 100 );
99.         }
100.         return arr;
101.     }
102.  
103.     private static void testFloatDiv( final float[] arr ) {
104.         for ( int i = 0; i < arr.length; i++ ) {
105.             arr[i] = arr[i] / arr[i];
106.         }
107.     }
108.  
109.     private static void testFloatMul( final float[] arr ) {
110.         for ( int i = 0; i < arr.length; i++ ) {
111.             arr[i] = arr[i] * arr[i];
112.         }
113.     }
114.  
115.     private static void testFloatSumm( final float[] arr ) {
116.         for ( int i = 0; i < arr.length; i++ ) {
117.             arr[i] = arr[i] + arr[i];
118.         }
119.     }
120.  
121.     private static void testFloatCompare( final float[] arr ) {
122.         int _ = 0;
123.         for ( int i = 0; i < arr.length; i++ ) {
124.             if ( arr[i] < 50 ) {
125.                 _++;//просто чтобы компилятор не выкинул сравнение
126.             }
127.         }
128.     }
129.  
130.  
131.     private static double[] createDoubleArray( final int size ) {
132.         final double[] arr = new double[size];
133.         for ( int i = 0; i < arr.length; i++ ) {
134.             arr[i] = Math.random() * 100;
135.         }
136.         return arr;
137.     }
138.  
139.     private static void testDoubleDiv( final double[] arr ) {
140.         for ( int i = 0; i < arr.length; i++ ) {
141.             arr[i] = arr[i] / arr[i];
142.         }
143.     }
144.  
145.     private static void testDoubleMul( final double[] arr ) {
146.         for ( int i = 0; i < arr.length; i++ ) {
147.             arr[i] = arr[i] * arr[i];
148.         }
149.     }
150.  
151.     private static void testDoubleSumm( final double[] arr ) {
152.         for ( int i = 0; i < arr.length; i++ ) {
153.             arr[i] = arr[i] + arr[i];
154.         }
155.     }
156.  
157.     private static void testDoubleCompare( final double[] arr ) {
158.         int _ = 0;
159.         for ( int i = 0; i < arr.length; i++ ) {
160.             if ( arr[i] < 50 ) {
161.                 _++;//просто чтобы компилятор не выкинул сравнение
162.             }
163.         }
164.     }
165. }

* This source code was highlighted with Source Code Highlighter.

тестировал под 1.6.0_14, с опциями -server -XX:+DoEscapeAnalysis -Xmx300m -Xms300m

В итоге, получается, что сложение/умножение примерно в два раза медленнее для double, деление -- примерно на 30%, сравнение мало отличается. Буду теперь для операций, не требующих особой точности, активно использовать float.

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

Например: есть множество движков игровой физики -- часть из них даже бесплатна, и портирована на джаву (JBullet, ODE4j). Но они заточены под максимально быструю симуляцию физики -- т.е. система должна выглядеть похоже на реальную. Основные задачи, которые они решают -- движение большого количества тел в простых силовых полях (чаще всего в равномерном поле тяжести), расчет столкновений, учет связей, диссипация. Самое сложное там -- как раз быстрый расчет столкновений и взаимопроникновений тел, связи и силы реакции. Все остальное делается по остаточному принципу -- как правило используется простейшая интеграционная схема Эйлера, которая, вообще говоря, расходится по энергии, но в присутствии трения становится более-менее стабильной, зато обеспечивает завидную скорость расчета. Демки JBullet, где в real-time просчитываются (и рендерятся!) до полусотни объектов довольно сложной геометрии (на java!) впечатляют. При этом попытка смоделировать идеальный биллиард (без диссипации) на этом движке проваливается -- схема эйлера "взрывается" уже через пару минут симуляции.

С другой стороны находятся чисто научные движки. Но они, как правило, заточены под моделирование какой-то конкретной системы. Т.е. программа с их использованием пишется для решения конкретной задачи или класса задач. Кроме того о real-time речь обычно не идет. Да и на джаву они не портированы -- джава, все-таки, пока еще не особо котируется в high perfomance computing

Нам же нужен движок во-первых универсальный -- т.е. единообразно моделировать большинство задач механики, просто накидал объектов и взаимодействий, и сказал "поехали!" -- и он сам со всем разобрался. Во-вторых, с производительностью, близкой к реальному времени для расчета хотя бы десятка сложных объектов, или до сотни простых (моделировать молекулярную механику -- чем больше, тем лучше). Нужно моделировать всяческие идеализации, используемые в механике: идеальные связи (нерастяжимые невесомые стержни, нити, блоки, цепи), идеальные столкновения. Причем все это должно моделироваться максимально "физически" -- т.е. с как можно лучшим сохранением основных инвариантов уравнений движения (энергии, импульса, момента импульса, и т.п.) на временах моделирования хотя бы порядка 10 минут. Такого пока нигде не видел.

Прикидываю, как буду писать сам.