计算科学 - 使用傅立叶变换求解泊松方程时如何避免被零除？ - 吾爱随笔录

我想尝试使用傅里叶变换来实现以下文章中的部分方法。

http://www.shodor.org/media/content/jocse/student_submissions/nocito2010/nocito2010_pdf

现在我正在使用 FFT 来近似每个点的梯度，并使用它来计算每个物体上的力。我目前得到奇怪的结果。当我只在系统中添加一个巨大的身体时，我得到以下作为我的潜力，这似乎是不正确的。我知道傅里叶变换的性质会导致环绕会扭曲球形，但我仍然不会期待这一点。

这是我的代码。我遇到的一个问题是，在论文中你应该除以|k|^2。然而，在某些时候，有一个零会导致问题。这就是为什么我将它限制在 0.1。不确定是否有更好的方法来处理它。

    #include <iostream>
#include <arrayfire.h>
#include <memory>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <cstring>

const int N = 2;
const int GRID_DIM = 1024;
const double CELL_SIZE = 100.0;
const double G = 1.0;
const double DT = 0.001;

struct Vec2
{
    float x, y;
    Vec2(float x = 0, float y = 0) :
        x(x), y(y) { }
};

int wrapIndex(int i, int i_max) 
{
    return ((i % i_max) + i_max) % i_max;
}


af::array makeGrid(const Vec2* p, const float* m)
{
    float* grid_h = new float[GRID_DIM * GRID_DIM];
    memset(grid_h, 0, sizeof(float) * GRID_DIM * GRID_DIM);

    for(int i = 0; i < N; i++)
    {
        int x = wrapIndex(roundf(p[i].x / CELL_SIZE), GRID_DIM);
        int y = wrapIndex(roundf(p[i].y / CELL_SIZE), GRID_DIM);
        grid_h[y * GRID_DIM + x] += m[i];
    }

    af::array grid(GRID_DIM, GRID_DIM, grid_h);
    delete[] grid_h;
    return grid;
}


int main()
{
    af::setBackend(AF_BACKEND_CUDA);
    af::info();
    Vec2* p = new Vec2[N];
    Vec2* v = new Vec2[N];
    float* m = new float[N];

    std::fill(v, v + N, 0.0f);
    std::fill(m, m + N, 1.0f);
//  for(int i = 0; i < N; i++)
//      p[i] = Vec2(rand() % GRID_DIM / 5 + GRID_DIM / 2, rand() % GRID_DIM / 5 + GRID_DIM / 2);

    //add one massive object at the center
    m[0] = 1000000.0;
    p[0] = Vec2(GRID_DIM * CELL_SIZE / 2, GRID_DIM * CELL_SIZE / 2);

    m[1] = 0.0;
    p[1] = Vec2(GRID_DIM * CELL_SIZE / 4, GRID_DIM * CELL_SIZE / 4);

    af::array rng = af::range(GRID_DIM);
    af::array k2 = 2.0 * af::Pi * af::select(rng > GRID_DIM / 2, (rng - GRID_DIM) , rng) / (GRID_DIM * CELL_SIZE);
    //af::array k2 = 2.0 * af::Pi * rng / (GRID_DIM * CELL_SIZE);
    k2(0) = 1.0;
    k2 = af::tile((k2 * k2).T(), GRID_DIM) + af::tile(k2 * k2, 1, GRID_DIM);

    af::Window wnd(GRID_DIM, GRID_DIM);

    //load the grid with bodies
    af::array grid = makeGrid(p, m);
    while(!wnd.close())
    {
        af::timer t = af::timer::start();

        //fft the density grid
        grid = af::fft2(grid);

        //apply the formula from the paper
        grid *= -G / (af::Pi * k2);
        //invert the fft to get the potentials
        grid = af::ifft2(grid);
        grid = af::abs(grid);

        //compute the gradients in the potential
        af::array dx, dy;
        af::grad(dx, dy, grid);

        float* dx_h = dx.host<float>();
        float* dy_h = dy.host<float>();

        //use the gradients to time step each of the bodies
        #pragma omp parallel for simd
        for(int i = 0; i < N; i++)
        {
            int xi = wrapIndex(roundf(p[i].x / CELL_SIZE), GRID_DIM);
            int yi = wrapIndex(roundf(p[i].y / CELL_SIZE), GRID_DIM);

            v[i].x += dx_h[yi * GRID_DIM + xi] * DT / m[i];
            v[i].y += dy_h[yi * GRID_DIM + xi] * DT / m[i];
            p[i].x += v[i].x * DT;
            p[i].y += v[i].y * DT;
        }
        af::freeHost(dx_h);
        af::freeHost(dy_h);

        //display the potential

        wnd.image(af::abs(grid) / af::max<float>(af::abs(grid)));
        //remake the new grid
        grid = makeGrid(p, m);

        std::cout << "dt = " << af::timer::stop(t) << std::endl;

    }

    delete[] p;
    delete[] v;
    delete[] m;
    return 0;
}

这是另一张不太饱和的单点潜力照片。我不知道为什么它不是近似圆形/径向对称的。