Я пытаюсь написать имитацию N-тела из сетки частиц.В таком моделировании потенциальное поле находится путем решения уравнения Пуассона с использованием преобразований Фурье.Я следил за выступлением Андрея Кравцова (http://astro.uchicago.edu/~andrey/talks/PM/pm.pdf),, но слайд 15 меня смутил. Пока что я назначил плотности для 3d-сетки по позициям частиц, и Фурье преобразовал сетку плотности. Следующий шаг -вычислите функцию Грина в пространстве Фурье и умножьте ее на сетку плотности, преобразованную Фурье, а затем примените обратное преобразование Фурье к реальному пространству, чтобы получить потенциальную сетку.и, в частности, вектор k-пространства.
Итак, для вычисления функции Грина в пространстве Фурье мне нужны оси Фурье, обычно называемые векторами k-пространства k_x, k_y, k_z. При использовании слайда это должно быть 2 * pi* (k, l, m) / N_g для компонентов k, l, m, где N_g - количество ячеек сетки. До сих пор я пробовал с этими компонентами, работающими с 0, + 1, + 2, ...,N_g. And -N_particle / 2, ..., + N_particle / 2 и несколько других итераций. Единственное, что дало разумные результатыts (можно увидеть кластер в срезе плотности, спроецированный на тот же срез потенциального поля), был с использованием numpy.fft.freq в Python для конкретных значений разрешения / интервала выборки.Однако любое выбранное мной разрешение (например, L / N_g, N_p / N_g, 2pi / N_g и т. Д.) Неправильно масштабировалось с размером блока L, числом ячеек сетки или числом частиц и больше не работало, например, для большего числаячейки сетки.
Мой вопрос:
Как мне определить мои векторы в k-пространстве (т. Е. Оси Фурье в обратном пространстве) для моделирования с помощью одного направления в одном поле?размер L, количество ячеек сетки N_g и количество частиц N_p?
Я должен добавить, что все положения и скорости частиц указаны в единицах кода, как определено в первых нескольких слайдах.
Минимальный рабочий пример:
#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
M = 30 #Number of particles in 1 direction
Mn = 90 #Number of grid cells in 1 direction
Lx = 10 #grid physical size
u = np.random.random(M*M*M)
v = np.random.random(M*M*M)
w = np.random.random(M*M*M)
#Have purposefully taken smaller cube, to show potential works
planex = M*u
planey = M*v
planez = M*w
#Create a new grid
grid = np.zeros([Mn,Mn,Mn], dtype='cfloat')
#cell center coordinates
x_c = np.floor(planex).astype(int)%Mn
y_c = np.floor(planey).astype(int)%Mn
z_c = np.floor(planez).astype(int)%Mn
#in terms of the average density of the universe, doesnt matter for the
#example
mass = 1.
#Update the grid
grid[z_c,y_c,x_c] += mass
fig = plt.figure()
ax = fig.add_subplot(111)
plt.imshow(grid[:,:,2].real)
plt.show()
#FFT the grid
grid = np.fft.fftn(grid)
#The resolution and the k-space vectors are the parts I am unsure about
resolution = np.pi*2/(M/Mn)
resolution = Lx/Mn
#Define the k-space vectors
k_x = np.fft.fftfreq(Mn, resolution)
k_y = np.fft.fftfreq(Mn, resolution)
k_z = np.fft.fftfreq(Mn, resolution)
kz, ky, kx = np.meshgrid(k_z, k_y, k_x)
Omega_0 = 0.27
a = 0.3
#Calculate Greens function
k_squared = np.sin(kz/2)**2 + np.sin(ky/2)**2 + np.sin(kx/2)**2
Greens = -3*Omega_0/8/a*np.divide(1, k_squared, where=k_squared!=0)
#Multiply the grids in Fourier space
grid = Greens*grid
#IFFT to real space
potentials = np.fft.ifftn(grid)
fig1 = plt.figure()
ax1 = fig1.add_subplot(111)
plt.imshow(potentials[:,:,0].real)
plt.show()
Большое значение дляразрешение делает скорости взрывными, малое значение и очень малые скорости.Так что делает правильное разрешение?
Я впервые спрашиваю о переполнении стека, пожалуйста, дайте мне знать, если я делаю что-то не так.
Бест, р.