From 556d9d840e7fbfeb73974249a5da2ac423a56200 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 17 May 2022 23:56:05 +0200 Subject: [PATCH 1/3] specifying V_e is more natural than F_i --- PRIVATE | 2 +- python/damask/_configmaterial.py | 6 +- python/tests/test_ConfigMaterial.py | 6 +- src/IO.f90 | 2 + src/material.f90 | 10 ++-- src/phase_mechanical.f90 | 89 +++++++++++++++-------------- src/rotations.f90 | 9 ++- 7 files changed, 64 insertions(+), 60 deletions(-) diff --git a/PRIVATE b/PRIVATE index 1766c855e..d2c229e6d 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 1766c855e54668d6225c9b22db536be7052605bc +Subproject commit d2c229e6dd81fb399b2c1b0658b750f7d309f175 diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py index 725567b7f..7a06a2c69 100644 --- a/python/damask/_configmaterial.py +++ b/python/damask/_configmaterial.py @@ -482,7 +482,7 @@ class ConfigMaterial(Config): """ N,n,shaped = 1,1,{} - map_dim = {'O':-1,'F_i':-2} + map_dim = {'O':-1,'V_e':-2} for k,v in kwargs.items(): shaped[k] = np.array(v) s = shaped[k].shape[:map_dim.get(k,None)] @@ -494,12 +494,12 @@ class ConfigMaterial(Config): if 'v' not in kwargs: shaped['v'] = np.broadcast_to(1/n,(N,n)) - map_shape = {'O':(N,n,4),'F_i':(N,n,3,3)} + map_shape = {'O':(N,n,4),'V_e':(N,n,3,3)} for k,v in shaped.items(): target = map_shape.get(k,(N,n)) obj = np.broadcast_to(v.reshape(util.shapeshifter(v.shape, target, mode = 'right')), target) for i in range(N): - if k in ['phase','O','v','F_i']: + if k in ['phase','O','v','V_e']: for j in range(n): mat[i]['constituents'][j][k] = obj[i,j].item() if isinstance(obj[i,j],np.generic) else obj[i,j] else: diff --git a/python/tests/test_ConfigMaterial.py b/python/tests/test_ConfigMaterial.py index d3ea986c5..f003e254e 100644 --- a/python/tests/test_ConfigMaterial.py +++ b/python/tests/test_ConfigMaterial.py @@ -113,15 +113,15 @@ class TestConfigMaterial: @pytest.mark.parametrize('N,n,kw',[ (1,1,{'phase':'Gold', 'O':[1,0,0,0], - 'F_i':np.eye(3), + 'V_e':np.eye(3), 'homogenization':'SX'}), (3,1,{'phase':'Gold', 'O':Rotation.from_random(3), - 'F_i':np.broadcast_to(np.eye(3),(3,3,3)), + 'V_e':np.broadcast_to(np.eye(3),(3,3,3)), 'homogenization':'SX'}), (2,3,{'phase':np.broadcast_to(['a','b','c'],(2,3)), 'O':Rotation.from_random((2,3)), - 'F_i':np.broadcast_to(np.eye(3),(2,3,3,3)), + 'V_e':np.broadcast_to(np.eye(3),(2,3,3,3)), 'homogenization':['SX','PX']}), ]) def test_material_add(self,kw,N,n): diff --git a/src/IO.f90 b/src/IO.f90 index 240c9e992..047d11add 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -428,6 +428,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'too many systems requested' case (146) msg = 'number of values does not match' + case (147) + msg = 'V_e needs to be symmetric' case (148) msg = 'Nconstituents mismatch between homogenization and material' diff --git a/src/material.f90 b/src/material.f90 index 2258c40a1..ce4cabea9 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -27,7 +27,7 @@ module material type(tRotationContainer), dimension(:), allocatable, public, protected :: material_O_0 - type(tTensorContainer), dimension(:), allocatable, public, protected :: material_F_i_0 + type(tTensorContainer), dimension(:), allocatable, public, protected :: material_V_e_0 integer, dimension(:), allocatable, public, protected :: & homogenization_Nconstituents !< number of grains in each homogenization @@ -133,7 +133,7 @@ subroutine parse() allocate(material_v(homogenization_maxNconstituents,discretization_Ncells),source=0.0_pReal) allocate(material_O_0(materials%length)) - allocate(material_F_i_0(materials%length)) + allocate(material_V_e_0(materials%length)) allocate(ho_of(materials%length)) allocate(ph_of(materials%length,homogenization_maxNconstituents),source=-1) @@ -154,7 +154,7 @@ subroutine parse() if (constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148) allocate(material_O_0(ma)%data(constituents%length)) - allocate(material_F_i_0(ma)%data(1:3,1:3,constituents%length)) + allocate(material_V_e_0(ma)%data(1:3,1:3,constituents%length)) do co = 1, constituents%length constituent => constituents%get(co) @@ -162,7 +162,9 @@ subroutine parse() ph_of(ma,co) = phases%getIndex(constituent%get_asString('phase')) call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) - material_F_i_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('F_i',defaultVal=math_I3,requiredShape=[3,3]) + material_V_e_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('V_e',defaultVal=math_I3,requiredShape=[3,3]) + if (any(dNeq(material_V_e_0(ma)%data(1:3,1:3,co),transpose(material_V_e_0(ma)%data(1:3,1:3,co))))) & + call IO_error(147) end do if (dNeq(sum(v_of(ma,:)),1.0_pReal,1.e-9_pReal)) call IO_error(153,ext_msg='constituent') diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index fc9e40e02..bddd4775e 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -210,7 +210,6 @@ module subroutine mechanical_init(phases) Nmembers class(tNode), pointer :: & num_crystallite, & - phase, & mech @@ -239,11 +238,8 @@ module subroutine mechanical_init(phases) allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_F(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers),source=0.0_pReal) @@ -265,23 +261,21 @@ module subroutine mechanical_init(phases) ma = discretization_materialAt((ce-1)/discretization_nIPs+1) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) ph = material_phaseID(co,ce) - phase_mechanical_Fi0(ph)%data(1:3,1:3,material_phaseEntry(co,ce)) = material_F_i_0(ma)%data(1:3,1:3,co) + en = material_phaseEntry(co,ce) + phase_mechanical_F(ph)%data(1:3,1:3,en) = math_I3 + phase_mechanical_Fp(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) + phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp(ph)%data(1:3,1:3,en) & + / math_det33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal) + phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(material_V_e_0(ma)%data(1:3,1:3,co), & + transpose(phase_mechanical_Fp(ph)%data(1:3,1:3,en))) + phase_mechanical_Fi(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%rotate(math_inv33(material_V_e_0(ma)%data(1:3,1:3,co))) enddo enddo do ph = 1, phases%length - do en = 1, count(material_phaseID == ph) - - phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_O_0(ph)%data(en)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) - phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) & - / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal) - phase_mechanical_F0(ph)%data(1:3,1:3,en) = math_I3 - phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), & - phase_mechanical_Fp0(ph)%data(1:3,1:3,en))) ! assuming that euler angles are given in internal strain free configuration - enddo - phase_mechanical_Fp(ph)%data = phase_mechanical_Fp0(ph)%data - phase_mechanical_Fi(ph)%data = phase_mechanical_Fi0(ph)%data - phase_mechanical_F(ph)%data = phase_mechanical_F0(ph)%data + phase_mechanical_F0(ph)%data = phase_mechanical_F(ph)%data + phase_mechanical_Fp0(ph)%data = phase_mechanical_Fp(ph)%data + phase_mechanical_Fi0(ph)%data = phase_mechanical_Fi(ph)%data enddo @@ -318,7 +312,6 @@ module subroutine mechanical_init(phases) end select - call eigen_init(phases) @@ -1244,6 +1237,9 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF) end function phase_mechanical_dPdF +!-------------------------------------------------------------------------------------------------- +!< @brief Write restart information to file. +!-------------------------------------------------------------------------------------------------- module subroutine mechanical_restartWrite(groupHandle,ph) integer(HID_T), intent(in) :: groupHandle @@ -1251,36 +1247,41 @@ module subroutine mechanical_restartWrite(groupHandle,ph) call HDF5_write(plasticState(ph)%state,groupHandle,'omega_plastic') - call HDF5_write(phase_mechanical_Fi(ph)%data,groupHandle,'F_i') - call HDF5_write(phase_mechanical_Li(ph)%data,groupHandle,'L_i') - call HDF5_write(phase_mechanical_Lp(ph)%data,groupHandle,'L_p') - call HDF5_write(phase_mechanical_Fp(ph)%data,groupHandle,'F_p') call HDF5_write(phase_mechanical_S(ph)%data,groupHandle,'S') call HDF5_write(phase_mechanical_F(ph)%data,groupHandle,'F') + call HDF5_write(phase_mechanical_Fp(ph)%data,groupHandle,'F_p') + call HDF5_write(phase_mechanical_Fi(ph)%data,groupHandle,'F_i') + call HDF5_write(phase_mechanical_Lp(ph)%data,groupHandle,'L_p') + call HDF5_write(phase_mechanical_Li(ph)%data,groupHandle,'L_i') end subroutine mechanical_restartWrite +!-------------------------------------------------------------------------------------------------- +!< @brief Read restart information from file. +!-------------------------------------------------------------------------------------------------- module subroutine mechanical_restartRead(groupHandle,ph) integer(HID_T), intent(in) :: groupHandle integer, intent(in) :: ph + integer :: en + call HDF5_read(plasticState(ph)%state0,groupHandle,'omega_plastic') - call HDF5_read(phase_mechanical_Fi0(ph)%data,groupHandle,'F_i') - call HDF5_read(phase_mechanical_Li0(ph)%data,groupHandle,'L_i') - call HDF5_read(phase_mechanical_Lp0(ph)%data,groupHandle,'L_p') - call HDF5_read(phase_mechanical_Fp0(ph)%data,groupHandle,'F_p') call HDF5_read(phase_mechanical_S0(ph)%data,groupHandle,'S') call HDF5_read(phase_mechanical_F0(ph)%data,groupHandle,'F') + call HDF5_read(phase_mechanical_Fp0(ph)%data,groupHandle,'F_p') + call HDF5_read(phase_mechanical_Fi0(ph)%data,groupHandle,'F_i') + call HDF5_read(phase_mechanical_Lp0(ph)%data,groupHandle,'L_p') + call HDF5_read(phase_mechanical_Li0(ph)%data,groupHandle,'L_i') end subroutine mechanical_restartRead -!---------------------------------------------------------------------------------------------- -!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics) -!---------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics). +!-------------------------------------------------------------------------------------------------- module function mechanical_S(ph,en) result(S) integer, intent(in) :: ph,en @@ -1292,9 +1293,9 @@ module function mechanical_S(ph,en) result(S) end function mechanical_S -!---------------------------------------------------------------------------------------------- -!< @brief Get plastic velocity gradient (for use by non-mech physics) -!---------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +!< @brief Get plastic velocity gradient (for use by non-mech physics). +!-------------------------------------------------------------------------------------------------- module function mechanical_L_p(ph,en) result(L_p) integer, intent(in) :: ph,en @@ -1306,9 +1307,9 @@ module function mechanical_L_p(ph,en) result(L_p) end function mechanical_L_p -!---------------------------------------------------------------------------------------------- -!< @brief Get elastic deformation gradient (for use by non-mech physics) -!---------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +!< @brief Get elastic deformation gradient (for use by non-mech physics). +!-------------------------------------------------------------------------------------------------- module function mechanical_F_e(ph,en) result(F_e) integer, intent(in) :: ph,en @@ -1320,9 +1321,9 @@ module function mechanical_F_e(ph,en) result(F_e) end function mechanical_F_e -!---------------------------------------------------------------------------------------------- -!< @brief Get second Piola-Kichhoff stress (for use by homogenization) -!---------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +!< @brief Get second Piola-Kichhoff stress (for use by homogenization). +!-------------------------------------------------------------------------------------------------- module function phase_P(co,ce) result(P) integer, intent(in) :: co, ce @@ -1334,9 +1335,9 @@ module function phase_P(co,ce) result(P) end function phase_P -!---------------------------------------------------------------------------------------------- -!< @brief Get deformation gradient (for use by homogenization) -!---------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +!< @brief Get deformation gradient (for use by homogenization). +!-------------------------------------------------------------------------------------------------- module function phase_F(co,ce) result(F) integer, intent(in) :: co, ce @@ -1348,9 +1349,9 @@ module function phase_F(co,ce) result(F) end function phase_F -!---------------------------------------------------------------------------------------------- -!< @brief Set deformation gradient (for use by homogenization) -!---------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +!< @brief Set deformation gradient (for use by homogenization). +!-------------------------------------------------------------------------------------------------- module subroutine phase_set_F(F,co,ce) real(pReal), dimension(3,3), intent(in) :: F diff --git a/src/rotations.f90 b/src/rotations.f90 index 7aa81fcab..3340bf2a0 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -323,17 +323,16 @@ pure function rotTensor2(self,T,active) result(tRot) logical :: passive + if (present(active)) then passive = .not. active else passive = .true. endif - if (passive) then - tRot = matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())) - else - tRot = matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()) - endif + tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), & + matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), & + passive) end function rotTensor2 From 491e2ec0b2b963353e09347c38383c78710477fa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 20 May 2022 06:30:07 +0200 Subject: [PATCH 2/3] avoid negative zero when not needed --- src/material.f90 | 2 +- src/rotations.f90 | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index ce4cabea9..eafc411cd 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -164,7 +164,7 @@ subroutine parse() call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4)) material_V_e_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('V_e',defaultVal=math_I3,requiredShape=[3,3]) if (any(dNeq(material_V_e_0(ma)%data(1:3,1:3,co),transpose(material_V_e_0(ma)%data(1:3,1:3,co))))) & - call IO_error(147) + call IO_error(147) end do if (dNeq(sum(v_of(ma,:)),1.0_pReal,1.e-9_pReal)) call IO_error(153,ext_msg='constituent') diff --git a/src/rotations.f90 b/src/rotations.f90 index 3340bf2a0..6257fe7e4 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -574,7 +574,7 @@ end function qu2cu !--------------------------------------------------------------------------------------------------- !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief convert rotation matrix to cubochoric +!> @brief convert rotation matrix to unit quaternion !> @details the original formulation (direct conversion) had (numerical?) issues !--------------------------------------------------------------------------------------------------- pure function om2qu(om) result(qu) @@ -601,14 +601,14 @@ pure function om2qu(om) result(qu) endif endif if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu - qu = qu*[1.0_pReal,P,P,P] + qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4))) end function om2qu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief orientation matrix to Euler angles +!> @brief convert orientation matrix to Euler angles !> @details Two step check for special cases to avoid invalid operations (not needed for python) !--------------------------------------------------------------------------------------------------- pure function om2eu(om) result(eu) @@ -1333,8 +1333,8 @@ pure function cu2ho(cu) result(ho) ! transform to sphere grid (inverse Lambert) ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero] c = sum(T**2) - s = PI * c/(24.0*XYZ(3)**2) - c = sqrt(PI) * c / sqrt(24.0_pReal) / XYZ(3) + s = c * PI/(24.0*XYZ(3)**2) + c = c * sqrt(PI/24.0_pReal) / XYZ(3) q = sqrt( 1.0 - s ) LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ] end if special From f42bb5d1755a77ea02a9d34a42fdff68b2746fef Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 20 May 2022 06:43:32 +0200 Subject: [PATCH 3/3] mitigate rounding errors --- src/phase_mechanical.f90 | 2 -- src/rotations.f90 | 2 ++ 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index bddd4775e..88b86a8d9 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -264,8 +264,6 @@ module subroutine mechanical_init(phases) en = material_phaseEntry(co,ce) phase_mechanical_F(ph)%data(1:3,1:3,en) = math_I3 phase_mechanical_Fp(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) - phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp(ph)%data(1:3,1:3,en) & - / math_det33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal) phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(material_V_e_0(ma)%data(1:3,1:3,co), & transpose(phase_mechanical_Fp(ph)%data(1:3,1:3,en))) phase_mechanical_Fi(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%rotate(math_inv33(material_V_e_0(ma)%data(1:3,1:3,co))) diff --git a/src/rotations.f90 b/src/rotations.f90 index 6257fe7e4..5367f2c07 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -449,6 +449,7 @@ pure function qu2om(qu) result(om) om(1,3) = 2.0_pReal*(qu(2)*qu(4)+qu(1)*qu(3)) if (sign(1.0_pReal,P) < 0.0_pReal) om = transpose(om) + om = om/math_det33(om)**(1.0_pReal/3.0_pReal) end function qu2om @@ -602,6 +603,7 @@ pure function om2qu(om) result(qu) endif if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4))) + qu = qu/norm2(qu) end function om2qu