diff --git a/src/fft/pfft_engine.cc b/src/fft/pfft_engine.cc index 01ed8dc..7e0dbc3 100644 --- a/src/fft/pfft_engine.cc +++ b/src/fft/pfft_engine.cc @@ -1,205 +1,205 @@ /** * @file pfft_engine.cc * * @author Lars Pastewka * * @date 06 Mar 2017 * * @brief implements the MPI-parallel pfft engine * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/ccoord_operations.hh" #include "fft/pfft_engine.hh" namespace muSpectre { template int PFFTEngine::nb_engines{0}; template PFFTEngine::PFFTEngine(Ccoord resolutions, Rcoord lengths, Communicator comm) :Parent{resolutions, lengths, comm} { if (!this->nb_engines) pfft_init(); this->nb_engines++; std::array narr; std::copy(this->domain_resolutions.begin(), this->domain_resolutions.end(), narr.begin()); ptrdiff_t res[DimS], loc[DimS], fres[DimS], floc[DimS]; this->workspace_size = pfft_local_size_many_dft_r2c(DimS, narr.data(), narr.data(), narr.data(), Field_t::nb_components, PFFT_DEFAULT_BLOCK, PFFT_DEFAULT_BLOCKS, this->comm.get_mpi_comm(), PFFT_TRANSPOSED_OUT, res, loc, fres, floc); std::copy(res, res+DimS, this->resolutions.begin()); std::copy(loc, loc+DimS, this->locations.begin()); std::copy(fres, fres+DimS, this->fourier_resolutions.begin()); std::copy(floc, floc+DimS, this->fourier_locations.begin()); for (auto & n: this->resolutions) { if (n == 0) { throw std::runtime_error("PFFT planning returned zero resolution. " "You may need to run on fewer processes."); } } for (auto & n: this->fourier_resolutions) { if (n == 0) { throw std::runtime_error("PFFT planning returned zero Fourier " "resolution. You may need to run on fewer " "processes."); } } - for (auto && pixel: CcoordOps::Pixels(this->fourier_resolutions, - this->fourier_locations)) { + for (auto && pixel: CcoordOps::Pixels(this->fourier_resolutions, + this->fourier_locations)) { this->work_space_container.add_pixel(pixel); } } /* ---------------------------------------------------------------------- */ template void PFFTEngine::initialise(FFT_PlanFlags plan_flags) { if (this->initialised) { throw std::runtime_error("double initialisation, will leak memory"); } // Initialize parent after local resolutions have been determined and // work space has been initialized Parent::initialise(plan_flags); this->real_workspace = pfft_alloc_real(2*this->workspace_size); // We need to check whether the workspace provided by our field is large // enough. PFFT may request a workspace size larger than the nominal size // of the complex buffer. if (long(this->work.size()*Field_t::nb_components) < this->workspace_size) { this->work.set_pad_size(this->workspace_size - Field_t::nb_components*this->work.size()); } unsigned int flags; switch (plan_flags) { case FFT_PlanFlags::estimate: { flags = PFFT_ESTIMATE; break; } case FFT_PlanFlags::measure: { flags = PFFT_MEASURE; break; } case FFT_PlanFlags::patient: { flags = PFFT_PATIENT; break; } default: throw std::runtime_error("unknown planner flag type"); break; } std::array narr; std::copy(this->domain_resolutions.begin(), this->domain_resolutions.end(), narr.begin()); Real * in = this->real_workspace; pfft_complex * out = reinterpret_cast(this->work.data()); this->plan_fft = pfft_plan_many_dft_r2c(DimS, narr.data(), narr.data(), narr.data(), Field_t::nb_components, PFFT_DEFAULT_BLOCKS, PFFT_DEFAULT_BLOCKS, in, out, this->comm.get_mpi_comm(), PFFT_FORWARD, - PFFT_TRANSPOSED_OUT | flags); + flags); if (this->plan_fft == nullptr) { throw std::runtime_error("r2c plan failed"); } pfft_complex * i_in = reinterpret_cast(this->work.data()); Real * i_out = this->real_workspace; this->plan_ifft = pfft_plan_many_dft_c2r(DimS, narr.data(), narr.data(), narr.data(), Field_t::nb_components, PFFT_DEFAULT_BLOCKS, PFFT_DEFAULT_BLOCKS, i_in, i_out, this->comm.get_mpi_comm(), PFFT_BACKWARD, - PFFT_TRANSPOSED_IN | flags); + flags); if (this->plan_ifft == nullptr) { throw std::runtime_error("c2r plan failed"); } this->initialised = true; } /* ---------------------------------------------------------------------- */ template PFFTEngine::~PFFTEngine() noexcept { if (this->real_workspace != nullptr) pfft_free(this->real_workspace); if (this->plan_fft != nullptr) pfft_destroy_plan(this->plan_fft); if (this->plan_ifft != nullptr) pfft_destroy_plan(this->plan_ifft); this->nb_engines--; if (!this->nb_engines) pfft_cleanup(); } /* ---------------------------------------------------------------------- */ template typename PFFTEngine::Workspace_t & PFFTEngine::fft (Field_t & field) { if (!this->plan_fft) { throw std::runtime_error("fft plan not allocated"); } if (field.size() != CcoordOps::get_size(this->resolutions)) { throw std::runtime_error("size mismatch"); } // Copy field data to workspace buffer. This is necessary because workspace // buffer is larger than field buffer. std::copy(field.data(), field.data()+Field_t::nb_components*field.size(), this->real_workspace); pfft_execute_dft_r2c( this->plan_fft, this->real_workspace, reinterpret_cast(this->work.data())); return this->work; } /* ---------------------------------------------------------------------- */ template void PFFTEngine::ifft (Field_t & field) const { if (!this->plan_ifft) { throw std::runtime_error("ifft plan not allocated"); } if (field.size() != CcoordOps::get_size(this->resolutions)) { throw std::runtime_error("size mismatch"); } pfft_execute_dft_c2r( this->plan_ifft, reinterpret_cast(this->work.data()), this->real_workspace); std::copy(this->real_workspace, this->real_workspace+Field_t::nb_components*field.size(), field.data()); } template class PFFTEngine; template class PFFTEngine; template class PFFTEngine; } // muSpectre diff --git a/tests/mpi_test_projection.hh b/tests/mpi_test_projection.hh index f0fddbb..0fc3130 100644 --- a/tests/mpi_test_projection.hh +++ b/tests/mpi_test_projection.hh @@ -1,88 +1,86 @@ /** * @file test_projection.hh * * @author Till Junge * * @date 16 Jan 2018 * * @brief common declarations for testing both the small and finite strain * projection operators * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "tests.hh" #include "mpi_context.hh" -#include "fft/fftwmpi_engine.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template struct Sizes { }; template<> struct Sizes { constexpr static Ccoord_t get_resolution() { return Ccoord_t{3, 5};} constexpr static Rcoord_t get_lengths() { return Rcoord_t{3.4, 5.8};} }; template<> struct Sizes { constexpr static Ccoord_t get_resolution() { return Ccoord_t{3, 5, 7};} constexpr static Rcoord_t get_lengths() { return Rcoord_t{3.4, 5.8, 6.7};} }; template struct Squares { }; template<> struct Squares { constexpr static Ccoord_t get_resolution() { return Ccoord_t{5, 5};} constexpr static Rcoord_t get_lengths() { return Rcoord_t{5, 5};} }; template<> struct Squares { constexpr static Ccoord_t get_resolution() { return Ccoord_t{7, 7, 7};} constexpr static Rcoord_t get_lengths() { return Rcoord_t{7, 7, 7};} }; /* ---------------------------------------------------------------------- */ - template + template struct ProjectionFixture { - using Engine = FFTWMPIEngine; using Parent = Proj; constexpr static Dim_t sdim{DimS}; constexpr static Dim_t mdim{DimM}; ProjectionFixture() :projector(std::make_unique(SizeGiver::get_resolution(), SizeGiver::get_lengths(), MPIContext::get_context().comm)){} Parent projector; }; } // muSpectre diff --git a/tests/mpi_test_projection_finite.cc b/tests/mpi_test_projection_finite.cc index bdb3500..c115b1a 100644 --- a/tests/mpi_test_projection_finite.cc +++ b/tests/mpi_test_projection_finite.cc @@ -1,126 +1,171 @@ /** * @file test_projection_finite.cc * * @author Till Junge * * @date 07 Dec 2017 * * @brief tests for standard finite strain projection operator * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "fft/projection_finite_strain.hh" #include "fft/projection_finite_strain_fast.hh" #include "fft/fft_utils.hh" #include "mpi_test_projection.hh" +#ifdef WITH_FFTWMPI +#include "fft/fftwmpi_engine.hh" +#endif +#ifdef WITH_PFFT +#include "fft/pfft_engine.hh" +#endif + #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(mpi_projection_finite_strain); /* ---------------------------------------------------------------------- */ using fixlist = boost::mpl::list< +#ifdef WITH_FFTWMPI + ProjectionFixture, + ProjectionFiniteStrain, + FFTWMPIEngine>, + ProjectionFixture, + ProjectionFiniteStrain, + FFTWMPIEngine>, + ProjectionFixture, + ProjectionFiniteStrain, + FFTWMPIEngine>, + ProjectionFixture, + ProjectionFiniteStrain, + FFTWMPIEngine>, + + ProjectionFixture, + ProjectionFiniteStrainFast, + FFTWMPIEngine>, + ProjectionFixture, + ProjectionFiniteStrainFast, + FFTWMPIEngine>, + ProjectionFixture, + ProjectionFiniteStrainFast, + FFTWMPIEngine>, + ProjectionFixture, + ProjectionFiniteStrainFast, + FFTWMPIEngine>, +#endif +#ifdef WITH_PFFT ProjectionFixture, - ProjectionFiniteStrain>, + ProjectionFiniteStrain, + PFFTEngine>, ProjectionFixture, - ProjectionFiniteStrain>, + ProjectionFiniteStrain, + PFFTEngine>, ProjectionFixture, - ProjectionFiniteStrain>, + ProjectionFiniteStrain, + PFFTEngine>, ProjectionFixture, - ProjectionFiniteStrain>, + ProjectionFiniteStrain, + PFFTEngine>, ProjectionFixture, - ProjectionFiniteStrainFast>, + ProjectionFiniteStrainFast, + PFFTEngine>, ProjectionFixture, - ProjectionFiniteStrainFast>, + ProjectionFiniteStrainFast, + PFFTEngine>, ProjectionFixture, - ProjectionFiniteStrainFast>, + ProjectionFiniteStrainFast, + PFFTEngine>, ProjectionFixture, - ProjectionFiniteStrainFast>>; + ProjectionFiniteStrainFast, + PFFTEngine> +#endif + >; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test, fix, fixlist, fix) { // create a gradient field with a zero mean gradient and verify // that the projection preserves it constexpr Dim_t dim{fix::sdim}, sdim{fix::sdim}, mdim{fix::mdim}; static_assert(dim == fix::mdim, "These tests assume that the material and spatial dimension are " "identical"); using Fields = GlobalFieldCollection; using FieldT = TensorField; using FieldMap = MatrixFieldMap; using Vector = Eigen::Matrix; Fields fields{}; FieldT & f_grad{make_field("gradient", fields)}; FieldT & f_var{make_field("working field", fields)}; FieldMap grad(f_grad); FieldMap var(f_var); fields.initialise(fix::projector.get_resolutions(), fix::projector.get_locations()); Vector k; for (Dim_t i = 0; i < dim; ++i) { // the wave vector has to be such that it leads to an integer // number of periods in each length of the domain k(i) = (i+1)*2*pi/fix::projector.get_lengths()[i]; } for (auto && tup: akantu::zip(fields, grad, var)) { auto & ccoord = std::get<0>(tup); auto & g = std::get<1>(tup); auto & v = std::get<2>(tup); Vector vec = CcoordOps::get_vector(ccoord, fix::projector.get_lengths()/ fix::projector.get_domain_resolutions()); g.row(0) = k.transpose() * cos(k.dot(vec)); v.row(0) = g.row(0); } fix::projector.initialise(FFT_PlanFlags::estimate); fix::projector.apply_projection(f_var); for (auto && tup: akantu::zip(fields, grad, var)) { auto & ccoord = std::get<0>(tup); auto & g = std::get<1>(tup); auto & v = std::get<2>(tup); Real error = (g-v).norm(); BOOST_CHECK_LT(error, tol); if (error >=tol) { Vector vec = CcoordOps::get_vector(ccoord, fix::projector.get_lengths()/ fix::projector.get_domain_resolutions()); std::cout << std::endl << "grad_ref :" << std::endl << g << std::endl; std::cout << std::endl << "grad_proj :" << std::endl << v << std::endl; std::cout << std::endl << "ccoord :" << std::endl << ccoord << std::endl; std::cout << std::endl << "vector :" << std::endl << vec.transpose() << std::endl; } } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre