Решатель уравнения Пуассона в трехмерном синусоидальном преобразовании типа r2r с использованием FFTW в Фортране

Question

Я написал следующий код на Фортране, чтобы решить уравнение Пуассона, используя преобразование типа FFTW типа r2r (от реального к реальному). В этом коде сначала я сделал нормальный FFTW, т.е. тип математической функции r2c (от реального к сложному) 27*sin(3x)*sin(3y)*sin(3z), а затем разделил его на 27 (3 * 3 + 3 * 3 + 3 * 3) для вычисления второй производной от ввода функция. Трехмерный график амплитуды входной функции показывает ее правильную амплитуду. амплитуда вдоль оси z относительно x-, y-координаты .
Затем обратное БПФ типа c2r (от сложного к реальному) восстанавливает входную функцию, но амплитуда теперь уменьшается до 1, как показано на трехмерном графике 2 . Это говорит о том, что мой код для решения уравнений Пуассона работает нормально в случае нормального FFTW.

Program Derivative

! To run this code: gcc dst_3d.c -c -std=c99 && gfortran derivative.f95 dst_3d.o -lfftw3 && ./a.out

implicit none

include "fftw3.f"

integer ( kind = 4 ), parameter :: Nx = 64
integer ( kind = 4 ), parameter :: Ny = 64
integer ( kind = 4 ), parameter :: Nz = 64

real ( kind = 8 ), parameter :: pi=3.14159265358979323846d0

integer ( kind = 4 ) i,j,k
real ( kind = 8 ) Lx,Ly,Lz,dx,dy,dz,kx,ky,kz
real ( kind = 8 ) x(Nx),y(Ny),z(Nz)

real ( kind = 8 ) in_dst(Nx,Ny,Nz),out_dst(Nx,Ny,Nz) ! DST
real ( kind = 8 ) in_k_dst(Nx,Ny,Nz),out_k_dst(Nx,Ny,Nz) ! DST

real ( kind = 8 ) in_dft(Nx,Ny,Nz),out_dft(Nx,Ny,Nz) ! DFT
complex ( kind = 8 ) in_k_dft(Nx/2+1,Ny,Nz),out_k_dft(Nx/2+1,Ny,Nz) ! DFT
integer ( kind = 8 ) plan_forward,plan_backward ! DFT

! System Size.
Lx = 2.0d0*pi; Ly = 2.0d0*pi; Lz = 2.0d0*pi 

! Grid Resolution.
dx = Lx/dfloat(Nx); dy = Ly/dfloat(Ny); dz = Lz/dfloat(Nz)

! =================================== INPUT ===========================================

! Initial Profile Details.
do i = 1, Nx
  x(i) = dfloat(i-1)*dx
  do j = 1, Ny
    y(j) = dfloat(j-1)*dy
    do k = 1, Nz
      z(k) = dfloat(k-1)*dz
      in_dst(i,j,k) = 27.0d0*sin(3.0d0*x(i))*sin(3.0d0*y(j))*sin(3.0d0*z(k))
      in_dft(i,j,k) = in_dst(i,j,k)
      write(10,*) x(i), y(j), z(k), in_dft(i,j,k)
    enddo
  enddo
enddo

! =================================== 3D DFT ===========================================

  call dfftw_plan_dft_r2c_3d_ (plan_forward, Nx, Ny, Nz, in_dft, in_k_dft, FFTW_ESTIMATE)
  call dfftw_execute_ (plan_forward)
  call dfftw_destroy_plan_ (plan_forward)

do i = 1, Nx/2+1
  do j = 1, Ny/2
    do k = 1, Nz/2
      kx = 2.0d0*pi*dfloat(i-1)/Lx
      ky = 2.0d0*pi*dfloat(j-1)/Ly
      kz = 2.0d0*pi*dfloat(k-1)/Lz
      out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo 
    do k = Nz/2+1, Nz
      kx = 2.0d0*pi*dfloat(i-1)/Lx
      ky = 2.0d0*pi*dfloat(j-1)/Ly
      kz = 2.0d0*pi*dfloat((k-1)-Nz)/Lz
      out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
  enddo
  do j = Ny/2+1, Ny
    do k = 1, Nz/2
      kx = 2.0d0*pi*dfloat(i-1)/Lx
      ky = 2.0d0*pi*dfloat((j-1)-Ny)/Ly
      kz = 2.0d0*pi*dfloat(k-1)/Lz
      out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
    do k = Nz/2+1, Nz
      kx = 2.0d0*pi*dfloat(i-1)/Lx
      ky = 2.0d0*pi*dfloat((j-1)-Ny)/Ly
      kz = 2.0d0*pi*dfloat((k-1)-Nz)/Lz
      out_k_dft(i,j,k) = in_k_dft(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
  enddo   
enddo

out_k_dft(1,1,1) = in_k_dft(1,1,1)

  call dfftw_plan_dft_c2r_3d_ (plan_backward, Nx, Ny, Nz, out_k_dft, out_dft, FFTW_ESTIMATE)
  call dfftw_execute_ (plan_backward)
  call dfftw_destroy_plan_ (plan_backward)

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    out_dft(i,j,k) = out_dft(i,j,k)/dfloat(Nx*Ny*Nz)
    write(20,*) x(i), y(j), z(k), out_dft(i,j,k)    
    enddo
  enddo   
enddo

! =================================== 3D DST ===========================================

  call Forward_FFT(Nx, Ny, Nz, in_dst, in_k_dst)

do k = 1, Nz
  do j = 1, Ny/2
    do i = 1, Nx/2
      kz = 2.0d0*pi*dfloat((k-1))/Lz
      ky = 2.0d0*pi*dfloat((j-1))/Ly
      kx = 2.0d0*pi*dfloat((i-1))/Lx
      out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo 
    do i = Nx/2+1, Nx
      kz = 2.0d0*pi*dfloat((k-1))/Lz
      ky = 2.0d0*pi*dfloat((j-1))/Ly
      kx = 2.0d0*pi*dfloat(Nx-(i-1))/Lx
      out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
  enddo  
  do j = Ny/2+1, Ny
    do i = 1, Nx/2
      kz = 2.0d0*pi*dfloat((k-1))/Lz
      ky = 2.0d0*pi*dfloat(Ny-(j-1))/Ly
      kx = 2.0d0*pi*dfloat((i-1))/Lx
      out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
    do i = Nx/2+1, Nx
      kz = 2.0d0*pi*dfloat((k-1))/Lz
      ky = 2.0d0*pi*dfloat(Ny-(j-1))/Ly
      kx = 2.0d0*pi*dfloat(Nx-(i-1))/Lx
      out_k_dst(i,j,k) = in_k_dst(i,j,k)/(kx*kx+ky*ky+kz*kz)
    enddo
  enddo   
enddo

out_k_dst(1,1,1) = in_k_dst(1,1,1)

  call Backward_FFT(Nx, Ny, Nz, out_k_dst, out_dst)  

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    out_dst(i,j,k) = out_dst(i,j,k)/dfloat(8*Nx*Ny*Nz)
    write(30,*) x(i), y(j), z(k), out_dst(i,j,k)
    enddo
  enddo   
enddo

end program Derivative

! ================================== FFTW SUBROUTINES 
====================================================

! ================================================================= !
! Wrapper Subroutine to call forward fftw c functions for 3d arrays !
! ================================================================= !

subroutine Forward_FFT(Nx, Ny, Nz, in, out)
implicit none
integer ( kind = 4 ) Nx,Ny,Nz,i,j,k
real ( kind = 8 ) in(Nx, Ny, Nz),out(Nx, Ny, Nz),dum(Nx*Ny*Nz)

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    dum(((i-1)*Ny+(j-1))*Nz+k) = in(i,j,k)
    enddo
  enddo
enddo

  call dst3f(Nx, Ny, Nz, dum)

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    out(i,j,k) = dum(((i-1)*Ny+(j-1))*Nz+k)
    enddo
  enddo
enddo

end subroutine

! ================================================================== !
! Wrapper Subroutine to call backward fftw c functions for 3d arrays !
! ================================================================== !

subroutine Backward_FFT(Nx, Ny, Nz, in, out)
implicit none
integer ( kind = 4 ) Nx,Ny,Nz,i,j,k
real ( kind = 8 ) in(Nx, Ny, Nz),out(Nx, Ny, Nz),dum(Nx*Ny*Nz)

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    dum(((i-1)*Ny+(j-1))*Nz+k) = in(i,j,k)
    enddo
  enddo
enddo

  call dst3b(Nx, Ny, Nz, dum)

do i = 1, Nx
  do j = 1, Ny
    do k = 1, Nz
    out(i,j,k) = dum(((i-1)*Ny+(j-1))*Nz+k)
    enddo
  enddo
enddo

end subroutine
! ==================================================================

Этот код использует C-оболочку ниже для вычисления прямого синусоидального преобразования FFTW и обратного синусоидального преобразования FFTW,

#include <fftw3.h>

int dst3f_(int *n0, int *n1, int *n2, double *in3cs)
{
    double *out3cs;
    out3cs = (double*) fftw_malloc(sizeof(double) * (*n0) * (*n1) * (*n2));
    fftw_plan p3cs;
    p3cs = fftw_plan_r2r_3d(*n0, *n1, *n2, in3cs, out3cs, FFTW_RODFT10, FFTW_RODFT10, FFTW_RODFT10, FFTW_ESTIMATE);

    fftw_execute(p3cs);
    fftw_destroy_plan(p3cs);

    for(int i0=0;i0<*n0;i0++) {
        for(int i1=0;i1<*n1;i1++) {
            for(int i2=0;i2<*n2;i2++) {
                in3cs[(i0*(*n1)+i1)*(*n2)+i2] = out3cs[(i0*(*n1)+i1)*(*n2)+i2];
            }
        }
    }

    return 0;
}

int dst3b_(int *n0, int *n1, int *n2, double *in3cs)
{
    double *out3cs;
    out3cs = (double*) fftw_malloc(sizeof(double) * (*n0) * (*n1) * (*n2));
    fftw_plan p3cs;
    p3cs = fftw_plan_r2r_3d(*n0, *n1, *n2, in3cs, out3cs, FFTW_RODFT01, FFTW_RODFT01, FFTW_RODFT01, FFTW_ESTIMATE);

    fftw_execute(p3cs);
    fftw_destroy_plan(p3cs);

    for(int i0=0;i0<*n0;i0++) {
        for(int i1=0;i1<*n1;i1++) {
            for(int i2=0;i2<*n2;i2++) {
                in3cs[(i0*(*n1)+i1)*(*n2)+i2] = out3cs[(i0*(*n1)+i1)*(*n2)+i2];
            }
        }
    }

    return 0;
}

Затем я попытался решить то же самое уравнение Пуассона, которое теперь использует синусоидальное преобразование FFTW, то есть тип r2r (от реального к реальному). Когда я делаю трехмерный график вывода, как показано здесь 3 , амплитуда теперь уменьшается до значения меньше 1. Я не смог выяснить, где моя ошибка в коде, из-за которой эта амплитуда равна уменьшается в случае синусоидального преобразования.

Answer 1

Использование реальных преобразований для решения уравнения Пуассона очень привлекательно, поскольку позволяет использовать различные граничные условия. Тем не менее, точки чувствительности не совсем соответствуют тем, которые рассматриваются для периодических граничных условий.

Для периодических граничных условий чувствительные точки находятся на регулярной сетке, скажем, 0,1 .., n-1. Если размер элементарной ячейки равен Lx, то расстояние между точками равно Lx / n. Конец истории.

Теперь давайте рассмотрим граничные условия, такие, что DST III должен использоваться для прямого преобразования, то есть флаг RODFT01. Как показано там , в документации FFTW сообщается, что:

FFTW_RODFT01 (DST-III): нечетное около j = -1 и четное около j = n-1.

Чувствительные точки все еще 0,1 .., n-1. Если длина Lx, расстояние по-прежнему равно dx = Lx / (n). Но входные и базовые функции DST III нечетны при j = -1 и даже при j = n-1. Это объясняет несоответствие величины :

in_dst(i,j,k) = 27.0d0*sin(3.0d0*x(i))*sin(3.0d0*y(j))*sin(3.0d0*z(k))

Действительно, этот входной сигнал не является нечетным для i = -1 и даже для i = n-1. Это странно вокруг i = 0 и i = n. Как следствие, могут быть полезны следующие действия:

• Убедиться, что рассматриваемое прямое преобразование является правильным. Для нечетно-нечетных граничных условий это может быть FFTW_RODFT00 (DST-I): odd around j=-1 and odd around j=n..
• Оцените ввод в правильных точках измерения. Для DST-I, используя n точек, dx = Lx / (n + 1) и x_i = dx * (i + 1) для i = 0,1, .., n-1, если функция нечетна в 0 и Lx .
• Вычислить частоты в преобразованном пространстве. Способ обратного преобразования пишет помогает найти эти частоты. Поскольку инверсией DST I является DST I, первая частота (k = 0) соответствует периоду 2Lx; вторая (k = 1) соответствует Lx, k-я частота соответствует периоду 2Lx / (k + 1).

Наконец, может потребоваться очистка выходных данных, поскольку преобразования FFTW не нормированы. Для DST-I, FFTW_RODFT00, это N=2(n+1). Совет Владимир Ф, несомненно, хороший. Действительно, хотя тестирование одной частоты идеально подходит для понимания и реализации алгоритма, финальный тест должен охватывать несколько частот, чтобы обеспечить надежность программы.