Я пишу псевдоспектральный код CFD на Фортране, который, по сути, шагает по времени в уравнениях Навье-Стокса в плоском слое.В моем случае это действительно трехмерный код, но проблему можно хорошо понять в 2d, поэтому я буду придерживаться этого случая.Геометрически мой слой 2d-плоскости ограничен y=0 и y=1 и является периодическим вдоль другого направления x.Не вдаваясь слишком много в сорняки, эффективная дискретизация заключается в разложении полей (например, скорости) на полиномах Чебышева вдоль y и мод Фурье вдоль x.Чебышевские многочлены по существу являются замаскированными косинусами на искаженной сетке.Уравнения Навье-Стокса имеют простую форму в спектральном пространстве, за исключением нелинейного члена.Таким образом, большинство вычислений выполняется в спектральном пространстве со случайными отклонениями в физическом пространстве для вычисления нелинейного члена: необходимо преобразовать двумерный массив комплексных коэффициентов Чебышева-Фурье в соответствующее двумерное поле (то есть массив вещественных значенийна сетке).Без других ограничений это преобразование относительно легко реализовать.Например, начиная с комплексных спектральных коэффициентов - назовем их c_in(NX, NY/2+1) - можно взять комплексное в реальное, дискретное преобразование Фурье вдоль x для каждого значения y, чтобы получить двумерный массив действительные чебышевские коэффициенты.Затем можно выполнить дискретное косинусное преобразование (FFTW_REDFT в FFTW) вдоль y для всех x и вуаля , наконец, получим реальное поле, r_out(NX,NY).
Корень всех проблем в том, что по некоторым причинам мне нужно , чтобы сначала вычислить DCT.Это проблема, потому что косинусные преобразования реализуются только для реальных данных в FFTW.Различные ограничения приводят к тому, что я не хочу разбивать свой сложный массив спектральных коэффициентов на два реальных массива для действительной и мнимой частей.Учитывая эти ограничения на структуру данных, возникает вопрос: как эффективно заставить FFTW вычислять несколько DCT по первому индексу массива комплексных чисел.
Пока мое решение состоит в использовании расширенного интерфейса plan_many_r2rчтобы определить преобразование, которое перепрыгнет через мнимые значения: я установлю idist в 2. В результате, если я использую этот план с указателем ptr2real_in, связанным с действительной частью c_in(1,1), косинусное преобразованиевсе реальные части вычисляются.Затем я повторяю выполнение плана с указателем ptr2imag_in, связанным с мнимой частью c_in(1,1).После этого вычисление комплексного и реального ДПФ во втором измерении становится простым.
Таким образом, суть этого подхода заключается в определении ptr2imag_in, который на самом деле является адресом памяти c_in(1,1), смещенным на размер впамять о C_double.Ниже приведен минимальный пример, который работает, но выглядит неуклюже для меня.В частности, я определяю указатель на мнимую часть сложного массива следующим образом
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 байтов.Как я мог это сделать?Сбой следующего кода:
cptr = C_loc(c_in(1,1))
cptr = cptr + 8
Call C_F_Pointer(cptr, ptr2imag_in, [2*NX, (NY/2+1)])
Ниже приведен полный минимальный пример взятия DCT, сопровождаемого комплексным и реальным ДПФ:
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