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

Merge branch 'new-structure-pt2' into development

This commit is contained in:
Martin Diehl 2021-01-08 13:18:31 +01:00
parent f9e69fa60b 7595fb7308
commit bf81131578
10 changed files with 778 additions and 604 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 45ef93dbfa3e0e6fa830914b3632e188c308a099 Subproject commit 591964dcf8521d95f6cccbfe840d462c430e63d9

30
src/CPFEM2.f90
View File

@ -74,9 +74,23 @@ end subroutine CPFEM_initAll
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_init subroutine CPFEM_init
integer(HID_T) :: fileHandle
character(len=pStringLen) :: fileName
print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
if (interface_restartInc > 0) call crystallite_restartRead
if (interface_restartInc > 0) then
print'(/,a,i0,a)', ' reading restart information of increment from file'; flush(IO_STDOUT)
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName)
call homogenization_restartRead(fileHandle)
call constitutive_restartRead(fileHandle)
call HDF5_closeFile(fileHandle)
endif
end subroutine CPFEM_init end subroutine CPFEM_init
@ -86,7 +100,19 @@ end subroutine CPFEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_restartWrite subroutine CPFEM_restartWrite
call crystallite_restartWrite integer(HID_T) :: fileHandle
character(len=pStringLen) :: fileName
print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT)
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'a')
call homogenization_restartWrite(fileHandle)
call constitutive_restartWrite(fileHandle)
call HDF5_closeFile(fileHandle)
end subroutine CPFEM_restartWrite end subroutine CPFEM_restartWrite

442
src/constitutive.f90
View File

@ -15,7 +15,6 @@ module constitutive
use discretization use discretization
use parallelization use parallelization
use HDF5_utilities use HDF5_utilities
use DAMASK_interface
use results use results
implicit none implicit none
@ -42,42 +41,14 @@ module constitutive
KINEMATICS_SLIPPLANE_OPENING_ID, & KINEMATICS_SLIPPLANE_OPENING_ID, &
KINEMATICS_THERMAL_EXPANSION_ID KINEMATICS_THERMAL_EXPANSION_ID
end enum end enum
real(pReal), dimension(:,:,:), allocatable :: &
crystallite_subdt !< substepped time increment of each grain
type(rotation), dimension(:,:,:), allocatable :: & type(rotation), dimension(:,:,:), allocatable :: &
crystallite_orientation !< current orientation crystallite_orientation !< current orientation
real(pReal), dimension(:,:,:,:,:), allocatable :: &
crystallite_F0, & !< def grad at start of FE inc
crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_subFp0,& !< plastic def grad at start of crystallite inc
crystallite_subFi0,& !< intermediate def grad at start of crystallite inc
crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
crystallite_partitionedLp0, & !< plastic velocity grad at start of homog inc
crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_partitionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc
real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
crystallite_P, & !< 1st Piola-Kirchhoff stress per grain
crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step)
crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
crystallite_partitionedF0, & !< def grad at start of homog inc
crystallite_F !< def grad to be reached at end of homog inc
type :: tTensorContainer type :: tTensorContainer
real(pReal), dimension(:,:,:), allocatable :: data real(pReal), dimension(:,:,:), allocatable :: data
end type end type
type(tTensorContainer), dimension(:), allocatable :: &
constitutive_mech_Fi, &
constitutive_mech_Fi0, &
constitutive_mech_partitionedFi0, &
constitutive_mech_Li, &
constitutive_mech_Li0, &
constitutive_mech_partitionedLi0, &
constitutive_mech_Fp, &
constitutive_mech_Fp0, &
constitutive_mech_partitionedFp0
type :: tNumerics type :: tNumerics
integer :: & integer :: &
iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp
@ -162,10 +133,6 @@ module constitutive
end subroutine damage_results end subroutine damage_results
module subroutine mech_restart_read(fileHandle)
integer(HID_T), intent(in) :: fileHandle
end subroutine mech_restart_read
module subroutine mech_initializeRestorationPoints(ph,me) module subroutine mech_initializeRestorationPoints(ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
end subroutine mech_initializeRestorationPoints end subroutine mech_initializeRestorationPoints
@ -185,6 +152,61 @@ module constitutive
includeL includeL
end subroutine mech_restore end subroutine mech_restore
module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
co, & !< counter in constituent loop
ip, & !< counter in integration point loop
el !< counter in element loop
real(pReal), dimension(3,3,3,3) :: dPdF
end function constitutive_mech_dPdF
module subroutine mech_restartWrite(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph
end subroutine mech_restartWrite
module subroutine mech_restartRead(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph
end subroutine mech_restartRead
module function constitutive_mech_getS(co,ip,el) result(S)
integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: S
end function constitutive_mech_getS
module function constitutive_mech_getLp(co,ip,el) result(Lp)
integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: Lp
end function constitutive_mech_getLp
module function constitutive_mech_getF(co,ip,el) result(F)
integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: F
end function constitutive_mech_getF
module function constitutive_mech_getF_e(co,ip,el) result(F_e)
integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: F_e
end function constitutive_mech_getF_e
module function constitutive_mech_getP(co,ip,el) result(P)
integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: P
end function constitutive_mech_getP
module function constitutive_thermal_T(co,ip,el) result(T)
integer, intent(in) :: co, ip, el
real(pReal) :: T
end function constitutive_thermal_T
module subroutine constitutive_mech_setF(F,co,ip,el)
real(pReal), dimension(3,3), intent(in) :: F
integer, intent(in) :: co, ip, el
end subroutine constitutive_mech_setF
! == cleaned:end =================================================================================== ! == cleaned:end ===================================================================================
module function crystallite_stress(dt,co,ip,el) result(converged_) module function crystallite_stress(dt,co,ip,el) result(converged_)
@ -238,15 +260,12 @@ module constitutive
dPhiDot_dPhi dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents end subroutine constitutive_damage_getRateAndItsTangents
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el) module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
T T
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
S, & !< current 2nd Piola Kitchoff stress vector
Lp !< plastic velocity gradient
real(pReal), intent(inout) :: & real(pReal), intent(inout) :: &
TDot, & TDot, &
dTDot_dT dTDot_dT
@ -342,13 +361,11 @@ module constitutive
end subroutine constitutive_plastic_LpAndItsTangents end subroutine constitutive_plastic_LpAndItsTangents
module subroutine constitutive_plastic_dependentState(F, co, ip, el) module subroutine constitutive_plastic_dependentState(co,ip,el)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point co, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(3,3) :: &
F !< elastic deformation gradient
end subroutine constitutive_plastic_dependentState end subroutine constitutive_plastic_dependentState
@ -390,12 +407,17 @@ module constitutive
converged, & converged, &
crystallite_init, & crystallite_init, &
crystallite_stress, & crystallite_stress, &
crystallite_stressTangent, & constitutive_mech_dPdF, &
crystallite_orientations, & crystallite_orientations, &
crystallite_push33ToRef, & crystallite_push33ToRef, &
crystallite_restartWrite, & constitutive_restartWrite, &
constitutive_restartRead, &
integrateSourceState, & integrateSourceState, &
crystallite_restartRead, & constitutive_mech_setF, &
constitutive_mech_getP, &
constitutive_mech_getLp, &
constitutive_mech_getF, &
constitutive_mech_getS, &
constitutive_initializeRestorationPoints, & constitutive_initializeRestorationPoints, &
constitutive_windForward, & constitutive_windForward, &
PLASTICITY_UNDEFINED_ID, & PLASTICITY_UNDEFINED_ID, &
@ -615,7 +637,7 @@ end subroutine constitutive_LiAndItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_damage_collectDotState(S, co, ip, el,ph,of) result(broken) function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID of integration point co, & !< component-ID of integration point
@ -623,8 +645,6 @@ function constitutive_damage_collectDotState(S, co, ip, el,ph,of) result(broken)
el, & !< element el, & !< element
ph, & ph, &
of of
real(pReal), intent(in), dimension(3,3) :: &
S !< 2nd Piola Kirchhoff stress (vector notation)
integer :: & integer :: &
so !< counter in source loop so !< counter in source loop
logical :: broken logical :: broken
@ -637,7 +657,7 @@ function constitutive_damage_collectDotState(S, co, ip, el,ph,of) result(broken)
sourceType: select case (phase_source(so,ph)) sourceType: select case (phase_source(so,ph))
case (SOURCE_damage_anisoBrittle_ID) sourceType case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState(S, co, ip, el) ! correct stress? call source_damage_anisoBrittle_dotState(constitutive_mech_getS(co,ip,el), co, ip, el) ! correct stress?
case (SOURCE_damage_isoDuctile_ID) sourceType case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_dotState(co, ip, el) call source_damage_isoDuctile_dotState(co, ip, el)
@ -748,7 +768,6 @@ subroutine constitutive_allocateState(state, &
allocate(state%atol (sizeState), source=0.0_pReal) allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal) allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal) allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal)
allocate(state%subState0 (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%state (sizeState,Nconstituents), source=0.0_pReal) allocate(state%state (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal) allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal)
@ -792,17 +811,14 @@ end subroutine constitutive_restore
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_forward subroutine constitutive_forward
integer :: i, j integer :: ph, so
crystallite_F0 = crystallite_F
crystallite_Lp0 = crystallite_Lp
crystallite_S0 = crystallite_S
call constitutive_mech_forward() call constitutive_mech_forward()
do i = 1, size(sourceState) do ph = 1, size(sourceState)
do j = 1,phase_Nsources(i) do so = 1,phase_Nsources(ph)
sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state sourceState(ph)%p(so)%state0 = sourceState(ph)%p(so)%state
enddo; enddo enddo; enddo
end subroutine constitutive_forward end subroutine constitutive_forward
@ -838,12 +854,12 @@ end subroutine constitutive_results
subroutine crystallite_init subroutine crystallite_init
integer :: & integer :: &
Nconstituents, &
ph, & ph, &
me, & me, &
co, & !< counter in integration point component loop co, & !< counter in integration point component loop
ip, & !< counter in integration point loop ip, & !< counter in integration point loop
el, & !< counter in element loop el, & !< counter in element loop
so, &
cMax, & !< maximum number of integration point components cMax, & !< maximum number of integration point components
iMax, & !< maximum number of integration points iMax, & !< maximum number of integration points
eMax !< maximum number of elements eMax !< maximum number of elements
@ -866,19 +882,6 @@ subroutine crystallite_init
iMax = discretization_nIPs iMax = discretization_nIPs
eMax = discretization_Nelems eMax = discretization_Nelems
allocate(crystallite_F(3,3,cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_S0, &
crystallite_F0,crystallite_Lp0, &
crystallite_partitionedS0, &
crystallite_partitionedF0,&
crystallite_partitionedLp0, &
crystallite_S,crystallite_P, &
crystallite_Fe,crystallite_Lp, &
crystallite_subFp0,crystallite_subFi0, &
source = crystallite_F)
allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_orientation(cMax,iMax,eMax)) allocate(crystallite_orientation(cMax,iMax,eMax))
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -914,27 +917,10 @@ subroutine crystallite_init
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(constitutive_mech_Fi(phases%length))
allocate(constitutive_mech_Fi0(phases%length))
allocate(constitutive_mech_partitionedFi0(phases%length))
allocate(constitutive_mech_Fp(phases%length))
allocate(constitutive_mech_Fp0(phases%length))
allocate(constitutive_mech_partitionedFp0(phases%length))
allocate(constitutive_mech_Li(phases%length))
allocate(constitutive_mech_Li0(phases%length))
allocate(constitutive_mech_partitionedLi0(phases%length))
do ph = 1, phases%length do ph = 1, phases%length
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs do so = 1, phase_Nsources(ph)
allocate(sourceState(ph)%p(so)%subState0,source=sourceState(ph)%p(so)%state0) ! ToDo: hack
allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents)) enddo
allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_partitionedFi0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_partitionedFp0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_partitionedLi0(ph)%data(3,3,Nconstituents))
enddo enddo
print'(a42,1x,i10)', ' # of elements: ', eMax print'(a42,1x,i10)', ' # of elements: ', eMax
@ -942,34 +928,6 @@ subroutine crystallite_init
print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax
flush(IO_STDOUT) flush(IO_STDOUT)
!$OMP PARALLEL DO PRIVATE(ph,me)
do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2)
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) &
/ math_det33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal)
constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = math_I3
crystallite_F0(1:3,1:3,co,ip,el) = math_I3
crystallite_Fe(1:3,1:3,co,ip,el) = math_inv33(matmul(constitutive_mech_Fi0(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
enddo
enddo; enddo
!$OMP END PARALLEL DO
crystallite_partitionedF0 = crystallite_F0
crystallite_F = crystallite_F0
!$OMP PARALLEL DO PRIVATE(ph,me) !$OMP PARALLEL DO PRIVATE(ph,me)
do el = 1, size(material_phaseMemberAt,3) do el = 1, size(material_phaseMemberAt,3)
@ -978,7 +936,7 @@ subroutine crystallite_init
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
call crystallite_orientations(co,ip,el) call crystallite_orientations(co,ip,el)
call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el),co,ip,el) ! update dependent state variables to be consistent with basic states call constitutive_plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
enddo enddo
enddo enddo
enddo enddo
@ -1003,15 +961,11 @@ subroutine constitutive_initializeRestorationPoints(ip,el)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp0(1:3,1:3,co,ip,el)
crystallite_partitionedF0(1:3,1:3,co,ip,el) = crystallite_F0(1:3,1:3,co,ip,el)
crystallite_partitionedS0(1:3,1:3,co,ip,el) = crystallite_S0(1:3,1:3,co,ip,el)
call mech_initializeRestorationPoints(ph,me) call mech_initializeRestorationPoints(ph,me)
do so = 1, phase_Nsources(material_phaseAt(co,el)) do so = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) = & sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state0(:,me)
sourceState(material_phaseAt(co,el))%p(so)%state0( :,material_phasememberAt(co,ip,el))
enddo enddo
enddo enddo
@ -1035,9 +989,6 @@ subroutine constitutive_windForward(ip,el)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
crystallite_partitionedF0 (1:3,1:3,co,ip,el) = crystallite_F (1:3,1:3,co,ip,el)
crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp(1:3,1:3,co,ip,el)
crystallite_partitionedS0 (1:3,1:3,co,ip,el) = crystallite_S (1:3,1:3,co,ip,el)
call constitutive_mech_windForward(ph,me) call constitutive_mech_windForward(ph,me)
do so = 1, phase_Nsources(material_phaseAt(co,el)) do so = 1, phase_Nsources(material_phaseAt(co,el))
@ -1048,136 +999,6 @@ subroutine constitutive_windForward(ip,el)
end subroutine constitutive_windForward end subroutine constitutive_windForward
!--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF).
!--------------------------------------------------------------------------------------------------
function crystallite_stressTangent(co,ip,el) result(dPdF)
real(pReal), dimension(3,3,3,3) :: dPdF
integer, intent(in) :: &
co, & !< counter in constituent loop
ip, & !< counter in integration point loop
el !< counter in element loop
integer :: &
o, &
p, ph, me
real(pReal), dimension(3,3) :: devNull, &
invSubFp0,invSubFi0,invFp,invFi, &
temp_33_1, temp_33_2, temp_33_3
real(pReal), dimension(3,3,3,3) :: dSdFe, &
dSdF, &
dSdFi, &
dLidS, & ! tangent in lattice configuration
dLidFi, &
dLpdS, &
dLpdFi, &
dFidS, &
dFpinvdF, &
rhs_3333, &
lhs_3333, &
temp_3333
real(pReal), dimension(9,9):: temp_99
logical :: error
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
crystallite_Fe(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el)
call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, &
crystallite_S (1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
co,ip,el)
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me))
invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el))
invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el))
if (sum(abs(dLidS)) < tol_math_check) then
dFidS = 0.0_pReal
else
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) &
+ crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p))
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) &
- crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
enddo; enddo
call math_invert(temp_99,error,math_3333to99(lhs_3333))
if (error) then
call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
ext_msg='inversion error in analytic tangent calculation')
dFidS = 0.0_pReal
else
dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
endif
dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
endif
call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
crystallite_S (1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
!--------------------------------------------------------------------------------------------------
! calculate dSdF
temp_33_1 = transpose(matmul(invFp,invFi))
temp_33_2 = matmul(crystallite_F(1:3,1:3,co,ip,el),invSubFp0)
temp_33_3 = matmul(matmul(crystallite_F(1:3,1:3,co,ip,el),invFp), invSubFi0)
do o=1,3; do p=1,3
rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
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 = crystallite_subdt(co,ip,el)*math_mul3333xx3333(dSdFe,temp_3333) &
+ math_mul3333xx3333(dSdFi,dFidS)
call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
if (error) then
call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
ext_msg='inversion error in analytic tangent calculation')
dSdF = rhs_3333
else
dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
endif
!--------------------------------------------------------------------------------------------------
! calculate dFpinvdF
temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
do o=1,3; do p=1,3
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(co,ip,el) &
* matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi))
enddo; enddo
!--------------------------------------------------------------------------------------------------
! assemble dPdF
temp_33_1 = matmul(crystallite_S(1:3,1:3,co,ip,el),transpose(invFp))
temp_33_2 = matmul(crystallite_F(1:3,1:3,co,ip,el),invFp)
temp_33_3 = matmul(temp_33_2,crystallite_S(1:3,1:3,co,ip,el))
dPdF = 0.0_pReal
do p=1,3
dPdF(p,1:3,p,1:3) = transpose(matmul(invFp,temp_33_1))
enddo
do o=1,3; do p=1,3
dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
+ matmul(matmul(crystallite_F(1:3,1:3,co,ip,el), &
dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)), &
transpose(invFp)) &
+ matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo
end function crystallite_stressTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates orientations !> @brief calculates orientations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1189,7 +1010,8 @@ subroutine crystallite_orientations(co,ip,el)
el !< counter in element loop el !< counter in element loop
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el)))) call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(&
constitutive_mech_getF_e(co,ip,el))))
if (plasticState(material_phaseAt(1,el))%nonlocal) & if (plasticState(material_phaseAt(1,el))%nonlocal) &
call plastic_nonlocal_updateCompatibility(crystallite_orientation, & call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
@ -1205,17 +1027,17 @@ end subroutine crystallite_orientations
function crystallite_push33ToRef(co,ip,el, tensor33) function crystallite_push33ToRef(co,ip,el, tensor33)
real(pReal), dimension(3,3), intent(in) :: tensor33 real(pReal), dimension(3,3), intent(in) :: tensor33
real(pReal), dimension(3,3) :: T
integer, intent(in):: & integer, intent(in):: &
el, & el, &
ip, & ip, &
co co
real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3) :: crystallite_push33ToRef
real(pReal), dimension(3,3) :: T
T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(constitutive_mech_getF(co,ip,el)))) ! ToDo: initial orientation correct?
T = matmul(material_orientation0(co,ip,el)%asMatrix(), & ! ToDo: initial orientation correct?
transpose(math_inv33(crystallite_F(1:3,1:3,co,ip,el))))
crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
end function crystallite_push33ToRef end function crystallite_push33ToRef
@ -1254,7 +1076,7 @@ function integrateSourceState(dt,co,ip,el) result(broken)
converged_ = .true. converged_ = .true.
broken = constitutive_thermal_collectDotState(ph,me) broken = constitutive_thermal_collectDotState(ph,me)
broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me) broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) return if(broken) return
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
@ -1272,7 +1094,7 @@ function integrateSourceState(dt,co,ip,el) result(broken)
enddo enddo
broken = constitutive_thermal_collectDotState(ph,me) broken = constitutive_thermal_collectDotState(ph,me)
broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me) broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) exit iteration if(broken) exit iteration
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
@ -1292,7 +1114,7 @@ function integrateSourceState(dt,co,ip,el) result(broken)
enddo enddo
if(converged_) then if(converged_) then
broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me) broken = constitutive_damage_deltaState(constitutive_mech_getF_e(co,ip,el),co,ip,el,ph,me)
exit iteration exit iteration
endif endif
@ -1345,90 +1167,60 @@ end function converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current restart information (Field and constitutive data) to file. !> @brief Write current restart information (Field and constitutive data) to file.
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively ! ToDo: Merge data into one file for MPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_restartWrite subroutine constitutive_restartWrite(fileHandle)
integer(HID_T), intent(in) :: fileHandle
integer(HID_T), dimension(2) :: groupHandle
integer :: ph integer :: ph
integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: fileName, datasetName
print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT)
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' groupHandle(1) = HDF5_addGroup(fileHandle,'phase')
fileHandle = HDF5_openFile(fileName,'a')
call HDF5_write(fileHandle,crystallite_F,'F') do ph = 1, size(material_name_phase)
call HDF5_write(fileHandle,crystallite_Lp, 'L_p')
call HDF5_write(fileHandle,crystallite_S, 'S') groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_phase(ph))
call mech_restartWrite(groupHandle(2),ph)
call HDF5_closeGroup(groupHandle(2))
groupHandle = HDF5_addGroup(fileHandle,'phase')
do ph = 1,size(material_name_phase)
write(datasetName,'(i0,a)') ph,'_omega'
call HDF5_write(groupHandle,plasticState(ph)%state,datasetName)
write(datasetName,'(i0,a)') ph,'_F_i'
call HDF5_write(groupHandle,constitutive_mech_Fi(ph)%data,datasetName)
write(datasetName,'(i0,a)') ph,'_L_i'
call HDF5_write(groupHandle,constitutive_mech_Li(ph)%data,datasetName)
write(datasetName,'(i0,a)') ph,'_F_p'
call HDF5_write(groupHandle,constitutive_mech_Fp(ph)%data,datasetName)
enddo enddo
call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_addGroup(fileHandle,'homogenization') call HDF5_closeGroup(groupHandle(1))
do ph = 1, size(material_name_homogenization)
write(datasetName,'(i0,a)') ph,'_omega'
call HDF5_write(groupHandle,homogState(ph)%state,datasetName)
enddo
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle) end subroutine constitutive_restartWrite
end subroutine crystallite_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Read data for restart !> @brief Read data for restart
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively ! ToDo: Merge data into one file for MPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_restartRead subroutine constitutive_restartRead(fileHandle)
integer(HID_T), intent(in) :: fileHandle
integer(HID_T), dimension(2) :: groupHandle
integer :: ph integer :: ph
integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: fileName, datasetName
print'(/,a,i0,a)', ' reading restart information of increment from file'
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5' groupHandle(1) = HDF5_openGroup(fileHandle,'phase')
fileHandle = HDF5_openFile(fileName)
call HDF5_read(fileHandle,crystallite_F0, 'F') do ph = 1, size(material_name_phase)
call HDF5_read(fileHandle,crystallite_Lp0,'L_p')
call HDF5_read(fileHandle,crystallite_S0, 'S') groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_phase(ph))
call mech_restartRead(groupHandle(2),ph)
call HDF5_closeGroup(groupHandle(2))
groupHandle = HDF5_openGroup(fileHandle,'phase')
do ph = 1,size(material_name_phase)
write(datasetName,'(i0,a)') ph,'_omega'
call HDF5_read(groupHandle,plasticState(ph)%state0,datasetName)
write(datasetName,'(i0,a)') ph,'_F_i'
call HDF5_read(groupHandle,constitutive_mech_Fi0(ph)%data,datasetName)
write(datasetName,'(i0,a)') ph,'_L_i'
call HDF5_read(groupHandle,constitutive_mech_Li0(ph)%data,datasetName)
write(datasetName,'(i0,a)') ph,'_F_p'
call HDF5_read(groupHandle,constitutive_mech_Fp0(ph)%data,datasetName)
enddo enddo
call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_openGroup(fileHandle,'homogenization') call HDF5_closeGroup(groupHandle(1))
do ph = 1,size(material_name_homogenization)
write(datasetName,'(i0,a)') ph,'_omega'
call HDF5_read(groupHandle,homogState(ph)%state0,datasetName)
enddo
call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle) end subroutine constitutive_restartRead
end subroutine crystallite_restartRead
end module constitutive end module constitutive

641
src/constitutive_mech.f90
View File

File diff suppressed because it is too large Load Diff

20
src/constitutive_plastic_nonlocal.f90
View File

@ -552,10 +552,8 @@ end function plastic_nonlocal_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure !> @brief calculates quantities characterizing the microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el) module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
real(pReal), dimension(3,3), intent(in) :: &
F
integer, intent(in) :: & integer, intent(in) :: &
instance, & instance, &
of, & of, &
@ -647,7 +645,7 @@ module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el)
ph = material_phaseAt(1,el) ph = material_phaseAt(1,el)
me = material_phaseMemberAt(1,ip,el) me = material_phaseMemberAt(1,ip,el)
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
invFe = matmul(constitutive_mech_Fp(ph)%data(1:3,1:3,me),math_inv33(F)) invFe = math_inv33(constitutive_mech_Fe(ph)%data(1:3,1:3,me))
rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
@ -976,13 +974,11 @@ end subroutine plastic_nonlocal_deltaState
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, & module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
instance,of,ip,el) instance,of,ip,el)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
F !< Deformation gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature, & !< temperature Temperature, & !< temperature
timestep !< substepped crystallite time increment timestep !< substepped crystallite time increment
@ -1149,7 +1145,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = rhoDotFlux(F,timestep, instance,of,ip,el) & rhoDot = rhoDotFlux(timestep, instance,of,ip,el) &
+ rhoDotMultiplication & + rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide & + rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation & + rhoDotAthermalAnnihilation &
@ -1178,10 +1174,8 @@ end subroutine plastic_nonlocal_dotState
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
function rhoDotFlux(F,timestep, instance,of,ip,el) function rhoDotFlux(timestep,instance,of,ip,el)
real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
F !< Deformation gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment timestep !< substepped crystallite time increment
integer, intent(in) :: & integer, intent(in) :: &
@ -1293,7 +1287,7 @@ function rhoDotFlux(F,timestep, instance,of,ip,el)
m(1:3,:,3) = -prm%slip_transverse m(1:3,:,3) = -prm%slip_transverse
m(1:3,:,4) = prm%slip_transverse m(1:3,:,4) = prm%slip_transverse
my_F = F(1:3,1:3,1,ip,el) my_F = constitutive_mech_F(ph)%data(1:3,1:3,of)
my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of))) my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of)))
neighbors: do n = 1,nIPneighbors neighbors: do n = 1,nIPneighbors
@ -1311,7 +1305,7 @@ function rhoDotFlux(F,timestep, instance,of,ip,el)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el) neighbor_F = constitutive_mech_F(np)%data(1:3,1:3,no)
neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no))) neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no)))
Favg = 0.5_pReal * (my_F + neighbor_F) Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average else ! if no neighbor, take my value as average

39
src/constitutive_thermal.f90
View File

@ -9,7 +9,7 @@ submodule(constitutive) constitutive_thermal
integer, intent(in) :: source_length integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: mySources
end function source_thermal_dissipation_init end function source_thermal_dissipation_init
module function source_thermal_externalheat_init(source_length) result(mySources) module function source_thermal_externalheat_init(source_length) result(mySources)
integer, intent(in) :: source_length integer, intent(in) :: source_length
logical, dimension(:,:), allocatable :: mySources logical, dimension(:,:), allocatable :: mySources
@ -55,8 +55,8 @@ module subroutine thermal_init
if(maxval(phase_Nsources) /= 0) then if(maxval(phase_Nsources) /= 0) then
where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID
where(source_thermal_externalheat_init(maxval(phase_Nsources))) phase_source = SOURCE_thermal_externalheat_ID where(source_thermal_externalheat_init(maxval(phase_Nsources))) phase_source = SOURCE_thermal_externalheat_ID
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!initialize kinematic mechanisms !initialize kinematic mechanisms
if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) & if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) &
@ -68,15 +68,13 @@ end subroutine thermal_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief calculates thermal dissipation rate !< @brief calculates thermal dissipation rate
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el) module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
T T !< plastic velocity gradient
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
S, & !< current 2nd Piola Kirchhoff stress
Lp !< plastic velocity gradient
real(pReal), intent(inout) :: & real(pReal), intent(inout) :: &
TDot, & TDot, &
dTDot_dT dTDot_dT
@ -84,6 +82,7 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T,
real(pReal) :: & real(pReal) :: &
my_Tdot, & my_Tdot, &
my_dTdot_dT my_dTdot_dT
real(pReal), dimension(3,3) :: Lp, S
integer :: & integer :: &
phase, & phase, &
homog, & homog, &
@ -101,10 +100,10 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T,
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase)) select case(phase_source(source,phase))
case (SOURCE_thermal_dissipation_ID) case (SOURCE_thermal_dissipation_ID)
Lp = constitutive_mech_getLp(grain,ip,el)
S = constitutive_mech_getS(grain,ip,el)
call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
S(1:3,1:3,grain,ip,el), & S, Lp, phase)
Lp(1:3,1:3,grain,ip,el), &
phase)
case (SOURCE_thermal_externalheat_ID) case (SOURCE_thermal_externalheat_ID)
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
@ -122,4 +121,22 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T,
end subroutine constitutive_thermal_getRateAndItsTangents end subroutine constitutive_thermal_getRateAndItsTangents
! getter for non-thermal (e.g. mech)
module function constitutive_thermal_T(co,ip,el) result(T)
integer, intent(in) :: co, ip, el
real(pReal) :: T
integer :: ho, tme
ho = material_homogenizationAt(el)
tme = material_homogenizationMemberAt(ip,el)
T = temperature(ho)%p(tme)
end function constitutive_thermal_T
end submodule constitutive_thermal end submodule constitutive_thermal

83
src/homogenization.f90
View File

@ -16,6 +16,7 @@ module homogenization
use thermal_conduction use thermal_conduction
use damage_none use damage_none
use damage_nonlocal use damage_nonlocal
use HDF5_utilities
use results use results
implicit none implicit none
@ -63,7 +64,8 @@ module homogenization
el !< element number el !< element number
end subroutine mech_partition end subroutine mech_partition
module subroutine mech_homogenize(ip,el) module subroutine mech_homogenize(dt,ip,el)
real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
@ -91,7 +93,9 @@ module homogenization
homogenization_init, & homogenization_init, &
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
homogenization_forward, & homogenization_forward, &
homogenization_results homogenization_results, &
homogenization_restartRead, &
homogenization_restartWrite
contains contains
@ -168,10 +172,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
converged = .false. ! pretend failed step ... converged = .false. ! pretend failed step ...
subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
if (homogState(ho)%sizeState > 0) & if (homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me)
homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me) if (damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me)
if (damageState(ho)%sizeState > 0) &
damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me)
cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog) cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog)
@ -184,10 +186,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
! wind forward grain starting point ! wind forward grain starting point
call constitutive_windForward(ip,el) call constitutive_windForward(ip,el)
if(homogState(ho)%sizeState > 0) & if(homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me)
homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me) if(damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me)
if(damageState(ho)%sizeState > 0) &
damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me)
endif steppingNeeded endif steppingNeeded
elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
@ -201,10 +201,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
call constitutive_restore(ip,el,subStep < 1.0_pReal) call constitutive_restore(ip,el,subStep < 1.0_pReal)
if(homogState(ho)%sizeState > 0) & if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me) if(damageState(ho)%sizeState > 0) damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me)
if(damageState(ho)%sizeState > 0) &
damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me)
endif endif
if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
@ -257,7 +255,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
do co = 1, myNgrains do co = 1, myNgrains
call crystallite_orientations(co,ip,el) call crystallite_orientations(co,ip,el)
enddo enddo
call mech_homogenize(ip,el) call mech_homogenize(dt,ip,el)
enddo IpLooping3 enddo IpLooping3
enddo elementLooping3 enddo elementLooping3
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -320,4 +318,59 @@ subroutine homogenization_forward
end subroutine homogenization_forward end subroutine homogenization_forward
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine homogenization_restartWrite(fileHandle)
integer(HID_T), intent(in) :: fileHandle
integer(HID_T), dimension(2) :: groupHandle
integer :: ho
groupHandle(1) = HDF5_addGroup(fileHandle,'homogenization')
do ho = 1, size(material_name_homogenization)
groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_homogenization(ho))
call HDF5_read(groupHandle(2),homogState(ho)%state,'omega') ! ToDo: should be done by mech
call HDF5_closeGroup(groupHandle(2))
enddo
call HDF5_closeGroup(groupHandle(1))
end subroutine homogenization_restartWrite
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine homogenization_restartRead(fileHandle)
integer(HID_T), intent(in) :: fileHandle
integer(HID_T), dimension(2) :: groupHandle
integer :: ho
groupHandle(1) = HDF5_openGroup(fileHandle,'homogenization')
do ho = 1, size(material_name_homogenization)
groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_homogenization(ho))
call HDF5_write(groupHandle(2),homogState(ho)%state,'omega') ! ToDo: should be done by mech
call HDF5_closeGroup(groupHandle(2))
enddo
call HDF5_closeGroup(groupHandle(1))
end subroutine homogenization_restartRead
end module homogenization end module homogenization

63
src/homogenization_mech.f90
View File

@ -52,12 +52,11 @@ submodule(homogenization) homogenization_mech
end subroutine mech_RGC_averageStressAndItsTangent end subroutine mech_RGC_averageStressAndItsTangent
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy) module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: & real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses P,& !< partitioned stresses
F,& !< partitioned deformation gradients F !< partitioned deformation gradients
F0 !< partitioned initial deformation gradients
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
real(pReal), dimension(3,3), intent(in) :: avgF !< average F real(pReal), dimension(3,3), intent(in) :: avgF !< average F
real(pReal), intent(in) :: dt !< time increment real(pReal), intent(in) :: dt !< time increment
@ -113,67 +112,73 @@ module subroutine mech_partition(subF,ip,el)
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
integer :: co
real(pReal) :: F(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
crystallite_F(1:3,1:3,1,ip,el) = subF F(1:3,1:3,1) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(& call mech_isostrain_partitionDeformation(F,subF)
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(& call mech_RGC_partitionDeformation(F,subF,ip,el)
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
subF,&
ip, &
el)
end select chosenHomogenization end select chosenHomogenization
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
call constitutive_mech_setF(F(1:3,1:3,co),co,ip,el)
enddo
end subroutine mech_partition end subroutine mech_partition
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Average P and dPdF from the individual constituents. !> @brief Average P and dPdF from the individual constituents.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_homogenize(ip,el) module subroutine mech_homogenize(dt,ip,el)
real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
integer :: co,ce integer :: co,ce
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
ce = (el-1)* discretization_nIPs + ip ce = (el-1)* discretization_nIPs + ip
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el) homogenization_P(1:3,1:3,ce) = constitutive_mech_getP(1,ip,el)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(1,ip,el) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = constitutive_mech_dPdF(dt,1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el)
Ps(:,:,co) = constitutive_mech_getP(co,ip,el)
enddo enddo
call mech_isostrain_averageStressAndItsTangent(& call mech_isostrain_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ce), & homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & Ps,dPdFs, &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el)
Ps(:,:,co) = constitutive_mech_getP(co,ip,el)
enddo enddo
call mech_RGC_averageStressAndItsTangent(& call mech_RGC_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ce), & homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & Ps,dPdFs, &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
end select chosenHomogenization end select chosenHomogenization
@ -198,21 +203,17 @@ module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
integer :: co integer :: co
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el) dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(subdt,co,ip,el)
Fs(:,:,co) = constitutive_mech_getF(co,ip,el)
Ps(:,:,co) = constitutive_mech_getP(co,ip,el)
enddo enddo
doneAndHappy = & doneAndHappy = mech_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el)
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),&
subF,&
subdt, &
dPdFs, &
ip, &
el)
else else
doneAndHappy = .true. doneAndHappy = .true.
endif endif

35
src/homogenization_mech_RGC.f90
View File

@ -24,9 +24,6 @@ submodule(homogenization:homogenization_mech) homogenization_mech_RGC
end type tParameters end type tParameters
type :: tRGCstate type :: tRGCstate
real(pReal), pointer, dimension(:) :: &
work, &
penaltyEnergy
real(pReal), pointer, dimension(:,:) :: & real(pReal), pointer, dimension(:,:) :: &
relaxationVector relaxationVector
end type tRGCstate end type tRGCstate
@ -170,8 +167,7 @@ module subroutine mech_RGC_init(num_homogMech)
nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) & nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) &
+ prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) & + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) &
+ prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1)) + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
sizeState = nIntFaceTot & sizeState = nIntFaceTot
+ size(['avg constitutive work ','average penalty energy'])
homogState(h)%sizeState = sizeState homogState(h)%sizeState = sizeState
allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal) allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal)
@ -180,8 +176,6 @@ module subroutine mech_RGC_init(num_homogMech)
stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:) stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:)
st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:) st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:)
stt%work => homogState(h)%state(nIntFaceTot+1,:)
stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:)
allocate(dst%volumeDiscrepancy( Nmaterialpoints), source=0.0_pReal) allocate(dst%volumeDiscrepancy( Nmaterialpoints), source=0.0_pReal)
allocate(dst%relaxationRate_avg( Nmaterialpoints), source=0.0_pReal) allocate(dst%relaxationRate_avg( Nmaterialpoints), source=0.0_pReal)
@ -243,12 +237,11 @@ end subroutine mech_RGC_partitionDeformation
!> @brief update the internal state of the homogenization scheme and tell whether "done" and !> @brief update the internal state of the homogenization scheme and tell whether "done" and
! "happy" with result ! "happy" with result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy) module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: & real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses P,& !< partitioned stresses
F,& !< partitioned deformation gradients F !< partitioned deformation gradients
F0 !< partitioned initial deformation gradients
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
real(pReal), dimension(3,3), intent(in) :: avgF !< average F real(pReal), dimension(3,3), intent(in) :: avgF !< average F
real(pReal), intent(in) :: dt !< time increment real(pReal), intent(in) :: dt !< time increment
@ -287,8 +280,8 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
allocate(resid(3*nIntFaceTot), source=0.0_pReal) allocate(resid(3*nIntFaceTot), source=0.0_pReal)
allocate(tract(nIntFaceTot,3), source=0.0_pReal) allocate(tract(nIntFaceTot,3), source=0.0_pReal)
relax = stt%relaxationVector(:,of) relax = stt%relaxationVector(:,of)
drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of) drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of)
@ -346,17 +339,6 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
if (residMax < num%rtol*stresMax .or. residMax < num%atol) then if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
doneAndHappy = .true. doneAndHappy = .true.
!--------------------------------------------------------------------------------------------------
! compute/update the state for postResult, i.e., all energy densities computed by time-integration
do iGrain = 1,product(prm%N_constituents)
do i = 1,3;do j = 1,3
stt%work(of) = stt%work(of) &
+ P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
stt%penaltyEnergy(of) = stt%penaltyEnergy(of) &
+ R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
enddo; enddo
enddo
dst%mismatch(1:3,of) = sum(NN,2)/real(nGrain,pReal) dst%mismatch(1:3,of) = sum(NN,2)/real(nGrain,pReal)
dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
dst%relaxationRate_max(of) = maxval(abs(drelax))/dt dst%relaxationRate_max(of) = maxval(abs(drelax))/dt
@ -523,7 +505,6 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
integer, dimension (3) :: iGrain3,iGNghb3,nGDim integer, dimension (3) :: iGrain3,iGNghb3,nGDim
real(pReal), dimension (3,3) :: gDef,nDef real(pReal), dimension (3,3) :: gDef,nDef
real(pReal), dimension (3) :: nVect,surfCorr real(pReal), dimension (3) :: nVect,surfCorr
real(pReal), dimension (2) :: Gmoduli
integer :: iGrain,iGNghb,iFace,i,j,k,l integer :: iGrain,iGNghb,iFace,i,j,k,l
real(pReal) :: muGrain,muGNghb,nDefNorm real(pReal) :: muGrain,muGNghb,nDefNorm
real(pReal), parameter :: & real(pReal), parameter :: &
@ -755,15 +736,9 @@ module subroutine mech_RGC_results(instance,group)
associate(stt => state(instance), dst => dependentState(instance), prm => param(instance)) associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
outputsLoop: do o = 1,size(prm%output) outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o))) select case(trim(prm%output(o)))
case('W')
call results_writeDataset(group,stt%work,trim(prm%output(o)), &
'work density','J/m³')
case('M') case('M')
call results_writeDataset(group,dst%mismatch,trim(prm%output(o)), & call results_writeDataset(group,dst%mismatch,trim(prm%output(o)), &
'average mismatch tensor','1') 'average mismatch tensor','1')
case('R')
call results_writeDataset(group,stt%penaltyEnergy,trim(prm%output(o)), &
'mismatch penalty density','J/m³')
case('Delta_V') case('Delta_V')
call results_writeDataset(group,dst%volumeDiscrepancy,trim(prm%output(o)), & call results_writeDataset(group,dst%volumeDiscrepancy,trim(prm%output(o)), &
'volume discrepancy','m³') 'volume discrepancy','m³')

27
src/thermal_conduction.f90
View File

@ -94,7 +94,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
dTdot_dT = 0.0_pReal dTdot_dT = 0.0_pReal
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el) call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el)
Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal)
@ -112,14 +112,16 @@ function thermal_conduction_getConductivity(ip,el)
el !< element number el !< element number
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
thermal_conduction_getConductivity thermal_conduction_getConductivity
integer :: & integer :: &
grain co
thermal_conduction_getConductivity = 0.0_pReal thermal_conduction_getConductivity = 0.0_pReal
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getConductivity = thermal_conduction_getConductivity + & thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el))) crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el)))
enddo enddo
thermal_conduction_getConductivity = thermal_conduction_getConductivity & thermal_conduction_getConductivity = thermal_conduction_getConductivity &
@ -138,14 +140,16 @@ function thermal_conduction_getSpecificHeat(ip,el)
el !< element number el !< element number
real(pReal) :: & real(pReal) :: &
thermal_conduction_getSpecificHeat thermal_conduction_getSpecificHeat
integer :: & integer :: &
grain co
thermal_conduction_getSpecificHeat = 0.0_pReal thermal_conduction_getSpecificHeat = 0.0_pReal
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
+ lattice_c_p(material_phaseAt(grain,el)) + lattice_c_p(material_phaseAt(co,el))
enddo enddo
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
@ -164,15 +168,16 @@ function thermal_conduction_getMassDensity(ip,el)
el !< element number el !< element number
real(pReal) :: & real(pReal) :: &
thermal_conduction_getMassDensity thermal_conduction_getMassDensity
integer :: & integer :: &
grain co
thermal_conduction_getMassDensity = 0.0_pReal thermal_conduction_getMassDensity = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
+ lattice_rho(material_phaseAt(grain,el)) + lattice_rho(material_phaseAt(co,el))
enddo enddo
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &