Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_fft.cpp
98f53bde7826features/netcdf-mpi

Last Change
Raw File

test_fft.cpp
View Options

/**
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "test.hh"
#include "grid.hh"
#include "grid_hermitian.hh"
#include "fft_plan_manager.hh"
#ifdef USE_PYTHON
#include <pybind11/pybind11.h>
#include <pybind11/embed.h>
#include <pybind11/numpy.h>
#endif
using namespace tamaas;
/* -------------------------------------------------------------------------- */
TEST(TestFFTInterface, FFT1D) {
constexpr UInt size = 20;
double data[size] = {0};
fftw_complex solution[size/2 + 1] = {{0}};
std::iota(std::begin(data), std::end(data), 0);
fftw_plan solution_plan = fftw_plan_dft_r2c_1d(size, data, solution, FFTW_ESTIMATE);
fftw_execute(solution_plan);
Grid<Real, 1> grid({size}, 1);
std::iota(grid.begin(), grid.end(), 0);
GridHermitian<Real, 1> result({size/2 + 1}, 1);
FFTPlanManager::get().createPlan(grid, result).forwardTransform();
#ifdef USE_CUDA
cudaDeviceSynchronize();
#endif
ASSERT_TRUE(compare(result, solution, AreComplexEqual())) << "1D FFTW transform failed";
fftw_destroy_plan(solution_plan);
FFTPlanManager::get().destroyPlan(grid, result);
}
/* -------------------------------------------------------------------------- */
TEST(TestFFTInterface, FFT2D) {
constexpr UInt size = 10;
double data[size * size] = {0};
fftw_complex solution[size * (size/2 + 1)] = {{0}};
std::iota(std::begin(data), std::end(data), 0);
fftw_plan solution_plan = fftw_plan_dft_r2c_2d(size, size, data, solution, FFTW_ESTIMATE);
fftw_execute(solution_plan);
Grid<Real, 2> grid({size, size}, 1);
std::iota(grid.begin(), grid.end(), 0);
GridHermitian<Real, 2> result({size, size/2 + 1}, 1);
FFTPlanManager::get().createPlan(grid, result).forwardTransform();
#ifdef USE_CUDA
cudaDeviceSynchronize();
#endif
ASSERT_TRUE(compare(result, solution, AreComplexEqual())) << "2D FFTW transform failed";
fftw_destroy_plan(solution_plan);
FFTPlanManager::get().destroyPlan(grid, result);
}
/* -------------------------------------------------------------------------- */
TEST(TestFFTInterface, FFT1D2Comp) {
constexpr UInt size = 20;
constexpr UInt size_c = size/2 + 1;
double data[2*size] = {0};
fftw_complex solution[2 * size_c] = {{0}};
fftw_complex solution_reordered[2 * size_c] = {{0}};
for (UInt i = 0; i < size; ++i) {
data[i] = 2*i;
data[size + i] = 2*i + 1;
}
for (UInt i = 0; i < 2; i++) {
fftw_plan solution_plan = fftw_plan_dft_r2c_1d(size, data + i * size,
solution + i * size_c,
FFTW_ESTIMATE);
fftw_execute(solution_plan);
fftw_destroy_plan(solution_plan);
}
for (UInt i = 0; i < size_c; ++i) {
for (UInt j = 0; j < 2; ++j) {
solution_reordered[2*i][j] = solution[i][j];
solution_reordered[2*i + 1][j] = solution[i + size_c][j];
}
}
Grid<Real, 1> grid({size}, 2);
std::iota(grid.begin(), grid.end(), 0);
GridHermitian<Real, 1> result({size_c}, 2);
FFTPlanManager::get().createPlan(grid, result).forwardTransform();
#ifdef USE_CUDA
cudaDeviceSynchronize();
#endif
ASSERT_TRUE(compare(result, solution_reordered, AreComplexEqual())) << \
"1D FFTW transform with 2 components failed";
FFTPlanManager::get().destroyPlan(grid, result);
}
/* -------------------------------------------------------------------------- */
TEST(TestFFTInterface, FFT2D3Comp) {
constexpr UInt size = 20;
constexpr UInt size_c = size/2 + 1;
double data[3 * size * size] = {0};
fftw_complex solution[3 * size * size_c] = {{0}};
fftw_complex solution_reordered[3 * size * size_c] = {{0}};
for (UInt i = 0; i < size * size; ++i) {
data[i] = 3*i;
data[size * size + i] = 3*i + 1;
data[2 * size * size + i] = 3*i + 2;
}
for (UInt i = 0; i < 3; i++) {
fftw_plan solution_plan = fftw_plan_dft_r2c_2d(size, size,
data + i * size * size,
solution + i * size * size_c,
FFTW_ESTIMATE);
fftw_execute(solution_plan);
fftw_destroy_plan(solution_plan);
}
for (UInt i = 0; i < size * size_c; ++i) {
for (UInt j = 0; j < 2; ++j) {
solution_reordered[3*i][j] = solution[i][j];
solution_reordered[3*i + 1][j] = solution[i + size * size_c][j];
solution_reordered[3*i + 2][j] = solution[i + 2 * size * size_c][j];
}
}
Grid<Real, 2> grid({size, size}, 3);
std::iota(grid.begin(), grid.end(), 0);
GridHermitian<Real, 2> result({size, size_c}, 3);
FFTPlanManager::get().createPlan(grid, result).forwardTransform();
#ifdef USE_CUDA
cudaDeviceSynchronize();
#endif
ASSERT_TRUE(compare(result, solution_reordered, AreComplexEqual())) << \
"2D FFTW transform with 3 components failed";
FFTPlanManager::get().destroyPlan(grid, result);
}
/* -------------------------------------------------------------------------- */
TEST(TestFFTInterface, FFTI1D2Comp) {
constexpr UInt size = 20;
Grid<Real, 1> grid({size}, 2);
std::iota(grid.begin(), grid.end(), 0);
GridHermitian<Real, 1> grid_hermitian({size/2 + 1}, 2);
Grid<Real, 1> result({size}, 2);
FFTPlanManager::get().createPlan(grid, grid_hermitian).forwardTransform();
FFTPlanManager::get().createPlan(result, grid_hermitian).backwardTransform();
#ifdef USE_CUDA
cudaDeviceSynchronize();
#endif
ASSERT_TRUE(compare(grid, result, AreFloatEqual())) << \
"1D FFTI transform with 2 components failed";
FFTPlanManager::get().destroyPlan(grid, grid_hermitian);
FFTPlanManager::get().destroyPlan(result, grid_hermitian);
}
/* -------------------------------------------------------------------------- */
TEST(TestFFTInterface, FFTI2D3Comp) {
constexpr UInt size = 20;
Grid<Real, 2> grid({size, size}, 3);
std::iota(grid.begin(), grid.end(), 0);
GridHermitian<Real, 2> grid_hermitian({size, size/2 + 1}, 3);
Grid<Real, 2> result({size, size}, 3);
FFTPlanManager::get().createPlan(grid, grid_hermitian).forwardTransform();
FFTPlanManager::get().createPlan(result, grid_hermitian).backwardTransform();
#ifdef USE_CUDA
cudaDeviceSynchronize();
#endif
ASSERT_TRUE(compare(grid, result, AreFloatEqual())) << \
"2D FFTI transform with 3 components failed";
FFTPlanManager::get().destroyPlan(grid, grid_hermitian);
FFTPlanManager::get().destroyPlan(result, grid_hermitian);
}
/* -------------------------------------------------------------------------- */
#ifdef USE_PYTHON
namespace py = pybind11;
TEST(TestFFTInterface, Frequencies1D) {
std::array<UInt, 1> sizes = {{10}};
auto freq = FFTransform<Real, 1>::computeFrequencies<false>(sizes);
py::module fftfreq = py::module::import("fftfreq");
std::vector<Real> reference(freq.dataSize());
py::array py_arg(reference.size(), reference.data(), py::none());
fftfreq.attr("frequencies1D")(py_arg);
ASSERT_TRUE(compare(reference, freq, AreFloatEqual()))
<< "Non hermitian frequencies are wrong";
auto hfreq = FFTransform<Real, 1>::computeFrequencies<true>(sizes);
std::iota(reference.begin(), reference.end(), 0);
ASSERT_TRUE(compare(reference, hfreq, AreFloatEqual()))
<< "Hermitian frequencies are wrong";
}
TEST(TestFFTInterface, Frequencies2D) {
std::array<UInt, 2> sizes = {{10, 10}};
auto freq = FFTransform<Real, 2>::computeFrequencies<false>(sizes);
py::module fftfreq = py::module::import("fftfreq");
std::vector<Real> reference(freq.dataSize());
py::array py_arg({10, 10, 2}, reference.data(), py::none());
fftfreq.attr("frequencies2D")(py_arg);
ASSERT_TRUE(compare(reference, freq, AreFloatEqual()))
<< "Non hermitian frequencies are wrong";
auto hfreq = FFTransform<Real, 2>::computeFrequencies<true>(sizes);
fftfreq.attr("hfrequencies2D")(py_arg);
ASSERT_TRUE(compare(reference, hfreq, AreFloatEqual()))
<< "Hermitian frequencies are wrong";
}
#endif
c4science · Help