Поиск причины аномально больших ускорений в моей симуляции N-Body на основе PyOpenCl - PullRequest
Новая
       28

Поиск причины аномально больших ускорений в моей симуляции N-Body на основе PyOpenCl

1 голос
/ 07 февраля 2020

Мы с другом создавали симуляцию N-тела на основе Python с целью моделирования солнечной системы или галактик. Наш код использует PyOpenCl для доступа к параллельным вычислениям OpenCl и PyQtGraphs для построения симуляции. Мы также используем библиотеку Astropy для инициализации объектов в нашем моделировании с правильными скоростями и позициями.

Большая часть кода, который я здесь показываю, находится в разделе cal c .py нашего моделирования. Именно здесь мы передаем нашим ядрам начальные позиции / скорости, извлекаем выходные данные и вводим их обратно в ядра.

Далее мы настраиваем код PyOpenCl. 

        self.cmdqueue = cl.CommandQueue(self.cntxt)
        self.bodytype = np.dtype([("velocity", cltypes.double3),("position", cltypes.double3),("mass", cltypes.double), ("acceleration", cltypes.double3), ("deltax", cltypes.double3), ("trace", cltypes.int)])
        self.bodytype, body_c = cl.tools.match_dtype_to_c_struct(self.cntxt.devices[0], "body", self.bodytype)
        self.bodytype = cl.tools.get_or_register_dtype("body", self.bodytype)

Затем мы создаем программу расчета, которая работает на ядрах. Мы даем каждому ядру один объект, а затем - используя позиции других объектов - он вычисляет новые позиции и скорости для этого объекта. Числовой метод, который мы используем, является интегратором Velocity Verlet (более подробную информацию вы можете найти здесь: https://en.wikipedia.org/wiki/Verlet_integration).

Именно тогда, когда мы впервые реализовали этот Veletity Verlet, наша программа начал выводить большие ускорения. До этого мы использовали метод Эйлера, который, казалось, работал нормально. Для моделирования Земли вокруг Солнца наша единица расстояния - Астрономические Единицы (AU), а наша единица времени - дни (d), 

__kernel void frst_prog(double delta, __global body *in,int len, __global body *outbodies,__global double3* outpositions, __global double3* outpositionsOld, __global double3* outaccelerations, __global double3* outaccelerationsOld)
             {
                // const double G = 6.67408 * pow(10.0, -11.0); // Normal G
                const double G = 1.4881361162906675 * pow(10.0, -34); // Solar System G
                //const double G = 1.42458 * pow(10.0, -22.0); // Galactic G
                int i  = get_global_id(0);
                body current = in[i];

                double3 vOld = current.velocity;
                double3 posOld = current.position;
                double3 oldAcc = findForces(len, in, current) * G;
                current.position = current.position + (vOld * delta) + (oldAcc * (delta * delta * 0.5));
                double3 acc = findForces(len, in, current) * G;
                current.velocity = vOld + ((acc + oldAcc) * (delta * 0.5));

                outbodies[i] = current;
                outpositions[i] = current.position;
                outpositionsOld[i] = posOld;
                outaccelerations[i] = acc;
                outaccelerationsOld[i] = oldAcc;
                }


            double3 findForces(int len,__global body *in, body current) {
                double3 forces["""  +  str(self.numThings) + """];
                for(int j = 0; j<len; j++){
                    body new = in[j];
                    double3 displacement = new.position - current.position;
                    double dist = length(displacement);

                    if(dist == 0.f){
                        //forces[j] = (double3)0;
                        continue;
                    }

                    forces[j] = ((new.mass)/(dist*dist))*(displacement/dist);
                    }

                double3 total = (double3)0;

                for(int j = 0; j<len; j++){
                    total += forces[j];
                }

                return total;
            }""").build()

Затем мы обработаем выходные данные и подготовим следующий расчет. Извиняюсь за длинный список входных данных в launch = self.prg.frst_prog(), некоторые из них используются для отладки программы. calcRepeats = 1 для процесса отладки. 

def update(self, delta, calcRepeats, realRunTime):
        delta = np.double(delta)
        for i in range(calcRepeats):
            launch = self.prg.frst_prog(self.cmdqueue, self.bodies.shape, None, delta, self.bodies_buf.data, np.int32(self.numThings), self.bodies_buf.data, self.outPosBuf, self.outPosOldBuf, self.outAccBuf, self.outAccOldBuf)
            launch.wait()
            realRunTime += delta

        cl.enqueue_copy(self.cmdqueue, self.outpositions, self.outPosBuf).wait()
        cl.enqueue_copy(self.cmdqueue, self.outpositionsOld, self.outPosOldBuf).wait()
        cl.enqueue_copy(self.cmdqueue, self.outaccelerations, self.outAccBuf).wait()
        cl.enqueue_copy(self.cmdqueue, self.outaccelerationsOld, self.outAccOldBuf).wait()
        return self.outpositions, self.outpositionsOld, self.bodies_buf, realRunTime, self.outaccelerations, self.outaccelerationsOld

Чтобы протестировать программу, я инициализировал симуляцию только с Землей и Солнцем. Когда я запускал программу, она рассчитывала смехотворно большие ускорения. «Старая позиция» относится к положению объекта перед расчетом. «Новая позиция» относится к позиции после расчета. «Старый a cc» относится к ускорению до изменения позиции. «New a cc» относится к ускорению, рассчитанному с использованием новой позиции. Как вы можете видеть на картинке, Старый pos и Новый pos идентичны и, хотя Old a cc возможен, New a cc слишком велик.

Распечатки командной строки в самом начале программы

Как мы можем найти причину этих больших ускорений?

...