Homec4science
Diffusion SP4E_Homework_Ashtari_Sieber (master) Browse

R9482/Homework3/src/test_fft.cc
77031402684emaster

Last Change
Raw File

test_fft.cc
View Options

#include "my_types.hh"
#include "fft.hh"
#include <gtest/gtest.h>
/* ------------------------------------------------------ */
TEST(FFT, transform) {
UInt N = 64;
Matrix<complex> m(N);
Real k = 2 * M_PI / N;
for (auto&& entry : index(m)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
val = cos(k * i);
}
Matrix<complex> res = FFT::transform(m);
for (auto&& entry : index(res)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
if (std::abs(val) > 1e-10)
std::cout << i << "," << j << " = " << val << std::endl;
if (i == 1 && j == 0)
ASSERT_NEAR(std::abs(val), 0.5, 1e-10);
else if (i == N - 1 && j == 0)
ASSERT_NEAR(std::abs(val), 0.5, 1e-10);
else
ASSERT_NEAR(std::abs(val), 0, 1e-10);
}
}
/* ------------------------------------------------------ */
TEST(FFT, inverse_transform) {
// Assuming that the first test will pass and
// FFT works, to test IFFT we transform a
// discrete fucntion to spectral space, then
// inverse transform to physical space again,
// and expect to end up with the same function.
// u(x,y)=sin(6.pi.x/L)+cos(4.pi.y/L)
UInt N = 64; // Grid size in x and y directions
Matrix<complex> physical_original(N); // Discrete function U in physical space
Matrix<complex> physical_inverse(N); // IFFT(FFT(U)) expected to equal U
Matrix<complex> spectral(N); // Discrete function U in spectral space
// Defining U as U(x,y)=sin(6.pi.x/L)+cos(4.pi.y/L)
Real k = 2 * M_PI / N;
for (auto&& entry : index(physical_original)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
val = sin(3.0 * k * j)+cos(2.0 * k * i);
}
// Transforming U to spectral space
spectral = FFT::transform(physical_original);
// Inverse transformin FFT(U) to physical space
physical_inverse = FFT::itransform(spectral);
// Testing equality
for(int i=0; i<N; i++){
for(int j=0; j<N; j++){
ASSERT_NEAR(std::abs(physical_original(i,j)-physical_inverse(i,j)), 0, 1e-11);
}
}
}
/* ------------------------------------------------------ */
TEST(FFT, compute_frequencies) {
// Assuming that the first two tests will pass and FFt and
// IFFT work, we transform an analytical 2D expression to
// spectral space, then calculate its derivative with respect
// to x and y in spectral space, then inverse transform to
// physical space again and expect to end up with derivative
// of the same function, found analytically.
// U(x,y)=sin(6.pi.x/L)+cos(4.pi.y/L)
// dU/dx=(+6.pi/L)cos(6.pi.x/L)
// dU/dy=(-4.pi/L)sin(4.pi.y/L)
// Assume L=4
UInt N = 64; // Grid size in x and y directions
Real L = 4.0; // Length of the domain in both x & y
Real dw = 1.0/L; // Frequency step in both x & y (dw=ds=1/L)
Matrix<complex> U(N); // Discrete function U
Matrix<complex> spectral_U(N); // Fourier coefficients of U
Matrix<complex> spectral_U_x(N); // Fourier coefficients of dU/dx (2.pi.I.w.FFT(U))
Matrix<complex> spectral_U_y(N); // Fourier coefficients of dU/dy (2.pi.I.s.FFT(U))
Matrix<complex> computed_U_x(N); // IFFT of FFT of dU/dx
Matrix<complex> computed_U_y(N); // IFFT of FFT of dU/dy
Matrix<complex> expected_U_x(N); // Analytical dU/dx
Matrix<complex> expected_U_y(N); // Analytical dU/dy
// Defining U as U(x,y)=sin(6.pi.x/L)+cos(4.pi.y/L)
for (auto&& entry : index(U)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
val = sin(6.0*M_PI*j/N)+cos(4.0*M_PI*i/N);
}
// Analytical dU/dx
for (auto&& entry : index(expected_U_x)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
val = (6.0*M_PI/L)*cos(6.0*M_PI*j/N);
}
// Analytical dU/dy
for (auto&& entry : index(expected_U_y)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
val = (-4.0*M_PI/L)*sin(4.0*M_PI*i/N);
}
// Matrix of wavenumbers
Matrix<std::complex<int>> Freq(N);
Freq = FFT::computeFrequencies(N);
// Transforming U to spectral space
spectral_U = FFT::transform(U);
// FFT(dU/dx)=(2.pi.I.w)FFT(U)
for (auto&& entry : index(spectral_U_x)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
val = (2.0*M_PI*_I_)*(Freq(i,j).imag()*dw)*spectral_U(i,j);
}
// FFT(dU/dy)=(2.pi.I.s)FFT(U)
for (auto&& entry : index(spectral_U_y)) {
int i = std::get<0>(entry);
int j = std::get<1>(entry);
auto& val = std::get<2>(entry);
val = (2.0*M_PI*_I_)*(Freq(i,j).real()*dw)*spectral_U(i,j);
}
// dU/d*=IFFT(FFT(dU/d*))
computed_U_x = FFT::itransform(spectral_U_x);
computed_U_y = FFT::itransform(spectral_U_y);
// Testing equality
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
ASSERT_NEAR(std::abs(expected_U_x(i,j)-computed_U_x(i,j)), 0, 1e-11);
ASSERT_NEAR(std::abs(expected_U_y(i,j)-computed_U_y(i,j)), 0, 1e-11);
}
}
}
/* ------------------------------------------------------ */
c4science · Help