From 6f65de27fcf718c98ba5a3275c5f86e63ebe5a34 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 7 Jan 2021 14:59:12 +0100 Subject: [PATCH] not used was only used for reporting (see v.2.0.0) --- examples/FEM/polyXtal/material.yaml | 2 +- src/constitutive.f90 | 218 +++++++++++++--------------- src/constitutive_mech.f90 | 148 +++++++++++-------- src/homogenization.f90 | 144 ++++++------------ src/homogenization_mech.f90 | 60 +++++++- src/homogenization_mech_RGC.f90 | 82 +++++------ src/lattice.f90 | 34 +++-- src/math.f90 | 7 +- src/prec.f90 | 32 ++-- 9 files changed, 365 insertions(+), 362 deletions(-) diff --git a/examples/FEM/polyXtal/material.yaml b/examples/FEM/polyXtal/material.yaml index c7d17657d..333073150 100644 --- a/examples/FEM/polyXtal/material.yaml +++ b/examples/FEM/polyXtal/material.yaml @@ -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: diff --git a/src/constitutive.f90 b/src/constitutive.f90 index e65ce864d..696611549 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -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) diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index de6f2ae9f..c48c59ec9 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -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 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index ebf5fd50d..52553b57b 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -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 diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index e4499e9b7..641e960fd 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -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. !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index a89008e96..04ec73845 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -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 !-------------------------------------------------------------------------------------------------- diff --git a/src/lattice.f90 b/src/lattice.f90 index 676232efe..6af135e4e 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -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 diff --git a/src/math.f90 b/src/math.f90 index 8005b5406..6b89a9923 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -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 diff --git a/src/prec.f90 b/src/prec.f90 index 95b1116cd..4d73462c4 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -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'