Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/tests/test_fft.cpp
9731ad93e058features/netcdf-mpi

Last Change
Raw File

test_fft.cpp
View Options

/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "fft_plan_manager.hh"
#include "grid.hh"
#include "grid_hermitian.hh"
#include "grid_view.hh"
#include "test.hh"
#ifdef USE_PYTHON
#include <pybind11/embed.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#endif
using namespace tamaas;
using fft = fftw::helper<Real>;
/* -------------------------------------------------------------------------- */
template <typename T>
struct span {
T* data;
std::size_t size;
const T* begin() const { return data; }
const T* end() const { return data + size; }
};
TEST(TestFFTInterface, FFT1D) {
constexpr UInt size = 100000;
Real* data = fft::alloc_real(size);
fft::complex* solution = fft::alloc_complex(size / 2 + 1);
std::iota(data, data + size, 0);
auto solution_plan =
fftw::plan_1d_forward(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, span<fft::complex>{solution, size / 2 + 1},
AreComplexEqual()))
<< "1D FFTW transform failed";
fftw::destroy(solution_plan);
FFTPlanManager::get().destroyPlan(grid, result);
}
/* -------------------------------------------------------------------------- */
TEST(TestFFTInterface, FFT2D) {
constexpr UInt size = 10;
Real* data = fft::alloc_real(size * size);
fft::complex* solution = fft::alloc_complex(size * (size / 2 + 1));
std::iota(data, data + (size * size), 0);
auto solution_plan =
fftw::plan_2d_forward(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
for (UInt i = 0; i < size * (size / 2+ 1); ++i) {
std::cout << result(i) << " " << solution[i][0] << ' ' << solution[i][1] << '\n';
}
ASSERT_TRUE(compare(result,
span<fft::complex>{solution, size * (size / 2 + 1)},
AreComplexEqual()))
<< "2D FFTW transform failed";
fftw::destroy(solution_plan);
FFTPlanManager::get().destroyPlan(grid, result);
}
/* -------------------------------------------------------------------------- */
TEST(TestFFTInterface, FFT1D2Comp) {
constexpr UInt size = 20;
/// 1D single component FFT should be working here
Grid<Real, 1> grid({size}, 2), data({size}, 1);
std::iota(grid.begin(), grid.end(), 0);
std::iota(data.begin(), data.end(), 0);
GridHermitian<Real, 1> result({size / 2 + 1}, 2), solution({size / 2 + 1}, 1);
FFTPlanManager::get().createPlan(grid, result).forwardTransform();
#ifdef USE_CUDA
cudaDeviceSynchronize();
#endif
auto& plan = FFTPlanManager::get().createPlan(data, solution);
std::iota(data.begin(), data.end(), 0);
data *= 2;
plan.forwardTransform();
ASSERT_TRUE(
compare(make_component_view(result, 0), solution, AreComplexEqual()))
<< "1D FFTW transform with 2 components failed on 1st component";
data += 1;
plan.forwardTransform();
ASSERT_TRUE(
compare(make_component_view(result, 1), solution, AreComplexEqual()))
<< "1D FFTW transform with 2 components failed on 2nd component";
FFTPlanManager::get().destroyPlan(grid, result);
}
/* -------------------------------------------------------------------------- */
TEST(TestFFTInterface, FFT2D3Comp) {
constexpr UInt size = 20;
/// 2D single component FFT should be working here
Grid<Real, 2> grid({size, size}, 3), data({size, size}, 1);
std::iota(grid.begin(), grid.end(), 0);
std::iota(data.begin(), data.end(), 0);
data *= 3;
GridHermitian<Real, 2> result({size, size / 2 + 1}, 3),
solution({size, size / 2 + 1}, 1);
FFTPlanManager::get().createPlan(grid, result).forwardTransform();
auto& plan = FFTPlanManager::get().createPlan(data, solution);
#ifdef USE_CUDA
cudaDeviceSynchronize();
#endif
for (UInt i = 0; i < 3; ++i) {
plan.forwardTransform();
ASSERT_TRUE(
compare(make_component_view(result, i), solution, AreComplexEqual()))
<< "2D FFTW transform with 3 components failed on " << i
<< "th component";
data += 1;
}
FFTPlanManager::get().destroyPlan(grid, result);
}
/* -------------------------------------------------------------------------- */
TEST(TestFFTInterface, FFT2DViewTransform) {
constexpr UInt size = 20;
Grid<Real, 2> data({size, size}, 1);
GridHermitian<Real, 2> solution({size, size / 2 + 1}, 1);
std::iota(std::begin(data), std::end(data), 0);
FFTPlanManager::get().createPlan(data, solution).forwardTransform();
Grid<Real, 2> grid({size, size}, 3);
auto view = make_component_view(grid, 1);
std::iota(view.begin(), view.end(), 0);
GridHermitian<Real, 2> result({size, size / 2 + 1}, 1);
FFTPlanManager::get().createPlan(view, result).forwardTransform();
ASSERT_TRUE(compare(result, solution, AreComplexEqual()))
<< "Fourier transform on component view fail";
}
/* -------------------------------------------------------------------------- */
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