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

not used 

was only used for reporting (see v.2.0.0)
This commit is contained in:
Martin Diehl 2021-01-07 14:59:12 +01:00 committed by Sharan Roongta
parent 29954763a2
commit 6f65de27fc
9 changed files with 365 additions and 362 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/FEM/polyXtal/material.yaml
View File

@ -5,8 +5,8 @@ homogenization:
phase:
Aluminum:
lattice: cF
mechanics:
lattice: cF
output: [F, P, F_e, F_p, L_p]
elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke}
plasticity:

218
src/constitutive.f90
View File

@ -48,7 +48,6 @@ module constitutive
crystallite_orientation !< current orientation
real(pReal), dimension(:,:,:,:,:), allocatable :: &
crystallite_F0, & !< def grad at start of FE inc
crystallite_subF, & !< def grad to be reached at end of crystallite 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
@ -60,9 +59,8 @@ module constitutive
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
real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
crystallite_partitionedF !< def grad to be reached at end of homog inc
crystallite_partitionedF0, & !< def grad at start of homog inc
crystallite_F !< def grad to be reached at end of homog inc
type :: tTensorContainer
real(pReal), dimension(:,:,:), allocatable :: data
@ -179,6 +177,14 @@ module constitutive
module subroutine constitutive_mech_forward
end subroutine constitutive_mech_forward
module subroutine mech_restore(ip,el,includeL)
integer, intent(in) :: &
ip, &
el
logical, intent(in) :: &
includeL
end subroutine mech_restore
! == cleaned:end ===================================================================================
module function crystallite_stress(dt,co,ip,el) result(converged_)
@ -392,8 +398,7 @@ module constitutive
crystallite_restartRead, &
constitutive_initializeRestorationPoints, &
constitutive_windForward, &
crystallite_restore, &
PLASTICITY_UNDEFINED_ID, &
PLASTICITY_UNDEFINED_ID, &
PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_PHENOPOWERLAW_ID, &
@ -734,20 +739,21 @@ subroutine constitutive_allocateState(state, &
sizeDotState, &
sizeDeltaState
state%sizeState = sizeState
state%sizeDotState = sizeDotState
state%sizeDeltaState = sizeDeltaState
state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (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%subState0 (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)
allocate(state%deltaState(sizeDeltaState,Nconstituents), source=0.0_pReal)
allocate(state%deltaState (sizeDeltaState,Nconstituents), source=0.0_pReal)
end subroutine constitutive_allocateState
@ -756,22 +762,27 @@ end subroutine constitutive_allocateState
!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
subroutine constitutive_restore(ip,el)
subroutine constitutive_restore(ip,el,includeL)
logical, intent(in) :: includeL
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer :: &
co, & !< constituent number
s
so
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
do s = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(material_phaseAt(co,el))%p(s)%state( :,material_phasememberAt(co,ip,el)) = &
sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phasememberAt(co,ip,el))
do so = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = &
sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el))
enddo
enddo
call mech_restore(ip,el,includeL)
end subroutine constitutive_restore
@ -783,7 +794,7 @@ subroutine constitutive_forward
integer :: i, j
crystallite_F0 = crystallite_partitionedF
crystallite_F0 = crystallite_F
crystallite_Lp0 = crystallite_Lp
crystallite_S0 = crystallite_S
@ -830,12 +841,13 @@ subroutine crystallite_init
Nconstituents, &
ph, &
me, &
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el, & !< counter in element loop
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el, & !< counter in element loop
cMax, & !< maximum number of integration point components
iMax, & !< maximum number of integration points
eMax !< maximum number of elements
class(tNode), pointer :: &
num_crystallite, &
@ -854,23 +866,21 @@ subroutine crystallite_init
iMax = discretization_nIPs
eMax = discretization_Nelems
allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal)
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_partitionedLp0, &
crystallite_S,crystallite_P, &
crystallite_Fe,crystallite_Lp, &
crystallite_subF, &
crystallite_subFp0,crystallite_subFi0, &
source = crystallite_partitionedF)
source = crystallite_F)
allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_orientation(cMax,iMax,eMax))
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal)
@ -933,8 +943,8 @@ subroutine crystallite_init
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))
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)
@ -953,12 +963,12 @@ subroutine crystallite_init
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
enddo
enddo; enddo
!$OMP END PARALLEL DO
crystallite_partitionedF0 = crystallite_F0
crystallite_partitionedF = crystallite_F0
crystallite_F = crystallite_F0
!$OMP PARALLEL DO PRIVATE(ph,me)
@ -978,9 +988,6 @@ subroutine crystallite_init
end subroutine crystallite_init
!--------------------------------------------------------------------------------------------------
!> @brief Backup data for homog cutback.
!--------------------------------------------------------------------------------------------------
@ -991,7 +998,7 @@ subroutine constitutive_initializeRestorationPoints(ip,el)
el !< element number
integer :: &
co, & !< constituent number
s,ph, me
so,ph, me
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
@ -1002,9 +1009,9 @@ subroutine constitutive_initializeRestorationPoints(ip,el)
call mech_initializeRestorationPoints(ph,me)
do s = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phasememberAt(co,ip,el)) = &
sourceState(material_phaseAt(co,el))%p(s)%state0( :,material_phasememberAt(co,ip,el))
do so = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) = &
sourceState(material_phaseAt(co,el))%p(so)%state0( :,material_phasememberAt(co,ip,el))
enddo
enddo
@ -1019,57 +1026,28 @@ subroutine constitutive_windForward(ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer :: &
co, & !< constituent number
s, ph, me
so, ph, me
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
crystallite_partitionedF0 (1:3,1:3,co,ip,el) = crystallite_partitionedF(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)
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)
do s = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(ph)%p(s)%partitionedState0(:,me) = sourceState(ph)%p(s)%state(:,me)
do so = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state(:,me)
enddo
enddo
end subroutine constitutive_windForward
!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restore(ip,el,includeL)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
logical, intent(in) :: &
includeL !< protect agains fake cutback
integer :: &
co, p, m !< constituent number
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
p = material_phaseAt(co,el)
m = material_phaseMemberAt(co,ip,el)
if (includeL) then
crystallite_Lp(1:3,1:3,co,ip,el) = crystallite_partitionedLp0(1:3,1:3,co,ip,el)
constitutive_mech_Li(p)%data(1:3,1:3,m) = constitutive_mech_partitionedLi0(p)%data(1:3,1:3,m)
endif ! maybe protecting everything from overwriting makes more sense
constitutive_mech_Fp(p)%data(1:3,1:3,m) = constitutive_mech_partitionedFp0(p)%data(1:3,1:3,m)
constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_partitionedFi0(p)%data(1:3,1:3,m)
crystallite_S (1:3,1:3,co,ip,el) = crystallite_partitionedS0 (1:3,1:3,co,ip,el)
plasticState (material_phaseAt(co,el))%state( :,material_phasememberAt(co,ip,el)) = &
plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phasememberAt(co,ip,el))
enddo
end subroutine crystallite_restore
!--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF).
!--------------------------------------------------------------------------------------------------
@ -1080,13 +1058,13 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
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, temp_33_4
temp_33_1, temp_33_2, temp_33_3
real(pReal), dimension(3,3,3,3) :: dSdFe, &
dSdF, &
dSdFi, &
@ -1102,6 +1080,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
real(pReal), dimension(9,9):: temp_99
logical :: error
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
@ -1149,8 +1128,8 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
!--------------------------------------------------------------------------------------------------
! calculate dSdF
temp_33_1 = transpose(matmul(invFp,invFi))
temp_33_2 = matmul(crystallite_subF(1:3,1:3,co,ip,el),invSubFp0)
temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,co,ip,el),invFp), invSubFi0)
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)
@ -1180,21 +1159,20 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
!--------------------------------------------------------------------------------------------------
! assemble dPdF
temp_33_1 = matmul(crystallite_S(1:3,1:3,co,ip,el),transpose(invFp))
temp_33_2 = matmul(invFp,temp_33_1)
temp_33_3 = matmul(crystallite_subF(1:3,1:3,co,ip,el),invFp)
temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,co,ip,el))
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(temp_33_2)
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_subF(1:3,1:3,co,ip,el), &
+ 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_3,dSdF(1:3,1:3,p,o)), &
+ matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)), &
transpose(invFp)) &
+ matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
+ matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo
end function crystallite_stressTangent
@ -1237,7 +1215,7 @@ function crystallite_push33ToRef(co,ip,el, tensor33)
T = matmul(material_orientation0(co,ip,el)%asMatrix(), & ! ToDo: initial orientation correct?
transpose(math_inv33(crystallite_subF(1:3,1:3,co,ip,el))))
transpose(math_inv33(crystallite_F(1:3,1:3,co,ip,el))))
crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
end function crystallite_push33ToRef
@ -1247,8 +1225,9 @@ end function crystallite_push33ToRef
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
function integrateSourceState(co,ip,el) result(broken)
function integrateSourceState(dt,co,ip,el) result(broken)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
@ -1281,8 +1260,7 @@ function integrateSourceState(co,ip,el) result(broken)
do so = 1, phase_Nsources(ph)
size_so(so) = sourceState(ph)%p(so)%sizeDotState
sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) &
+ sourceState(ph)%p(so)%dotState (1:size_so(so),me) &
* crystallite_subdt(co,ip,el)
+ sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt
source_dotState(1:size_so(so),2,so) = 0.0_pReal
enddo
@ -1304,8 +1282,8 @@ function integrateSourceState(co,ip,el) result(broken)
sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta &
+ source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta)
r(1:size_so(so)) = sourceState(ph)%p(so)%state (1:size_so(so),me) &
- sourceState(ph)%p(so)%subState0(1:size_so(so),me) &
- sourceState(ph)%p(so)%dotState (1:size_so(so),me) * crystallite_subdt(co,ip,el)
- sourceState(ph)%p(so)%subState0(1:size_so(so),me) &
- sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt
sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) &
- r(1:size_so(so))
converged_ = converged_ .and. converged(r(1:size_so(so)), &
@ -1371,7 +1349,7 @@ end function converged
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restartWrite
integer :: i
integer :: ph
integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: fileName, datasetName
@ -1380,27 +1358,27 @@ subroutine crystallite_restartWrite
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'a')
call HDF5_write(fileHandle,crystallite_partitionedF,'F')
call HDF5_write(fileHandle,crystallite_F,'F')
call HDF5_write(fileHandle,crystallite_Lp, 'L_p')
call HDF5_write(fileHandle,crystallite_S, 'S')
groupHandle = HDF5_addGroup(fileHandle,'phase')
do i = 1,size(material_name_phase)
write(datasetName,'(i0,a)') i,'_omega'
call HDF5_write(groupHandle,plasticState(i)%state,datasetName)
write(datasetName,'(i0,a)') i,'_F_i'
call HDF5_write(groupHandle,constitutive_mech_Fi(i)%data,datasetName)
write(datasetName,'(i0,a)') i,'_L_i'
call HDF5_write(groupHandle,constitutive_mech_Li(i)%data,datasetName)
write(datasetName,'(i0,a)') i,'_F_p'
call HDF5_write(groupHandle,constitutive_mech_Fp(i)%data,datasetName)
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
call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_addGroup(fileHandle,'homogenization')
do i = 1, size(material_name_homogenization)
write(datasetName,'(i0,a)') i,'_omega'
call HDF5_write(groupHandle,homogState(i)%state,datasetName)
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)
@ -1415,7 +1393,7 @@ end subroutine crystallite_restartWrite
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restartRead
integer :: i
integer :: ph
integer(HID_T) :: fileHandle, groupHandle
character(len=pStringLen) :: fileName, datasetName
@ -1429,22 +1407,22 @@ subroutine crystallite_restartRead
call HDF5_read(fileHandle,crystallite_S0, 'S')
groupHandle = HDF5_openGroup(fileHandle,'phase')
do i = 1,size(material_name_phase)
write(datasetName,'(i0,a)') i,'_omega'
call HDF5_read(groupHandle,plasticState(i)%state0,datasetName)
write(datasetName,'(i0,a)') i,'_F_i'
call HDF5_read(groupHandle,constitutive_mech_Fi0(i)%data,datasetName)
write(datasetName,'(i0,a)') i,'_L_i'
call HDF5_read(groupHandle,constitutive_mech_Li0(i)%data,datasetName)
write(datasetName,'(i0,a)') i,'_F_p'
call HDF5_read(groupHandle,constitutive_mech_Fp0(i)%data,datasetName)
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
call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_openGroup(fileHandle,'homogenization')
do i = 1,size(material_name_homogenization)
write(datasetName,'(i0,a)') i,'_omega'
call HDF5_read(groupHandle,homogState(i)%state0,datasetName)
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)

148
src/constitutive_mech.f90
View File

@ -800,7 +800,7 @@ function integrateStress(F,Delta_t,co,ip,el) result(broken)
broken = .true.
call constitutive_plastic_dependentState(crystallite_partitionedF(1:3,1:3,co,ip,el),co,ip,el)
call constitutive_plastic_dependentState(crystallite_F(1:3,1:3,co,ip,el),co,ip,el)
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
@ -959,19 +959,21 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: &
broken
integer :: &
NiterationState, & !< number of iterations in state loop
ph, &
me, &
size_pl
sizeDotState
real(pReal) :: &
zeta
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: &
r ! state residuum
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: &
plastic_dotState
logical :: &
broken
dotState
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
@ -979,15 +981,15 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
if(broken) return
size_pl = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%subState0(1:size_pl,me) &
+ plasticState(ph)%dotState (1:size_pl,me) * Delta_t
plastic_dotState(1:size_pl,2) = 0.0_pReal
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) &
+ plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t
dotState(1:sizeDotState,2) = 0.0_pReal
iteration: do NiterationState = 1, num%nState
if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1)
plastic_dotState(1:size_pl,1) = plasticState(ph)%dotState(:,me)
if(nIterationState > 1) dotState(1:sizeDotState,2) = dotState(1:sizeDotState,1)
dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,me)
broken = integrateStress(F,Delta_t,co,ip,el)
if(broken) exit iteration
@ -995,16 +997,16 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
if(broken) exit iteration
zeta = damper(plasticState(ph)%dotState(:,me),plastic_dotState(1:size_pl,1),&
plastic_dotState(1:size_pl,2))
zeta = damper(plasticState(ph)%dotState(:,me),dotState(1:sizeDotState,1),&
dotState(1:sizeDotState,2))
plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) * zeta &
+ plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta)
r(1:size_pl) = plasticState(ph)%state (1:size_pl,me) &
- plasticState(ph)%subState0(1:size_pl,me) &
- plasticState(ph)%dotState (1:size_pl,me) * Delta_t
plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%state(1:size_pl,me) &
- r(1:size_pl)
if (converged(r(1:size_pl),plasticState(ph)%state(1:size_pl,me),plasticState(ph)%atol(1:size_pl))) then
+ dotState(1:sizeDotState,1) * (1.0_pReal - zeta)
r(1:sizeDotState) = plasticState(ph)%state (1:sizeDotState,me) &
- plasticState(ph)%subState0(1:sizeDotState,me) &
- plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) &
- r(1:sizeDotState)
if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then
broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me)
exit iteration
@ -1012,6 +1014,7 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
enddo iteration
contains
!--------------------------------------------------------------------------------------------------
@ -1048,12 +1051,14 @@ function integrateStateEuler(F_0,F,Delta_t,co,ip,el) result(broken)
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: &
broken
integer :: &
ph, &
me, &
sizeDotState
logical :: &
broken
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
@ -1085,13 +1090,13 @@ function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken)
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: &
broken
integer :: &
ph, &
me, &
sizeDotState
logical :: &
broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic
@ -1105,7 +1110,7 @@ function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken)
residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,me) * 0.5_pReal * Delta_t
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) &
+ plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t
+ plasticState(ph)%dotstate(1:sizeDotState,me) * Delta_t
broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me)
@ -1145,6 +1150,7 @@ function integrateStateRK4(F_0,F,Delta_t,co,ip,el) result(broken)
real(pReal), dimension(4), parameter :: &
B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal]
broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C)
end function integrateStateRK4
@ -1178,6 +1184,7 @@ function integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) result(broken)
[2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,&
13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal]
broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
end function integrateStateRKCK45
@ -1215,18 +1222,18 @@ function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken)
broken = mech_collectDotState(Delta_t,co,ip,el,ph,me)
if(broken) return
sizeDotState = plasticState(ph)%sizeDotState
do stage = 1, size(A,1)
sizeDotState = plasticState(ph)%sizeDotState
plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,me)
plasticState(ph)%dotState(:,me) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1)
do n = 2, stage
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) &
+ A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
+ A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
enddo
sizeDotState = plasticState(ph)%sizeDotState
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) &
+ plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t
@ -1239,7 +1246,6 @@ function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken)
enddo
if(broken) return
sizeDotState = plasticState(ph)%sizeDotState
plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,me)
plasticState(ph)%dotState(:,me) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B)
@ -1282,7 +1288,7 @@ subroutine crystallite_results(group,ph)
select case (output_constituent(ph)%label(ou))
case('F')
selected_tensors = select_tensors(crystallite_partitionedF,ph)
selected_tensors = select_tensors(crystallite_F,ph)
call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),&
'deformation gradient','1')
case('F_e')
@ -1482,25 +1488,24 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
formerSubStep
integer :: &
NiterationCrystallite, & ! number of iterations in crystallite loop
s, ph, me
so, ph, me
logical :: todo
real(pReal) :: subFrac,subStep
real(pReal), dimension(3,3) :: &
subLp0, & !< plastic velocity grad at start of crystallite inc
subLi0, & !< intermediate velocity grad at start of crystallite inc
subF0
subF0, &
subF
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
subLi0 = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me)
subLp0 = crystallite_partitionedLp0(1:3,1:3,co,ip,el)
plasticState (material_phaseAt(co,el))%subState0( :,material_phaseMemberAt(co,ip,el)) = &
plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phaseMemberAt(co,ip,el))
do s = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(material_phaseAt(co,el))%p(s)%subState0( :,material_phaseMemberAt(co,ip,el)) = &
sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phaseMemberAt(co,ip,el))
plasticState(ph)%subState0(:,me) = plasticState(ph)%partitionedState0(:,me)
do so = 1, phase_Nsources(ph)
sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%partitionedState0(:,me)
enddo
crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)
crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)
@ -1525,16 +1530,14 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
todo = subStep > 0.0_pReal ! still time left to integrate on?
if (todo) then
subF0 = crystallite_subF(1:3,1:3,co,ip,el)
subF0 = subF
subLp0 = crystallite_Lp (1:3,1:3,co,ip,el)
subLi0 = constitutive_mech_Li(ph)%data(1:3,1:3,me)
crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_Fp(ph)%data(1:3,1:3,me)
crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_Fi(ph)%data(1:3,1:3,me)
plasticState( material_phaseAt(co,el))%subState0(:,material_phaseMemberAt(co,ip,el)) &
= plasticState(material_phaseAt(co,el))%state( :,material_phaseMemberAt(co,ip,el))
do s = 1, phase_Nsources(material_phaseAt(co,el))
sourceState( material_phaseAt(co,el))%p(s)%subState0(:,material_phaseMemberAt(co,ip,el)) &
= sourceState(material_phaseAt(co,el))%p(s)%state( :,material_phaseMemberAt(co,ip,el))
plasticState(ph)%subState0(:,me) = plasticState(ph)%state(:,me)
do so = 1, phase_Nsources(ph)
sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%state(:,me)
enddo
endif
!--------------------------------------------------------------------------------------------------
@ -1548,11 +1551,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
crystallite_Lp (1:3,1:3,co,ip,el) = subLp0
constitutive_mech_Li(ph)%data(1:3,1:3,me) = subLi0
endif
plasticState (material_phaseAt(co,el))%state( :,material_phaseMemberAt(co,ip,el)) &
= plasticState(material_phaseAt(co,el))%subState0(:,material_phaseMemberAt(co,ip,el))
do s = 1, phase_Nsources(material_phaseAt(co,el))
sourceState( material_phaseAt(co,el))%p(s)%state( :,material_phaseMemberAt(co,ip,el)) &
= sourceState(material_phaseAt(co,el))%p(s)%subState0(:,material_phaseMemberAt(co,ip,el))
plasticState(ph)%state(:,me) = plasticState(ph)%subState0(:,me)
do so = 1, phase_Nsources(ph)
sourceState(ph)%p(so)%state(:,me) = sourceState(ph)%p(so)%subState0(:,me)
enddo
todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair)
@ -1561,21 +1562,50 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
!--------------------------------------------------------------------------------------------------
! prepare for integration
if (todo) then
crystallite_subF(1:3,1:3,co,ip,el) = subF0 &
+ subStep *( crystallite_partitionedF (1:3,1:3,co,ip,el) &
-crystallite_partitionedF0(1:3,1:3,co,ip,el))
crystallite_Fe(1:3,1:3,co,ip,el) = matmul(crystallite_subF(1:3,1:3,co,ip,el), &
math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
subF = subF0 &
+ subStep * (crystallite_F(1:3,1:3,co,ip,el) - crystallite_partitionedF0(1:3,1:3,co,ip,el))
crystallite_Fe(1:3,1:3,co,ip,el) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
crystallite_subdt(co,ip,el) = subStep * dt
converged_ = .not. integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),&
crystallite_subdt(co,ip,el),co,ip,el)
converged_ = converged_ .and. .not. integrateSourceState(co,ip,el)
converged_ = .not. integrateState(subF0,subF,subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el)
endif
enddo cutbackLooping
end function crystallite_stress
!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
module subroutine mech_restore(ip,el,includeL)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
logical, intent(in) :: &
includeL !< protect agains fake cutback
integer :: &
co, p, m !< constituent number
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
p = material_phaseAt(co,el)
m = material_phaseMemberAt(co,ip,el)
if (includeL) then
crystallite_Lp(1:3,1:3,co,ip,el) = crystallite_partitionedLp0(1:3,1:3,co,ip,el)
constitutive_mech_Li(p)%data(1:3,1:3,m) = constitutive_mech_partitionedLi0(p)%data(1:3,1:3,m)
endif ! maybe protecting everything from overwriting makes more sense
constitutive_mech_Fp(p)%data(1:3,1:3,m) = constitutive_mech_partitionedFp0(p)%data(1:3,1:3,m)
constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_partitionedFi0(p)%data(1:3,1:3,m)
crystallite_S (1:3,1:3,co,ip,el) = crystallite_partitionedS0 (1:3,1:3,co,ip,el)
plasticState (material_phaseAt(co,el))%state( :,material_phasememberAt(co,ip,el)) = &
plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phasememberAt(co,ip,el))
enddo
end subroutine mech_restore
end submodule constitutive_mech

144
src/homogenization.f90
View File

@ -70,29 +70,22 @@ module homogenization
end subroutine mech_homogenize
module subroutine mech_results(group_base,h)
character(len=*), intent(in) :: group_base
integer, intent(in) :: h
end subroutine mech_results
! -------- ToDo ---------------------------------------------------------
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
logical, dimension(2) :: mech_RGC_updateState
real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses
F,& !< partitioned deformation gradients
F0 !< partitioned initial deformation gradients
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
real(pReal), intent(in) :: dt !< time increment
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
end function mech_RGC_updateState
module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
real(pReal), intent(in) :: &
subdt !< current time step
real(pReal), intent(in), dimension(3,3) :: &
subF
integer, intent(in) :: &
ip, & !< integration point
el !< element number
logical, dimension(2) :: doneAndHappy
end function mech_updateState
end interface
! -----------------------------------------------------------------------
public :: &
homogenization_init, &
@ -148,11 +141,10 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
real(pReal), intent(in) :: dt !< time increment
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
integer :: &
NiterationHomog, &
NiterationMPstate, &
ip, & !< integration point number
el, & !< element number
myNgrains, co, ce, ho
myNgrains, co, ce, ho, me
real(pReal) :: &
subFrac, &
subStep
@ -162,12 +154,12 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
doneAndHappy
!$OMP PARALLEL DO PRIVATE(ce,ho,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,doneAndHappy)
!$OMP PARALLEL DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho)
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
me = material_homogenizationMemberAt(ip,el)
!--------------------------------------------------------------------------------------------------
! initialize restoration points
call constitutive_initializeRestorationPoints(ip,el)
@ -177,15 +169,10 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
if (homogState(ho)%sizeState > 0) &
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(ho)%State0( :,material_homogenizationMemberAt(ip,el))
homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me)
if (damageState(ho)%sizeState > 0) &
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(ho)%State0( :,material_homogenizationMemberAt(ip,el))
damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me)
NiterationHomog = 0
cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog)
if (converged) then
@ -198,33 +185,26 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
call constitutive_windForward(ip,el)
if(homogState(ho)%sizeState > 0) &
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
homogState(ho)%State (:,material_homogenizationMemberAt(ip,el))
homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me)
if(damageState(ho)%sizeState > 0) &
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
damageState(ho)%State (:,material_homogenizationMemberAt(ip,el))
damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me)
endif steppingNeeded
else
if ( (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
num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep
! cutback makes no sense
if (.not. terminallyIll) & ! so first signals terminally ill...
print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
terminallyIll = .true. ! ...and kills all others
else ! cutback makes sense
subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback
if (.not. terminallyIll) & ! so first signals terminally ill...
print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
terminallyIll = .true. ! ...and kills all others
else ! cutback makes sense
subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback
call crystallite_restore(ip,el,subStep < 1.0_pReal)
call constitutive_restore(ip,el)
call constitutive_restore(ip,el,subStep < 1.0_pReal)
if(homogState(ho)%sizeState > 0) &
homogState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = &
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el))
if(damageState(ho)%sizeState > 0) &
damageState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = &
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el))
endif
if(homogState(ho)%sizeState > 0) &
homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
if(damageState(ho)%sizeState > 0) &
damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me)
endif
if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
@ -253,18 +233,16 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
doneAndHappy = [.true.,.false.]
else
ce = (el-1)*discretization_nIPs + ip
doneAndHappy = updateState(dt*subStep, &
homogenization_F0(1:3,1:3,ce) &
+ (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) &
doneAndHappy = mech_updateState(dt*subStep, &
homogenization_F0(1:3,1:3,ce) &
+ (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) &
*(subStep+subFrac), &
ip,el)
ip,el)
converged = all(doneAndHappy)
endif
endif
enddo convergenceLooping
NiterationHomog = NiterationHomog + 1
enddo cutBackLooping
enddo
enddo
@ -290,74 +268,35 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
end subroutine materialpoint_stressAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
!> "happy" with result
!--------------------------------------------------------------------------------------------------
function updateState(subdt,subF,ip,el)
real(pReal), intent(in) :: &
subdt !< current time step
real(pReal), intent(in), dimension(3,3) :: &
subF
integer, intent(in) :: &
ip, & !< integration point
el !< element number
integer :: c
logical, dimension(2) :: updateState
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
updateState = .true.
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c=1,homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo
updateState = &
updateState .and. &
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
crystallite_partitionedF(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)
end select chosenHomogenization
end function updateState
!--------------------------------------------------------------------------------------------------
!> @brief writes homogenization results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine homogenization_results
use material, only: &
material_homogenization_type => homogenization_type
integer :: p
integer :: ho
character(len=:), allocatable :: group_base,group
call results_closeGroup(results_addGroup('current/homogenization/'))
do p=1,size(material_name_homogenization)
group_base = 'current/homogenization/'//trim(material_name_homogenization(p))
do ho=1,size(material_name_homogenization)
group_base = 'current/homogenization/'//trim(material_name_homogenization(ho))
call results_closeGroup(results_addGroup(group_base))
call mech_results(group_base,p)
call mech_results(group_base,ho)
group = trim(group_base)//'/damage'
call results_closeGroup(results_addGroup(group))
select case(damage_type(p))
select case(damage_type(ho))
case(DAMAGE_NONLOCAL_ID)
call damage_nonlocal_results(p,group)
call damage_nonlocal_results(ho,group)
end select
group = trim(group_base)//'/thermal'
call results_closeGroup(results_addGroup(group))
select case(thermal_type(p))
select case(thermal_type(ho))
case(THERMAL_CONDUCTION_ID)
call thermal_conduction_results(p,group)
call thermal_conduction_results(ho,group)
end select
enddo
@ -373,6 +312,7 @@ subroutine homogenization_forward
integer :: ho
do ho = 1, size(material_name_homogenization)
homogState (ho)%state0 = homogState (ho)%state
damageState(ho)%state0 = damageState(ho)%state

60
src/homogenization_mech.f90
View File

@ -52,6 +52,21 @@ submodule(homogenization) homogenization_mech
end subroutine mech_RGC_averageStressAndItsTangent
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses
F,& !< partitioned deformation gradients
F0 !< partitioned initial deformation gradients
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
real(pReal), intent(in) :: dt !< time increment
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
end function mech_RGC_updateState
module subroutine mech_RGC_results(instance,group)
integer, intent(in) :: instance !< homogenization instance
character(len=*), intent(in) :: group !< group name in HDF5 file
@ -101,16 +116,16 @@ module subroutine mech_partition(subF,ip,el)
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
crystallite_partitionedF(1:3,1:3,1,ip,el) = subF
crystallite_F(1:3,1:3,1,ip,el) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(&
crystallite_partitionedF(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), &
subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(&
crystallite_partitionedF(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), &
subF,&
ip, &
el)
@ -166,6 +181,45 @@ module subroutine mech_homogenize(ip,el)
end subroutine mech_homogenize
!--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
!> "happy" with result
!--------------------------------------------------------------------------------------------------
module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
real(pReal), intent(in) :: &
subdt !< current time step
real(pReal), intent(in), dimension(3,3) :: &
subF
integer, intent(in) :: &
ip, & !< integration point
el !< element number
logical, dimension(2) :: doneAndHappy
integer :: co
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
enddo
doneAndHappy = &
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
doneAndHappy = .true.
endif
end function mech_updateState
!--------------------------------------------------------------------------------------------------
!> @brief Write results to file.
!--------------------------------------------------------------------------------------------------

82
src/homogenization_mech_RGC.f90
View File

@ -8,6 +8,7 @@
!--------------------------------------------------------------------------------------------------
submodule(homogenization:homogenization_mech) homogenization_mech_RGC
use rotations
use lattice
type :: tParameters
integer, dimension(:), allocatable :: &
@ -242,7 +243,18 @@ end subroutine mech_RGC_partitionDeformation
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
! "happy" with result
!--------------------------------------------------------------------------------------------------
module procedure mech_RGC_updateState
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses
F,& !< partitioned deformation gradients
F0 !< partitioned initial deformation gradients
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
real(pReal), intent(in) :: dt !< time increment
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
integer, dimension(4) :: intFaceN,intFaceP,faceID
integer, dimension(3) :: nGDim,iGr3N,iGr3P
@ -256,7 +268,7 @@ module procedure mech_RGC_updateState
real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
zeroTimeStep: if(dEq0(dt)) then
mech_RGC_updateState = .true. ! pretend everything is fine and return
doneAndHappy = .true. ! pretend everything is fine and return
return
endif zeroTimeStep
@ -327,12 +339,12 @@ module procedure mech_RGC_updateState
stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress
residMax = maxval(abs(tract)) ! get the maximum of the residual
mech_RGC_updateState = .false.
doneAndHappy = .false.
!--------------------------------------------------------------------------------------------------
! If convergence reached => done and happy
if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
mech_RGC_updateState = .true.
doneAndHappy = .true.
!--------------------------------------------------------------------------------------------------
! compute/update the state for postResult, i.e., all energy densities computed by time-integration
@ -354,7 +366,7 @@ module procedure mech_RGC_updateState
!--------------------------------------------------------------------------------------------------
! if residual blows-up => done but unhappy
elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then ! try to restart when residual blows up exceeding maximum bound
mech_RGC_updateState = [.true.,.false.] ! with direct cut-back
doneAndHappy = [.true.,.false.] ! with direct cut-back
return
endif
@ -484,7 +496,7 @@ module procedure mech_RGC_updateState
enddo; enddo
stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration
if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
mech_RGC_updateState = [.true.,.false.]
doneAndHappy = [.true.,.false.]
!$OMP CRITICAL (write2out)
print'(a,i3,a,i3,a)',' RGC_updateState: ip ',ip,' | el ',el,' enforces cutback'
print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax))
@ -513,8 +525,10 @@ module procedure mech_RGC_updateState
real(pReal), dimension (3) :: nVect,surfCorr
real(pReal), dimension (2) :: Gmoduli
integer :: iGrain,iGNghb,iFace,i,j,k,l
real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
real(pReal), parameter :: nDefToler = 1.0e-10_pReal
real(pReal) :: muGrain,muGNghb,nDefNorm
real(pReal), parameter :: &
nDefToler = 1.0e-10_pReal, &
b = 2.5e-10_pReal ! Length of Burgers vector
nGDim = param(instance)%N_constituents
rPen = 0.0_pReal
@ -532,9 +546,7 @@ module procedure mech_RGC_updateState
!-----------------------------------------------------------------------------------------------
! computing the mismatch and penalty stress tensor of all grains
grainLoop: do iGrain = 1,product(prm%N_constituents)
Gmoduli = equivalentModuli(iGrain,ip,el)
muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain
bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector
muGrain = equivalentMu(iGrain,ip,el)
iGrain3 = grain1to3(iGrain,prm%N_constituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position
interfaceLoop: do iFace = 1,6
@ -546,9 +558,7 @@ module procedure mech_RGC_updateState
where(iGNghb3 < 1) iGNghb3 = nGDim
where(iGNghb3 >nGDim) iGNghb3 = 1
iGNghb = grain3to1(iGNghb3,prm%N_constituents) ! get the ID of the neighboring grain
Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor
muGNghb = Gmoduli(1)
bgGNghb = Gmoduli(2)
muGNghb = equivalentMu(iGNghb,ip,el)
gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! difference/jump in deformation gradeint across the neighbor
!-------------------------------------------------------------------------------------------
@ -568,7 +578,7 @@ module procedure mech_RGC_updateState
!-------------------------------------------------------------------------------------------
! compute the stress penalty of all interfaces
do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xi_alpha &
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*b + muGNghb*b)*prm%xi_alpha &
*surfCorr(abs(intFace(1)))/prm%D_alpha(abs(intFace(1))) &
*cosh(prm%c_alpha*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
@ -655,44 +665,26 @@ module procedure mech_RGC_updateState
end function surfaceCorrection
!--------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
!--------------------------------------------------------------------------------------------------
function equivalentModuli(grainID,ip,el)
real(pReal), dimension(2) :: equivalentModuli
!-------------------------------------------------------------------------------------------------
real(pReal) function equivalentMu(grainID,ip,el)
integer, intent(in) :: &
grainID,&
ip, & !< integration point number
el !< element number
real(pReal), dimension(6,6) :: elasTens
real(pReal) :: &
cEquiv_11, &
cEquiv_12, &
cEquiv_44
elasTens = constitutive_homogenizedC(grainID,ip,el)
!----------------------------------------------------------------------------------------------
! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005)
cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal
cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + &
elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal
cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal
equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44
!----------------------------------------------------------------------------------------------
! obtain the length of Burgers vector (could be model dependend)
equivalentModuli(2) = 2.5e-10_pReal
end function equivalentModuli
!--------------------------------------------------------------------------------------------------
equivalentMu = lattice_equivalent_mu(constitutive_homogenizedC(grainID,ip,el),'voigt')
end function equivalentMu
!-------------------------------------------------------------------------------------------------
!> @brief calculating the grain deformation gradient (the same with
! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
!--------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
subroutine grainDeformation(F, avgF, instance, of)
real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain
@ -707,7 +699,7 @@ module procedure mech_RGC_updateState
integer, dimension(3) :: iGrain3
integer :: iGrain,iFace,i,j
!-------------------------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
associate(prm => param(instance))
@ -729,7 +721,7 @@ module procedure mech_RGC_updateState
end subroutine grainDeformation
end procedure mech_RGC_updateState
end function mech_RGC_updateState
!--------------------------------------------------------------------------------------------------

34
src/lattice.f90
View File

@ -421,6 +421,8 @@ module lattice
lattice_BCT_ID, &
lattice_HEX_ID, &
lattice_ORT_ID, &
lattice_equivalent_nu, &
lattice_equivalent_mu, &
lattice_applyLatticeSymmetry33, &
lattice_SchmidMatrix_slip, &
lattice_SchmidMatrix_twin, &
@ -508,8 +510,8 @@ subroutine lattice_init
lattice_C66(1:6,1:6,p) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,p),phase%get_asString('lattice'))
lattice_mu(p) = equivalent_mu(lattice_C66(1:6,1:6,p),'voigt')
lattice_nu(p) = equivalent_nu(lattice_C66(1:6,1:6,p),'voigt')
lattice_nu(p) = lattice_equivalent_nu(lattice_C66(1:6,1:6,p),'voigt')
lattice_mu(p) = lattice_equivalent_mu(lattice_C66(1:6,1:6,p),'voigt')
lattice_C66(1:6,1:6,p) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,p))) ! Literature data is in Voigt notation
do i = 1, 6
@ -2188,15 +2190,16 @@ end function getlabels
!> @brief Equivalent Poisson's ratio (ν)
!> @details https://doi.org/10.1143/JPSJ.20.635
!--------------------------------------------------------------------------------------------------
function equivalent_nu(C,assumption) result(nu)
function lattice_equivalent_nu(C,assumption) result(nu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
real(pReal) :: K, mu, nu
logical :: error
real(pReal), dimension(6,6) :: S
if (IO_lc(assumption) == 'voigt') then
K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) &
/ 9.0_pReal
@ -2210,25 +2213,26 @@ function equivalent_nu(C,assumption) result(nu)
K = 0.0_pReal
endif
mu = equivalent_mu(C,assumption)
mu = lattice_equivalent_mu(C,assumption)
nu = (1.5_pReal*K -mu)/(3.0_pReal*K+mu)
end function equivalent_nu
end function lattice_equivalent_nu
!--------------------------------------------------------------------------------------------------
!> @brief Equivalent shear modulus (μ)
!> @details https://doi.org/10.1143/JPSJ.20.635
!--------------------------------------------------------------------------------------------------
function equivalent_mu(C,assumption) result(mu)
function lattice_equivalent_mu(C,assumption) result(mu)
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
real(pReal) :: mu
logical :: error
real(pReal), dimension(6,6) :: S
if (IO_lc(assumption) == 'voigt') then
mu = (1.0_pReal*(C(1,1)+C(2,2)+C(3,3)) -1.0_pReal*(C(1,2)+C(2,3)+C(1,3)) +3.0_pReal*(C(4,4)+C(5,5)+C(6,6))) &
/ 15.0_pReal
@ -2242,7 +2246,7 @@ function equivalent_mu(C,assumption) result(mu)
mu = 0.0_pReal
endif
end function equivalent_mu
end function lattice_equivalent_mu
!--------------------------------------------------------------------------------------------------
@ -2266,14 +2270,14 @@ subroutine selfTest
call random_number(C)
C(1,1) = C(1,1) + 1.0_pReal
C = applyLatticeSymmetryC66(C,'aP')
if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
if(dNeq(C(6,6),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
if(dNeq(C(6,6),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
lambda = C(1,2)
if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'voigt')),equivalent_nu(C,'voigt'),1.0e-12_pReal)) &
error stop 'equivalent_nu/voigt'
if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'reuss')),equivalent_nu(C,'reuss'),1.0e-12_pReal)) &
error stop 'equivalent_nu/reuss'
if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), &
lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt'
if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), &
lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss'
end subroutine selfTest

7
src/math.f90
View File

@ -279,9 +279,12 @@ real(pReal) pure function math_LeviCivita(i,j,k)
integer, intent(in) :: i,j,k
if (all([i,j,k] == [1,2,3]) .or. all([i,j,k] == [2,3,1]) .or. all([i,j,k] == [3,1,2])) then
integer :: o
if (any([(all(cshift([i,j,k],o) == [1,2,3]),o=0,2)])) then
math_LeviCivita = +1.0_pReal
elseif (all([i,j,k] == [3,2,1]) .or. all([i,j,k] == [2,1,3]) .or. all([i,j,k] == [1,3,2])) then
elseif (any([(all(cshift([i,j,k],o) == [3,2,1]),o=0,2)])) then
math_LeviCivita = -1.0_pReal
else
math_LeviCivita = 0.0_pReal

32
src/prec.f90
View File

@ -108,8 +108,10 @@ logical elemental pure function dEq(a,b,tol)
real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal) :: eps
if (present(tol)) then
eps = tol
else
@ -132,11 +134,8 @@ logical elemental pure function dNeq(a,b,tol)
real(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
if (present(tol)) then
dNeq = .not. dEq(a,b,tol)
else
dNeq = .not. dEq(a,b)
endif
dNeq = .not. dEq(a,b,tol)
end function dNeq
@ -151,8 +150,10 @@ logical elemental pure function dEq0(a,tol)
real(pReal), intent(in) :: a
real(pReal), intent(in), optional :: tol
real(pReal) :: eps
if (present(tol)) then
eps = tol
else
@ -175,11 +176,8 @@ logical elemental pure function dNeq0(a,tol)
real(pReal), intent(in) :: a
real(pReal), intent(in), optional :: tol
if (present(tol)) then
dNeq0 = .not. dEq0(a,tol)
else
dNeq0 = .not. dEq0(a)
endif
dNeq0 = .not. dEq0(a,tol)
end function dNeq0
@ -195,8 +193,10 @@ logical elemental pure function cEq(a,b,tol)
complex(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
real(pReal) :: eps
if (present(tol)) then
eps = tol
else
@ -220,11 +220,8 @@ logical elemental pure function cNeq(a,b,tol)
complex(pReal), intent(in) :: a,b
real(pReal), intent(in), optional :: tol
if (present(tol)) then
cNeq = .not. cEq(a,b,tol)
else
cNeq = .not. cEq(a,b)
endif
cNeq = .not. cEq(a,b,tol)
end function cNeq
@ -238,6 +235,7 @@ pure function prec_bytesToC_FLOAT(bytes)
real(C_FLOAT), dimension(size(bytes,kind=pI64)/(storage_size(0._C_FLOAT,pI64)/8_pI64)) :: &
prec_bytesToC_FLOAT
prec_bytesToC_FLOAT = transfer(bytes,prec_bytesToC_FLOAT,size(prec_bytesToC_FLOAT))
end function prec_bytesToC_FLOAT
@ -252,6 +250,7 @@ pure function prec_bytesToC_DOUBLE(bytes)
real(C_DOUBLE), dimension(size(bytes,kind=pI64)/(storage_size(0._C_DOUBLE,pI64)/8_pI64)) :: &
prec_bytesToC_DOUBLE
prec_bytesToC_DOUBLE = transfer(bytes,prec_bytesToC_DOUBLE,size(prec_bytesToC_DOUBLE))
end function prec_bytesToC_DOUBLE
@ -266,6 +265,7 @@ pure function prec_bytesToC_INT32_T(bytes)
integer(C_INT32_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT32_T,pI64)/8_pI64)) :: &
prec_bytesToC_INT32_T
prec_bytesToC_INT32_T = transfer(bytes,prec_bytesToC_INT32_T,size(prec_bytesToC_INT32_T))
end function prec_bytesToC_INT32_T
@ -280,6 +280,7 @@ pure function prec_bytesToC_INT64_T(bytes)
integer(C_INT64_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT64_T,pI64)/8_pI64)) :: &
prec_bytesToC_INT64_T
prec_bytesToC_INT64_T = transfer(bytes,prec_bytesToC_INT64_T,size(prec_bytesToC_INT64_T))
end function prec_bytesToC_INT64_T
@ -295,6 +296,7 @@ subroutine selfTest
integer(pInt), dimension(1) :: i
real(pReal), dimension(2) :: r
realloc_lhs_test = [1,2]
if (any(realloc_lhs_test/=[1,2])) error stop 'LHS allocation'