- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
我正在学习如何使用 Fortran 处理 FFTW 包。为了生成一个易于验证的示例,我计算了一个二维平面的功率谱,我用两个不同的叠加波填充它。这样,我就可以确切地知道功率谱中的峰值在哪里。
由于 FFTW 文档针对 C 的内容更为广泛,我首先在 C 中实现了该算法,这给了我相当满意的结果: “lambda1”和“lambda2”是已知的预期波长。蓝线是得到的功率谱。
我不知道去哪里寻找可能的错误。代码执行顺利。有人可以帮忙吗?
这是 C 代码,编译时使用:gcc stackexchange.c -o a.out -I/home/myname/.local/include -lfftw3 -lm -g -Wall
(gcc 5.4 .0)
#include <fftw3.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
int main(void)
{
// parameters
int Nx = 200;
int Ny = 100;
int nsamples = 200;
float pi = 3.1415926;
float physical_length_x = 20;
float physical_length_y = 10;
float lambda1 = 0.5;
float lambda2 = 0.7;
float dx = physical_length_x/Nx;
float dy = physical_length_y/Ny;
float dkx = 2*pi/physical_length_x;
float dky = 2*pi/physical_length_y;
// counters/iterators
int ind, i, j, ix, iy, ik;
// power spectra stuff
float *Pk, *distances_k, d, kmax;
float *Pk_field;
// FFTW vars and arrays
fftw_complex *in, *out;
fftw_plan my_plan;
// allocate arrays for input/output
in = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * Nx * Ny);
out = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * Nx * Ny);
// Create Plan
int n[2]; n[0] = Nx; n[1] = Ny;
my_plan = fftw_plan_dft(2, n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
// fill up array with a wave
for (i=0; i<Nx; i++){
for (j=0; j<Ny; j++){
ind = i*Ny + j;
in[ind][0] = cos(2.0*pi/lambda1*i*dx) + sin(2.0*pi/lambda2*j*dy);
}
}
// execute fft
fftw_execute(my_plan); /* repeat as needed */
// Calculate power spectrum: P(kx, ky) = |F[delta(x,y)]|^2
Pk_field = malloc(sizeof(float)*Nx*Ny);
for (i=0; i<Nx; i++){
for (j=0; j<Ny; j++){
ind = i*Ny+j;
Pk_field[ind] = out[ind][0]*out[ind][0] + out[ind][1]*out[ind][1];
}
}
Pk = calloc(nsamples, sizeof(float));
distances_k = malloc(nsamples*sizeof(float));
kmax = sqrt(pow((Nx/2+1)*dkx, 2) + pow((Ny/2+1)*dky, 2));
for(i=0; i<nsamples; i++){
distances_k[i]= 1.0001*i/nsamples*kmax; // add a little more to make sure kmax will fit
}
// histogrammize P(|k|)
for (i=0; i<Nx; i++){
if (i<Nx/2+1)
{ ix = i; }
else
{ ix = -Nx+i; }
for (j=0; j<Ny; j++){
if (j<Ny/2+1)
{ iy = j; }
else
{ iy = -Ny+j; }
ind = i*Ny + j;
d = sqrt(pow(ix*dkx,2)+pow(iy*dky,2));
for(ik=0; ik<nsamples; ik++){
if (d<=distances_k[ik] || ik==nsamples){
break;
}
}
Pk[ik] += Pk_field[ind];
}
}
//-----------------------------------
// write arrays to file.
// can plot them with plot_fftw.py
//-----------------------------------
FILE *filep;
filep = fopen("./fftw_output_2d_complex.txt", "w");
for (i=0; i<nsamples; i++){
fprintf(filep, "%f\t%f\n", distances_k[i], Pk[i]);
}
fclose(filep);
printf("Finished! Written results to ./fftw_output_2d_complex.txt\n");
//----------------------
// deallocate arrays
//----------------------
fftw_destroy_plan(my_plan);
fftw_free(in); fftw_free(out);
return 0;
}
这是 Fortran 代码,编译时使用:gfortran stackexchange.f03 -o a.out -L/home/my_name/.local/lib/-I/home/my_name/.local/include/-lfftw3 -g -Wall
program use_fftw
use,intrinsic :: iso_c_binding
implicit none
include 'fftw3.f03'
integer, parameter :: dp=kind(1.0d0)
integer, parameter :: Nx = 200
integer, parameter :: Ny = 100
integer, parameter :: nsamples = 200
real(dp), parameter :: pi = 3.1415926d0
real(dp), parameter :: physical_length_x = 20.d0
real(dp), parameter :: physical_length_y = 10.d0
real(dp), parameter :: lambda1 = 0.5d0
real(dp), parameter :: lambda2 = 0.7d0
real(dp), parameter :: dx = physical_length_x/real(Nx,dp)
real(dp), parameter :: dy = physical_length_y/real(Ny,dp)
real(dp), parameter :: dkx = 2.d0 * pi / physical_length_x
real(dp), parameter :: dky = 2.d0 * pi / physical_length_y
integer :: i, j, ix, iy, ik
real(dp):: kmax, d
complex(dp), allocatable, dimension(:,:) :: arr_in, arr_out, Pk_field
real(dp), allocatable, dimension(:) :: Pk, distances_k
integer*8 :: my_plan
integer :: n(2) = (/Nx, Ny/)
allocate(arr_in( 1:Nx,1:Ny))
allocate(arr_out(1:Nx,1:Ny))
! call dfftw_plan_dft_2d(my_plan, Nx, Ny, arr_in, arr_out, FFTW_FORWARD, FFTW_ESTIMATE) ! doesn't help either
call dfftw_plan_dft(my_plan, 2, n, arr_in, arr_out, FFTW_FORWARD, FFTW_ESTIMATE)
! fill up wave
do i = 1, Nx
do j = 1, Ny
arr_in(i,j) =cmplx(cos(2.0*pi/lambda1*i*dx)+sin(2.0*pi/lambda2*j*dy) , 0.d0, kind=dp)
enddo
enddo
! execute fft
call dfftw_execute_dft(my_plan, arr_in, arr_out)
allocate(Pk_field(1:Nx, 1:Ny))
allocate(distances_k(1:nsamples))
allocate(Pk(1:nsamples))
Pk=0
distances_k=0
! Get bin distances
kmax = sqrt(((Nx/2+1)*dkx)**2+((Ny/2+1)*dky)**2)*1.001
do i = 1, nsamples
distances_k(i) = i*kmax/nsamples
enddo
! Compute P(k) field, distances field
do i = 1, Nx
do j = 1, Ny
Pk_field(i,j) = arr_out(i,j)*conjg(arr_out(i,j))
if (i<=Nx/2+1) then
ix = i
else
ix = -Nx+i
endif
if (j<=Ny/2+1) then
iy = j
else
iy = -Ny+j
endif
d = sqrt((ix*dkx)**2+(iy*dky)**2)
do ik=1, nsamples
if (d<=distances_k(ik) .or. ik==nsamples) exit
enddo
Pk(ik) = Pk(ik)+real(Pk_field(i,j))
enddo
enddo
! write file
open(unit=666,file='./fftw_output_2d_complex.txt', form='formatted')
do i = 1, nsamples
write(666, '(2E14.5,x)') distances_k(i), Pk(i)
enddo
close(666)
deallocate(arr_in, arr_out, Pk, Pk_field, distances_k)
call dfftw_destroy_plan(my_plan)
write(*,*) "Finished! Written results to fftw_output_2d_complex.txt"
end program use_fftw
为了您的方便,这是我用来绘制结果的简短 python 脚本:
#!/usr/bin/python3
#====================================
# Plots the results of the FFTW
# example programs.
#====================================
import numpy as np
import matplotlib.pyplot as plt
from sys import argv
from time import sleep
errormessage="""
I require an argument: Which output file to plot.
Usage: ./plot_fftw.py <case>
options for case:
1 fftw_1d_complex.txt
2 fftw_2d_complex.txt
3 fftw_1d_real.txt
4 fftw_2d_real.txt
Please select a case: """
#----------------------
# Hardcoded stuff
#----------------------
file_dict={}
file_dict['1'] = ('fftw_output_1d_complex.txt', '1d complex fftw')
file_dict['2'] = ('fftw_output_2d_complex.txt', '2d complex fftw')
file_dict['3'] = ('fftw_output_1d_real.txt', '1d real fftw')
file_dict['4'] = ('fftw_output_2d_real.txt', '2d real fftw')
lambda1=0.5
lambda2=0.7
#------------------------
# Get case from cmdline
#------------------------
case = ''
def enforce_integer():
global case
while True:
case = input(errormessage)
try:
int(case)
break
except ValueError:
print("\n\n!!! Error: Case must be an integer !!!\n\n")
sleep(2)
if len(argv) != 2:
enforce_integer()
else:
try:
int(argv[1])
case = argv[1]
except ValueError:
enforce_integer()
filename,title=file_dict[case]
#-------------------------------
# Read and plot data
#-------------------------------
k, Pk = np.loadtxt(filename, dtype=float, unpack=True)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title("Power spectrum for "+title)
ax.set_xlabel("k")
ax.set_ylabel("P(k)")
# ax.plot(k, Pk, label='power spectrum')
ax.semilogx(k[k>0], Pk[k>0], label='power spectrum') # ignore negative k
ax.plot([2*np.pi/lambda1]*2, [Pk.min()-1, Pk.max()+1], label='expected lambda1')
ax.plot([2*np.pi/lambda2]*2, [Pk.min()-1, Pk.max()+1], label='expected lambda2')
ax.legend()
plt.show()
最佳答案
在直方图的计算中,循环的第一个索引对应于平均值 i。 e. “零频率”。
在C中,当i=0,j=0,ix=0,iy=0,d=0。那是正确的。对于fortran中数组的同一项,索引i=1,则ix=1,d不再为null。
那是零频率,但由于 C 和 Fortran 之间的索引差异,整个直方图被轻微污染。特别是,修改可能以不同方式影响正频率和负频率,从而触发图中观察到的双峰。
您能否尝试将 Fortran 中 ix 和 iy 的计算修改为:
if (i-1<Nx/2+1) then
ix = i-1
else
ix = -Nx+i-1
endif
if (j-1<Ny/2+1) then
iy = j-1
else
iy = -Ny+j-1
endif
关于c - fortran 2d-FFTW 与 C 2d-FFTW 结果不一致,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/50726456/
是的,我知道..,这不是想象的...这是一个真正的 Fortran 问题。 以前的版本是指 Fortran 2003、95、90,甚至 77。 我所说的“向后兼容”是指可以轻松运行为 2008 年以前
我有一个程序,它的变量中有一个值。一旦确定了该值,我想调用另一个程序并使用该变量的值来确定在新程序中的位置。有人知道该怎么做吗? 最佳答案 如果您有 Fortran 2008 编译器,您将拥有标准子例
namelist 是一种有用的 fortran 结构,可以从文件中快速初始化变量。 namelist 有一个名称并包含一组具有已知类型的变量。这使得它类似于 type 结构。 通常情况下,给程序或子例
我正在遍历索引,我正在检查我是否不在第一个循环交互和另一个条件中。如果第一个条件是 .False.,我不想评估第二个条件。 do i = 1, n if ( i /= 1 .and. var(
Fortran 2003 具有用于数组连接的方括号语法,Intel fortran 编译器也支持它。我在这里为矩阵连接写了一个简单的代码: program matrix implicit none r
我正在尝试通过重载类型名称来制作自定义数据类型构造函数。但是,在进行调用时,将调用默认构造函数。我不明白我做错了什么。 这是有问题的代码片段。 module test type, pu
我的最终目标是在 Fortran 中有一个通用的映射函数,即一个接受任意类型 A 的数组和一个 A->B 类型的函数的函数,将此函数应用于给定数组的所有元素并返回一个B 类型的数组。我无法用数组实现它
我正在学习 Fortran,在使用格式编写时发现了一些奇怪的东西(我使用的是 Fortran onlinegdb) Program Hello real, dimension(3,2):: array
Fortran 中的INTERFACE 语句是否使其成为正式实现multiple dispatch 的编程语言? ? (我问是因为所链接的维基百科文章在其看似全面的支持相关范式的示例编程语言列表中并未
我可以使用 Fortran 95 编译器编译 Fortran 90 文件吗? Fortran 95 似乎有很多,但 Fortran 90 没有。 最佳答案 这个可以: NAGWare f95 Comp
嗨,我在 Fortran 中对二维离散化问题强加边界条件时遇到了麻烦。我的离散化网格是一个二维正方形,在 x,y 方向上从 -L 到 L。 我想强加这样的边界条件, 在 x=L 的边界线上,指定了函数
Fortran 是否有与 C assert 等效的标准函数/关键字? ? 我找不到 assert我在Fortran2003标准中提到过。我发现了一些如何使用预处理器的方法,但是在这个 answer建议
我有一系列的作业,使用“;”将它们分配给同一个ike。分开statemnts,但我收到此错误: 1.0;磅(1,9) 1个 错误:(1)处无法分类的陈述 在文件LJ.F90:223中 如果每个语句都在
我正在使用 gfortran -std=f2008。我有一个函数,它返回一个包含可分配数组的派生类型。该函数在返回之前调用allocate()。似乎在分配数组的函数返回之后,数组会自动释放一段时间,并
我制作了这个小型测试程序来“证明”在编译之前(或者如果你让它们可分配),你不能在不指定它们的大小的情况下使用向量。我的观点失败了。我期待本地向量“num”会失败。程序在执行程序之前无法知道它的大小。大
出于优化原因,Fortran 强制子例程或函数的虚拟参数不是别名,即它们不指向相同的内存位置。 我想知道相同的约束是否适用于函数的返回值。 换句话说,对于给定的 myfunc 函数: function
我已经在Fortran 90中编写了一个相当大的程序。它已经运行了一段时间了,但是今天我尝试将其提高一个档次并增加问题的大小(这是研究非标准的有限元求解器,如果那样的话)。可以帮助任何人...)现在,
在 C 和 C++ 中,有许多操作会导致未定义的行为,即允许编译器做任何它想做的事情的情况。 Examples包括在释放变量后使用它,释放变量两次和取消引用空指针。 Fortran 是否也有未定义的行
通常我使用fortran进行数值分析,然后使用matlab、R和python进行后期和前期工作。 我发现 matlab、R 和 python 在终端中提供了命令提示符,以便您可以运行脚本以及从命令行立
在 Fortran 中将变量设置为 +Infinity 的最安全方法是什么?目前我正在使用: program test implicit none print *,infinity() con
我是一名优秀的程序员,十分优秀!