223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/development' into MatrixInteraction_clean

This commit is contained in:
Martin Diehl 2021-05-09 08:04:15 +02:00
parent 0d66320182 62203bcddb
commit 693c2f4e3f
10 changed files with 181 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 6ce67b4cd86df6402b5e43e354e1c77f685aa469
Subproject commit 27add707f520816fc05e81d2b00769cd60216010

2
VERSION
View File

@ -1 +1 @@
v3.0.0-alpha3-51-g38ca85b58
v3.0.0-alpha3-73-g43db3bd8d

7
examples/config/phase/thermal/Al.yaml Normal file
View File

@ -0,0 +1,7 @@
references:
- www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html
- www.engineeringtoolbox.com/specific-heat-metals-d_152.html
c_p: 0.91e3
K_11: 236.0
K_22: 236.0
K_33: 236.0

7
examples/config/phase/thermal/Steel-0.5C.yaml Normal file
View File

@ -0,0 +1,7 @@
references:
- www.engineeringtoolbox.com/thermal-conductivity-metals-d_858.html
- www.engineeringtoolbox.com/specific-heat-metals-d_152.html
c_p: 0.49e3
K_11: 54.0
K_22: 54.0
K_33: 54.0

10
python/damask/seeds.py
View File

@ -3,8 +3,8 @@
from scipy import spatial as _spatial
import numpy as _np
from . import util
from . import grid_filters
from . import util as _util
from . import grid_filters as _grid_filters
def from_random(size,N_seeds,cells=None,rng_seed=None):
@ -34,7 +34,7 @@ def from_random(size,N_seeds,cells=None,rng_seed=None):
if cells is None:
coords = rng.random((N_seeds,3)) * size
else:
grid_coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F')
grid_coords = _grid_filters.coordinates0_point(cells,size).reshape(-1,3,order='F')
coords = grid_coords[rng.choice(_np.prod(cells),N_seeds, replace=False)] \
+ _np.broadcast_to(size/cells,(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25) # wobble without leaving cells
@ -73,7 +73,7 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed=
s = 1
i = 0
progress = util._ProgressBar(N_seeds+1,'',50)
progress = _util._ProgressBar(N_seeds+1,'',50)
while s < N_seeds:
candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3))
tree = _spatial.cKDTree(coords[:s],boxsize=size) if periodic else \
@ -120,7 +120,7 @@ def from_grid(grid,selection=None,invert=False,average=False,periodic=True):
material = grid.material.reshape((-1,1),order='F')
mask = _np.full(grid.cells.prod(),True,dtype=bool) if selection is None else \
_np.isin(material,selection,invert=invert).flatten()
coords = grid_filters.coordinates0_point(grid.cells,grid.size).reshape(-1,3,order='F')
coords = _grid_filters.coordinates0_point(grid.cells,grid.size).reshape(-1,3,order='F')
if not average:
return (coords[mask],material[mask])

1
src/commercialFEM_fileList.f90
View File

@ -21,6 +21,7 @@
#include "lattice.f90"
#include "phase.f90"
#include "phase_mechanical.f90"
#include "phase_mechanical_elastic.f90"
#include "phase_mechanical_plastic.f90"
#include "phase_mechanical_plastic_none.f90"
#include "phase_mechanical_plastic_isotropic.f90"

3
src/grid/grid_thermal_spectral.f90
View File

@ -90,7 +90,8 @@ subroutine grid_thermal_spectral_init(T_0)
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type ngmres',ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type newtonls -thermal_snes_mf &
&-thermal_snes_ksp_ew -thermal_ksp_type fgmres',ierr)
CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr)
CHKERRQ(ierr)

5
src/phase.f90
View File

@ -58,10 +58,6 @@ module phase
grain
end type tDebugOptions
integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler)
phase_elasticityInstance, &
phase_NstiffnessDegradations
logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler)
phase_localPlasticity !< flags phases with local constitutive law
@ -298,7 +294,6 @@ module phase
end interface
type(tDebugOptions) :: debugConstitutive
#if __INTEL_COMPILER >= 1900
public :: &

131
src/phase_mechanical.f90
View File

@ -5,10 +5,6 @@ submodule(phase) mechanical
enum, bind(c); enumerator :: &
ELASTICITY_UNDEFINED_ID, &
ELASTICITY_HOOKE_ID, &
STIFFNESS_DEGRADATION_UNDEFINED_ID, &
STIFFNESS_DEGRADATION_DAMAGE_ID, &
PLASTICITY_UNDEFINED_ID, &
PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, &
@ -23,11 +19,6 @@ submodule(phase) mechanical
KINEMATICS_THERMAL_EXPANSION_ID
end enum
integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: &
phase_elasticity !< elasticity of each phase
integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: &
phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase
type(tTensorContainer), dimension(:), allocatable :: &
! current value
phase_mechanical_Fe, &
@ -57,9 +48,27 @@ submodule(phase) mechanical
class(tNode), pointer :: phases
end subroutine eigendeformation_init
module subroutine elastic_init(phases)
class(tNode), pointer :: phases
end subroutine elastic_init
module subroutine plastic_init
end subroutine plastic_init
module subroutine phase_hooke_SandItsTangents(S,dS_dFe,dS_dFi,Fe,Fi,ph,en)
integer, intent(in) :: &
ph, &
en
real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real(pReal), intent(out), dimension(3,3,3,3) :: &
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
end subroutine phase_hooke_SandItsTangents
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
@ -73,7 +82,6 @@ submodule(phase) mechanical
end subroutine plastic_isotropic_LiAndItsTangent
module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
@ -198,17 +206,11 @@ module subroutine mechanical_init(materials,phases)
constituents, &
constituent, &
phase, &
mech, &
elastic, &
stiffDegradation
mech
print'(/,a)', ' <<<+- phase:mechanical init -+>>>'
!-------------------------------------------------------------------------------------------------
! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC?
allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
allocate(phase_elasticityInstance(phases%length), source = 0)
allocate(phase_NstiffnessDegradations(phases%length),source=0)
allocate(output_constituent(phases%length))
allocate(phase_mechanical_Fe(phases%length))
@ -253,32 +255,8 @@ module subroutine mechanical_init(materials,phases)
#else
output_constituent(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
#endif
elastic => mech%get('elastic')
if (IO_lc(elastic%get_asString('type')) == 'hooke') then ! accept small letter h for the moment
phase_elasticity(ph) = ELASTICITY_HOOKE_ID
else
call IO_error(200,ext_msg=elastic%get_asString('type'))
endif
stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms
phase_NstiffnessDegradations(ph) = stiffDegradation%length
enddo
allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), &
source=STIFFNESS_DEGRADATION_undefined_ID)
if(maxVal(phase_NstiffnessDegradations)/=0) then
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanical')
stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList)
do stiffDegradationCtr = 1, stiffDegradation%length
if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') &
phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID
enddo
enddo
endif
do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2)
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
material => materials%get(discretization_materialAt(el))
@ -306,6 +284,9 @@ module subroutine mechanical_init(materials,phases)
enddo; enddo
! initialize elasticity
call elastic_init(phases)
! initialize plasticity
allocate(plasticState(phases%length))
allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID)
@ -313,9 +294,6 @@ module subroutine mechanical_init(materials,phases)
call plastic_init()
do ph = 1, phases%length
phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph))
enddo
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -348,51 +326,6 @@ module subroutine mechanical_init(materials,phases)
end subroutine mechanical_init
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermediate deformation gradients using Hooke's law
!--------------------------------------------------------------------------------------------------
subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, ph, en)
integer, intent(in) :: &
ph, &
en
real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real(pReal), intent(out), dimension(3,3,3,3) :: &
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C
integer :: &
d, & !< counter in degradation loop
i, j
C = math_66toSym3333(phase_homogenizedC(ph,en))
DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph)
degradationType: select case(phase_stiffnessDegradation(d,ph))
case (STIFFNESS_DEGRADATION_damage_ID) degradationType
C = C * damage_phi(ph,en)**2
end select degradationType
enddo DegradationLoop
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
do i =1, 3;do j=1,3
dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
enddo; enddo
end subroutine phase_hooke_SandItsTangents
module subroutine mechanical_results(group,ph)
character(len=*), intent(in) :: group
@ -1082,26 +1015,6 @@ module subroutine mechanical_forward()
end subroutine mechanical_forward
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent
!--------------------------------------------------------------------------------------------------
module function phase_homogenizedC(ph,en) result(C)
real(pReal), dimension(6,6) :: C
integer, intent(in) :: ph, en
plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,en)
case default plasticType
C = lattice_C66(1:6,1:6,ph)
end select plasticType
end function phase_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------

135
src/phase_mechanical_elastic.f90 Normal file
View File

@ -0,0 +1,135 @@
submodule(phase:mechanical) elastic
enum, bind(c); enumerator :: &
ELASTICITY_UNDEFINED_ID, &
ELASTICITY_HOOKE_ID, &
STIFFNESS_DEGRADATION_UNDEFINED_ID, &
STIFFNESS_DEGRADATION_DAMAGE_ID
end enum
integer, dimension(:), allocatable :: &
phase_NstiffnessDegradations
integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: &
phase_elasticity !< elasticity of each phase
integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: &
phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase
contains
module subroutine elastic_init(phases)
class(tNode), pointer :: &
phases
integer :: &
ph, &
stiffDegradationCtr
class(tNode), pointer :: &
phase, &
mech, &
elastic, &
stiffDegradation
print'(/,a)', ' <<<+- phase:mechanical:elastic init -+>>>'
allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
allocate(phase_NstiffnessDegradations(phases%length),source=0)
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanical')
elastic => mech%get('elastic')
if(IO_lc(elastic%get_asString('type')) == 'hooke') then ! accept small letter h for the moment
phase_elasticity(ph) = ELASTICITY_HOOKE_ID
else
call IO_error(200,ext_msg=elastic%get_asString('type'))
endif
stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms
phase_NstiffnessDegradations(ph) = stiffDegradation%length
enddo
allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), &
source=STIFFNESS_DEGRADATION_undefined_ID)
if(maxVal(phase_NstiffnessDegradations)/=0) then
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanical')
stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList)
do stiffDegradationCtr = 1, stiffDegradation%length
if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') &
phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID
enddo
enddo
endif
end subroutine elastic_init
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermediate deformation gradients using Hooke's law
!--------------------------------------------------------------------------------------------------
module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, ph, en)
integer, intent(in) :: &
ph, &
en
real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real(pReal), intent(out), dimension(3,3,3,3) :: &
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C
integer :: &
d, & !< counter in degradation loop
i, j
C = math_66toSym3333(phase_homogenizedC(ph,en))
DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph)
degradationType: select case(phase_stiffnessDegradation(d,ph))
case (STIFFNESS_DEGRADATION_damage_ID) degradationType
C = C * damage_phi(ph,en)**2
end select degradationType
enddo DegradationLoop
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
do i =1, 3;do j=1,3
dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
enddo; enddo
end subroutine phase_hooke_SandItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix
!> ToDo: homogenizedC66 would be more consistent
!--------------------------------------------------------------------------------------------------
module function phase_homogenizedC(ph,en) result(C)
real(pReal), dimension(6,6) :: C
integer, intent(in) :: ph, en
plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,en)
case default plasticType
C = lattice_C66(1:6,1:6,ph)
end select plasticType
end function phase_homogenizedC
end submodule elastic