Previous Entry Share Next Entry
Численное решение задачи N тел
evatutin
Давно хотелось сделать для параллельного программирования что-то наглядное и хорошо распараллеливаемое, решил попробовать задачу n тел. Выбранная постановка задачи довольно сильно упрощена, однако в работе нет претензии на научную новизну, есть на наглядный пример, который, я надеюсь, будет интересен студентам. Допущения:


  1. рассмотрение двумерного пространства (под N-мерное код элементарно переделывается копипастом в случае чего);
  2. моделирование по классическим законам (законы Ньютона, не учитываем эффекты ОТО);
  3. моделирование движения точечных объектов (не учитываем форму, вращение);
  4. фиксированный шаг интегрирования (и это в жесткой задаче — закидали бы меня специалисты тухлыми помидорами :) );
  5. абсолютно неупругое столкновение тел при сближении ближе порога (опять таки тут есть над чем подискутировать);
  6. не учитывается потеря энергии на излучение гравитационных волн.


В математической постановке решено не решать численно систему дифуров (например, методом Рунге-Кутты), а, зная шаг по времени , ограничиться расчетом приращений для скоростей

и координат

Сила, действующая на тело, считается как сумма сил притяжения других тел

ускорение — по второму закону Ньютона

а сами силы — по закону всемирного тяготения

Пространство у нас Евклидово, поэтому

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

заметить, что расстояния от головы до хвоста такое же, как и от хвоста до головы).
В результате проделанных алгоритмических оптимизаций (сокр. АО, а по большей части математических) можно получить более менее быстро работающий высокоуровневый вариант расчетного кода. Далее, чтобы убедиться в том, что код, генерируемый компилятором Delphi для сопроцессора, из рук вон плох (кто сомневается — в дизассемблер либо ко мне на лекции), можно перейти к микроархитектурным оптимизациям (сокр. МО), написать расчетный код вручную под FPU и получить прибавку в скорости.
А далее в игру вступают векторные расширения. Наиболее правильно моделировать динамику системы с использованием вещественных чисел повышенной точности (Extended или long double), но их векторные расширения, к сожалению или к счастью, не поддерживают. Значит возьмем числа двойной точности (Double или double) и для простоты реализуем скалярный вариант SSE2-кода (с командами вроде ADDSD, MULSD и т.д.). Несмотря на то, что вроде бы вычисления по сравнению с FPU-кодом те же самые, у меня на Core 2 Duo E6300 он работает чуть быстрее (на процессорах от AMD может быть по другому). Если получился скалярный SSE2-код, получится и векторный, надо лишь придумать, как (по или по ) загружать данные на обработку парами (команды ADDPD, MULPD и т.д. работают с двумя double'ами), чтобы не выйти за пределы массивов. В данном случае можно воспользоваться диагональными элементами , главное не забыть про возможность деления на ноль, т.к. само от себя тело находится на нулевом расстоянии.
Ну а нельзя ли теперь понизить точность вещественных вычислений до одинарной? Понимая весь риск накопления численных ошибок, можно. В таком случае в векторном SSE-коде (скалярный мне было писать лень, т.к. выигрыша в скорости он не даст) можно вести обработку не парами, как в SSE2, а четверками операндов (ADDPS, MULPS и т.д.), что также положительно сказывается на скорости расчета. Правда тут диагональными элементами уже не ограничиться и придется доопределять число тел до кратного 4 телами нулевой массы, которые не будут никого притягивать. В итоге выглядит это как-то так:

  { Изменение скорости }
  { fParticlesArr[I].V.Y := fParticlesArr[I].V.Y + fAy[I] * dt }
  movq  xmm2, qword ptr [edi+ecx*8]  // xmm2 = fAy[i]
  mulsd xmm2, xmm0                   // xmm2 = fAy[i]*dt
  addsd xmm2, qword ptr [ebx+$18]    // xmm2 = fAy[i]*dt + Vy_i = Vy_i(new)
  movq  qword ptr [ebx+$18], xmm2

  { Изменение координат }
  { fParticlesArr[I].Coords.Y := fParticlesArr[I].Coords.Y + fParticlesArr[I].V.Y * dt }
  mulsd xmm2, xmm0                   // xmm2 = Vy_i(new)*dt
  addsd xmm2, qword ptr [ebx+8]      // xmm2 = Vy_i(new)*dt + y_i
  movq  qword ptr [ebx+8], xmm2


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



Разумеется это не все, т.к. еще есть AVX, многоядерность (программируемая через thread'ы или OpenMP), видеокарты с CUDA/OpenCL, кластеры с MPI. Все это можно попробовать использовать в задаче, но чуть позже, свободного времени много не бывает.
В программе имеется возможность выбора численных параметров моделирования, слежения за траекториями группы тел и различные сценарии моделирования, на которые у меня и моих друзей хватило фантазии (спасибо Максим Игоревич, Владимир Алексеевич и др.).



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



Если будут интересные предложения — welcome в обсуждение! Только сразу оговорюсь, что наш Млечный путь пока моделировать не буду, т.к. там очень много тел и без дальнейшей оптимизации там делать нечего. На данный момент порог для моделирования — несколько тысяч тел (хотя если терпения и вычислительного времени много...). Да и с темным гало надо что-то делать... Одним словом, приятного обучения!

Рабочий бинарник, 516 КБ

  • 1

Guten Tag

(Anonymous)
Mit freundlichen Grüßen!

а можно посмотрель на программу ..?

screen 21
dim as double j=230e9 'window
dim shared as double t,r
dim shared as double Rdif
dim shared as integer n,m
window(-j,-j*.75)-(j,j*.75)
dim shared as double x (9):dim shared as double y (9)
dim shared as double vy (9,9):dim shared as double vx (9,9)
dim shared as double mc (9)
dim shared as double dx,dy ,a,df,dfy,dfx,dt
dt=100
mc(1)=1.9891e30 'Sun /static place
y(1)=0

mc(2)=5e24 'EArth
y(2)=150e9
vx(2,1)=29783

mc(5)=1e22 'EArth moon
y(5)=150e9+384e6
vx(5,1)=29783+940

mc(3)=6.4185e23 'MARS 6.4185e23
y(3)=-225e9
vx(3,1)=-24130

mc(4)=1.8986e27 'YOUPiteR 1.8986e27
y(4)=-778e9
vx(4,1)=13070

mc(6)=1.5e22 'Ganimed
y(6)=-778e9+1070e6
vx(6,1)=13070+10880


sub vals ()
dim as double dx,dy ,a,df,dfy,dfx,dt,xm,ym
dx=x(m)-x(n):dy=y(m)-y(n)
df=6.67e-11*mc(m)/(dx*dx+dy*dy)
a=atn(dy/abs(dx)):dfy=df*sin(a):dfx=df*cos(a)*sgn(dx)
dt=100
vy(n,m)=vy(n,m)+dfy*dt:vx(n,m)=vx(n,m)+dfx*dt
y(n)=y(n)+vy(n,m)*dt:x(n)=x(n)+vx(n,m)*dt

end sub


for t=0 to 31556925.9*20 step 1
for n=2 to 6
for m=1 to 6
if n=m then goto label1
pset(x(1),y(1)),2
pset(x(3),y(3)),4
pset(x(4),y(4)),6
pset(x(5),y(5)),7
vals () ' step numeric intr
pset(x(2),y(2)),13
pset(x(6),y(6)),15

label1:

next m
next n

Rdif= 0.001 *Sqr ( ((y(5)-y(2)) ^2 ) + ( (x(5)-x(2))^2) ) 'earth -moon distance
locate 1.1
print using "########### KM earth -moon dist";Rdif

Rdif= 0.001 *Sqr ( ((y(6)-y(4)) ^2 ) + ( (x(6)-x(4))^2) ) 'yupiter- ganimed distance
locate 2.1
print using "########### KM yupiter- ganimed dist ";Rdif


next t

end
есть что то похожее на это? чтобы можно было самому редактировать парметры и давать некоторое искусственное ускорение и проверьте ..почему то улуны круговая скорость ~940 а должна быт 1030
для проверки внес ганимед вокруг юпитера так там все как по теории

sub vals ()
dim as double dx,dy ,a,df,dfy,dfx,dt,xm,ym
dx=x(m)-x(n):dy=y(m)-y(n)
df=6.67e-11*mc(m)/(dx*dx+dy*dy)
a=atn(dy/abs(dx)):dfy=df*sin(a):dfx=df*cos(a)*sgn(dx)
dt=1000
vy(n,m)=vy(n,m)+dfy*dt:vx(n,m)=vx(n,m)+dfx*dt
y(n)=y(n)+vy(n,m)*dt:x(n)=x(n)+vx(n,m)*dt

end sub

ка обойтись без синусов и косинусов при разложении силы на составляющие .. кста я эту прогу написал еще в 1998 где то)) тогда еще под турбобейсик

  • 1
?

Log in

No account? Create an account