diff --git a/CMakeLists.txt b/CMakeLists.txt index 3848e4d..96b48d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,186 +1,176 @@ # ============================================================================= # 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() -<<<<<<< HEAD # 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) -======= -#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) ->>>>>>> feat/store_Lame_consts_per_pixle_T1920 add_definitions(-DNO_EXPERIMENTAL) endif() endif() -<<<<<<< HEAD - -======= ->>>>>>> feat/store_Lame_consts_per_pixle_T1920 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) include_directories(SYSTEM ${PFFT_INCLUDES}) 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_suite1 ${MPIEXEC_EXECUTABLE} ${MPIEXEC_PREFLAGS} ${MPIEXEC_NUMPROC_FLAG} 1 mpi_main_test_suite --report_level=detailed --build_info=TRUE) - add_test(mpi_main_test_suite2 ${MPIEXEC_EXECUTABLE} ${MPIEXEC_PREFLAGS} ${MPIEXEC_NUMPROC_FLAG} 2 mpi_main_test_suite --report_level=detailed --build_info=TRUE) + add_test(mpi_main_test_suite1 ${MPIEXEC} ${MPIEXEC_PREFLAGS} ${MPIEXEC_NUMPROC_FLAG} 1 mpi_main_test_suite --report_level=detailed --build_info=TRUE) + add_test(mpi_main_test_suite2 ${MPIEXEC} ${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/language_bindings/python/muSpectre/__init__.py b/language_bindings/python/muSpectre/__init__.py index 175f148..4c94d3c 100644 --- a/language_bindings/python/muSpectre/__init__.py +++ b/language_bindings/python/muSpectre/__init__.py @@ -1,87 +1,93 @@ # # @file __init__.py # # @author Lars Pastewka # # @date 21 Mar 2018 # # @brief Main entry point for muSpectre Python module # # 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. # -import mpi4py +try: + from mpi4py import MPI +except ImportError: + MPI = None import _muSpectre from _muSpectre import (Formulation, get_domain_ccoord, get_domain_index, get_hermitian_sizes, material, solvers) import muSpectre.fft _factories = {'fftw': ('CellFactory', False), 'fftwmpi': ('FFTWMPICellFactory', True), 'pfft': ('PFFTCellFactory', True), 'p3dfft': ('P3DFFTCellFactory', True)} def Cell(resolutions, lengths, formulation=Formulation.finite_strain, fft='fftw', communicator=None): """ Instantiate a muSpectre Cell class. Parameters ---------- resolutions: list Grid resolutions in the Cartesian directions. lengths: list Physical size of the cell in the Cartesian directions. formulation: Formulation Formulation for strains and stresses used by the solver. Options are `Formulation.finite_strain` and `Formulation.small_strain`. Finite strain formulation is the default. fft: string FFT engine to use. Options are 'fftw', 'fftwmpi', 'pfft' and 'p3dfft'. Default is 'fftw'. communicator: mpi4py communicator mpi4py communicator object passed to parallel FFT engines. Note that the default 'fftw' engine does not support parallel execution. Returns ------- cell: object Return a muSpectre Cell object. """ try: factory_name, is_parallel = _factories[fft] except KeyError: raise KeyError("Unknown FFT engine '{}'.".format(fft)) try: factory = _muSpectre.__dict__[factory_name] except KeyError: raise KeyError("FFT engine '{}' has not been compiled into the " "muSpectre library.".format(fft)) if is_parallel: + if MPI is None: + raise RuntimeError('Parallel solver requested but mpi4py could' + ' not be imported.') if communicator is None: - communicator = mpi4py.MPI.COMM_SELF + communicator = MPI.COMM_SELF return factory(resolutions, lengths, formulation, - mpi4py.MPI._handleof(communicator)) + MPI._handleof(communicator)) else: if communicator is not None: raise ValueError("FFT engine '{}' does not support parallel " "execution.".format(fft)) return factory(resolutions, lengths, formulation) diff --git a/language_bindings/python/muSpectre/fft.py b/language_bindings/python/muSpectre/fft.py index 09e5a79..3fdc5f1 100644 --- a/language_bindings/python/muSpectre/fft.py +++ b/language_bindings/python/muSpectre/fft.py @@ -1,89 +1,95 @@ # # @file fft.py # # @author Lars Pastewka # # @date 27 Mar 2018 # # @brief Wrapper for muSpectre's FFT engines # # 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. # -import mpi4py +try: + from mpi4py import MPI +except ImportError: + MPI = None import _muSpectre _factories = {'fftw': ('FFTW_2d', 'FFTW_3d', False), 'fftwmpi': ('FFTWMPI_2d', 'FFTWMPI_3d', True), 'pfft': ('PFFT_2d', 'PFFT_3d', True), 'p3dfft': ('P3DFFT_2d', 'P3DFFT_3d', True)} def FFT(resolutions, lengths, fft='fftw', communicator=None): """ Instantiate a muSpectre FFT class. Parameters ---------- resolutions: list Grid resolutions in the Cartesian directions. lengths: list Physical size of the cell in the Cartesian directions. fft: string FFT engine to use. Options are 'fftw', 'fftwmpi', 'pfft' and 'p3dfft'. Default is 'fftw'. communicator: mpi4py communicator mpi4py communicator object passed to parallel FFT engines. Note that the default 'fftw' engine does not support parallel execution. Returns ------- cell: object Return a muSpectre Cell object. """ if len(resolutions) != len(lengths): raise ValueError("'resolutions' and 'lengths' must have identical " "lengths.") try: factory_name_2d, factory_name_3d, is_parallel = _factories[fft] except KeyError: raise KeyError("Unknown FFT engine '{}'.".format(fft)) if len(resolutions) == 2: factory_name = factory_name_2d elif len(resolutions) == 3: factory_name = factory_name_3d else: raise ValueError('{}-d transforms are not supported' .format(len(resolutions))) try: factory = _muSpectre.fft.__dict__[factory_name] except KeyError: raise KeyError("FFT engine '{}' has not been compiled into the " "muSpectre library.".format(fft)) if is_parallel: + if MPI is None: + raise RuntimeError('Parallel solver requested but mpi4py could' + ' not be imported.') if communicator is None: - communicator = mpi4py.MPI.COMM_SELF + communicator = MPI.COMM_SELF return factory(resolutions, lengths, - mpi4py.MPI._handleof(communicator)) + MPI._handleof(communicator)) else: if communicator is not None: raise ValueError("FFT engine '{}' does not support parallel " "execution.".format(fft)) return factory(resolutions, lengths) diff --git a/src/fft/fftw_engine.cc b/src/fft/fftw_engine.cc index f500b0f..0bfabc5 100644 --- a/src/fft/fftw_engine.cc +++ b/src/fft/fftw_engine.cc @@ -1,158 +1,153 @@ /** * @file fftw_engine.cc * * @author Till Junge * * @date 03 Dec 2017 * * @brief implements the fftw 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/fftw_engine.hh" namespace muSpectre { template FFTWEngine::FFTWEngine(Ccoord resolutions, Rcoord lengths, Communicator comm) :Parent{resolutions, lengths, comm} { for (auto && pixel: CcoordOps::Pixels(this->fourier_resolutions)) { this->work_space_container.add_pixel(pixel); } } /* ---------------------------------------------------------------------- */ template void FFTWEngine::initialise(FFT_PlanFlags plan_flags) { - if (this->comm.size() > 1) { - throw std::runtime_error("Attempting to run in parallel, but standard " - "FFTW engine does not support parallel " - "execution."); - } if (this->initialised) { throw std::runtime_error("double initialisation, will leak memory"); } Parent::initialise(plan_flags); const int & rank = DimS; std::array narr; const int * const n = &narr[0]; std::copy(this->resolutions.begin(), this->resolutions.end(), narr.begin()); int howmany = Field_t::nb_components; //temporary buffer for plan size_t alloc_size = CcoordOps::get_size(this->resolutions) *howmany; Real * r_work_space = fftw_alloc_real(alloc_size); Real * in = r_work_space; const int * const inembed = nullptr;//nembed are tricky: they refer to physical layout int istride = howmany; int idist = 1; fftw_complex * out = reinterpret_cast(this->work.data()); const int * const onembed = nullptr; int ostride = howmany; int odist = idist; unsigned int flags; switch (plan_flags) { case FFT_PlanFlags::estimate: { flags = FFTW_ESTIMATE; break; } case FFT_PlanFlags::measure: { flags = FFTW_MEASURE; break; } case FFT_PlanFlags::patient: { flags = FFTW_PATIENT; break; } default: throw std::runtime_error("unknown planner flag type"); break; } this->plan_fft = fftw_plan_many_dft_r2c(rank, n, howmany, in, inembed, istride, idist, out, onembed, ostride, odist, FFTW_PRESERVE_INPUT | flags); if (this->plan_fft == nullptr) { throw std::runtime_error("Plan failed"); } fftw_complex * i_in = reinterpret_cast(this->work.data()); Real * i_out = r_work_space; this->plan_ifft = fftw_plan_many_dft_c2r(rank, n, howmany, i_in, inembed, istride, idist, i_out, onembed, ostride, odist, flags); if (this->plan_ifft == nullptr) { throw std::runtime_error("Plan failed"); } fftw_free(r_work_space); this->initialised = true; } /* ---------------------------------------------------------------------- */ template FFTWEngine::~FFTWEngine() noexcept { fftw_destroy_plan(this->plan_fft); fftw_destroy_plan(this->plan_ifft); // TODO: We cannot run fftw_cleanup since subsequent FFTW calls will fail // but multiple FFT engines can be active at the same time. //fftw_cleanup(); } /* ---------------------------------------------------------------------- */ template typename FFTWEngine::Workspace_t & FFTWEngine::fft (Field_t & field) { if (this->plan_fft == nullptr) { throw std::runtime_error("fft plan not initialised"); } if (field.size() != CcoordOps::get_size(this->resolutions)) { throw std::runtime_error("size mismatch"); } fftw_execute_dft_r2c(this->plan_fft, field.data(), reinterpret_cast(this->work.data())); return this->work; } /* ---------------------------------------------------------------------- */ template void FFTWEngine::ifft (Field_t & field) const { if (this->plan_ifft == nullptr) { throw std::runtime_error("ifft plan not initialised"); } if (field.size() != CcoordOps::get_size(this->resolutions)) { throw std::runtime_error("size mismatch"); } fftw_execute_dft_c2r(this->plan_ifft, reinterpret_cast(this->work.data()), field.data()); } template class FFTWEngine; template class FFTWEngine; template class FFTWEngine; } // muSpectre diff --git a/tests/mpi_test_fft_engine.cc b/tests/mpi_test_fft_engine.cc index ea0a1b2..21d1f56 100644 --- a/tests/mpi_test_fft_engine.cc +++ b/tests/mpi_test_fft_engine.cc @@ -1,178 +1,174 @@ /** * @file mpi_test_fft_engine.cc * * @author Lars Pastewka * * @date 06 Mar 2017 * * @brief tests for MPI-parallel fft engine implementations * * 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. */ #define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS #define BOOST_MPL_LIMIT_LIST_SIZE 50 #include #include "tests.hh" #include "mpi_context.hh" #include "fft/fftw_engine.hh" #ifdef WITH_FFTWMPI #include "fft/fftwmpi_engine.hh" #endif #ifdef WITH_PFFT #include "fft/pfft_engine.hh" #endif #include "common/ccoord_operations.hh" #include "common/field_collection.hh" #include "common/field_map.hh" #include "common/iterators.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(mpi_fft_engine); /* ---------------------------------------------------------------------- */ template struct FFTW_fixture { constexpr static Dim_t box_resolution{resolution}; constexpr static Dim_t serial_engine{serial}; constexpr static Real box_length{4.5}; constexpr static Dim_t sdim{Engine::sdim}; constexpr static Dim_t mdim{Engine::mdim}; constexpr static Ccoord_t res() { return CcoordOps::get_cube(box_resolution); } FFTW_fixture(): engine(res(), CcoordOps::get_cube(box_length), MPIContext::get_context().comm) {} Engine engine; }; template struct FFTW_fixture_python_segfault{ constexpr static Dim_t serial_engine{false}; constexpr static Dim_t dim{twoD}; constexpr static Dim_t sdim{twoD}; constexpr static Dim_t mdim{twoD}; constexpr static Ccoord_t res() {return {6, 4};} FFTW_fixture_python_segfault(): engine{res(), {3., 3}, MPIContext::get_context().comm} {} Engine engine; }; using fixlist = boost::mpl::list< #ifdef WITH_FFTWMPI FFTW_fixture, 3>, FFTW_fixture, 3>, FFTW_fixture, 3>, FFTW_fixture, 4>, FFTW_fixture, 4>, FFTW_fixture, 4>, FFTW_fixture_python_segfault>, #endif #ifdef WITH_PFFT FFTW_fixture, 3>, FFTW_fixture, 3>, FFTW_fixture, 3>, FFTW_fixture, 4>, FFTW_fixture, 4>, FFTW_fixture, 4>, FFTW_fixture_python_segfault>, #endif FFTW_fixture, 3, true>>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Constructor_test, Fix, fixlist, Fix) { Communicator &comm = MPIContext::get_context().comm; if (Fix::serial_engine && comm.size() > 1) { - BOOST_CHECK_THROW(Fix::engine.initialise(FFT_PlanFlags::estimate), - std::runtime_error); return; } else { BOOST_CHECK_NO_THROW(Fix::engine.initialise(FFT_PlanFlags::estimate)); } BOOST_CHECK_EQUAL(comm.sum(Fix::engine.size()), CcoordOps::get_size(Fix::res())); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(fft_test, Fix, fixlist, Fix) { if (Fix::serial_engine && Fix::engine.get_communicator().size() > 1) { - // expect initialization failure - BOOST_CHECK_THROW(Fix::engine.initialise(FFT_PlanFlags::estimate), - std::runtime_error); + // dont test serial engies in parallel return; } else { Fix::engine.initialise(FFT_PlanFlags::estimate); } constexpr Dim_t order{2}; using FC_t = GlobalFieldCollection; FC_t fc; auto & input{make_field>("input", fc)}; auto & ref {make_field>("reference", fc)}; auto & result{make_field>("result", fc)}; fc.initialise(Fix::engine.get_resolutions(), Fix::engine.get_locations()); using map_t = MatrixFieldMap; map_t inmap{input}; auto refmap{map_t{ref}}; auto resultmap{map_t{result}}; size_t cntr{0}; for (auto tup: akantu::zip(inmap, refmap)) { cntr++; auto & in_{std::get<0>(tup)}; auto & ref_{std::get<1>(tup)}; in_.setRandom(); ref_ = in_; } auto & complex_field = Fix::engine.fft(input); using cmap_t = MatrixFieldMap, Complex, Fix::mdim, Fix::mdim>; cmap_t complex_map(complex_field); if (Fix::engine.get_locations() == CcoordOps::get_cube(0)) { // Check that 0,0 location has no imaginary part. Real error = complex_map[0].imag().norm(); BOOST_CHECK_LT(error, tol); } /* make sure, the engine has not modified input (which is unfortunately const-casted internally, hence this test) */ for (auto && tup: akantu::zip(inmap, refmap)) { Real error{(std::get<0>(tup) - std::get<1>(tup)).norm()}; BOOST_CHECK_LT(error, tol); } /* make sure that the ifft of fft returns the original*/ Fix::engine.ifft(result); for (auto && tup: akantu::zip(resultmap, refmap)) { Real error{(std::get<0>(tup)*Fix::engine.normalisation() - std::get<1>(tup)).norm()}; BOOST_CHECK_LT(error, tol); if (error > tol) { std::cout << std::get<0>(tup).array()/std::get<1>(tup).array() << std::endl << std::endl; } } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/mpi_test_projection_finite.cc b/tests/mpi_test_projection_finite.cc index b20016c..5811186 100644 --- a/tests/mpi_test_projection_finite.cc +++ b/tests/mpi_test_projection_finite.cc @@ -1,189 +1,185 @@ /** * @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. */ #define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS #define BOOST_MPL_LIMIT_LIST_SIZE 50 #include "fft/projection_finite_strain.hh" #include "fft/projection_finite_strain_fast.hh" #include "fft/fft_utils.hh" #include "mpi_test_projection.hh" #include "fft/fftw_engine.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, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrain, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrain, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrain, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrainFast, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrainFast, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrainFast, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrainFast, PFFTEngine>, #endif ProjectionFixture, ProjectionFiniteStrain, FFTWEngine, false> >; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { if (fix::is_parallel || fix::projector.get_communicator().size() == 1) { BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); } - else { - BOOST_CHECK_THROW(fix::projector.initialise(FFT_PlanFlags::estimate), - std::runtime_error); - } } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test, fix, fixlist, fix) { if (!fix::is_parallel || fix::projector.get_communicator().size() > 1) { return; } // 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 diff --git a/tests/mpi_test_projection_small.cc b/tests/mpi_test_projection_small.cc index cc078b4..b745f5e 100644 --- a/tests/mpi_test_projection_small.cc +++ b/tests/mpi_test_projection_small.cc @@ -1,173 +1,170 @@ /** * @file test_projection_small.cc * * @author Till Junge * * @date 16 Jan 2018 * * @brief tests for standard small strain projection operator * * 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. */ #define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS #define BOOST_MPL_LIMIT_LIST_SIZE 50 #include "fft/projection_small_strain.hh" #include "mpi_test_projection.hh" #include "fft/fft_utils.hh" #include "fft/fftw_engine.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_small_strain); using fixlist = boost::mpl::list< #ifdef WITH_FFTWMPI ProjectionFixture, ProjectionSmallStrain, FFTWMPIEngine>, ProjectionFixture, ProjectionSmallStrain, FFTWMPIEngine>, ProjectionFixture, ProjectionSmallStrain, FFTWMPIEngine>, ProjectionFixture, ProjectionSmallStrain, FFTWMPIEngine>, #endif #ifdef WITH_PFFT ProjectionFixture, ProjectionSmallStrain, PFFTEngine>, ProjectionFixture, ProjectionSmallStrain, PFFTEngine>, ProjectionFixture, ProjectionSmallStrain, PFFTEngine>, ProjectionFixture, ProjectionSmallStrain, PFFTEngine>, #endif ProjectionFixture, ProjectionSmallStrain, FFTWEngine, false> >; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { if (fix::is_parallel || fix::projector.get_communicator().size() == 1) { BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); - } else { - BOOST_CHECK_THROW(fix::projector.initialise(FFT_PlanFlags::estimate), - std::runtime_error); } } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test, fix, fixlist, fix) { if (!fix::is_parallel || fix::projector.get_communicator().size() > 1) { return; } // 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("strain", 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_resolutions()); g.row(0) << k.transpose() * cos(k.dot(vec)); // We need to add I to the term, because this field has a net // zero gradient, which leads to a net -I strain g = 0.5*((g-g.Identity()).transpose() + (g-g.Identity())).eval()+g.Identity(); v = g; } fix::projector.initialise(FFT_PlanFlags::estimate); fix::projector.apply_projection(f_var); constexpr bool verbose{false}; 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) || verbose) { 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; std::cout << "means:" << std::endl << ":" << std::endl << grad.mean() << std::endl << ":" << std::endl << var.mean(); } } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/python_material_linear_elastic4_test.py b/tests/python_material_linear_elastic4_test.py index d0a660b..70a2e4e 100644 --- a/tests/python_material_linear_elastic4_test.py +++ b/tests/python_material_linear_elastic4_test.py @@ -1,77 +1,77 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ @file python_material_linear_elastic4_test.py @author Richard Leute @date 27 Mar 2018 @brief description @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. """ import unittest import numpy as np from python_test_imports import µ class MaterialLinearElastic4_Check(unittest.TestCase): """ Check the implementation of storing the first and second Lame constant in each cell. Assign the same Youngs modulus and Poisson ratio to each cell, from which the two Lame constants are internally computed. Then calculate the stress and compare the result with stress=2*mu*Del0 (Hooke law for small symmetric strains). """ def setUp(self): self.resolution = [7,7] self.lengths = [2.3, 3.9] self.formulation = µ.Formulation.small_strain - self.sys = µ.CellFactory(self.resolution, - self.lengths, - self.formulation) + self.sys = µ.Cell(self.resolution, + self.lengths, + self.formulation) self.mat = µ.material.MaterialLinearElastic4_2d.make( self.sys, "material") def test_solver(self): Youngs_modulus = 10. Poisson_ratio = 0.3 for i, pixel in enumerate(self.sys): self.mat.add_pixel(pixel, Youngs_modulus, Poisson_ratio) self.sys.initialise() tol = 1e-6 Del0 = np.array([[0, 0.025], [0.025, 0]]) maxiter = 100 verbose = 1 solver=µ.solvers.SolverCG(self.sys, tol, maxiter, verbose) r = µ.solvers.newton_cg(self.sys, Del0, solver, tol, tol, verbose) #compare the computed stress with the trivial by hand computed stress mu = (Youngs_modulus/(2*(1+Poisson_ratio))) stress = 2*mu*Del0 self.assertLess(np.linalg.norm(r.stress-stress.reshape(-1,1)), 1e-8)