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

descriptive names

This commit is contained in:
Martin Diehl 2021-07-17 11:50:21 +02:00
parent e6d25294d3
commit f9edeb40a5
13 changed files with 83 additions and 96 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/CPFEM.f90
View File

@ -190,9 +190,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
else validCalculation
if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP])
call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP])
if (.not. terminallyIll) &
call materialpoint_stressAndItsTangent2(dt,[ip,ip],[elCP,elCP])
call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP])
terminalIllness: if (terminallyIll) then

2
src/grid/grid_mech_FEM.f90
View File

@ -261,7 +261,7 @@ subroutine grid_mechanical_FEM_init
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(F)
call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
F, & ! target F

2
src/grid/grid_mech_spectral_basic.f90
View File

@ -211,7 +211,7 @@ subroutine grid_mechanical_spectral_basic_init
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F

2
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -237,7 +237,7 @@ subroutine grid_mechanical_spectral_polarisation_init
F_tau_lastInc = 2.0_pReal*F_lastInc
endif restartRead
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F

6
src/grid/spectral_utilities.f90
View File

@ -815,11 +815,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field
call homogenization_mechanical_response(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field
if (.not. terminallyIll) &
call materialpoint_stressAndItsTangent3(timeinc,[1,1],[1,product(grid(1:2))*grid3])
call homogenization_thermal_response(timeinc,[1,1],[1,product(grid(1:2))*grid3])
if (.not. terminallyIll) &
call materialpoint_stressAndItsTangent2(timeinc,[1,1],[1,product(grid(1:2))*grid3])
call homogenization_mechanical_response2(timeinc,[1,1],[1,product(grid(1:2))*grid3])
P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt

66
src/homogenization.f90
View File

@ -101,8 +101,8 @@ module homogenization
integer, intent(in) :: ce
end subroutine damage_partition
module subroutine mechanical_homogenize(dt,ce)
real(pReal), intent(in) :: dt
module subroutine mechanical_homogenize(Delta_t,ce)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
ce !< cell
end subroutine mechanical_homogenize
@ -178,9 +178,9 @@ module homogenization
public :: &
homogenization_init, &
materialpoint_stressAndItsTangent, &
materialpoint_stressAndItsTangent2, &
materialpoint_stressAndItsTangent3, &
homogenization_mechanical_response, &
homogenization_mechanical_response2, &
homogenization_thermal_response, &
homogenization_mu_T, &
homogenization_K_T, &
homogenization_f_T, &
@ -231,15 +231,15 @@ end subroutine homogenization_init
!--------------------------------------------------------------------------------------------------
!> @brief
!--------------------------------------------------------------------------------------------------
subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem)
subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
real(pReal), intent(in) :: dt !< time increment
real(pReal), intent(in) :: Delta_t !< time increment
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
integer :: &
NiterationMPstate, &
ip, & !< integration point number
el, & !< element number
co, ce, ho, en, ph
co, ce, ho, en
logical :: &
converged
logical, dimension(2) :: &
@ -268,13 +268,10 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
NiterationMPstate = NiterationMPstate + 1
call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
converged = .true.
do co = 1, homogenization_Nconstituents(ho)
converged = converged .and. phase_mechanical_constitutive(dt,co,ip,el)
enddo
converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
if (converged) then
doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ce)
doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce)
converged = all(doneAndHappy)
else
doneAndHappy = [.true.,.false.]
@ -282,34 +279,30 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
enddo convergenceLooping
if (.not. converged) then
if (.not. terminallyIll) print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
endif
enddo
enddo
!$OMP END PARALLEL DO
end subroutine materialpoint_stressAndItsTangent
end subroutine homogenization_mechanical_response
!--------------------------------------------------------------------------------------------------
!> @brief
!--------------------------------------------------------------------------------------------------
subroutine materialpoint_stressAndItsTangent3(dt,FEsolving_execIP,FEsolving_execElem)
subroutine homogenization_thermal_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
real(pReal), intent(in) :: dt !< time increment
real(pReal), intent(in) :: Delta_t !< time increment
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
integer :: &
NiterationMPstate, &
ip, & !< integration point number
el, & !< element number
co, ce, ho, en, ph
logical :: &
converged
logical, dimension(2) :: &
doneAndHappy
co, ce, ho
!$OMP PARALLEL DO PRIVATE(ho,ph,ce)
!$OMP PARALLEL DO PRIVATE(ho,ce)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
if (terminallyIll) continue
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
@ -317,36 +310,29 @@ subroutine materialpoint_stressAndItsTangent3(dt,FEsolving_execIP,FEsolving_exec
ho = material_homogenizationID(ce)
call thermal_partition(ce)
do co = 1, homogenization_Nconstituents(ho)
ph = material_phaseID(co,ce)
if (.not. phase_thermal_constitutive(dt,ph,material_phaseMemberAt(co,ip,el))) then
if (.not. terminallyIll) & ! so first signals terminally ill...
print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
terminallyIll = .true. ! ...and kills all others
if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
endif
enddo
enddo
enddo
!$OMP END PARALLEL DO
end subroutine materialpoint_stressAndItsTangent3
end subroutine homogenization_thermal_response
!--------------------------------------------------------------------------------------------------
!> @brief
!--------------------------------------------------------------------------------------------------
subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_execElem)
subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolving_execElem)
real(pReal), intent(in) :: dt !< time increment
real(pReal), intent(in) :: Delta_t !< time increment
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
integer :: &
NiterationMPstate, &
ip, & !< integration point number
el, & !< element number
co, ce, ho, en, ph
logical :: &
converged
logical, dimension(2) :: &
doneAndHappy
co, ce, ho
!$OMP PARALLEL DO PRIVATE(ho,ce)
@ -357,13 +343,13 @@ subroutine materialpoint_stressAndItsTangent2(dt,FEsolving_execIP,FEsolving_exec
do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el)
enddo
call mechanical_homogenize(dt,ce)
call mechanical_homogenize(Delta_t,ce)
enddo IpLooping3
enddo elementLooping3
!$OMP END PARALLEL DO
end subroutine materialpoint_stressAndItsTangent2
end subroutine homogenization_mechanical_response2
!--------------------------------------------------------------------------------------------------

8
src/homogenization_mechanical.f90
View File

@ -123,21 +123,21 @@ end subroutine mechanical_partition
!--------------------------------------------------------------------------------------------------
!> @brief Average P and dPdF from the individual constituents.
!--------------------------------------------------------------------------------------------------
module subroutine mechanical_homogenize(dt,ce)
module subroutine mechanical_homogenize(Delta_t,ce)
real(pReal), intent(in) :: dt
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: ce
integer :: co
homogenization_P(1:3,1:3,ce) = phase_P(1,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) &
+ phase_P(co,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
+ phase_mechanical_dPdF(dt,co,ce)
+ phase_mechanical_dPdF(Delta_t,co,ce)
enddo
homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) &

4
src/mesh/FEM_utilities.f90
View File

@ -162,9 +162,9 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
print'(/,a)', ' ... evaluating constitutive response ......................................'
call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field
call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field
if (.not. terminallyIll) &
call materialpoint_stressAndItsTangent2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])
call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])
cutBack = .false.
P_av = sum(homogenization_P,dim=3) * wgt

14
src/phase.f90
View File

@ -113,8 +113,8 @@ module phase
end subroutine damage_restore
module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
real(pReal), intent(in) :: dt
module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
co, & !< counter in constituent loop
ce
@ -221,21 +221,21 @@ module phase
end function phase_thermal_constitutive
module function integrateDamageState(dt,co,ce) result(broken)
real(pReal), intent(in) :: dt
module function integrateDamageState(Delta_t,co,ce) result(broken)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
ce, &
co
logical :: broken
end function integrateDamageState
module function phase_mechanical_constitutive(dt,co,ip,el) result(converged_)
real(pReal), intent(in) :: dt
module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co, ip, el
logical :: converged_
end function phase_mechanical_constitutive
!ToDo: Try to merge the all stiffness functions
!ToDo: Merge all the stiffness functions
module function phase_homogenizedC(ph,en) result(C)
integer, intent(in) :: ph, en
real(pReal), dimension(6,6) :: C

8
src/phase_damage.f90
View File

@ -197,9 +197,9 @@ end function phase_f_phi
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
module function integrateDamageState(dt,co,ce) result(broken)
module function integrateDamageState(Delta_t,co,ce) result(broken)
real(pReal), intent(in) :: dt
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
ce, &
co
@ -232,7 +232,7 @@ module function integrateDamageState(dt,co,ce) result(broken)
size_so = damageState(ph)%sizeDotState
damageState(ph)%state(1:size_so,me) = damageState(ph)%state0 (1:size_so,me) &
+ damageState(ph)%dotState(1:size_so,me) * dt
+ damageState(ph)%dotState(1:size_so,me) * Delta_t
source_dotState(1:size_so,2) = 0.0_pReal
iteration: do NiterationState = 1, num%nState
@ -249,7 +249,7 @@ module function integrateDamageState(dt,co,ce) result(broken)
+ source_dotState(1:size_so,1)* (1.0_pReal - zeta)
r(1:size_so) = damageState(ph)%state (1:size_so,me) &
- damageState(ph)%State0 (1:size_so,me) &
- damageState(ph)%dotState(1:size_so,me) * dt
- damageState(ph)%dotState(1:size_so,me) * Delta_t
damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so)
converged_ = converged_ .and. converged(r(1:size_so), &
damageState(ph)%state(1:size_so,me), &

18
src/phase_damage_anisobrittle.f90
View File

@ -40,7 +40,7 @@ module function anisobrittle_init() result(mySources)
phase, &
sources, &
src
integer :: Nmembers,p
integer :: Nmembers,ph
integer, dimension(:), allocatable :: N_cl
character(len=pStringLen) :: extmsg = ''
@ -56,12 +56,12 @@ module function anisobrittle_init() result(mySources)
allocate(param(phases%length))
do p = 1, phases%length
if(mySources(p)) then
phase => phases%get(p)
do ph = 1, phases%length
if(mySources(ph)) then
phase => phases%get(ph)
sources => phase%get('damage')
associate(prm => param(p))
associate(prm => param(ph))
src => sources%get(1)
N_cl = src%get_as1dInt('N_cl',defaultVal=emptyIntArray)
@ -92,10 +92,10 @@ module function anisobrittle_init() result(mySources)
if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
Nmembers = count(material_phaseID==p)
call phase_allocateState(damageState(p),Nmembers,1,1,0)
damageState(p)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
Nmembers = count(material_phaseID==ph)
call phase_allocateState(damageState(ph),Nmembers,1,1,0)
damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
end associate

27
src/phase_mechanical.f90
View File

@ -44,9 +44,9 @@ submodule(phase) mechanical
interface
module subroutine eigendeformation_init(phases)
module subroutine eigen_init(phases)
class(tNode), pointer :: phases
end subroutine eigendeformation_init
end subroutine eigen_init
module subroutine elastic_init(phases)
class(tNode), pointer :: phases
@ -302,7 +302,7 @@ module subroutine mechanical_init(materials,phases)
end select
call eigendeformation_init(phases)
call eigen_init(phases)
end subroutine mechanical_init
@ -976,9 +976,9 @@ end subroutine mechanical_forward
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
module function phase_mechanical_constitutive(dt,co,ip,el) result(converged_)
module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
real(pReal), intent(in) :: dt
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
co, &
ip, &
@ -1054,12 +1054,12 @@ module function phase_mechanical_constitutive(dt,co,ip,el) result(converged_)
if (todo) then
subF = subF0 &
+ subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el)
endif
enddo cutbackLooping
converged_ = converged_ .and. .not. integrateDamageState(dt,co,(el-1)*discretization_nIPs + ip)
converged_ = converged_ .and. .not. integrateDamageState(Delta_t,co,(el-1)*discretization_nIPs + ip)
end function phase_mechanical_constitutive
@ -1094,12 +1094,13 @@ module subroutine mechanical_restore(ce,includeL)
end subroutine mechanical_restore
!--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF).
!--------------------------------------------------------------------------------------------------
module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
real(pReal), intent(in) :: dt
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
co, & !< counter in constituent loop
ce
@ -1149,11 +1150,11 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
do o=1,3; do p=1,3
lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
+ matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * dt
+ matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * Delta_t
lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
+ invFi*invFi(p,o)
rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
- matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * dt
- matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * Delta_t
enddo; enddo
call math_invert(temp_99,error,math_3333to99(lhs_3333))
if (error) then
@ -1182,7 +1183,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
+ matmul(temp_33_3,dLidS(1:3,1:3,p,o))
enddo; enddo
lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * dt &
lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * Delta_t &
+ math_mul3333xx3333(dSdFi,dFidS)
call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
@ -1198,7 +1199,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
! calculate dFpinvdF
temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
do o=1,3; do p=1,3
dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * dt
dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * Delta_t
enddo; enddo
!--------------------------------------------------------------------------------------------------

4
src/phase_mechanical_eigen.f90
View File

@ -32,7 +32,7 @@ submodule(phase:mechanical) eigen
contains
module subroutine eigendeformation_init(phases)
module subroutine eigen_init(phases)
class(tNode), pointer :: &
phases
@ -68,7 +68,7 @@ module subroutine eigendeformation_init(phases)
where(damage_anisobrittle_init()) model_damage = KINEMATICS_cleavage_opening_ID
end subroutine eigendeformation_init
end subroutine eigen_init
!--------------------------------------------------------------------------------------------------