Migrate code to numpy

[noooway]

Migration to numpy from lists and Vec3d classes should accelerate the code.
Besides, pyCUDA and pyOpenCL libraries can transfer numpy arrays to GPU semiautomatically.

[noooway]

Что есть:

В нескольких местах в программе встречаются величины, которые удобно задавать векторами в трехмерном пространстве. Это импульсы и координаты частиц, а также компоненты полей в узлах PIC-сетки или внешних полей. С этими величинами делается довольно много арифметики (см. класс ParticleToMeshMap и . boris_update_momentum в Particle. Сейчас трехмерные вектора задаются объектами класса Vec3d. Это осталось от переписывания C++ версии. Думаю, что было бы лучше сделать реализацию этих классов на numpy - арифметика ускорится; также numpy это де-факто стандарт. Кроме того, есть примеры использования GPU библиотек в питоне ( см. https://github.com/fizmat/compare-gpu-lib или https://colab.research.google.com/drive/1Pt8SRp-3zwuAE2DZot3odG-8LgrpoHTg ). Эти библиотеки поддерживают загрузку/выгрузку numpy-массивов в/из память видеокарты. Это будет проще, чем возиться с массивами объектов.

Что имеет смысл сделать:

Заменить массивы Vec3d на numpy-массивы. Возможный вариант - использовать dtype из numpy.
Например, определить vec3d = numpy.dtype( float, float, float ). Вместо класса Vec3d оставить только модуль, где определить нужные функции - они будут обертками для работы с numpy-массивами.

Т.е. заменить

from math import sqrt

class Vec3d():
    def __init__(self, x, y, z):
        self.x = x
        self.y = y
        self.z = z

    @classmethod
    def zero(cls):
        return cls(0.0, 0.0, 0.0)

    def length(self):
        return sqrt(self.x ** 2 + self.y ** 2 + self.z ** 2)

    def negate(self):
        return Vec3d(-self.x, -self.y, -self.z)

    .....

на что-то типа

from math import sqrt

vec3d_dtype = numpy.dtype( [ ('x', float), ('y', float), ('z', float) ] ) # или как-то так

def new(x, y, z):
    return np.array( [x, y, z], dtype=vec3d_dtype )  # или как-то так

def zero():
    return np.zero(dtype=vec3d_dtype)  # или как-то так

def length(vec):
    return np.linalg.norm(vec)

    .....

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

[noooway]

Что касается класса Particle, для него также можно определить отдельный dtype.
Можно сделать particle = np.dtype( [ ('id',int) , ('charge', float), ('mass', float), ('position', vec3d), ('momentum', vec3d ) ] ) . Но у меня осталось впечатление, что numpy не всегда очевидно работает с dtypами, содержащими другие dtypы. Поэтому лучше сделать particle dtype плоским: particle = np.dtype( [ ('id',int) , ('charge', float), ('mass', float), ('position_x', float), ('position_y', float ), ..... ] ) .

Numpy вроде позволяет довольно гибко индексировать массивы, поэтому останется возможность работать с компонентами как с векторами, что-то типа: Vec3d.norm( particles_array[10][ ('position_x', 'position_y', 'position_z') ] ) .

Также нужно будет заменить классе Particle на модуль с функциями и аккуратно обновить их использование в остальной программе.

[fizmat]

Создал класс ParticleArray, которй хранит все поля частиц как однотипные массивы. Должно заметно ускорить симуляцию в питоне. Надо бы добавить тесты производительности... С dtype как бы не возникло проблем с передачей данных на графический ускоритель.
#54