diff --git a/CMakeLists.txt b/CMakeLists.txt index 4692a90..0364085 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,167 +1,175 @@ # ============================================================================= # file CMakeLists.txt # # @author Till Junge # # @date 08 Jan 2018 # # @brief Main configuration file # # @section LICENSE # # 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. # ============================================================================= cmake_minimum_required(VERSION 3.5.0) project(µSpectre) set(CMAKE_CXX_STANDARD 14) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) set(BUILD_SHARED_LIBS ON) add_compile_options(-Wall -Wextra -Weffc++) enable_testing() set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) include(muspectreTools) if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang") # using Clang add_compile_options(-Wno-missing-braces) elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") # using GCC add_compile_options(-Wno-non-virtual-dtor) string( TOLOWER "${CMAKE_BUILD_TYPE}" build_type ) if ("release" STREQUAL "${build_type}" ) add_compile_options(-march=native) endif() if ("debug" STREQUAL "${build_type}" ) add_compile_options(-O0) endif() elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") # using Intel C++ elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") # using Visual Studio C++ endif() +# Do not trust old gcc. the std::optional has memory bugs +if(${CMAKE_COMPILER_IS_GNUCC}) + if(${CMAKE_CXX_COMPILER_VERSION} VERSION_LESS 6.0.0) + add_definitions(-DNO_EXPERIMENTAL) + endif() +endif() + + find_package(Boost COMPONENTS unit_test_framework REQUIRED) add_external_package(Eigen3 VERSION 3.3 CONFIG) add_external_package(pybind11 VERSION 2.2 CONFIG) include_directories( ${CMAKE_SOURCE_DIR}/tests ${CMAKE_SOURCE_DIR}/src ${CMAKE_SOURCE_DIR} ) if(APPLE) include_directories(${CMAKE_INSTALL_PREFIX}/include ${Boost_INCLUDE_DIRS}) endif() find_package(MPI) find_package(FFTWMPI) find_package(PFFT) if(MPI_FOUND) add_definitions(-DWITH_MPI) include_directories(SYSTEM ${MPI_C_INCLUDE_PATH}) endif() if (FFTWMPI_FOUND) add_definitions(-DWITH_FFTWMPI) endif() if (PFFT_FOUND) add_definitions(-DWITH_PFFT) endif() if (P3DFFT_FOUND) add_definitions(-DWITH_P3DFFT) endif() #build tests file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/test_*.cc") add_executable(main_test_suite tests/main_test_suite.cc ${TEST_SRCS}) target_link_libraries(main_test_suite ${Boost_LIBRARIES} muSpectre) add_test(main_test_suite main_test_suite --report_level=detailed --build_info=TRUE) if(MPI_FOUND) file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/mpi_test_*.cc") add_executable(mpi_main_test_suite tests/mpi_main_test_suite.cc ${TEST_SRCS}) target_link_libraries(mpi_main_test_suite ${Boost_LIBRARIES} muSpectre) add_test(mpi_main_test_suite ${MPIEXEC_EXECUTABLE} ${MPIEXEC_PREFLAGS} ${MPIEXEC_NUMPROC_FLAG} 2 mpi_main_test_suite --report_level=detailed --build_info=TRUE) endif() # compile the library add_compile_options( -Werror) add_subdirectory( ${CMAKE_SOURCE_DIR}/src/ ) add_subdirectory( ${CMAKE_SOURCE_DIR}/language_bindings/ ) add_subdirectory( ${CMAKE_SOURCE_DIR}/doc/ ) #compile executables file( GLOB binaries "${CMAKE_SOURCE_DIR}/bin/*.cc") foreach(binaryfile ${binaries}) get_filename_component(binaryname ${binaryfile} NAME_WE) add_executable(${binaryname} ${binaryfile}) target_link_libraries(${binaryname} ${Boost_LIBRARIES} muSpectre) endforeach(binaryfile ${binaries}) #or copy them file (GLOB pybins "${CMAKE_SOURCE_DIR}/bin/*.py") foreach(pybin ${pybins}) get_filename_component(binaryname ${pybin} NAME_WE) configure_file( ${pybin} "${CMAKE_BINARY_DIR}/${binaryname}.py" COPYONLY) endforeach(pybin ${pybins}) # compile benchmarks file(GLOB benchmarks "${CMAKE_SOURCE_DIR}/benchmarks/benchmark*cc") foreach(benchmark ${benchmarks}) get_filename_component(benchmark_name ${benchmark} NAME_WE) add_executable(${benchmark_name} ${benchmark}) target_link_libraries(${benchmark_name} ${BOOST_LIBRARIES} muSpectre) endforeach(benchmark ${benchmark}) # copy python test #build tests file( GLOB PY_TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/python_*.py") foreach(pytest ${PY_TEST_SRCS}) get_filename_component(pytest_name ${pytest} NAME) configure_file( ${pytest} "${CMAKE_BINARY_DIR}/${pytest_name}" COPYONLY) endforeach(pytest ${PY_TEST_SRCS}) add_test(python_binding_test python_binding_tests.py) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e9dea52..dff9aa7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,101 +1,106 @@ # ============================================================================= # file CMakeLists.txt # # @author Till Junge # # @date 08 Jan 2018 # # @brief Configuration for libmuSpectre # # @section LICENSE # # 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_directories( common materials cell ) set (muSpectre_SRC ) add_subdirectory(common) add_subdirectory(materials) add_subdirectory(fft) add_subdirectory(cell) add_subdirectory(solver) find_package(MPI) find_package(FFTW REQUIRED) find_package(FFTWMPI) find_package(PFFT) # The following checks whether std::optional exists and replaces it by # boost::optional if necessary include(CheckCXXSourceCompiles) check_cxx_source_compiles( "#include int main() { std::experimental::optional A{}; }" HAS_STD_OPTIONAL) +add_definitions(-DBAR) if( NOT HAS_STD_OPTIONAL) add_definitions(-DNO_EXPERIMENTAL) endif() file(GLOB_RECURSE _headers RELATIVE "${CMAKE_CURRENT_SOURCE_DIR}" "*.hh") list(APPEND muSpectre_SRC ${_headers}) add_library(muSpectre ${muSpectre_SRC}) + set(muSpectre_INCLUDES ${FFTW_INCLUDES}) set(muSpectre_LIBRARIES Eigen3::Eigen ${FFTW_LIBRARIES}) if(FFTWMPI_FOUND) set(muSpectre_INCLUDES ${muSpectre_INCLUDES} ${FFTWMPI_INCLUDES}) set(muSpectre_LIBRARIES ${muSpectre_LIBRARIES} ${FFTWMPI_LIBRARIES}) endif() if(PFFT_FOUND) set(muSpectre_INCLUDES ${muSpectre_INCLUDES} ${PFFT_INCLUDES}) set(muSpectre_LIBRARIES ${muSpectre_LIBRARIES} ${PFFT_LIBRARIES}) endif() if(FFTWMPI_FOUND OR PFFT_FOUND) set(muSpectre_INCLUDES ${muSpectre_INCLUDES} ${MPI_INCLUDES}) set(muSpectre_LIBRARIES ${muSpectre_LIBRARIES} ${MPI_LIBRARIES}) endif() target_include_directories(muSpectre INTERFACE ${muSpectre_INCLUDES}) target_link_libraries(muSpectre ${muSpectre_LIBRARIES}) +target_link_libraries(muSpectre Eigen3::Eigen ${FFTW_LIBRARIES}) +target_compile_definitions(muSpectre PRIVATE ${-DNO_EXPERIMENTAL}) + set_property(TARGET muSpectre PROPERTY PUBLIC_HEADER ${_headers}) install(TARGETS muSpectre LIBRARY DESTINATION lib PUBLIC_HEADER DESTINATION include) \ No newline at end of file diff --git a/src/fft/pfft_engine.cc b/src/fft/pfft_engine.cc index b827223..3b617bd 100644 --- a/src/fft/pfft_engine.cc +++ b/src/fft/pfft_engine.cc @@ -1,242 +1,239 @@ /** * @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} + :Parent{resolutions, lengths, comm}, mpi_comm{comm.get_mpi_comm()} { if (!this->nb_engines) pfft_init(); this->nb_engines++; int size{comm.size()}; int dimx{size}; int dimy{1}; //if (DimS > 2) { if (false) { dimy = int(sqrt(size)); while ((size/dimy)*dimy != size) dimy--; dimx = size/dimy; } std::cout << "PFFT process mesh: " << dimx << "x" << dimy << std::endl; //if (DimS > 2) { if (false) { if (pfft_create_procmesh_2d(this->comm.get_mpi_comm(), dimx, dimy, &this->mpi_comm)) { throw std::runtime_error("Failed to create 2d PFFT process mesh."); } } - else { - this->mpi_comm = this->comm.get_mpi_comm(); - } // FIXME! Code presently always use a one-dimensional process mesh but // should use a two-dimensional process mesh for three-dimensional data. 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->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 (int i = 0; i < DimS-1; i++) { for (int i = 0; i < 1; i++) { std::swap(this->fourier_resolutions[i], this->fourier_resolutions[i+1]); std::swap(this->fourier_locations[i], this->fourier_locations[i+1]); } 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: std::conditional_t< DimS==2, CcoordOps::Pixels, 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->mpi_comm, PFFT_FORWARD, PFFT_TRANSPOSED_OUT | 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->mpi_comm, PFFT_BACKWARD, PFFT_TRANSPOSED_IN | 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); if (this->mpi_comm != this->comm.get_mpi_comm()) { MPI_Comm_free(&this->mpi_comm); } 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/src/fft/projection_base.hh b/src/fft/projection_base.hh index c76fb19..57a7768 100644 --- a/src/fft/projection_base.hh +++ b/src/fft/projection_base.hh @@ -1,156 +1,169 @@ /** * @file projection_base.hh * * @author Till Junge * * @date 03 Dec 2017 * * @brief Base class for Projection operators * * 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. */ #ifndef PROJECTION_BASE_H #define PROJECTION_BASE_H #include "common/common.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "fft/fft_engine_base.hh" #include namespace muSpectre { + /** + * base class for projection related exceptions + */ + class ProjectionError: public std::runtime_error { + public: + //! constructor + explicit ProjectionError(const std::string& what) + :std::runtime_error(what){} + //! constructor + explicit ProjectionError(const char * what) + :std::runtime_error(what){} + }; + template struct Projection_traits { }; /** * defines the interface which must be implemented by projection operators */ template class ProjectionBase { public: //! type of fft_engine used using FFTEngine = FFTEngineBase; //! reference to fft engine is safely managed through a `std::unique_ptr` using FFTEngine_ptr = std::unique_ptr; //! cell coordinates type using Ccoord = typename FFTEngine::Ccoord; //! spatial coordinates type using Rcoord = typename FFTEngine::Rcoord; //! global FieldCollection using GFieldCollection_t = typename FFTEngine::GFieldCollection_t; //! local FieldCollection (for Fourier-space pixels) using LFieldCollection_t = typename FFTEngine::LFieldCollection_t; //! Field type on which to apply the projection using Field_t = typename FFTEngine::Field_t; /** * iterator over all pixels. This is taken from the FFT engine, * because depending on the real-to-complex FFT employed, only * roughly half of the pixels are present in Fourier space * (because of the hermitian nature of the transform) */ using iterator = typename FFTEngine::iterator; //! Default constructor ProjectionBase() = delete; //! Constructor with cell sizes ProjectionBase(FFTEngine_ptr engine, Formulation form); //! Copy constructor ProjectionBase(const ProjectionBase &other) = delete; //! Move constructor ProjectionBase(ProjectionBase &&other) = default; //! Destructor virtual ~ProjectionBase() = default; //! Copy assignment operator ProjectionBase& operator=(const ProjectionBase &other) = delete; //! Move assignment operator ProjectionBase& operator=(ProjectionBase &&other) = default; //! initialises the fft engine (plan the transform) virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate); //! apply the projection operator to a field virtual void apply_projection(Field_t & field) = 0; //! returns the process-local resolutions of the cell const Ccoord & get_resolutions() const { return this->fft_engine->get_resolutions();} //! returns the process-local locations of the cell const Ccoord & get_locations() const { return this->fft_engine->get_locations();} //! returns the resolutions of the cell const Ccoord & get_domain_resolutions() const { return this->fft_engine->get_domain_resolutions();} //! returns the physical sizes of the cell const Rcoord & get_lengths() const { return this->fft_engine->get_lengths();} /** * return the `muSpectre::Formulation` that is used in solving * this cell. This allows tho check whether a projection is * compatible with the chosen formulation */ const Formulation & get_formulation() const {return this->form;} //! return the raw projection operator. This is mainly intended //! for maintenance and debugging and should never be required in //! regular use virtual Eigen::Map get_operator() = 0; //! return the communicator object const Communicator & get_communicator() const { return this->fft_engine->get_communicator(); } protected: //! handle on the fft_engine used FFTEngine_ptr fft_engine; /** * formulation this projection can be applied to (determines * whether the projection enforces gradients, small strain tensor * or symmetric smal strain tensor */ const Formulation form; /** * A local `muSpectre::FieldCollection` to store the projection * operator per k-space point. This is a local rather than a * global collection, since the pixels considered depend on the * FFT implementation. See * http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data * for an example */ LFieldCollection_t & projection_container{}; private: }; } // muSpectre #endif /* PROJECTION_BASE_H */ diff --git a/src/fft/projection_finite_strain.cc b/src/fft/projection_finite_strain.cc index edef674..bc020c5 100644 --- a/src/fft/projection_finite_strain.cc +++ b/src/fft/projection_finite_strain.cc @@ -1,85 +1,92 @@ /** * @file projection_finite_strain.cc * * @author Till Junge * * @date 05 Dec 2017 * * @brief implementation of 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/fftw_engine.hh" #include "fft/fft_utils.hh" #include "common/field_map.hh" #include "common/tensor_algebra.hh" #include "common/iterators.hh" #include "Eigen/Dense" namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionFiniteStrain:: ProjectionFiniteStrain(FFTEngine_ptr engine) :Parent{std::move(engine), Formulation::finite_strain} - {} + { + for (auto res: this->fft_engine->get_resolutions()) { + if (res % 2 == 0) { + throw ProjectionError + ("Only an odd number of gridpoints in each direction is supported"); + } + } + } /* ---------------------------------------------------------------------- */ template void ProjectionFiniteStrain:: initialise(FFT_PlanFlags flags) { Parent::initialise(flags); FFT_freqs fft_freqs(this->fft_engine->get_domain_resolutions(), this->fft_engine->get_lengths()); for (auto && tup: akantu::zip(*this->fft_engine, this->Ghat)) { const auto & ccoord = std::get<0> (tup); auto & G = std::get<1>(tup); auto xi = fft_freqs.get_unit_xi(ccoord); //! this is simplifiable using Curnier's Méthodes numériques, 6.69(c) G = Matrices::outer_under(Matrices::I2(), xi*xi.transpose()); // The commented block below corresponds to the original // definition of the operator in de Geus et // al. (https://doi.org/10.1016/j.cma.2016.12.032). However, // they use a bizarre definition of the double contraction // between fourth-order and second-order tensors that has a // built-in transpose operation (i.e., C = A:B <-> AᵢⱼₖₗBₗₖ = // Cᵢⱼ , note the inverted ₗₖ instead of ₖₗ), here, we define // the double contraction without the transposition. As a // result, the Projection operator produces the transpose of de // Geus's // for (Dim_t im = 0; im < DimS; ++im) { // for (Dim_t j = 0; j < DimS; ++j) { // for (Dim_t l = 0; l < DimS; ++l) { // get(G, im, j, l, im) = xi(j)*xi(l); // } // } // } } if (this->get_locations() == Ccoord{}) { this->Ghat[0].setZero(); } } template class ProjectionFiniteStrain; template class ProjectionFiniteStrain; } // muSpectre diff --git a/src/fft/projection_finite_strain_fast.cc b/src/fft/projection_finite_strain_fast.cc index f1ef432..89b2cf2 100644 --- a/src/fft/projection_finite_strain_fast.cc +++ b/src/fft/projection_finite_strain_fast.cc @@ -1,85 +1,92 @@ /** * @file projection_finite_strain_fast.cc * * @author Till Junge * * @date 12 Dec 2017 * * @brief implementation for fast projection in finite strain * * 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_fast.hh" #include "fft/fft_utils.hh" #include "common/tensor_algebra.hh" #include "common/iterators.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionFiniteStrainFast:: ProjectionFiniteStrainFast(FFTEngine_ptr engine) :Parent{std::move(engine), Formulation::finite_strain}, xiField{make_field("Projection Operator", this->projection_container)}, xis(xiField) - {} + { + for (auto res: this->fft_engine->get_resolutions()) { + if (res % 2 == 0) { + throw ProjectionError + ("Only an odd number of gridpoints in each direction is supported"); + } + } + } /* ---------------------------------------------------------------------- */ template void ProjectionFiniteStrainFast:: initialise(FFT_PlanFlags flags) { Parent::initialise(flags); FFT_freqs fft_freqs(this->fft_engine->get_domain_resolutions(), this->fft_engine->get_lengths()); for (auto && tup: akantu::zip(*this->fft_engine, this->xis)) { const auto & ccoord = std::get<0> (tup); auto & xi = std::get<1>(tup); xi = fft_freqs.get_unit_xi(ccoord); } if (this->get_locations() == Ccoord{}) { this->xis[0].setZero(); } } /* ---------------------------------------------------------------------- */ template void ProjectionFiniteStrainFast::apply_projection(Field_t & field) { Grad_map field_map{this->fft_engine->fft(field)}; Real factor = this->fft_engine->normalisation(); for (auto && tup: akantu::zip(this->xis, field_map)) { auto & xi{std::get<0>(tup)}; auto & f{std::get<1>(tup)}; f = factor * ((f*xi).eval()*xi.transpose()); } this->fft_engine->ifft(field); } /* ---------------------------------------------------------------------- */ template Eigen::Map ProjectionFiniteStrainFast:: get_operator() { return this->xiField.dyn_eigen(); } template class ProjectionFiniteStrainFast; template class ProjectionFiniteStrainFast; } // muSpectre diff --git a/src/fft/projection_small_strain.cc b/src/fft/projection_small_strain.cc index 91303fd..f045ce6 100644 --- a/src/fft/projection_small_strain.cc +++ b/src/fft/projection_small_strain.cc @@ -1,76 +1,83 @@ /** * @file projection_small_strain.cc * * @author Till Junge * * @date 14 Jan 2018 * * @brief Implementation for ProjectionSmallStrain * * 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 "fft/projection_small_strain.hh" #include "fft/fft_utils.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionSmallStrain:: ProjectionSmallStrain(FFTEngine_ptr engine) : Parent{std::move(engine), Formulation::small_strain} - {} + { + for (auto res: this->fft_engine->get_resolutions()) { + if (res % 2 == 0) { + throw ProjectionError + ("Only an odd number of gridpoints in each direction is supported"); + } + } + } /* ---------------------------------------------------------------------- */ template void ProjectionSmallStrain::initialise(FFT_PlanFlags flags) { Parent::initialise(flags); FFT_freqs fft_freqs(this->fft_engine->get_domain_resolutions(), this->fft_engine->get_lengths()); for (auto && tup: akantu::zip(*this->fft_engine, this->Ghat)) { const auto & ccoord = std::get<0> (tup); auto & G = std::get<1>(tup); auto xi = fft_freqs.get_unit_xi(ccoord); auto kron = [](const Dim_t i, const Dim_t j) -> Real{ return (i==j) ? 1. : 0.; }; for (Dim_t i{0}; i < DimS; ++i) { for (Dim_t j{0}; j < DimS; ++j) { for (Dim_t l{0}; l < DimS; ++l) { for (Dim_t m{0}; m < DimS; ++m ) { Real & g = get(G, i, j, l, m); g = 0.5* (xi(i) * kron(j, l) * xi(m) + xi(i) * kron(j, m) * xi(l) + xi(j) * kron(i, l) * xi(m) + xi(j) * kron(i, m) * xi(l)) - xi(i)*xi(j)*xi(l)*xi(m); } } } } } if (this->get_locations() == Ccoord{}) { this->Ghat[0].setZero(); } } template class ProjectionSmallStrain; template class ProjectionSmallStrain; } // muSpectre diff --git a/tests/test_projection_finite.cc b/tests/test_projection_finite.cc index 5bd227e..60700af 100644 --- a/tests/test_projection_finite.cc +++ b/tests/test_projection_finite.cc @@ -1,128 +1,138 @@ /** * @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 "test_projection.hh" #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(projection_finite_strain); /* ---------------------------------------------------------------------- */ using fixlist = boost::mpl::list< ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrainFast>, ProjectionFixture, ProjectionFiniteStrainFast>, ProjectionFixture, ProjectionFiniteStrainFast>, ProjectionFixture, ProjectionFiniteStrainFast>>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); } + /* ---------------------------------------------------------------------- */ + BOOST_AUTO_TEST_CASE(even_grid_test) { + using Engine = FFTWEngine; + using proj = ProjectionFiniteStrainFast; + auto engine = std::make_unique(Ccoord_t{2, 2}, + Rcoord_t{4.3, 4.3}); + BOOST_CHECK_THROW(proj (std::move(engine)), + std::runtime_error); + } + /* ---------------------------------------------------------------------- */ 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()); FFT_freqs freqs{fix::projector.get_domain_resolutions(), fix::projector.get_lengths()}; 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_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); Vector vec = CcoordOps::get_vector(ccoord, fix::projector.get_lengths()/ fix::projector.get_resolutions()); Real error = (g-v).norm(); BOOST_CHECK_LT(error, tol); if (error >=tol) { 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