diff --git a/cmake/muspectreTools.cmake b/cmake/muspectreTools.cmake index 1ab055a..753bcae 100644 --- a/cmake/muspectreTools.cmake +++ b/cmake/muspectreTools.cmake @@ -1,263 +1,262 @@ #============================================================================== # file muspectreTools.cmake # # @author Nicolas Richart # # @date 11 Jan 2018 # # @brief some tool to help to do stuff with cmake in µSpectre # # @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 µSpectre; see the file COPYING. If not, write to the # Free Software Foundation, Inc., 59 Temple Place - Suite 330, # Boston, MA 02111-1307, USA. # #[=[.rst: µSpectreTools ------------- This module provide some helper functions for µSpectre :: download_external_project(project_name URL BACKEND [TAG ] [THIRD_PARTY_SRC_DIR ] [NO_UPDATE] ) Downloads the external project based on the uri ``THIRD_PARTY_SRC_DIR `` Specifies where to download the source ``NO_UPDATE`` Does not try to update existing download :: mark_as_advanced_prefix(prefix) Marks as advanced all variables whoes names start with prefix :: add_external_dep(package [IGNORE_SYSTEM] [VERSION ] [...]) Tries to find the external package and if not found but a local ${package}.cmake files exists, it includes it. The extra arguments are passed to find_package() ``IGNORE_SYSTEM`` does not do the find_package ``VERSION`` Specifies the required minimum version ]=] # function(download_external_project project_name) include(CMakeParseArguments) set(_dep_flags NO_UPDATE) set(_dep_one_variables URL TAG BACKEND THIRD_PARTY_SRC_DIR ) set(_dep_multi_variables) cmake_parse_arguments(_dep_args "${_dep_flags}" "${_dep_one_variables}" "${_dep_multi_variables}" ${ARGN} ) if(NOT _dep_args_URL) message(FATAL_ERROR "You have to provide a URL for the project ${project_name}") endif() if(NOT _dep_args_BACKEND) message(FATAL_ERROR "You have to provide a backend to download ${project_name}") endif() if(_dep_args_TAG) set(_ep_tag "${_dep_args_BACKEND}_TAG ${_dep_args_TAG}") endif() set(_src_dir ${PROJECT_SOURCE_DIR}/third-party/${project_name}) if (_dep_args_THIRD_PARTY_SRC_DIR) set(_src_dir ${_dep_args_THIRD_PARTY_SRC_DIR}/${project_name}) endif() if(EXISTS ${_src_dir}/.DOWNLOAD_SUCCESS AND _dep_args_NO_UPDATE) return() endif() set(_working_dir ${PROJECT_BINARY_DIR}/third-party/${project_name}-download) file(WRITE ${_working_dir}/CMakeLists.txt " cmake_minimum_required(VERSION 3.1) project(${project_name}-download NONE) include(ExternalProject) ExternalProject_Add(${project_name} SOURCE_DIR ${_src_dir} BINARY_DIR ${_working_dir} ${_dep_args_BACKEND}_REPOSITORY ${_dep_args_URL} ${_ep_tag} CONFIGURE_COMMAND \"\" BUILD_COMMAND \"\" INSTALL_COMMAND \"\" TEST_COMMAND \"\" ) ") message(STATUS "Downloading ${project_name} ${_dep_args_GIT_TAG}") execute_process(COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" . RESULT_VARIABLE _result WORKING_DIRECTORY ${_working_dir} OUTPUT_FILE ${_working_dir}/configure-out.log ERROR_FILE ${_working_dir}/configure-error.log) if(_result) message(FATAL_ERROR "Something went wrong (${_result}) during the download" " process of ${project_name} check the file" " ${_working_dir}/configure-error.log for more details:") file(STRINGS "${_working_dir}/configure-error.log" ERROR_MSG) message("${ERROR_MSG}") endif() execute_process(COMMAND "${CMAKE_COMMAND}" --build . RESULT_VARIABLE _result WORKING_DIRECTORY ${_working_dir} OUTPUT_FILE ${_working_dir}/build-out.log ERROR_FILE ${_working_dir}/build-error.log) if(_result) message(FATAL_ERROR "Something went wrong (${_result}) during the download" " process of ${project_name} check the file" " ${_working_dir}/build-error.log for more details") endif() file(WRITE ${_src_dir}/.DOWNLOAD_SUCCESS "") endfunction() # ------------------------------------------------------------------------------ function(mark_as_advanced_prefix prefix) get_property(_list DIRECTORY PROPERTY VARIABLES) foreach(_var ${_list}) if(${_var} MATCHES "^${prefix}.*") mark_as_advanced(${_var}) endif() endforeach() endfunction() # ------------------------------------------------------------------------------ function(add_external_package package) include(CMakeParseArguments) set(_cmake_includes ${PROJECT_SOURCE_DIR}/cmake) set(_${package}_external_dir ${PROJECT_BINARY_DIR}/external) set(_aep_flags IGNORE_SYSTEM ) set(_aep_one_variables VERSION ) set(_aep_multi_variables) cmake_parse_arguments(_aep_args "${_aep_flags}" "${_aep_one_variables}" "${_aep_multi_variables}" ${ARGN} ) if(_aep_args_VERSION) set(_${package}_version ${_aep_args_VERSION}) endif() if(NOT EXISTS ${_cmake_includes}/${package}.cmake) set(_required REQUIRED) endif() if(NOT _aep_args_IGNORE_SYSTEM) find_package(${package} ${_${package}_version} ${_required} ${_aep_UNPARSED_ARGUMENTS} QUIET) if(${package}_FOUND AND NOT ${package}_FOUND_EXTERNAL) return() endif() endif() if(EXISTS ${_cmake_includes}/${package}.cmake) include(${_cmake_includes}/${package}.cmake) endif() endfunction() function(muSpectre_add_test test_name) include(CMakeParseArguments) set(_mat_flags ) set(_mat_one_variables TYPE MPI_NB_PROCS ) set(_mat_multi_variables) cmake_parse_arguments(_mat_args "${_mat_flags}" "${_mat_one_variables}" "${_mat_multi_variables}" ${ARGN} ) if ("${_mat_args_TYPE}" STREQUAL "BOOST") elseif("${_mat_args_TYPE}" STREQUAL "PYTHON") else () message (SEND_ERROR "Can only handle types 'BOOST' and 'PYTHON'") endif ("${_mat_args_TYPE}" STREQUAL "BOOST") set(_exe ${_mat_args_UNPARSED_ARGUMENTS}) if (${RUNNING_IN_CI}) if ("${_mat_args_TYPE}" STREQUAL "BOOST") LIST(APPEND _exe "--logger=JUNIT,all,test_results_${test_name}.xml") - LIST(APPEND cat test_results_${test_name}.xml) elseif("${_mat_args_TYPE}" STREQUAL "PYTHON") set(_exe ${PYTHON_EXECUTABLE} -m pytest --junitxml test_results_${test_name}.xml ${_exe}) endif ("${_mat_args_TYPE}" STREQUAL "BOOST") else () if("${_mat_args_TYPE}" STREQUAL "PYTHON") set(_exe ${PYTHON_EXECUTABLE} ${_exe}) endif ("${_mat_args_TYPE}" STREQUAL "PYTHON") endif (${RUNNING_IN_CI}) if(${_mat_args_MPI_NB_PROCS}) set(_exe ${MPIEXEC_EXECUTABLE} ${MPIEXEC_PREFLAGS} ${MPIEXEC_NUMPROC_FLAG} ${_mat_args_MPI_NB_PROCS} ${_exe}) endif(${_mat_args_MPI_NB_PROCS}) add_test(${test_name} ${_exe}) endfunction() diff --git a/tests/python_exact_reference_elastic_test.py b/tests/python_exact_reference_elastic_test.py index f723e34..f595a54 100644 --- a/tests/python_exact_reference_elastic_test.py +++ b/tests/python_exact_reference_elastic_test.py @@ -1,434 +1,434 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ @file python_exact_reference_test.py @author Till Junge @date 18 Jun 2018 @brief Tests exactness of each iterate with respect to python reference implementation from GooseFFT for elasticity 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 numpy.linalg import norm from python_test_imports import µ import scipy.sparse.linalg as sp import itertools np.set_printoptions(linewidth=180) comparator_nb_cols=9 # ----------------------------------- GRID ------------------------------------ ndim = 3 # number of dimensions N = 3 # number of voxels (assumed equal for all directions) Nx = Ny = Nz = N def deserialise_t4(t4): turnaround = np.arange(ndim**2).reshape(ndim,ndim).T.reshape(-1) retval = np.zeros([ndim*ndim, ndim*ndim]) for i,j in itertools.product(range(ndim**2), repeat=2): retval[i,j] = t4[:ndim, :ndim, :ndim, :ndim, 0,0].reshape(ndim**2, ndim**2)[turnaround[i], turnaround[j]] pass return retval def scalar_to_goose(s_msp): s_goose = np.zeros((Nx, Ny, Nz)) for i in range(Nx): for j in range(Ny): for k in range(Nz): s_goose[i,j,k] = s_msp[Nz*Ny*i + Nz*j + k] pass pass return s_goose def t2_to_goose(t2_msp): t2_goose = np.zeros((ndim, ndim, Nx, Ny, Nz)) for i in range(Nx): for j in range(Ny): for k in range(Nz): t2_goose[:,:,i,j,k] = t2_msp[:, Nz*Ny*i + Nz*j + k].reshape(ndim, ndim).T pass pass return t2_goose def t2_vec_to_goose(t2_msp_vec): return t2_to_goose(t2_msp_vec.reshape(ndim*ndim, Nx*Ny*Nz)).reshape(-1) def scalar_vec_to_goose(s_msp_vec): return scalar_to_goose(s_msp_vec.reshape(Nx*Ny*N)).reshape(-1) def t4_to_goose(t4_msp, right_transposed=True): t4_goose = np.zeros((ndim, ndim, ndim, ndim, Nx, Ny, Nz)) turnaround = np.arange(ndim**2).reshape(ndim,ndim).T.reshape(-1) for i in range(Nx): for j in range(Ny): for k in range(Nz): tmp = t4_msp[:, Nz*Ny*i + Nz*j + k].reshape(ndim**2, ndim**2) goose_view = t4_goose[:,:,:,:,i,j,k].reshape(ndim**2, ndim**2) for a,b in itertools.product(range(ndim**2), repeat=2): a_id = a if right_transposed else turnaround[a] goose_view[a,b] = tmp[a_id, turnaround[b]] pass pass return t4_goose def t4_vec_to_goose(t4_msp_vec): return t4_to_goose(t4_msp_vec.reshape(ndim**4, Nx*Ny*Nz)).reshape(-1) def t2_from_goose(t2_goose): nb_pix = Nx*Ny*Nz t2_msp = np.zeros((ndim**2, nb_pix), order='F') for i in range(Nx): for j in range(Ny): for k in range(Nz): view = t2_msp[:, i + Nx*j + Nx*Ny*k].reshape(ndim, ndim).T view = t2goose[:,:,i,j,k].T pass pass return t2_msp # ---------------------- PROJECTION, TENSORS, OPERATIONS ---------------------- # tensor operations/products: np.einsum enables index notation, avoiding loops # e.g. ddot42 performs $C_ij = A_ijkl B_lk$ for the entire grid trans2 = lambda A2 : np.einsum('ijxyz ->jixyz ',A2 ) ddot42 = lambda A4,B2: np.einsum('ijklxyz,lkxyz ->ijxyz ',A4,B2) ddot44 = lambda A4,B4: np.einsum('ijklxyz,lkmnxyz->ijmnxyz',A4,B4) dot22 = lambda A2,B2: np.einsum('ijxyz ,jkxyz ->ikxyz ',A2,B2) dot24 = lambda A2,B4: np.einsum('ijxyz ,jkmnxyz->ikmnxyz',A2,B4) dot42 = lambda A4,B2: np.einsum('ijklxyz,lmxyz ->ijkmxyz',A4,B2) dyad22 = lambda A2,B2: np.einsum('ijxyz ,klxyz ->ijklxyz',A2,B2) # identity tensor [single tensor] i = np.eye(ndim) # identity tensors [grid of tensors] I = np.einsum('ij,xyz' , i ,np.ones([N,N,N])) I4 = np.einsum('ijkl,xyz->ijklxyz',np.einsum('il,jk',i,i),np.ones([N,N,N])) I4rt = np.einsum('ijkl,xyz->ijklxyz',np.einsum('ik,jl',i,i),np.ones([N,N,N])) I4s = (I4+I4rt)/2. II = dyad22(I,I) # projection operator [grid of tensors] # NB can be vectorized (faster, less readable), see: "elasto-plasticity.py" # - support function / look-up list / zero initialize delta = lambda i,j: np.float(i==j) # Dirac delta function freq = np.arange(-(N-1)/2.,+(N+1)/2.) # coordinate axis -> freq. axis Ghat4 = np.zeros([ndim,ndim,ndim,ndim,N,N,N]) # zero initialize # - compute for i,j,l,m in itertools.product(range(ndim),repeat=4): for x,y,z in itertools.product(range(N), repeat=3): q = np.array([freq[x], freq[y], freq[z]]) # frequency vector if not q.dot(q) == 0: # zero freq. -> mean Ghat4[i,j,l,m,x,y,z] = delta(i,m)*q[j]*q[l]/(q.dot(q)) # (inverse) Fourier transform (for each tensor component in each direction) fft = lambda x : np.fft.fftshift(np.fft.fftn (np.fft.ifftshift(x),[N,N,N])) ifft = lambda x : np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(x),[N,N,N])) # functions for the projection 'G', and the product 'G : K^LT : (delta F)^T' G = lambda A2 : np.real( ifft( ddot42(Ghat4,fft(A2)) ) ).reshape(-1) K_dF = lambda dFm: trans2(ddot42(K4,trans2(dFm.reshape(ndim,ndim,N,N,N)))) G_K_dF = lambda dFm: G(K_dF(dFm)) # ------------------- PROBLEM DEFINITION / CONSTITIVE MODEL ------------------- # phase indicator: cubical inclusion of volume fraction (9**3)/(31**3) phase = np.zeros([N,N,N]); phase[:2,:2,:2] = 1. # material parameters + function to convert to grid of scalars param = lambda M0,M1: M0*np.ones([N,N,N])*(1.-phase)+M1*np.ones([N,N,N])*phase K = param(0.833,8.33) # bulk modulus [grid of scalars] mu = param(0.386,3.86) # shear modulus [grid of scalars] # constitutive model: grid of "F" -> grid of "P", "K4" [grid of tensors] def constitutive(F): C4 = K*II+2.*mu*(I4s-1./3.*II) S = ddot42(C4,.5*(dot22(trans2(F),F)-I)) P = dot22(F,S) K4 = dot24(S,I4)+ddot44(ddot44(I4rt,dot42(dot24(F,C4),trans2(F))),I4rt) return P,K4 F = np.array(I,copy=True) P,K4 = constitutive(F) class Counter(object): def __init__(self): self.count = self.reset() def reset(self): self.count = 0 return self.count def get(self): return self.count def __call__(self, dummy): self.count += 1 class LinearElastic_Check(unittest.TestCase): def t2_comparator(self, µT2, gT2): err_sum = 0. err_max = 0. for counter, (i, j, k) in enumerate(self.rve): print((i,j,k)) µ_arr = µT2[:, counter].reshape(ndim, ndim).T g_arr = gT2[:,:,i,j,k] self.assertEqual(Nz*Ny*i+Nz*j + k, counter) print("µSpectre:") print(µ_arr) print("Goose:") print(g_arr) print(µ_arr-g_arr) err = norm(µ_arr-g_arr) print("error norm for pixel {} = {}".format((i, j, k), err)) err_sum += err err_max = max(err_max, err) pass print("∑(err) = {}, max(err) = {}".format (err_sum, err_max)) return err_sum def t4_comparator(self, µT4, gT4, right_transposed=True): """ right_transposed: in de Geus's notation, e.g., stiffness tensors have the last two dimensions inverted """ err_sum = 0. err_max = 0. errs = dict() turnaround = np.arange(ndim**2).reshape(ndim,ndim).T.reshape(-1) def zero_repr(arr): arrcopy = arr.copy() arrcopy[abs(arr)<1e-13] = 0. return arrcopy for counter, (i, j, k) in enumerate(self.rve): µ_arr_tmp = µT4[:, counter].reshape(ndim**2, ndim**2).T µ_arr = np.empty((ndim**2, ndim**2)) for a,b in itertools.product(range(ndim**2), repeat=2): a = a if right_transposed else turnaround[a] µ_arr[a,b] = µ_arr_tmp[a, turnaround[b]] g_arr = gT4[:,:,:,:,i,j,k].reshape(ndim**2, ndim**2) self.assertEqual(Nz*Ny*i+Nz*j + k, counter) print("µSpectre:") print(zero_repr(µ_arr[:, :comparator_nb_cols])) print("Goose:") print(zero_repr(g_arr[:, :comparator_nb_cols])) print("Diff") print(zero_repr((µ_arr-g_arr)[:, :comparator_nb_cols])) err = norm(µ_arr-g_arr)/norm(g_arr) print("error norm for pixel {} = {}".format((i, j, k), err)) err_sum += err errs[(i,j,k)] = err err_max = max(err_max, err) print("count {:>2}: err_norm = {:.5f}, err_sum = {:.5f}".format( counter, err, err_sum)) pass print("∑(err) = {}, max(err) = {}".format (err_sum, err_max)) return err_sum, errs def setUp(self): #---------------------------- µSpectre init ----------------------------------- resolution = list(phase.shape) dim = len(resolution) self.dim=dim center = np.array([r//2 for r in resolution]) incl = resolution[0]//5 ## Domain dimensions lengths = [float(r) for r in resolution] ## formulation (small_strain or finite_strain) formulation = µ.Formulation.finite_strain ## build a computational domain self.rve = µ.Cell(resolution, lengths, formulation) def get_E_nu(bulk, shear): Young = 9*bulk*shear/(3*bulk + shear) Poisson = Young/(2*shear) - 1 return Young, Poisson mat = µ.material.MaterialLinearElastic1_3d E, nu = get_E_nu(.833, .386) hard = mat.make(self.rve, 'hard', 10*E, nu) soft = mat.make(self.rve, 'soft', E, nu) for pixel in self.rve: if pixel[0] < 2 and pixel[1] < 2 and pixel[2] < 2: hard.add_pixel(pixel) else: soft.add_pixel(pixel) def test_solve(self): before_cg_tol = 1e-11 cg_tol = 1e-11 after_cg_tol = 1e-9 newton_tol = 1e-4 # ----------------------------- NEWTON ITERATIONS --------------------- # initialize deformation gradient, and stress/stiffness [tensor grid] global K4, P, F F = np.array(I,copy=True) F2 = np.array(I,copy=True)*1.1 P2,K42 = constitutive(F2) P,K4 = constitutive(F) self.rve.set_uniform_strain(np.array(np.eye(ndim))) µF = self.rve.get_strain() self.assertLess(norm(t2_vec_to_goose(µF) - F.reshape(-1))/norm(F), before_cg_tol) # set macroscopic loading DbarF = np.zeros([ndim,ndim,N,N,N]); DbarF[0,1] += 1.0 # initial residual: distribute "barF" over grid using "K4" b = -G_K_dF(DbarF) F += DbarF Fn = np.linalg.norm(F) iiter = 0 # µSpectre inits µbarF = np.zeros_like(µF) µbarF[ndim, :] += 1. µF2 = µF.copy()*1.1 µP2, µK2 = self.rve.evaluate_stress_tangent(µF2) err = norm(t2_vec_to_goose(µP2) - P2.reshape(-1))/norm(P2) if not (err < before_cg_tol): self.t2_comparator(µP2, µK2) self.assertLess(err, before_cg_tol) self.rve.set_uniform_strain(np.array(np.eye(ndim))) µP, µK = self.rve.evaluate_stress_tangent(µF) err = norm(t2_vec_to_goose(µP) - P.reshape(-1)) if not (err < before_cg_tol): print(µF) self.t2_comparator(µP, P) self.assertLess(err, before_cg_tol) err = norm(t4_vec_to_goose(µK) - K4.reshape(-1))/norm(K4) if not (err < before_cg_tol): print ("err = {}".format(err)) self.assertLess(err, before_cg_tol) µF += µbarF µFn = norm(µF) self.assertLess(norm(t2_vec_to_goose(µF) - F.reshape(-1))/norm(F), before_cg_tol) µG_K_dF = lambda x: self.rve.directional_stiffness(x.reshape(µF.shape)).reshape(-1) µG = lambda x: self.rve.project(x).reshape(-1) µb = -µG_K_dF(µbarF) err = (norm(t2_vec_to_goose(µb.reshape(µF.shape)) - b) / norm(b)) if not (err < before_cg_tol): print("|µb| = {}".format(norm(µb))) print("|b| = {}".format(norm(b))) print("total error = {}".format(err)) self.t2_comparator(µb.reshape(µF.shape), b.reshape(F.shape)) self.assertLess(err, before_cg_tol) # iterate as long as the iterative update does not vanish while True: # solve linear system using CG g_counter = Counter() dFm,_ = sp.cg(tol=cg_tol, A = sp.LinearOperator(shape=(F.size,F.size), matvec=G_K_dF,dtype='float'), b = b, callback=g_counter ) µ_counter = Counter() µdFm,_ = sp.cg(tol=cg_tol, A = sp.LinearOperator(shape=(F.size,F.size), matvec=µG_K_dF,dtype='float'), b = µb, callback=µ_counter) err = g_counter.get()-µ_counter.get() if err != 0: - print("n_iter(g) = {}, n_iter(µ) = {}".format(g_counter.get(), + print("n_iter(g) = {}, n_iter(mu) = {}".format(g_counter.get(), µ_counter.get())) pass #self.assertEqual(g_counter.get(), µ_counter.get()) # in the last iteration, the increment is essentially # zero, so we don't care about relative error anymore err = norm(t2_vec_to_goose(µdFm) - dFm)/norm(dFm) if norm(dFm)/Fn > newton_tol and norm(µdFm)/Fn > newton_tol: if not (err < after_cg_tol): self.t2_comparator(µdFm.reshape(µF.shape), dFm.reshape(F.shape)) print("µdFm.shape = {}".format(µdFm.shape)) print("|µdFm| = {}".format(norm(µdFm))) print("|dFm| = {}".format(norm(dFm))) print("|µdFm - dFm| = {}".format(norm(µdFm-dFm))) print("AssertionWarning: {} is not less than {}".format(err, after_cg_tol)) self.assertLess(err, after_cg_tol) # update DOFs (array -> tens.grid) F += dFm.reshape(ndim,ndim,N,N,N) µF += µdFm.reshape(µF.shape) # new residual stress and tangent P,K4 = constitutive(F) µP, µK = self.rve.evaluate_stress_tangent(µF) err = norm(t2_vec_to_goose(µP) - P.reshape(-1))/norm(P) self.assertLess(err, before_cg_tol) err = norm(t4_vec_to_goose(µK) - K4.reshape(-1))/norm(K4) if not (err < before_cg_tol): print ("err = {}".format(err)) self.t4_comparator(µK, K4) self.assertLess(err, before_cg_tol) # convert res.stress to residual b = -G(P) µb = -µG(µP) # in the last iteration, the rhs is essentianly zero, # leading to large relative errors, which are ok. So we # want either the relative error for the rhs to be small, # or their absolute error to be small compared to unity diff_norm = norm(t2_vec_to_goose(µb) - b.reshape(-1)) err = diff_norm/norm(b) if not ((err < after_cg_tol) or (diff_norm < before_cg_tol)): self.t2_comparator(µb.reshape(µF.shape), b.reshape(F.shape)) print("|µb| = {}".format(norm(µb))) print("|b| = {}".format(norm(b))) print("err = {}".format(err)) print("|µb-b| = {}".format(norm(t2_vec_to_goose(µb) - b.reshape(-1)))) print("AssertionWarning: {} is not less than {}".format(err, before_cg_tol)) self.assertTrue((err < after_cg_tol) or (diff_norm < after_cg_tol)) # print residual to the screen print('Goose: %10.15e'%(np.linalg.norm(dFm)/Fn)) - print('µSpectre: %10.15e'%(np.linalg.norm(µdFm)/µFn)) + print('muSpectre: %10.15e'%(np.linalg.norm(µdFm)/µFn)) if np.linalg.norm(dFm)/Fn0: break # check convergence iiter += 1 if __name__ == '__main__': unittest.main() diff --git a/tests/python_exact_reference_plastic_test.py b/tests/python_exact_reference_plastic_test.py index 17206b0..1532ad5 100644 --- a/tests/python_exact_reference_plastic_test.py +++ b/tests/python_exact_reference_plastic_test.py @@ -1,690 +1,690 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ file python_exact_reference_plastic_test.py @author Till Junge @date 22 Jun 2018 @brief Tests exactness of each iterate with respect to python reference implementation from GooseFFT for plasticity 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 import numpy.linalg as linalg from python_test_imports import µ # turn of warning for zero division # (which occurs in the linearization of the logarithmic strain) np.seterr(divide='ignore', invalid='ignore') import scipy.sparse.linalg as sp import itertools import sys from python_exact_reference_elastic_test import deserialise_t4, t2_to_goose from python_exact_reference_elastic_test import t2_vec_to_goose from python_exact_reference_elastic_test import scalar_vec_to_goose from python_exact_reference_elastic_test import t4_to_goose from python_exact_reference_elastic_test import t4_vec_to_goose from python_exact_reference_elastic_test import t2_from_goose from python_exact_reference_elastic_test import Counter import python_exact_reference_elastic_test as elastic_ref from material_hyper_elasto_plastic1 import PK1_fun_3d # ----------------------------------- GRID ------------------------------------ from python_exact_reference_elastic_test import ndim, N, Nx, Ny, Nz shape = [Nx, Ny, Nz] standalone_dyad22 = lambda A2,B2: np.einsum('ij ,kl ->ijkl',A2,B2) standalone_dyad11 = lambda A1,B1: np.einsum('i ,j ->ij ',A1,B1) standalone_dot22 = lambda A2,B2: np.einsum('ij ,jk ->ik ',A2,B2) standalone_dot24 = lambda A2,B4: np.einsum('ij ,jkmn->ikmn',A2,B4) standalone_dot42 = lambda A4,B2: np.einsum('ijkl,lm ->ijkm',A4,B2) standalone_inv2 = np.linalg.inv standalone_ddot22 = lambda A2,B2: np.einsum('ij ,ji -> ',A2,B2) standalone_ddot42 = lambda A4,B2: np.einsum('ijkl,lk ->ij ',A4,B2) standalone_ddot44 = lambda A4,B4: np.einsum('ijkl,lkmn->ijmn',A4,B4) def constitutive_standalone(K, mu, H, tauy0, F,F_t,be_t,ep_t, dim): I = np.eye(dim) II = standalone_dyad22(I,I) I4 = np.einsum('il,jk',I,I) I4rt = np.einsum('ik,jl',I,I) I4s = (I4+I4rt)/2. def ln2(A2): vals,vecs = np.linalg.eig(A2) return sum( [np.log(vals[i])*standalone_dyad11(vecs[:,i],vecs[:,i]) for i in range(dim)]) def exp2(A2): vals,vecs = np.linalg.eig(A2) return sum( [np.exp(vals[i])*standalone_dyad11(vecs[:,i],vecs[:,i]) for i in range(dim)]) # function to compute linearization of the logarithmic Finger tensor def dln2_d2(A2): vals,vecs = np.linalg.eig(A2) K4 = np.zeros([dim, dim, dim, dim]) for m, n in itertools.product(range(dim),repeat=2): if vals[n]==vals[m]: gc = (1.0/vals[m]) else: gc = (np.log(vals[n])-np.log(vals[m]))/(vals[n]-vals[m]) K4 += gc*standalone_dyad22(standalone_dyad11(vecs[:,m],vecs[:,n]),standalone_dyad11(vecs[:,m],vecs[:,n])) return K4 # elastic stiffness tensor C4e = K*II+2.*mu*(I4s-1./3.*II) # trial state Fdelta = standalone_dot22(F,standalone_inv2(F_t)) be_s = standalone_dot22(Fdelta,standalone_dot22(be_t,Fdelta.T)) lnbe_s = ln2(be_s) tau_s = standalone_ddot42(C4e,lnbe_s)/2. taum_s = standalone_ddot22(tau_s,I)/3. taud_s = tau_s-taum_s*I taueq_s = np.sqrt(3./2.*standalone_ddot22(taud_s,taud_s)) div = np.where(taueq_s < 1e-12, np.ones_like(taueq_s), taueq_s) N_s = 3./2.*taud_s/div phi_s = taueq_s-(tauy0+H*ep_t) phi_s = 1./2.*(phi_s+np.abs(phi_s)) # return map dgamma = phi_s/(H+3.*mu) ep = ep_t + dgamma tau = tau_s -2.*dgamma*N_s*mu lnbe = lnbe_s-2.*dgamma*N_s be = exp2(lnbe) P = standalone_dot22(tau,standalone_inv2(F).T) # consistent tangent operator a0 = dgamma*mu/taueq_s a1 = mu/(H+3.*mu) C4ep = (((K-2./3.*mu)/2.+a0*mu)*II+(1.-3.*a0)*mu* I4s+2.*mu*(a0-a1)*standalone_dyad22(N_s,N_s)) dlnbe4_s = dln2_d2(be_s) dbe4_s = 2.*standalone_dot42(I4s,be_s) #K4a = ((C4e/2.)*(phi_s<=0.).astype(np.float)+ # C4ep*(phi_s>0.).astype(np.float)) K4a = np.where(phi_s<=0, C4e/2., C4ep) K4b = standalone_ddot44(K4a,standalone_ddot44(dlnbe4_s,dbe4_s)) K4c = standalone_dot42(-I4rt,tau)+K4b K4 = standalone_dot42(standalone_dot24(standalone_inv2(F),K4c),standalone_inv2(F).T) return P,tau,K4,be,ep, dlnbe4_s, dbe4_s, K4a, K4b, K4c # ----------------------------- TENSOR OPERATIONS ----------------------------- # tensor operations / products: np.einsum enables index notation, avoiding loops # e.g. ddot42 performs $C_ij = A_ijkl B_lk$ for the entire grid trans2 = lambda A2 : np.einsum('ijxyz ->jixyz ',A2 ) ddot22 = lambda A2,B2: np.einsum('ijxyz ,jixyz ->xyz ',A2,B2) ddot42 = lambda A4,B2: np.einsum('ijklxyz,lkxyz ->ijxyz ',A4,B2) ddot44 = lambda A4,B4: np.einsum('ijklxyz,lkmnxyz->ijmnxyz',A4,B4) dot11 = lambda A1,B1: np.einsum('ixyz ,ixyz ->xyz ',A1,B1) dot22 = lambda A2,B2: np.einsum('ijxyz ,jkxyz ->ikxyz ',A2,B2) dot24 = lambda A2,B4: np.einsum('ijxyz ,jkmnxyz->ikmnxyz',A2,B4) dot42 = lambda A4,B2: np.einsum('ijklxyz,lmxyz ->ijkmxyz',A4,B2) dyad22 = lambda A2,B2: np.einsum('ijxyz ,klxyz ->ijklxyz',A2,B2) dyad11 = lambda A1,B1: np.einsum('ixyz ,jxyz ->ijxyz ',A1,B1) # eigenvalue decomposition of 2nd-order tensor: return in convention i,j,x,y,z # NB requires to swap default order of NumPy (in in/output) def eig2(A2): swap1i = lambda A1: np.einsum('xyzi ->ixyz ',A1) swap2 = lambda A2: np.einsum('ijxyz->xyzij',A2) swap2i = lambda A2: np.einsum('xyzij->ijxyz',A2) vals,vecs = np.linalg.eig(swap2(A2)) vals = swap1i(vals) vecs = swap2i(vecs) return vals,vecs # logarithm of grid of 2nd-order tensors def ln2(A2): vals,vecs = eig2(A2) return sum([np.log(vals[i])*dyad11(vecs[:,i],vecs[:,i]) for i in range(3)]) # exponent of grid of 2nd-order tensors def exp2(A2): vals,vecs = eig2(A2) return sum([np.exp(vals[i])*dyad11(vecs[:,i],vecs[:,i]) for i in range(3)]) # determinant of grid of 2nd-order tensors def det2(A2): return (A2[0,0]*A2[1,1]*A2[2,2]+A2[0,1]*A2[1,2]*A2[2,0]+A2[0,2]*A2[1,0]*A2[2,1])-\ (A2[0,2]*A2[1,1]*A2[2,0]+A2[0,1]*A2[1,0]*A2[2,2]+A2[0,0]*A2[1,2]*A2[2,1]) # inverse of grid of 2nd-order tensors def inv2(A2): A2det = det2(A2) A2inv = np.empty([3,3,Nx,Ny,Nz]) A2inv[0,0] = (A2[1,1]*A2[2,2]-A2[1,2]*A2[2,1])/A2det A2inv[0,1] = (A2[0,2]*A2[2,1]-A2[0,1]*A2[2,2])/A2det A2inv[0,2] = (A2[0,1]*A2[1,2]-A2[0,2]*A2[1,1])/A2det A2inv[1,0] = (A2[1,2]*A2[2,0]-A2[1,0]*A2[2,2])/A2det A2inv[1,1] = (A2[0,0]*A2[2,2]-A2[0,2]*A2[2,0])/A2det A2inv[1,2] = (A2[0,2]*A2[1,0]-A2[0,0]*A2[1,2])/A2det A2inv[2,0] = (A2[1,0]*A2[2,1]-A2[1,1]*A2[2,0])/A2det A2inv[2,1] = (A2[0,1]*A2[2,0]-A2[0,0]*A2[2,1])/A2det A2inv[2,2] = (A2[0,0]*A2[1,1]-A2[0,1]*A2[1,0])/A2det return A2inv # ------------------------ INITIATE (IDENTITY) TENSORS ------------------------ # identity tensor (single tensor) i = np.eye(3) # identity tensors (grid) I = np.einsum('ij,xyz' , i ,np.ones([Nx,Ny,Nz])) I4 = np.einsum('ijkl,xyz->ijklxyz',np.einsum('il,jk',i,i),np.ones([Nx,Ny,Nz])) I4rt = np.einsum('ijkl,xyz->ijklxyz',np.einsum('ik,jl',i,i),np.ones([Nx,Ny,Nz])) I4s = (I4+I4rt)/2. II = dyad22(I,I) # ------------------------------------ FFT ------------------------------------ # projection operator (only for non-zero frequency, associated with the mean) # NB: vectorized version of "hyper-elasticity.py" # - allocate / support function Ghat4 = np.zeros([3,3,3,3,Nx,Ny,Nz]) # projection operator x = np.zeros([3 ,Nx,Ny,Nz],dtype='int64') # position vectors q = np.zeros([3 ,Nx,Ny,Nz],dtype='int64') # frequency vectors delta = lambda i,j: np.float(i==j) # Dirac delta function # - set "x" as position vector of all grid-points [grid of vector-components] x[0],x[1],x[2] = np.mgrid[:Nx,:Ny,:Nz] # - convert positions "x" to frequencies "q" [grid of vector-components] for i in range(3): freq = np.arange(-(shape[i]-1)/2,+(shape[i]+1)/2,dtype='int64') q[i] = freq[x[i]] # - compute "Q = ||q||", and "norm = 1/Q" being zero for the mean (Q==0) # NB: avoid zero division q = q.astype(np.float) Q = dot11(q,q) Z = Q==0 Q[Z] = 1. norm = 1./Q norm[Z] = 0. # - set projection operator [grid of tensors] for i, j, l, m in itertools.product(range(3), repeat=4): Ghat4[i,j,l,m] = norm*delta(i,m)*q[j]*q[l] # (inverse) Fourier transform (for each tensor component in each direction) fft = lambda x: np.fft.fftshift(np.fft.fftn (np.fft.ifftshift(x),[Nx,Ny,Nz])) ifft = lambda x: np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(x),[Nx,Ny,Nz])) # functions for the projection 'G', and the product 'G : K^LT : (delta F)^T' G = lambda A2 : np.real( ifft( ddot42(Ghat4,fft(A2)) ) ).reshape(-1) K_dF = lambda dFm: trans2(ddot42(K4,trans2(dFm.reshape(3,3,Nx,Ny,Nz)))) G_K_dF = lambda dFm: G(K_dF(dFm)) # --------------------------- CONSTITUTIVE RESPONSE --------------------------- # constitutive response to a certain loading and history # NB: completely uncoupled from the FFT-solver, but implemented as a regular # grid of quadrature points, to have an efficient code; # each point is completely independent, just evaluated at the same time def constitutive(F,F_t,be_t,ep_t): # function to compute linearization of the logarithmic Finger tensor def dln2_d2(A2): vals,vecs = eig2(A2) K4 = np.zeros([3,3,3,3,Nx,Ny,Nz]) for m, n in itertools.product(range(3),repeat=2): gc = (np.log(vals[n])-np.log(vals[m]))/(vals[n]-vals[m]) gc[vals[n]==vals[m]] = (1.0/vals[m])[vals[n]==vals[m]] K4 += gc*dyad22(dyad11(vecs[:,m],vecs[:,n]),dyad11(vecs[:,m],vecs[:,n])) return K4 # elastic stiffness tensor C4e = K*II+2.*mu*(I4s-1./3.*II) # trial state Fdelta = dot22(F,inv2(F_t)) be_s = dot22(Fdelta,dot22(be_t,trans2(Fdelta))) lnbe_s = ln2(be_s) tau_s = ddot42(C4e,lnbe_s)/2. taum_s = ddot22(tau_s,I)/3. taud_s = tau_s-taum_s*I taueq_s = np.sqrt(3./2.*ddot22(taud_s,taud_s)) div = np.where(taueq_s < 1e-12, np.ones_like(taueq_s), taueq_s) N_s = 3./2.*taud_s/div phi_s = taueq_s-(tauy0+H*ep_t) phi_s = 1./2.*(phi_s+np.abs(phi_s)) # return map dgamma = phi_s/(H+3.*mu) ep = ep_t + dgamma tau = tau_s -2.*dgamma*N_s*mu lnbe = lnbe_s-2.*dgamma*N_s be = exp2(lnbe) P = dot22(tau,trans2(inv2(F))) # consistent tangent operator a0 = dgamma*mu/taueq_s a1 = mu/(H+3.*mu) C4ep = ((K-2./3.*mu)/2.+a0*mu)*II+(1.-3.*a0)*mu*I4s+2.*mu*(a0-a1)*dyad22(N_s,N_s) dlnbe4_s = dln2_d2(be_s) dbe4_s = 2.*dot42(I4s,be_s) #K4a = (C4e/2.)*(phi_s<=0.).astype(np.float)+C4ep*(phi_s>0.).astype(np.float) K4a = np.where(phi_s<=0, C4e/2., C4ep) K4b = ddot44(K4a,ddot44(dlnbe4_s,dbe4_s)) K4c = dot42(-I4rt,tau)+K4b K4 = dot42(dot24(inv2(F),K4c),trans2(inv2(F))) return P,K4,be,ep, dlnbe4_s, dbe4_s, K4a, K4b, K4c # phase indicator: square inclusion of volume fraction (3*3*15)/(11*13*15) phase = np.zeros([Nx,Ny,Nz]); phase[0,0,0] = 1. # function to convert material parameters to grid of scalars param = lambda M0,M1: M0*np.ones([Nx,Ny,Nz])*(1.-phase)+\ M1*np.ones([Nx,Ny,Nz])* phase # material parameters K = param(0.833,0.833) # bulk modulus Kmat = K mu = param(0.386,0.386) # shear modulus H = param(0.004,0.008) # hardening modulus tauy0 = param(0.003,0.006) # initial yield stress # ---------------------------------- LOADING ---------------------------------- # stress, deformation gradient, plastic strain, elastic Finger tensor # NB "_t" signifies that it concerns the value at the previous increment ep_t = np.zeros([ Nx,Ny,Nz]) P = np.zeros([3,3,Nx,Ny,Nz]) F = np.array(I,copy=True) F_t = np.array(I,copy=True) be_t = np.array(I,copy=True) # initialize macroscopic incremental loading ninc = 50 lam = 0.0 barF = np.array(I,copy=True) barF_t = np.array(I,copy=True) # initial tangent operator: the elastic tangent K4 = K*II+2.*mu*(I4s-1./3.*II) class ElastoPlastic_Check(unittest.TestCase): t2_comparator = elastic_ref.LinearElastic_Check.t2_comparator t4_comparator = elastic_ref.LinearElastic_Check.t4_comparator def scalar_comparator(self, µ, g): err_sum = 0. err_max = 0. for counter, (i, j, k) in enumerate(self.rve): print((i,j,k)) µ_arr = µ[counter] g_arr = g[i,j,k] self.assertEqual(Nz*Ny*i+Nz*j + k, counter) print("µSpectre:") print(µ_arr) print("Goose:") print(g_arr) print(µ_arr-g_arr) err = linalg.norm(µ_arr-g_arr) print("error norm = {}".format(err)) err_sum += err err_max = max(err_max, err) pass print("∑(err) = {}, max(err) = {}".format (err_sum, err_max)) return err_sum def setUp(self): #---------------------------- µSpectre init ----------------------------------- resolution = list(phase.shape) dim = len(resolution) self.dim=dim center = np.array([r//2 for r in resolution]) incl = resolution[0]//5 ## Domain dimensions lengths = [float(r) for r in resolution] ## formulation (small_strain or finite_strain) formulation = µ.Formulation.finite_strain ## build a computational domain self.rve = µ.Cell(resolution, lengths, formulation) def get_E_nu(bulk, shear): Young = 9*bulk*shear/(3*bulk + shear) Poisson = Young/(2*shear) - 1 return Young, Poisson mat = µ.material.MaterialHyperElastoPlastic1_3d E, nu = get_E_nu(.833, .386) H = 0.004 tauy0 = .003 self.hard = mat.make(self.rve, 'hard', E, nu, 2*tauy0, 2*H) self.soft = mat.make(self.rve, 'soft', E, nu, tauy0, H) for pixel in self.rve: if pixel[0] == 0 and pixel[1] == 0 and pixel[2] == 0: self.hard.add_pixel(pixel) else: self.soft.add_pixel(pixel) pass pass return def test_solve(self): strict_tol = 1e-11 cg_tol = 1e-11 after_cg_tol = 1e-10 newton_tol = 1e-5 self.rve.set_uniform_strain(np.array(np.eye(ndim))) µF = self.rve.get_strain() µF_t = µF.copy() µbarF_t = µF.copy() # incremental deformation for inc in range(1,ninc): print('=============================') print('inc: {0:d}'.format(inc)) # set macroscopic deformation gradient (pure-shear) global lam, F, F_t, barF_t, K4, be_t, ep_t lam += 0.2/float(ninc) barF = np.array(I,copy=True) barF[0,0] = (1.+lam) barF[1,1] = 1./(1.+lam) def rel_error_scalar(µ, g, tol, do_assert=True): err = (linalg.norm(scalar_vec_to_goose(µ) - g.reshape(-1)) / linalg.norm(g)) if not (err < tol): self.scalar_comparator(µ.reshape(-1), g) if do_assert: self.assertLess(err, tol) else: print("AssertionWarning: {} is not less than {}".format( err, tol)) pass return err def rel_error_t2(µ, g, tol, do_assert=True): err = linalg.norm(t2_vec_to_goose(µ) - g.reshape(-1)) / linalg.norm(g) if not (err < tol): self.t2_comparator(µ.reshape(µF.shape), g.reshape(F.shape)) if do_assert: self.assertLess(err, tol) else: print("AssertionWarning: {} is not less than {}".format( err, tol)) pass return err def rel_error_t4(µ, g, tol, right_transposed=True, do_assert=True, pixel_tol=1e-4): err = linalg.norm(t4_vec_to_goose(µ) - g.reshape(-1)) / linalg.norm(g) errors = None if not (err < tol): err_sum, errors = self.t4_comparator(µ.reshape(µK.shape), g.reshape(K4.shape), right_transposed) if do_assert: self.assertLess(err, tol) else: print("AssertionWarning: {} is not less than {}".format( err, tol)) pass return err, errors def abs_error_t2(µ, g, tol, do_assert=True): ref_norm = linalg.norm(g) if ref_norm > 1: return rel_error_t2(µ, g, tol, do_assert) else: err = linalg.norm(t2_vec_to_goose(µ) - g.reshape(-1)) if not (err < tol): self.t2_comparator(µ.reshape(µF.shape), g.reshape(F.shape)) if do_assert: self.assertLess(err, tol) else: print("AssertionWarning: {} is not less than {}".format( err, tol)) return err rel_error_t2(µF, F, strict_tol) # store normalization Fn = np.linalg.norm(F) # first iteration residual: distribute "barF" over grid using "K4" b = -G_K_dF(barF-barF_t) F += barF-barF_t # parameters for Newton iterations: normalization and iteration counter Fn = np.linalg.norm(F) iiter = 0 # µSpectre inits µbarF = np.zeros_like(µF) µbarF[0, :] = 1. + lam µbarF[ndim + 1, :] = 1./(1. + lam) µbarF[-1, :] = 1. rel_error_t2(µbarF, barF, strict_tol) if inc == 1: µP, µK = self.rve.evaluate_stress_tangent(µF) rel_error_t4(µK, K4, strict_tol) µF += µbarF - µbarF_t rel_error_t2(µF, F, strict_tol) µFn = linalg.norm(µF) self.assertLess((µFn-Fn)/Fn, strict_tol) µG_K_dF = lambda x: self.rve.directional_stiffness(x.reshape(µF.shape)).reshape(-1) µG = lambda x: self.rve.project(x).reshape(-1) µb = -µG_K_dF(µbarF-µbarF_t) abs_error_t2(µb, b, strict_tol) # because of identical elastic properties, µb has got to # be zero before plasticity kicks in print("inc = {}".format(inc)) if inc == 1: self.assertLess(linalg.norm(µb), strict_tol) #print(self.hard.list_fields()) #print(self.hard.collection.statefield_names) #sys.exit() global_be_t = self.rve.get_globalised_current_real_array( "Previous left Cauchy-Green deformation bₑᵗ") # iterate as long as the iterative update does not vanish while True: # solve linear system using the Conjugate Gradient iterative solver g_counter = Counter() dFm,_ = sp.cg(tol=cg_tol, atol=1e-10, A = sp.LinearOperator(shape=(F.size,F.size),matvec=G_K_dF,dtype='float'), b = b, callback=g_counter ) µ_counter = Counter() µdFm,_ = sp.cg(tol=cg_tol, atol=1e-10, A = sp.LinearOperator(shape=(F.size,F.size), matvec=µG_K_dF,dtype='float'), b = µb, callback=µ_counter) err = g_counter.get()-µ_counter.get() if err != 0: print("n_iter(g) = {}, n_iter(µ) = {}".format(g_counter.get(), µ_counter.get())) print("AssertionWarning: {} != {}".format(g_counter.get(), µ_counter.get())) #self.assertEqual(g_counter.get(), µ_counter.get()) try: err = abs_error_t2(µdFm, dFm, after_cg_tol, do_assert=True) except Exception as err: raise err # add solution of linear system to DOFs F += dFm.reshape(3,3,Nx,Ny,Nz) µF += µdFm.reshape(µF.shape) err = rel_error_t2(µF, F, strict_tol, do_assert=True) # compute residual stress and tangent, convert to residual P,K4,be,ep,dln, dbe4_s, K4a, K4b, K4c = constitutive(F,F_t,be_t,ep_t) µP, µK = self.rve.evaluate_stress_tangent(µF) err = rel_error_t2(µP, P, strict_tol) µbe = self.rve.get_globalised_current_real_array( "Previous left Cauchy-Green deformation bₑᵗ") err = rel_error_t2(µbe, be, strict_tol) µep = self.rve.get_globalised_current_real_array( "cumulated plastic flow εₚ") err = rel_error_scalar(µep, ep, strict_tol) #µK4b = self.rve.get_globalised_internal_real_array( # "debug-T4") #err = rel_error_t4(µK4c, K4c, strict_tol, do_assert=True) err, errors = rel_error_t4(µK, K4, strict_tol, do_assert=False) if not err < strict_tol: def t2_disp(name, µ, g, i,j,k, index): g_v = g[:,:,i,j,k].copy() µ_v = µ[:,index].reshape(3,3).T print("{}_g =\n{}".format(name, g_v)) print("{}_µ =\n{}".format(name, µ_v)) print("{}_err = {}".format(name, np.linalg.norm(g_v-µ_v))) return g_v, µ_v def t0_disp(name, µ, g, i,j,k, index): g_v = g[i,j,k].copy() µ_v = µ[:,index] print("{}_g =\n{}".format(name, g_v)) print("{}_µ =\n{}".format(name, µ_v)) print("{}_err = {}".format(name, abs(g_v-µ_v))) return g_v, µ_v for pixel, error in errors.items(): for index, (ir,jr,kr) in enumerate(self.rve): i,j,k = pixel if (i,j,k) == (ir,jr,kr): break if error > 1e-4: i,j,k = pixel print("error for pixel {} ({}) = {}".format( pixel, index, error)) t2_disp("F", µF, F, i, j, k, index) t2_disp("be", µbe, be, i, j, k, index) F_t_n, dummy = t2_disp("F_t", µF_t, F_t, i,j,k,index) print("ep.shape = {}".format(µep.shape)) t0_disp("ep", µep, ep, i, j, k, index) K_comp_g = K4[:,:,:,:,i,j,k].reshape(9,9) K_comp_µ = µK[:, index].reshape( 3,3,3,3).transpose( 1,0,3,2).reshape(9,9) F_n = F[:,:,i,j,k] be_n = be_t[:,:,i,j,k] ep_n = ep_t[i,j,k] response = constitutive_standalone(Kmat[pixel], mu[pixel], H[pixel], tauy0[pixel], F_n, F_t_n, be_n, ep_n, 3) P_gn, K_gn = response[0], response[2] P_µn, K_µn = PK1_fun_3d(Kmat[pixel], mu[pixel], H[pixel], tauy0[pixel], F_n, F_t_n, be_n, ep_n) K_µn = K_µn.reshape(3,3,3,3).transpose(1,0,3,2).reshape(9,9) print("P_µn =\n{}".format(P_µn)) print("P_gn =\n{}".format(P_gn)) P_comp_gn, P_comp_µn = t2_disp( "P_comp", µP, P, i, j, k, index) print("|P_µn - P_comp_µn| = {}".format( np.linalg.norm(P_µn-P_comp_µn))) print("|P_gn - P_comp_gn| = {}".format( np.linalg.norm(P_gn-P_comp_gn))) print("|P_gn - P_µn| = {}".format( np.linalg.norm(P_gn-P_µn))) print("|P_comp_gn - P_comp_µn| = {}".format( np.linalg.norm(P_comp_gn-P_comp_µn))) print() K_gn.shape = 9, 9 print("K_µn.shape = {}".format(K_µn.shape)) print("K_gn.shape = {}".format(K_gn.shape)) print("K_comp_g =\n{}".format(K_comp_g)) print("K_comp_µ =\n{}".format(K_comp_µ)) print("K_µn=\n{}".format(K_µn)) print("K_gn=\n{}".format(K_gn)) print("|K_µn - K_comp_µ| = {}".format( np.linalg.norm(K_µn-K_comp_µ))) print("|K_gn - K_comp_g| = {}".format( np.linalg.norm(K_gn-K_comp_g))) print("|K_gn - K_µn| = {}".format( np.linalg.norm(K_gn-K_µn))) print("|K_comp_g - K_comp_µ| = {}".format( np.linalg.norm(K_comp_g-K_comp_µ))) break raise AssertionError( "at iiter = {}, inc = {}, caught this: '{}'".format( iiter, inc, err)) b = -G(P) µb = -µG(µP) err = abs_error_t2(µb, b, strict_tol)#after_cg_tol) #print("inc, iiter, err: {}".format((inc, iiter, err))) #print() # check for convergence, print convergence info to screen #print('{0:10.2e}'.format(np.linalg.norm(dFm)/Fn)) - print('Goose: rel_residual {:10.15e}, |rhs|: {:10.15e}'.format( + print('Goose: rel_residual {:10.15e}, |rhs|: {:10.15e}'.format( np.linalg.norm(dFm)/Fn, linalg.norm(b))) - print('µSpectre:rel_residual {:10.15e}, |rhs|: {:10.15e}'.format( + print('muSpectre:rel_residual {:10.15e}, |rhs|: {:10.15e}'.format( np.linalg.norm(µdFm)/µFn, linalg.norm(µb))) if np.linalg.norm(dFm)/Fn<1.e-5 and iiter>0: break # update Newton iteration counter print("reached end of iiter = {}".format(iiter)) iiter += 1 # end-of-increment: update history barF_t = np.array(barF,copy=True) µbarF_t[:] = µbarF F_t = np.array(F ,copy=True) be_t = np.array(be ,copy=True) ep_t = np.array(ep ,copy=True) µF_t[:] = µF self.rve.save_history_variables() if __name__ == '__main__': unittest.main()