Merge branch 'development' into 53-separate-mesh-for-different-solvers-3

This commit is contained in:
Martin Diehl 2019-01-28 12:53:09 +01:00
commit 9260af2eaa
12 changed files with 1935 additions and 2703 deletions

@ -1 +1 @@
Subproject commit 683bf0074f3fa079989b51f5a67aa593b7577f0b
Subproject commit beb9682fff7d4d6c65aba12ffd04c7441dc6ba6b

View File

@ -1 +1 @@
v2.0.2-1453-g8e56f0d0
v2.0.2-1540-ge2582a8d

View File

@ -1236,6 +1236,10 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'zero entry on stiffness diagonal'
case (136_pInt)
msg = 'zero entry on stiffness diagonal for transformed phase'
case (137_pInt)
msg = 'not defined for lattice structure'
case (138_pInt)
msg = 'not enough interaction parameters given'
!--------------------------------------------------------------------------------------------------
! errors related to the parsing of material.config

View File

@ -550,7 +550,7 @@ end function getString
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!--------------------------------------------------------------------------------------------------
function getFloats(this,key,defaultVal,requiredShape,requiredSize)
function getFloats(this,key,defaultVal,requiredSize)
use IO, only: &
IO_error, &
IO_stringValue, &
@ -561,7 +561,6 @@ function getFloats(this,key,defaultVal,requiredShape,requiredSize)
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
real(pReal), dimension(:), intent(in), optional :: defaultVal
integer(pInt), dimension(:), intent(in), optional :: requiredShape ! not useful (is always 1D array)
integer(pInt), intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item
integer(pInt) :: i
@ -601,7 +600,7 @@ end function getFloats
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all
!! values from the last occurrence. If key is not found exits with error unless default is given.
!--------------------------------------------------------------------------------------------------
function getInts(this,key,defaultVal,requiredShape,requiredSize)
function getInts(this,key,defaultVal,requiredSize)
use IO, only: &
IO_error, &
IO_stringValue, &
@ -611,8 +610,7 @@ function getInts(this,key,defaultVal,requiredShape,requiredSize)
integer(pInt), dimension(:), allocatable :: getInts
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
integer(pInt), dimension(:), intent(in), optional :: defaultVal, &
requiredShape ! not useful (is always 1D array)
integer(pInt), dimension(:), intent(in), optional :: defaultVal
integer(pInt), intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item
integer(pInt) :: i
@ -653,7 +651,7 @@ end function getInts
!! values from the last occurrence. If key is not found exits with error unless default is given.
!! If raw is true, the the complete string is returned, otherwise the individual chunks are returned
!--------------------------------------------------------------------------------------------------
function getStrings(this,key,defaultVal,requiredShape,raw)
function getStrings(this,key,defaultVal,raw)
use IO, only: &
IO_error, &
IO_StringValue
@ -663,7 +661,6 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
character(len=65536),dimension(:), intent(in), optional :: defaultVal
integer(pInt), dimension(:), intent(in), optional :: requiredShape
logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item
character(len=65536) :: str

View File

@ -151,7 +151,7 @@ subroutine constitutive_init()
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init
if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init
if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init
if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then
call plastic_nonlocal_init(FILEUNIT)
@ -365,7 +365,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
use plastic_nonlocal, only: &
plastic_nonlocal_microstructure
use plastic_dislotwin, only: &
plastic_dislotwin_microstructure
plastic_dislotwin_dependentState
use plastic_disloUCLA, only: &
plastic_disloUCLA_dependentState
@ -389,7 +389,9 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_microstructure(temperature(ho)%p(tme),ipc,ip,el)
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
@ -409,9 +411,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
pReal
use math, only: &
math_mul33x33, &
math_Mandel6to33, &
math_Mandel33to6, &
math_Plain99to3333
math_6toSym33, &
math_sym33to6, &
math_99to3333
use material, only: &
phasememberAt, &
phase_plasticity, &
@ -470,7 +472,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
S = math_Mandel6to33(S6)
S = math_6toSym33(S6)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -495,9 +497,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp99, math_sym33to6(Mp), &
temperature(ho)%p(tme),ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
dLp_dMp = math_99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_DISLOTWIN_ID) plasticityType
of = phasememberAt(ipc,ip,el)
@ -540,7 +542,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
math_inv33, &
math_det33, &
math_mul33x33, &
math_Mandel6to33
math_6toSym33
use material, only: &
phasememberAt, &
phase_plasticity, &
@ -597,7 +599,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
case (PLASTICITY_isotropic_ID) plasticityType
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, math_Mandel6to33(S6),instance,of)
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, math_6toSym33(S6),instance,of)
case default plasticityType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
@ -716,7 +718,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip
use math, only : &
math_mul33x33, &
math_mul3333xx33, &
math_Mandel66to3333, &
math_66toSym3333, &
math_I3
use material, only: &
material_phase, &
@ -749,7 +751,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip
i, j
ho = material_homogenizationAt(el)
C = math_Mandel66to3333(constitutive_homogenizedC(ipc,ip,el))
C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el))
DegradationLoop: do d = 1_pInt, phase_NstiffnessDegradations(material_phase(ipc,ip,el))
degradationType: select case(phase_stiffnessDegradation(d,material_phase(ipc,ip,el)))
@ -784,8 +786,8 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
debug_levelBasic
use math, only: &
math_mul33x33, &
math_Mandel6to33, &
math_Mandel33to6, &
math_6toSym33, &
math_sym33to6, &
math_mul33x33
use mesh, only: &
mesh_NcpElems, &
@ -854,13 +856,13 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
integer(pInt) :: &
ho, & !< homogenization
tme, & !< thermal member position
s, & !< counter in source loop
s, & !< counter in source loop
instance, of
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_6toSym33(S6))
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -890,7 +892,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState (math_Mandel33to6(Mp),FeArray,FpArray,temperature(ho)%p(tme), &
call plastic_nonlocal_dotState (math_sym33to6(Mp),FeArray,FpArray,temperature(ho)%p(tme), &
subdt,subfracArray,ip,el)
end select plasticityType
@ -920,7 +922,7 @@ end subroutine constitutive_collectDotState
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
use prec, only: &
pReal, &
pLongInt
@ -929,8 +931,7 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
debug_constitutive, &
debug_levelBasic
use math, only: &
math_Mandel6to33, &
math_Mandel33to6, &
math_sym33to6, &
math_mul33x33
use material, only: &
phasememberAt, &
@ -954,18 +955,17 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola Kirchhoff stress
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), dimension(3,3) :: &
Mp
integer(pInt) :: &
s, & !< counter in source loop
i, &
instance, of
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -975,13 +975,13 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
call plastic_kinehardening_deltaState(Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(math_Mandel33to6(Mp),ip,el)
call plastic_nonlocal_deltaState(math_sym33to6(Mp),ip,el)
end select plasticityType
sourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, &
@ -1001,7 +1001,7 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
use prec, only: &
pReal
use math, only: &
math_Mandel6to33, &
math_6toSym33, &
math_mul33x33
use mesh, only: &
mesh_NcpElems, &
@ -1076,7 +1076,7 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
constitutive_postResults = 0.0_pReal
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_6toSym33(S6))
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)

View File

@ -1654,29 +1654,7 @@ subroutine integrateStateFPI()
!$OMP END PARALLEL DO
call update_dependentState
!$OMP PARALLEL
! --- STRESS INTEGRATION ---
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo before stress integration'
#endif
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
crystallite_todo(g,i,e) = integrateStress(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ... then all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
call update_stress(1.0_pReal)
call update_dotState(1.0_pReal)
!$OMP PARALLEL
! --- UPDATE STATE ---
@ -1874,17 +1852,6 @@ end subroutine integrateStateFPI
subroutine integrateStateEuler()
use, intrinsic :: &
IEEE_arithmetic
#ifdef DEBUG
use debug, only: &
debug_e, &
debug_i, &
debug_g, &
debug_level, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective
#endif
use mesh, only: &
mesh_element, &
mesh_NcpElems
@ -1896,7 +1863,6 @@ subroutine integrateStateEuler()
constitutive_microstructure
implicit none
integer(pInt) :: &
e, & ! element index in element loop
i, & ! integration point index in ip loop
@ -1921,63 +1887,10 @@ eIter = FEsolving_execElem(1:2)
call update_dotState(1.0_pReal)
call update_State(1.0_pReal)
!$OMP PARALLEL
! --- STATE JUMP ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
crystallite_todo(g,i,e) = stateJump(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
call update_dependentState
!$OMP PARALLEL
! --- STRESS INTEGRATION ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
crystallite_todo(g,i,e) = integrateStress(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
! --- SET CONVERGENCE FLAG ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
crystallite_converged(g,i,e) = crystallite_todo(g,i,e) .or. crystallite_converged(g,i,e) ! if still "to do" then converged per definitionem
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
call update_deltaState
call update_dependentState
call update_stress(1.0_pReal)
call setConvergenceFlag
! --- CHECK NON-LOCAL CONVERGENCE ---
@ -2109,74 +2022,13 @@ subroutine integrateStateAdaptiveEuler()
endif
enddo; enddo; enddo
!$OMP ENDDO
! --- STATE JUMP ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then
crystallite_todo(g,i,e) = stateJump(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
call update_deltaState
call update_dependentState
! --- STRESS INTEGRATION (EULER INTEGRATION) ---
!$OMP PARALLEL DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then
crystallite_todo(g,i,e) = integrateStress(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
call update_dotState(1.0_pReal)
call update_stress(1.0_pReal)
call update_dotState(1.0_pReal)
!$OMP PARALLEL
!$OMP DO PRIVATE(p,c,NaN)
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do mySource = 1_pInt, phase_Nsources(p)
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
else ! if broken local...
crystallite_todo(g,i,e) = .false. ! ... skip this one next time
endif
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
! --- ERROR ESTIMATE FOR STATE (HEUN METHOD) ---
!$OMP SINGLE
@ -2365,46 +2217,9 @@ subroutine integrateStateRK4()
!$OMP END PARALLEL
call update_state(TIMESTEPFRACTION(n))
!$OMP PARALLEL
! --- state jump ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then
crystallite_todo(g,i,e) = stateJump(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
call update_deltaState
call update_dependentState
!$OMP PARALLEL
! --- stress integration ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then
crystallite_todo(g,i,e) = integrateStress(g,i,e,timeStepFraction(n)) ! fraction of original times step
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
call update_stress(TIMESTEPFRACTION(n))
! --- dot state and RK dot state---
@ -2414,14 +2229,7 @@ subroutine integrateStateRK4()
enddo
! --- SET CONVERGENCE FLAG ---
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
crystallite_converged(g,i,e) = crystallite_todo(g,i,e) .or. crystallite_converged(g,i,e) ! if still "to do" then converged per definitionem
enddo; enddo; enddo
call setConvergenceFlag
! --- CHECK NONLOCAL CONVERGENCE ---
@ -2584,51 +2392,10 @@ subroutine integrateStateRKCK45()
!$OMP END PARALLEL
call update_state(1.0_pReal) !MD: 1.0 correct?
!$OMP PARALLEL
! --- state jump ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then
crystallite_todo(g,i,e) = stateJump(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
call update_dependentState
!$OMP PARALLEL
! --- stress integration ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then
crystallite_todo(g,i,e) = integrateStress(g,i,e,C(stage)) ! fraction of original time step
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
call update_dotState(C(stage))
call update_deltaState
call update_dependentState
call update_stress(C(stage))
call update_dotState(C(stage))
enddo
@ -2735,56 +2502,12 @@ subroutine integrateStateRKCK45()
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
! --- STATE JUMP ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then
crystallite_todo(g,i,e) = stateJump(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
call update_deltaState
call update_dependentState
!$OMP PARALLEL
!--------------------------------------------------------------------------------------------------
! --- FINAL STRESS INTEGRATION STEP IF RESIDUUM BELOW TOLERANCE ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e)) then
crystallite_todo(g,i,e) = integrateStress(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!--------------------------------------------------------------------------------------------------
! --- SET CONVERGENCE FLAG ---
!$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
crystallite_converged(g,i,e) = crystallite_todo(g,i,e) .or. crystallite_converged(g,i,e) ! if still "to do" then converged per definition
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
call update_stress(1.0_pReal)
call setConvergenceFlag
! --- nonlocal convergence check ---
@ -2798,6 +2521,67 @@ subroutine integrateStateRKCK45()
end subroutine integrateStateRKCK45
!--------------------------------------------------------------------------------------------------
!> @brief Sets convergence flag based on "todo": every point that survived the integration (todo is
! still .true. is considered as converged
!> @details: For explicitEuler, RK4 and RKCK45, adaptive Euler and FPI have their on criteria
!--------------------------------------------------------------------------------------------------
subroutine setConvergenceFlag()
implicit none
integer(pInt) :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g !< grain index in grain loop
!OMP DO PARALLEL PRIVATE(i,g)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
g = 1:homogenization_Ngrains(mesh_element(3,e)))
crystallite_converged(g,i,e) = crystallite_todo(g,i,e) .or. crystallite_converged(g,i,e) ! if still "to do" then converged per definition
end forall; enddo
!OMP END DO PARALLEL
end subroutine setConvergenceFlag
!--------------------------------------------------------------------------------------------------
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
!--------------------------------------------------------------------------------------------------
subroutine update_stress(timeFraction)
use material, only: &
plasticState, &
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt
implicit none
real(pReal), intent(in) :: &
timeFraction
integer(pInt) :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e))
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
crystallite_todo(g,i,e) = integrateStress(g,i,e,timeFraction)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
end subroutine update_stress
!--------------------------------------------------------------------------------------------------
!> @brief tbd
!--------------------------------------------------------------------------------------------------
@ -2824,6 +2608,7 @@ subroutine update_dependentState()
end subroutine update_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
!--------------------------------------------------------------------------------------------------
@ -2886,7 +2671,6 @@ subroutine update_dotState(timeFraction)
constitutive_collectDotState
implicit none
real(pReal), intent(in) :: &
timeFraction
integer(pInt) :: &
@ -2897,15 +2681,17 @@ subroutine update_dotState(timeFraction)
c, &
s
logical :: &
NaN
NaN, &
nonlocalStop
nonlocalStop = .false.
!$OMP PARALLEL
!$OMP DO PRIVATE (p,c,NaN)
!$OMP PARALLEL DO PRIVATE (p,c,NaN)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e))
!$OMP FLUSH(crystallite_todo)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
!$OMP FLUSH(nonlocalStop)
if (nonlocalStop .or. (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e))) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), &
@ -2918,20 +2704,94 @@ subroutine update_dotState(timeFraction)
enddo
if (NaN) then
crystallite_todo(g,i,e) = .false. ! this one done (and broken)
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
!$OMP END CRITICAL (checkTodo)
endif
if (.not. crystallite_localPlasticity(g,i,e)) nonlocalStop = .True.
endif
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
!$OMP END PARALLEL DO
if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity
end subroutine update_DotState
subroutine update_deltaState
use, intrinsic :: &
IEEE_arithmetic
use prec, only: &
dNeq0
use material, only: &
plasticState, &
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt
use constitutive, only: &
constitutive_collectDeltaState
use math, only: &
math_6toSym33
implicit none
integer(pInt) :: &
e, & !< element index in element loop
i, & !< integration point index in ip loop
g, & !< grain index in grain loop
p, &
mySize, &
myOffset, &
mySource, &
c, &
s
logical :: &
NaN, &
nonlocalStop
nonlocalStop = .false.
!$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize,mySource,NaN)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e))
!$OMP FLUSH(nonlocalStop)
if (nonlocalStop .or. (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e))) then
call constitutive_collectDeltaState(math_6toSym33(crystallite_Tstar_v(1:6,g,i,e)), &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fi(1:3,1:3,g,i,e), &
g,i,e)
p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
myOffset = plasticState(p)%offsetDeltaState
mySize = plasticState(p)%sizeDeltaState
NaN = any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))
if (.not. NaN) then
plasticState(p)%state(myOffset + 1_pInt: myOffset + mySize,c) = &
plasticState(p)%state(myOffset + 1_pInt: myOffset + mySize,c) + &
plasticState(p)%deltaState(1:mySize,c)
do mySource = 1_pInt, phase_Nsources(p)
myOffset = sourceState(p)%p(mySource)%offsetDeltaState
mySize = sourceState(p)%p(mySource)%sizeDeltaState
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))
if (.not. NaN) then
sourceState(p)%p(mySource)%state(myOffset + 1_pInt:myOffset +mySize,c) = &
sourceState(p)%p(mySource)%state(myOffset + 1_pInt:myOffset +mySize,c) + &
sourceState(p)%p(mySource)%deltaState(1:mySize,c)
endif
enddo
endif
crystallite_todo(g,i,e) = .not. NaN
if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken
crystallite_converged(g,i,e) = .false.
if (.not. crystallite_localPlasticity(g,i,e)) nonlocalStop = .true.
endif
endif
enddo; enddo; enddo
!$OMP END PARALLEL DO
if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity
end subroutine update_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief calculates a jump in the state according to the current state and the current stress
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
@ -2958,6 +2818,8 @@ logical function stateJump(ipc,ip,el)
phaseAt, phasememberAt
use constitutive, only: &
constitutive_collectDeltaState
use math, only: &
math_6toSym33
implicit none
integer(pInt), intent(in):: &
@ -2969,57 +2831,50 @@ logical function stateJump(ipc,ip,el)
c, &
p, &
mySource, &
myOffsetPlasticDeltaState, &
myOffsetSourceDeltaState, &
mySizePlasticDeltaState, &
mySizeSourceDeltaState
myOffset, &
mySize
c = phasememberAt(ipc,ip,el)
p = phaseAt(ipc,ip,el)
call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), &
call constitutive_collectDeltaState(math_6toSym33(crystallite_Tstar_v(1:6,ipc,ip,el)), &
crystallite_Fe(1:3,1:3,ipc,ip,el), &
crystallite_Fi(1:3,1:3,ipc,ip,el), &
ipc,ip,el)
myOffsetPlasticDeltaState = plasticState(p)%offsetDeltaState
mySizePlasticDeltaState = plasticState(p)%sizeDeltaState
myOffset = plasticState(p)%offsetDeltaState
mySize = plasticState(p)%sizeDeltaState
if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)))) then ! NaN occured in deltaState
if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then ! NaN occured in deltaState
stateJump = .false.
return
endif
plasticState(p)%state(myOffsetPlasticDeltaState + 1_pInt : &
myOffsetPlasticDeltaState + mySizePlasticDeltaState,c) = &
plasticState(p)%state(myOffsetPlasticDeltaState + 1_pInt : &
myOffsetPlasticDeltaState + mySizePlasticDeltaState,c) + &
plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)
plasticState(p)%state(myOffset + 1_pInt:myOffset + mySize,c) = &
plasticState(p)%state(myOffset + 1_pInt:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c)
do mySource = 1_pInt, phase_Nsources(p)
myOffsetSourceDeltaState = sourceState(p)%p(mySource)%offsetDeltaState
mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState
if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)))) then ! NaN occured in deltaState
myOffset = sourceState(p)%p(mySource)%offsetDeltaState
mySize = sourceState(p)%p(mySource)%sizeDeltaState
if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then ! NaN occured in deltaState
stateJump = .false.
return
endif
sourceState(p)%p(mySource)%state(myOffsetSourceDeltaState + 1_pInt : &
myOffsetSourceDeltaState + mySizeSourceDeltaState,c) = &
sourceState(p)%p(mySource)%state(myOffsetSourceDeltaState + 1_pInt : &
myOffsetSourceDeltaState + mySizeSourceDeltaState,c) + &
sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)
sourceState(p)%p(mySource)%state(myOffset + 1_pInt: myOffset + mySize,c) = &
sourceState(p)%p(mySource)%state(myOffset + 1_pInt: myOffset + mySize,c) + &
sourceState(p)%p(mySource)%deltaState(1:mySize,c)
enddo
#ifdef DEBUG
if (any(dNeq0(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c))) &
if (any(dNeq0(plasticState(p)%deltaState(1:mySize,c))) &
.and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3, /)') '<< CRYST >> update state at el ip ipc ',el,ip,ipc
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> deltaState', plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> deltaState', plasticState(p)%deltaState(1:mySize,c)
write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state', &
plasticState(p)%state(myOffsetPlasticDeltaState + 1_pInt : &
myOffsetPlasticDeltaState + mySizePlasticDeltaState,c)
plasticState(p)%state(myOffset + 1_pInt : &
myOffset + mySize,c)
endif
#endif

File diff suppressed because it is too large Load Diff

View File

@ -918,7 +918,8 @@ end subroutine material_parseTexture
!--------------------------------------------------------------------------------------------------
!> @brief allocates the plastic state of a phase
!--------------------------------------------------------------------------------------------------
subroutine material_allocatePlasticState(phase,NofMyPhase,sizeState,sizeDotState,sizeDeltaState,&
subroutine material_allocatePlasticState(phase,NofMyPhase,&
sizeState,sizeDotState,sizeDeltaState,&
Nslip,Ntwin,Ntrans)
use numerics, only: &
numerics_integrator2 => numerics_integrator ! compatibility hack
@ -936,9 +937,10 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,sizeState,sizeDotState
integer(pInt) :: numerics_integrator ! compatibility hack
numerics_integrator = numerics_integrator2(1) ! compatibility hack
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
plasticState(phase)%Nslip = Nslip
plasticState(phase)%Ntwin = Ntwin
plasticState(phase)%Ntrans= Ntrans

View File

@ -28,8 +28,7 @@ module plastic_disloUCLA
shearrate_ID, &
accumulatedshear_ID, &
mfp_ID, &
thresholdstress_ID, &
dipoledistance_ID
thresholdstress_ID
end enum
type, private :: tParameters
@ -73,7 +72,7 @@ module plastic_disloUCLA
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output
logical :: &
dipoleformation
dipoleFormation !< flag indicating consideration of dipole formation
end type !< container type for internal constitutive parameters
type, private :: tDisloUCLAState
@ -93,7 +92,7 @@ module plastic_disloUCLA
!--------------------------------------------------------------------------------------------------
! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param
type(tDisloUCLAState ), allocatable, dimension(:), private :: &
type(tDisloUCLAState), allocatable, dimension(:), private :: &
dotState, &
state
type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState
@ -127,7 +126,6 @@ subroutine plastic_disloUCLA_init()
debug_constitutive,&
debug_levelBasic
use math, only: &
math_mul3x3, &
math_expand
use IO, only: &
IO_error, &
@ -148,8 +146,6 @@ subroutine plastic_disloUCLA_init()
implicit none
integer(pInt) :: &
index_myFamily, index_otherFamily, &
f,j,k,o, &
Ninstance, &
p, i, &
NipcMyPhase, &
@ -164,7 +160,6 @@ subroutine plastic_disloUCLA_init()
outputID
character(len=pStringLen) :: &
structure = '',&
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
@ -197,8 +192,6 @@ subroutine plastic_disloUCLA_init()
dst => dependentState(phase_plasticityInstance(p)), &
config => config_phase(p))
structure = config%getString('lattice_structure')
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%mu = lattice_mu(p)
@ -213,36 +206,41 @@ subroutine plastic_disloUCLA_init()
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
else
prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
structure(1:3))
prm%rho0 = config%getFloats('rhoedge0', requiredShape=shape(prm%Nslip))
prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredShape=shape(prm%Nslip))
prm%v0 = config%getFloats('v0', requiredShape=shape(prm%Nslip))
prm%burgers = config%getFloats('slipburgers', requiredShape=shape(prm%Nslip))
prm%H0kp = config%getFloats('qedge', requiredShape=shape(prm%Nslip))
config%getString('lattice_structure'))
prm%forestProjectionEdge = lattice_forestProjection(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%clambda = config%getFloats('clambdaslip', requiredShape=shape(prm%Nslip))
prm%tau_Peierls = config%getFloats('tau_peierls', requiredShape=shape(prm%Nslip)) ! ToDo: Deprecated
prm%p = config%getFloats('p_slip', requiredShape=shape(prm%Nslip), &
prm%rho0 = config%getFloats('rhoedge0', requiredSize=size(prm%Nslip))
prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%Nslip))
prm%v0 = config%getFloats('v0', requiredSize=size(prm%Nslip))
prm%burgers = config%getFloats('slipburgers', requiredSize=size(prm%Nslip))
prm%H0kp = config%getFloats('qedge', requiredSize=size(prm%Nslip))
prm%clambda = config%getFloats('clambdaslip', requiredSize=size(prm%Nslip))
prm%tau_Peierls = config%getFloats('tau_peierls', requiredSize=size(prm%Nslip)) ! ToDo: Deprecated
prm%p = config%getFloats('p_slip', requiredSize=size(prm%Nslip), &
defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%q = config%getFloats('q_slip', requiredShape=shape(prm%Nslip), &
prm%q = config%getFloats('q_slip', requiredSize=size(prm%Nslip), &
defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%kink_height = config%getFloats('kink_height', requiredShape=shape(prm%Nslip))
prm%w = config%getFloats('kink_width', requiredShape=shape(prm%Nslip))
prm%omega = config%getFloats('omega', requiredShape=shape(prm%Nslip))
prm%B = config%getFloats('friction_coeff', requiredShape=shape(prm%Nslip))
prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%Nslip))
prm%w = config%getFloats('kink_width', requiredSize=size(prm%Nslip))
prm%omega = config%getFloats('omega', requiredSize=size(prm%Nslip))
prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%Nslip))
prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! ToDo: Deprecated
prm%grainSize = config%getFloat('grainsize')
@ -250,7 +248,7 @@ subroutine plastic_disloUCLA_init()
prm%Qsd = config%getFloat('qsd')
prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers**3.0_pReal
prm%minDipDistance = config%getFloat('cedgedipmindistance') * prm%burgers
prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-key
prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key
! expand: family => system
prm%rho0 = math_expand(prm%rho0, prm%Nslip)
@ -313,8 +311,6 @@ subroutine plastic_disloUCLA_init()
outputID = merge(mfp_ID,undefined_ID,prm%totalNslip>0_pInt)
case ('threshold_stress','threshold_stress_slip')
outputID = merge(thresholdstress_ID,undefined_ID,prm%totalNslip>0_pInt)
case ('edge_dipole_distance')
outputID = merge(dipoleDistance_ID,undefined_ID,prm%totalNslip>0_pInt)
end select
@ -336,24 +332,6 @@ subroutine plastic_disloUCLA_init()
prm%totalNslip,0_pInt,0_pInt)
plasticState(p)%sizePostResults = sum(plastic_disloUCLA_sizePostResult(:,phase_plasticityInstance(p)))
allocate(prm%forestProjectionEdge(prm%totalNslip,prm%totalNslip),source = 0.0_pReal)
i = 0_pInt
mySlipFamilies: do f = 1_pInt,size(prm%Nslip,1)
index_myFamily = sum(prm%Nslip(1:f-1_pInt))
slipSystemsLoop: do j = 1_pInt,prm%Nslip(f)
i = i + 1_pInt
do o = 1_pInt, size(prm%Nslip,1)
index_otherFamily = sum(prm%Nslip(1:o-1_pInt))
do k = 1_pInt,prm%Nslip(o) ! loop over (active) systems in other family (slip)
prm%forestProjectionEdge(index_myFamily+j,index_otherFamily+k) = &
abs(math_mul3x3(lattice_sn(:,sum(lattice_NslipSystem(1:f-1,p))+j,p), &
lattice_st(:,sum(lattice_NslipSystem(1:o-1,p))+k,p)))
enddo; enddo
enddo slipSystemsLoop
enddo mySlipFamilies
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1_pInt
@ -374,7 +352,7 @@ subroutine plastic_disloUCLA_init()
endIndex = endIndex + prm%totalNslip
stt%accshear=>plasticState(p)%state(startIndex:endIndex,:)
dot%accshear=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = 1e6_pReal !ToDo: better make optional parameter
plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:)
@ -579,16 +557,6 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe
postResults(c+1_pInt:c+prm%totalNslip) = dst%mfp(1_pInt:prm%totalNslip, of)
case (thresholdstress_ID)
postResults(c+1_pInt:c+prm%totalNslip) = dst%threshold_stress(1_pInt:prm%totalNslip,of)
case (dipoleDistance_ID) ! ToDo: Discuss required changes with Franz
do i = 1_pInt, prm%totalNslip
if (dNeq0(abs(math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))))) then
postResults(c+i) = (3.0_pReal*prm%mu*prm%burgers(i)) &
/ (16.0_pReal*pi*abs(math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))))
else
postResults(c+i) = huge(1.0_pReal)
endif
postResults(c+i)=min(postResults(c+i),dst%mfp(i,of))
enddo
end select

File diff suppressed because it is too large Load Diff

View File

@ -151,7 +151,6 @@ subroutine plastic_kinehardening_init
outputID
character(len=pStringLen) :: &
structure = '',&
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
@ -187,8 +186,6 @@ subroutine plastic_kinehardening_init
endif
#endif
structure = config%getString('lattice_structure')
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
@ -203,28 +200,29 @@ subroutine plastic_kinehardening_init
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
else
prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid
endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
structure(1:3))
config%getString('lattice_structure'))
prm%crss0 = config%getFloats('crss0', requiredShape=shape(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredShape=shape(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredShape=shape(prm%Nslip))
prm%theta0 = config%getFloats('theta0', requiredShape=shape(prm%Nslip))
prm%theta1 = config%getFloats('theta1', requiredShape=shape(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredShape=shape(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredShape=shape(prm%Nslip))
prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip))
prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip))
prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip))
prm%gdot0 = config%getFloat('gdot0')
prm%n = config%getFloat('n_slip')
@ -302,7 +300,6 @@ subroutine plastic_kinehardening_init
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState, &
prm%totalNslip,0_pInt,0_pInt)
plasticState(p)%sizePostResults = sum(plastic_kinehardening_sizePostResult(:,phase_plasticityInstance(p)))
plasticState(p)%offsetDeltaState = sizeDotState
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and aTolState

View File

@ -153,7 +153,6 @@ subroutine plastic_phenopowerlaw_init
outputID
character(len=pStringLen) :: &
structure = '',&
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
@ -181,8 +180,6 @@ subroutine plastic_phenopowerlaw_init
stt => state(phase_plasticityInstance(p)), &
config => config_phase(p))
structure = config%getString('lattice_structure')
!--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined
prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal)
@ -204,30 +201,31 @@ subroutine plastic_phenopowerlaw_init
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
else
prm%nonSchmid_pos = prm%Schmid_slip
prm%nonSchmid_neg = prm%Schmid_slip
prm%nonSchmid_pos = prm%Schmid_slip
prm%nonSchmid_neg = prm%Schmid_slip
endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
structure(1:3))
config%getString('lattice_structure'))
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), &
defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), &
defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip')
prm%a_slip = config%getFloat('a_slip')
prm%h0_SlipSlip = config%getFloat('h0_slipslip')
prm%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip')
prm%a_slip = config%getFloat('a_slip')
prm%h0_SlipSlip = config%getFloat('h0_slipslip')
! expand: family => system
prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip)
@ -250,12 +248,12 @@ subroutine plastic_phenopowerlaw_init
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%totalNtwin = sum(prm%Ntwin)
twinActive: if (prm%totalNtwin > 0_pInt) then
prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,structure(1:3),&
prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,&
config%getFloats('interaction_twintwin'), &
structure(1:3))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),&
config%getString('lattice_structure'))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))
@ -282,10 +280,10 @@ subroutine plastic_phenopowerlaw_init
slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then
prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,&
config%getFloats('interaction_sliptwin'), &
structure(1:3))
config%getString('lattice_structure'))
prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,&
config%getFloats('interaction_twinslip'), &
structure(1:3))
config%getString('lattice_structure'))
else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0