在Fortran中使用FFTW进行复杂数组的DCT:如何指向虚部数组? [英] DCT of complex arrays with FFTW in Fortran: How to point to the imaginary part array?

问题描述

我正在用Fortran编写伪光谱CFD代码，它本质上是平面层中Navier-Stokes方程的时间步长.在我的情况下，这确实是3d代码，但是在2d中可以很好地理解该问题，因此我会坚持这种情况.在几何上，我的2d平面层由y=0和y=1界定，并且沿着另一个方向x是周期性的.无需过多研究杂草，有效的离散化是分解沿y的切比雪夫多项式上的场(例如速度)和沿x的傅立叶模式. Chebyshev多项式本质上是伪装在扭曲网格上的余弦. Navier-Stokes方程在光谱空间中具有简单形式，非线性项除外.因此，大多数计算是在频谱空间中进行的，偶尔会游览物理空间以计算非线性项:需要将复杂的Chebyshev-Fourier系数的2d数组转换为相应的2d字段(即，实数值的数组)在网格上).没有其他限制，此转换相对容易实现.例如，从复数频谱系数开始-我们称它们为c_in(NX, NY/2+1)-对于y的每个值，可以沿x进行复数到实数，离散傅立叶变换，以获得的二维数组真正的切比雪夫系数.然后，可以对所有x和 voila 沿y执行离散余弦变换(FFTW中的FFTW_REDFT)，最后得到真实场r_out(NX,NY).

I am writing a pseudo-spectral CFD code in Fortran, which is essentially a time-stepper of the Navier-Stokes equations in a plane layer. It is really a 3d code in my case, but the problem can be very well understood in 2d, so I will stick to this case. Geometrically, my 2d plane layer is bounded by y=0 and y=1, and is periodic along the other direction x. Without going too much into the weeds, an efficient discretisation is to decompose fields (e.g., velocity) on Chebyshev polynomials along y and Fourier modes along x. Chebyshev polynomials are essentially cosines in disguise on a distorded grid. The Navier-Stokes equations have a simple form in spectral space, with the exception of the nonlinear term. Therefore, most of the computations are carried on in spectral space, with occasional excursions to physical space to compute the nonlinear term: one needs to transform a 2d array of complex Chebyshev-Fourier coefficients to the corresponding 2d field (i.e., array of real values on a grid). Without other constraint, this transform is relatively easy to implement. For instance, starting from complex spectral coefficients -- let us call them c_in(NX, NY/2+1) -- one may take a complex to real, dicrete Fourier transform along x for each value of y to obtain a 2d array of real Chebyshev coefficients. Then, one may perform a discrete cosine transform (FFTW_REDFT in FFTW) along y for all x and voila, one finally gets the real field, r_out(NX,NY).

所有麻烦的根源在于，出于某些原因，我需要首先计算DCT.这是一个问题，因为余弦变换仅在FFTW中用于实数数据.各种约束导致我不想将频谱系数的复杂阵列分为实部和虚部两个实数阵列.考虑到这些对数据结构的约束，问题是:如何有效地使FFTW沿着复数数组的第一个索引计算多个DCT.

The root of all troubles is that for some reasons, I need to compute the DCT first. This is a problem because cosine transforms are only implemented for real data in FFTW. Various constraints result in me not wanting to split my complex array of spectral coefficients into two real arrays for the real and imaginary parts. Given these constraint on data structure, the question is: how to efficiently get FFTW compute several DCTs along the first index of an array of complex numbers.

到目前为止，我的解决方案包括使用plan_many_r2r高级界面定义一个跳跃虚拟值的转换:我将idist设置为2.结果，如果我将此计划与指针一起使用与c_in(1,1)的实部相关联的ptr2real_in，将计算所有实部的余弦变换.然后，我用与c_in(1,1)的虚构部分相关的指针ptr2imag_in重复执行计划.之后，沿第二维计算复杂到实际的DFT很容易.

So far, my solution consists in using the plan_many_r2r advanced interface to define a transform that leap-frogs over the imaginary values: I set idist to 2. As a result, if I use this plan with a pointer ptr2real_in associated with the real part of c_in(1,1), a cosine transform of all the real parts is computed. Then I repeat the execution of the plan with a pointer ptr2imag_in associated with the imaginary part of c_in(1,1). After that, computing a complex to real DFT along the second dimension is easy.

因此，此方法的关键是定义ptr2imag_in，它实际上是c_in(1,1)的存储器地址，移位了C_double的内存大小.我在下面提供了一个可行的最小示例，但对我来说却显得笨拙.特别是，我以这种方式定义了指向复杂数组虚部的指针

So the crux of this approach is to define ptr2imag_in, which is really the memory address of c_in(1,1) shifted by the size in memory of a C_double. I include a minimal example below which works, but looks clumsy to me. In particular, I define the pointer to the imaginary part of the complex array in this way

cptr = C_loc(c_in(1,1))
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])
cptr = C_loc(ptr2imag_in(2,1))
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])

在我看来，我要做的就是将cptr移8个字节.我该怎么办?以下代码失败:

It seems to me that all I need to do would be to shift cptr by 8 bytes. How could I do that? The following code fails:

cptr = C_loc(c_in(1,1))
cptr = cptr + 8 
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])

---

下面是采用DCT，然后进行复数转换为实数DFT的完整最小示例:

---

The full minimal example for taking a DCT followed by a complex to real DFT is below:

Program monkeying_with_dct 
Use, Intrinsic :: ISO_C_BINDING
Implicit None
include 'fftw3.f03'

Integer, Parameter :: dp = C_Double
Complex(C_double), Pointer :: c_in (:,:)
Complex(C_double), Pointer :: c_out(:,:)
Real(C_Double),    Pointer :: r_out(:,:)
Real(C_Double),    Pointer :: ptr2real_in (:,:)
Real(C_Double),    Pointer :: ptr2real_out(:,:)
Real(C_Double),    Pointer :: ptr2imag_in (:,:)
Real(C_Double),    Pointer :: ptr2imag_out(:,:)
Type(C_ptr) :: cptr
Type(C_ptr) :: plan_IDCT
Type(C_ptr) :: plan_C2R 
Type(C_ptr) :: pdum


Integer, Parameter :: NX = 512
Integer, Parameter :: NY = 1024

print *,'... allocating memory ...'
pdum = fftw_alloc_complex(int((NY/2+1)*NX, C_size_T))
Call C_F_Pointer(pdum, c_in , [NX, NY/2+1])
pdum = fftw_alloc_complex(int((NY/2+1)*NX, C_size_T))
Call C_F_Pointer(pdum, c_out, [NX, NY/2+1])
pdum = fftw_alloc_real(int(NY*NX, C_size_T))
Call C_F_Pointer(pdum, r_out, [NX, NY])

print *,'... initializing data ...'
c_in      = Cmplx(0._dp, 0._dp,  Kind=dp)
c_in(2,3) = Cmplx(1._dp, 0.5_dp, Kind=dp) 

print *, '... defining a pointer to the real part of input complex data ...'                
cptr = C_loc(c_in(1,1))
Call C_F_Pointer(cptr, ptr2real_in, [2*NX, (NY/2+1)])
print *, '... defining a pointer to the imag part of input complex data ...'                
cptr = C_loc(c_in(1,1))
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])
cptr = C_loc(ptr2imag_in(2,1))
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])

print *, '... defining a pointer to the real part of transformed complex data ...'                
cptr = C_loc(c_out(1,1))
Call C_F_Pointer(cptr, ptr2real_out, [2*NX, (NY/2+1)])
print *, '... defining a pointer to the imag part of transformed complex data ...'                
cptr = C_loc(c_out(1,1))
Call C_F_Pointer(cptr, ptr2imag_out, [2*NX, (NY/2+1)])
cptr = C_loc(ptr2imag_out(2,1))
Call C_F_Pointer(cptr, ptr2imag_out, [2*NX, (NY/2+1)])


print*, '... planning IDCT ...'
plan_IDCT = fftw_plan_many_r2r(1, [NX], (NY/2+1), &
                        ptr2real_in,  [2*NX], 2, 2*NX, &
                        ptr2real_out, [2*NX], 2, 2*NX, &
                        [FFTW_REDFT01]  , FFTW_MEASURE) 

print*, '... planning C2R DFT ...' 
plan_C2R   = fftw_plan_many_dft_c2r(1, [NY], NX, & 
                        c_out, [NY/2+1], NX, 1, & 
                        r_out, [NY], NX, 1, & 
                        FFTW_MEASURE)  

print*, '... DCT of the real part ...' 
Call fftw_execute_r2r(plan_IDCT, ptr2real_in, ptr2real_out) 
print*, '... DCT of the imaginary part ...' 
Call fftw_execute_r2r(plan_IDCT, ptr2imag_in, ptr2imag_out) 
print*, '... DFT Complex to real ...' 
Call fftw_execute_dft_c2r(plan_C2R, c_out,r_out) 


End Program monkeying_with_dct 

推荐答案

如果需要的话，可以一直transfer()指向整数的指针，进行移位并transfer()返回.

One can alwas transfer() the pointer to integer, do the shift and transfer() back, if that is what you want.

cptr = transfer( transfer(cptr, 1_c_intptr_t) + c_sizeof(1._c_double) , cptr)

或只可以调用一个小的C函数，以更好的控制方式执行指针算术.但是我不确定这真的是您所需要的.

or one can just call a small C function that does the pointer arithmetic in a better controlled way. But I am not sure it is really what you need.

在Fortran 2008中，人们应该只能使用%im语法，但是对编译器的支持还不够好，gfortran根本不支持它..

In Fortran 2008 one should be able to just use the %im syntax, but the compiler support isn't good yet, gfortran does not support it at all..

这篇关于在Fortran中使用FFTW进行复杂数组的DCT:如何指向虚部数组?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！