diff --git a/PRIVATE b/PRIVATE index 6ce67b4cd..27add707f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 6ce67b4cd86df6402b5e43e354e1c77f685aa469 +Subproject commit 27add707f520816fc05e81d2b00769cd60216010 diff --git a/VERSION b/VERSION index eb33473c3..4e38c952a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha3-51-g38ca85b58 +v3.0.0-alpha3-73-g43db3bd8d diff --git a/examples/config/phase/thermal/Al.yaml b/examples/config/phase/thermal/Al.yaml new file mode 100644 index 000000000..8837f30af --- /dev/null +++ b/examples/config/phase/thermal/Al.yaml @@ -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 diff --git a/examples/config/phase/thermal/Steel-0.5C.yaml b/examples/config/phase/thermal/Steel-0.5C.yaml new file mode 100644 index 000000000..861cbbee5 --- /dev/null +++ b/examples/config/phase/thermal/Steel-0.5C.yaml @@ -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 diff --git a/python/damask/seeds.py b/python/damask/seeds.py index 126f138aa..b2e91a352 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -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]) diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 5ddc7219e..b0bb96402 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -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" diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 1ccb3b901..5a1745668 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -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) diff --git a/src/phase.f90 b/src/phase.f90 index 93743d64a..927ba6b2b 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -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 :: & diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 298ff57f3..29e9e8ada 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -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) !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 new file mode 100644 index 000000000..9c18bfa7f --- /dev/null +++ b/src/phase_mechanical_elastic.f90 @@ -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