Open noooway opened 5 years ago
Что есть:
В нескольких местах в программе встречаются величины, которые удобно задавать векторами в трехмерном пространстве. Это импульсы и координаты частиц, а также компоненты полей в узлах 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 и соответствующие операторы тоже уже перегружены, поэтому ими заниматься не придется. Но надо будет оставить функции для векторного и скалярного произведения и некоторые другие.
Что касается класса 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
на модуль с функциями и аккуратно обновить их использование в остальной программе.
Создал класс ParticleArray, которй хранит все поля частиц как однотипные массивы. Должно заметно ускорить симуляцию в питоне. Надо бы добавить тесты производительности... С dtype как бы не возникло проблем с передачей данных на графический ускоритель.
Migration to
numpy
from lists and Vec3d classes should accelerate the code. Besides, pyCUDA and pyOpenCL libraries can transfer numpy arrays to GPU semiautomatically.