diff --git a/CMakeLists.txt b/CMakeLists.txt index 27e22f8..ab2e834 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,217 +1,219 @@ # ============================================================================= # 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.0.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) set(MUSPECTRE_PYTHON_MAJOR_VERSION 3) add_compile_options(-Wall -Wextra -Weffc++) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) set(MAKE_DOC_TARGET "OFF" CACHE BOOL "If on, a target dev_doc (which builds the documentation) is added") set(MAKE_TESTS "ON" CACHE BOOL "If on, several ctest targets will be built automatically") set(MAKE_EXAMPLES "ON" CACHE BOOL "If on, the executables in the bin folder will be compiled") set(MAKE_BENCHMARKS "ON" CACHE BOOL "If on, the benchmarks will be compiled") set(MPI_PARALLEL "OFF" CACHE BOOL "If on, MPI-parallel solvers become available") set(RUNNING_IN_CI "OFF" CACHE INTERNAL "changes output format for tests") if(${MAKE_TESTS}) enable_testing() find_package(Boost COMPONENTS unit_test_framework REQUIRED) endif(${MAKE_TESTS}) if(${MPI_PARALLEL}) add_definitions(-DWITH_MPI) find_package(MPI) if (NOT ${MPI_FOUND}) message(SEND_ERROR "You chose MPI but CMake cannot find the MPI package") endif(NOT ${MPI_FOUND}) endif(${MPI_PARALLEL}) 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) add_compile_options(-march=native) string( TOLOWER "${CMAKE_BUILD_TYPE}" build_type ) if (("relwithdebinfo" STREQUAL "${build_type}") OR ("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() add_external_package(Eigen3 VERSION 3.3.0 CONFIG) add_external_package(pybind11 VERSION 2.2 CONFIG) include_directories( ${CMAKE_SOURCE_DIR}/src ${CMAKE_SOURCE_DIR} ) if(APPLE) include_directories(${CMAKE_INSTALL_PREFIX}/include ${Boost_INCLUDE_DIRS}) endif() #build tests (these are before we add -Werror to the compile options) if (${MAKE_TESTS}) ############################################################################## # 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) + muSpectre_add_test(main_test_suite TYPE BOOST main_test_suite --report_level=detailed) + add_executable(cp_test tests/main_test_suite tests/test_material_crystal_plasticity_finite.cc) target_link_libraries(cp_test ${Boost_LIBRARIES} muSpectre) + muSpectre_add_test(cp_test TYPE BOOST cp_test --report_level=detailed) - muSpectre_add_test(main_test_suite TYPE BOOST main_test_suite --report_level=detailed) ############################################################################## # copy python test 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}) find_package(PythonInterp ${MUSPECTRE_PYTHON_MAJOR_VERSION} REQUIRED) muSpectre_add_test(python_binding_test TYPE PYTHON python_binding_tests.py) if(${MPI_PARALLEL}) ############################################################################ # add MPI tests 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) muSpectre_add_test(mpi_main_test_suite1 TYPE BOOST MPI_NB_PROCS 1 mpi_main_test_suite --report_level=detailed) muSpectre_add_test(mpi_main_test_suite2 TYPE BOOST MPI_NB_PROCS 2 mpi_main_test_suite --report_level=detailed) muSpectre_add_test(python_mpi_binding_test1 TYPE PYTHON MPI_NB_PROCS 1 python_mpi_binding_tests.py) muSpectre_add_test(python_mpi_binding_test2 TYPE PYTHON MPI_NB_PROCS 2 python_mpi_binding_tests.py) endif(${MPI_PARALLEL}) endif(${MAKE_TESTS}) ################################################################################ # compile the library add_compile_options( -Werror) add_subdirectory( ${CMAKE_SOURCE_DIR}/src/ ) add_subdirectory( ${CMAKE_SOURCE_DIR}/language_bindings/ ) if (${MAKE_DOC_TARGET}) add_subdirectory( ${CMAKE_SOURCE_DIR}/doc/ ) endif() ################################################################################ if (${MAKE_EXAMPLES}) #compile executables set(binaries ${CMAKE_SOURCE_DIR}/bin/demonstrator1.cc ${CMAKE_SOURCE_DIR}/bin/demonstrator_dynamic_solve.cc ${CMAKE_SOURCE_DIR}/bin/demonstrator2.cc ${CMAKE_SOURCE_DIR}/bin/hyper-elasticity.cc ${CMAKE_SOURCE_DIR}/bin/small_case.cc) if (${MPI_PARALLEL}) set (binaries ${binaries} ${CMAKE_SOURCE_DIR}/bin/demonstrator_mpi.cc ) endif (${MPI_PARALLEL}) 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}) endif (${MAKE_EXAMPLES}) ################################################################################ # compile benchmarks if(${MAKE_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}) endif(${MAKE_BENCHMARKS}) diff --git a/tests/test_material_crystal_plasticity_finite.cc b/tests/test_material_crystal_plasticity_finite.cc index 050b6cd..862926a 100644 --- a/tests/test_material_crystal_plasticity_finite.cc +++ b/tests/test_material_crystal_plasticity_finite.cc @@ -1,422 +1,445 @@ /** * file test_material_crystal_plasticity_finite.cc * * @author Till Junge * * @date 22 May 2018 * * @brief Tests for the basic crystal plasticity material, * MaterialCrystalPlasticityFinite * * @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 "tests.hh" #include "materials/material_crystal_plasticity_finite.hh" #include "materials/material_linear_elastic1.hh" #include "materials/materials_toolbox.hh" #include "common/tensor_algebra.hh" #include "solver/solvers.hh" #include "solver/solver_cg.hh" #include "cell/cell_factory.hh" #include "common/iterators.hh" #include #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(crystal_plasticity); struct CPDefaultParams { Real bulk_m{175e9}; Real shear_m{120e9}; Real gamma_dot0{10e-2}; Real m_par{.1}; Real tau_y0{200e6}; Real h0{0e9}; Real delta_tau_y{100e6}; Real a_par{0}; Real q_n{1.4}; - Real delta_t{1e-4}; + Real delta_t{1e-3}; }; struct CPDefaultParamsElast: public CPDefaultParams { CPDefaultParamsElast(): CPDefaultParams{} { this->tau_y0 = 1e200; } }; template struct CrystalPlastFixture: public CPDefaultParams { using Mat_t = Material; constexpr static Dim_t get_DimS() {return Mat_t::sdim();} constexpr static Dim_t get_DimM() {return Mat_t::mdim();} constexpr static Dim_t get_NbSlip() {return Mat_t::get_NbSlip();} CrystalPlastFixture(): CPDefaultParams{}, mat(name, bulk_m, shear_m, gamma_dot0, m_par, tau_y0, h0, delta_tau_y, a_par, q_n, slip0, normals, delta_t) {} std::string name{"mateerial"}; typename Mat_t::SlipVecs slip0; typename Mat_t::SlipVecs normals; Mat_t mat; }; /** * material with ridiculously high yield strength to test elastic sanity */ template struct Crystal_non_plastFixture: public CPDefaultParamsElast { using Mat_t = Material; constexpr static Dim_t get_DimS() {return Mat_t::sdim();} constexpr static Dim_t get_DimM() {return Mat_t::mdim();} constexpr static Dim_t get_NbSlip() {return Mat_t::get_NbSlip();} Crystal_non_plastFixture(): CPDefaultParamsElast{}, slip0{Mat_t::SlipVecs::Zero()}, normals{Mat_t::SlipVecs::Zero()}, mat(name, bulk_m, shear_m, gamma_dot0, m_par, tau_y0, h0, delta_tau_y, a_par, q_n, slip0, normals, delta_t) { this->slip0.col(0).setConstant(1.); this->normals.col(1).setConstant(1.); } std::string name{"material"}; typename Mat_t::SlipVecs slip0; typename Mat_t::SlipVecs normals; Mat_t mat; }; using mat_list = boost::mpl::list< CrystalPlastFixture>, CrystalPlastFixture>, CrystalPlastFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, Fix, mat_list, Fix) { } using elast_mat_list = boost::mpl::list< Crystal_non_plastFixture>, Crystal_non_plastFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(elastic_test_no_loop, Fix, elast_mat_list, Fix) { constexpr Dim_t Dim{Fix::get_DimS()}; using Ccoord = Ccoord_t; using Euler_t = Eigen::Array; using T2_t = Eigen::Matrix; using T4_t = T4Mat; using Hooke = typename MatTB::Hooke; auto & mat{Fix::mat}; Euler_t angles = 2*pi*Euler_t::Random(); constexpr Ccoord pixel{0}; using MatElast_t = MaterialLinearElastic1; Real Young{MatTB::convert_elastic_modulus (Fix::bulk_m, Fix::shear_m)}; Real Poisson{MatTB::convert_elastic_modulus (Fix::bulk_m, Fix::shear_m)}; MatElast_t mat_ref("hard", Young, Poisson); mat.add_pixel(pixel, angles); T2_t F{T2_t::Identity()}; F(0,1) += .01; T2_t E{.5*(F.transpose()*F - T2_t::Identity())}; const Real lambda{MatTB::convert_elastic_modulus (Fix::shear_m, Fix::bulk_m)}; T2_t stress_ref{Hooke::evaluate_stress(lambda, Fix::shear_m, E)}; // for ∂E/∂F, see Curnier T4_t C{Hooke::compute_C_T4(lambda, Fix::shear_m)}; T4_t tangent_ref{ Matrices::ddot(C, (Matrices::outer_under(F.transpose(), T2_t::Identity())))}; auto & internals{mat.get_internals()}; mat.initialise(); auto & Fp_map = std::get<0>(internals); auto & gamma_dot_map = std::get<1>(internals); auto & tau_y_map = std::get<2>(internals); auto & Euler_map = std::get<3>(internals); auto & dummy_gamma_dot_map = std::get<4>(internals); auto & dummy_tau_inc_map = std::get<5>(internals); T2_t stress = mat.evaluate_stress(F, *Fp_map.begin(), *gamma_dot_map.begin(), *tau_y_map.begin(), *Euler_map.begin(), *dummy_gamma_dot_map.begin(), *dummy_tau_inc_map.begin()); Real error{(stress-stress_ref).norm()/stress_ref.norm()}; BOOST_CHECK_LT(error, tol); if (not(error < tol)) { std::cout << "stress_ref =" << std::endl << stress_ref << std::endl; std::cout << "stress =" << std::endl << stress << std::endl; } auto stress_tgt = mat.evaluate_stress_tangent(F, *Fp_map.begin(), *gamma_dot_map.begin(), *tau_y_map.begin(), *Euler_map.begin(), *dummy_gamma_dot_map.begin(), *dummy_tau_inc_map.begin()); auto stress_tgt_ref = mat_ref.evaluate_stress_tangent(E); T2_t & stress_2{get<0>(stress_tgt)}; T4_t & tangent{get<1>(stress_tgt)}; T2_t stress_2_lin{get<0>(stress_tgt_ref)}; T4_t C_lin{get<1>(stress_tgt_ref)}; T2_t I2{T2_t::Identity()}; T4_t tangent_lin{C_lin*Matrices::outer_under(F.transpose(), I2)}; error = (stress_2-stress_ref).norm()/stress_ref.norm(); BOOST_CHECK_LT(error, tol); if (not (error < tol)) { std::cout << "stress_ref =" << std::endl << stress_ref << std::endl; std::cout << "stress_2 =" << std::endl << stress_2 << std::endl; } error = (stress_2 - stress_2_lin).norm() / stress_2_lin.norm(); BOOST_CHECK_LT(error, tol); if (not (error < tol)) { std::cout << "stress_2 from linear elastic =" << std::endl << stress_2_lin << std::endl; std::cout << "stress_2 =" << std::endl << stress_2 << std::endl; } error = (tangent-tangent_ref).norm()/tangent_ref.norm(); BOOST_CHECK_LT(error, tol); if (not (error < tol)) { std::cout << "tangent_ref =" << std::endl << tangent_ref << std::endl; std::cout << "tangent =" << std::endl << tangent << std::endl; } error = (tangent_lin-tangent_ref).norm()/tangent_ref.norm(); BOOST_CHECK_LT(error, tol); if (not (error < tol)) { std::cout << "tangent_ref =" << std::endl << tangent_ref << std::endl; std::cout << "tangent from linear elastic =" << std::endl << tangent << std::endl; } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(elastic_bimaterial, Fix, elast_mat_list, Fix) { constexpr Dim_t Dim{Fix::get_DimS()}; using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; constexpr Ccoord resolutions{CcoordOps::get_cube(3)}; constexpr Rcoord lengths{CcoordOps::get_cube(1.)}; constexpr Formulation form{Formulation::finite_strain}; // number of layers in the hard material constexpr Uint nb_lays{1}; static_assert(nb_lays < resolutions[0], "the number or layers in the hard material must be smaller " "than the total number of layers in dimension 0"); auto sys{make_cell(resolutions, lengths, form)}; auto sys_ref{make_cell(resolutions, lengths, form)}; using Mat_t = typename Fix::Mat_t; Mat_t & hard = Mat_t::make(sys, "hard", 2*Fix::bulk_m, 2*Fix::shear_m, Fix::gamma_dot0, Fix::m_par, Fix::tau_y0, Fix::h0, Fix::delta_tau_y, Fix::a_par, Fix::q_n, Fix::slip0, Fix::normals, Fix::delta_t); Mat_t & soft = Mat_t::make(sys, "soft", Fix::bulk_m, Fix::shear_m, Fix::gamma_dot0, Fix::m_par, Fix::tau_y0, Fix::h0, Fix::delta_tau_y, Fix::a_par, Fix::q_n, Fix::slip0, Fix::normals, Fix::delta_t); using MatElast_t = MaterialLinearElastic1; Real Young{MatTB::convert_elastic_modulus (Fix::bulk_m, Fix::shear_m)}; Real Poisson{MatTB::convert_elastic_modulus (Fix::bulk_m, Fix::shear_m)}; MatElast_t & hard_ref = MatElast_t::make(sys_ref, "hard", 2*Young, Poisson); MatElast_t & soft_ref = MatElast_t::make(sys_ref, "soft", Young, Poisson); using Euler_t = Eigen::Array; Euler_t angles = Euler_t::Zero(); for (const auto & pixel: sys) { if (pixel[0] < Dim_t(nb_lays)) { hard.add_pixel(pixel, angles); hard_ref.add_pixel(pixel); } else { soft.add_pixel(pixel, angles); soft_ref.add_pixel(pixel); } } using Grad_t = Eigen::MatrixXd; Grad_t F_bar{Grad_t::Identity(Dim, Dim) + Grad_t::Ones(Dim, Dim)*.1}; constexpr Real cg_tol{1e-4}, newton_tol{1e-4}, equil_tol{1e-4}; constexpr Uint maxiter{Dim*10}; constexpr Dim_t verbose{0}; SolverCG solver{sys, cg_tol, maxiter, bool(verbose)}; auto result = newton_cg(sys, F_bar, solver, newton_tol, equil_tol); BOOST_CHECK(result.success); if (not result.success) { std::cout << result.message << std::endl; } SolverCG solver_ref{sys_ref, cg_tol, maxiter, bool(verbose)}; auto result_ref = newton_cg(sys_ref, F_bar, solver_ref, newton_tol, equil_tol); Real error = (result.stress - result_ref.stress).matrix().norm() / result_ref.stress.matrix().norm(); BOOST_CHECK_LT(error, tol); } //----------------------------------------------------------------------------// BOOST_AUTO_TEST_CASE(stress_strain_single_slip_system) { constexpr Dim_t Dim{twoD}; constexpr Dim_t NbSlip{1}; using Ccoord = Ccoord_t; using Mat_t = MaterialCrystalPlasticityFinite; using Euler_t = Eigen::Array; using T2_t = Eigen::Matrix; //using T4_t = T4Mat; using SlipVecs = Eigen::Matrix; CPDefaultParams params{}; SlipVecs normals{}; normals << 0, 1; SlipVecs slip_dirs{}; slip_dirs << 1, 0; - Mat_t mat("material", - params.bulk_m, - params.shear_m, - params.gamma_dot0, - params.m_par, - params.tau_y0, - params.h0, - params.delta_tau_y, - params.a_par, - params.q_n, - slip_dirs, - normals, - params.delta_t); Euler_t angles = Euler_t::Zero(); - mat.add_pixel(Ccoord{}, angles); - - auto & internals{mat.get_internals()}; - - mat.initialise(); - - auto & Fp_map = std::get<0>(internals); - auto & gamma_dot_map = std::get<1>(internals); - auto & tau_y_map = std::get<2>(internals); - auto & Euler_map = std::get<3>(internals); - auto & dummy_gamma_dot_map = std::get<4>(internals); - auto & dummy_tau_inc_map = std::get<5>(internals); + //! check for multiple time steps + constexpr std::array time_steps {1e-2, 1e-3, 1e-4}; + //! imposed strain increments + constexpr Real Delta_gamma{1e-4}; + //! maximum total strain + constexpr Real gamma_max_no_hard{1e-2}; + //! maximum total strain + constexpr Real gamma_max_hard{1e-1}; + + //! loose tolerance to check for saturation values (slowly reached asymptotes + constexpr Real saturation_tol{1e-3}; + + + auto stress_strain_runner = [&](auto h0, auto a_par, auto saturation_stress, + auto gamma_max) { + for (auto & dt: time_steps) { + Mat_t mat("material", + params.bulk_m, + params.shear_m, + Delta_gamma/dt, + params.m_par, + params.tau_y0, + h0, + params.delta_tau_y, + a_par, + params.q_n, + slip_dirs, + normals, + dt); + + mat.add_pixel(Ccoord{}, angles); + + auto & internals{mat.get_internals()}; + + mat.initialise(); + + auto & Fp_map = std::get<0>(internals); + auto & gamma_dot_map = std::get<1>(internals); + auto & tau_y_map = std::get<2>(internals); + auto & Euler_map = std::get<3>(internals); + auto & dummy_gamma_dot_map = std::get<4>(internals); + auto & dummy_tau_inc_map = std::get<5>(internals); + + T2_t F{T2_t::Identity()}; + + T2_t stress; + for (int i{}; i < gamma_max/Delta_gamma; ++i) { + F(0, 1) += Delta_gamma; + + stress = mat.evaluate_stress(F, + *Fp_map.begin(), + *gamma_dot_map.begin(), + *tau_y_map.begin(), + *Euler_map.begin(), + *dummy_gamma_dot_map.begin(), + *dummy_tau_inc_map.begin()); + mat.save_history_variables(); + } + auto error {std::abs(stress(0, 1)-saturation_stress)/params.tau_y0}; + BOOST_CHECK_LT(error, saturation_tol); + if (not(error < saturation_tol)) { + std::cout << "τ_xy = " << stress(0,1) << ", but shoud be " << saturation_stress << std::endl; + } + } + }; - T2_t F{T2_t::Identity()}; - Real shear_incr = 2e-3; - - for (int i{}; i < 25; ++i) { - F(0, 1) += shear_incr; - - T2_t stress = mat.evaluate_stress(F, - *Fp_map.begin(), - *gamma_dot_map.begin(), - *tau_y_map.begin(), - *Euler_map.begin(), - *dummy_gamma_dot_map.begin(), - *dummy_tau_inc_map.begin()); - mat.save_history_variables(); - - //std::cout << "for F =\n" << F <