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

:Merge branch '53-separate-mesh-for-different-solvers-3' of magit1.mpie.de:damask/DAMASK into 53-separate-mesh-for-different-solvers-3

This commit is contained in:
navyanthkusam 2019-01-28 14:25:16 +01:00
parent 5101a3796f 9260af2eaa
commit 848a81fd39
15 changed files with 1970 additions and 2708 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

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

2
VERSION
View File

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

2
src/CMakeLists.txt
View File

@ -62,7 +62,7 @@ add_dependencies(MATH FEsolving)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATH>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATH>)
add_library(MESH_BASE OBJECT "mesh_base.f90") add_library(MESH_BASE OBJECT "mesh_base.f90")
add_dependencies(MESH_BASE MATH) add_dependencies(MESH_BASE MATH ELEMENT)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH_BASE>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH_BASE>)
# SPECTRAL solver and FEM solver use different mesh files # SPECTRAL solver and FEM solver use different mesh files

4
src/FEM_zoo.f90
View File

@ -9,11 +9,11 @@ module FEM_Zoo
private private
integer(pInt), parameter, public:: & integer(pInt), parameter, public:: &
maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary) maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary)
real(pReal), dimension(2,3), private, protected :: & real(pReal), dimension(2,3), private, parameter :: &
triangle = reshape([-1.0_pReal, -1.0_pReal, & triangle = reshape([-1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, & 1.0_pReal, -1.0_pReal, &
-1.0_pReal, 1.0_pReal], shape=[2,3]) -1.0_pReal, 1.0_pReal], shape=[2,3])
real(pReal), dimension(3,4), private, protected :: & real(pReal), dimension(3,4), private, parameter :: &
tetrahedron = reshape([-1.0_pReal, -1.0_pReal, -1.0_pReal, & tetrahedron = reshape([-1.0_pReal, -1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, -1.0_pReal, & 1.0_pReal, -1.0_pReal, -1.0_pReal, &
-1.0_pReal, 1.0_pReal, -1.0_pReal, & -1.0_pReal, 1.0_pReal, -1.0_pReal, &

4
src/IO.f90
View File

@ -1236,6 +1236,10 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'zero entry on stiffness diagonal' msg = 'zero entry on stiffness diagonal'
case (136_pInt) case (136_pInt)
msg = 'zero entry on stiffness diagonal for transformed phase' 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 ! errors related to the parsing of material.config

11
src/config.f90
View File

@ -550,7 +550,7 @@ end function getString
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all !> @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. !! 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: & use IO, only: &
IO_error, & IO_error, &
IO_stringValue, & IO_stringValue, &
@ -561,7 +561,6 @@ function getFloats(this,key,defaultVal,requiredShape,requiredSize)
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
real(pReal), dimension(:), intent(in), optional :: defaultVal 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 integer(pInt), intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
integer(pInt) :: i integer(pInt) :: i
@ -601,7 +600,7 @@ end function getFloats
!> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all !> @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. !! 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: & use IO, only: &
IO_error, & IO_error, &
IO_stringValue, & IO_stringValue, &
@ -611,8 +610,7 @@ function getInts(this,key,defaultVal,requiredShape,requiredSize)
integer(pInt), dimension(:), allocatable :: getInts integer(pInt), dimension(:), allocatable :: getInts
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
integer(pInt), dimension(:), intent(in), optional :: defaultVal, & integer(pInt), dimension(:), intent(in), optional :: defaultVal
requiredShape ! not useful (is always 1D array)
integer(pInt), intent(in), optional :: requiredSize integer(pInt), intent(in), optional :: requiredSize
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
integer(pInt) :: i 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. !! 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 !! 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: & use IO, only: &
IO_error, & IO_error, &
IO_StringValue IO_StringValue
@ -663,7 +661,6 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
class(tPartitionedStringList), target, intent(in) :: this class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key character(len=*), intent(in) :: key
character(len=65536),dimension(:), intent(in), optional :: defaultVal character(len=65536),dimension(:), intent(in), optional :: defaultVal
integer(pInt), dimension(:), intent(in), optional :: requiredShape
logical, intent(in), optional :: raw logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
character(len=65536) :: str character(len=65536) :: str

58
src/constitutive.f90
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_ISOTROPIC_ID)) call plastic_isotropic_init
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_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_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_DISLOUCLA_ID)) call plastic_disloucla_init
if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then
call plastic_nonlocal_init(FILEUNIT) call plastic_nonlocal_init(FILEUNIT)
@ -365,7 +365,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
use plastic_nonlocal, only: & use plastic_nonlocal, only: &
plastic_nonlocal_microstructure plastic_nonlocal_microstructure
use plastic_dislotwin, only: & use plastic_dislotwin, only: &
plastic_dislotwin_microstructure plastic_dislotwin_dependentState
use plastic_disloUCLA, only: & use plastic_disloUCLA, only: &
plastic_disloUCLA_dependentState 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))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType 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 case (PLASTICITY_DISLOUCLA_ID) plasticityType
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(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 pReal
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_Mandel6to33, & math_6toSym33, &
math_Mandel33to6, & math_sym33to6, &
math_Plain99to3333 math_99to3333
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticity, & phase_plasticity, &
@ -470,7 +472,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
S = math_Mandel6to33(S6) S = math_6toSym33(S6)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) 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) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType 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) 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 case (PLASTICITY_DISLOTWIN_ID) plasticityType
of = phasememberAt(ipc,ip,el) 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_inv33, &
math_det33, & math_det33, &
math_mul33x33, & math_mul33x33, &
math_Mandel6to33 math_6toSym33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticity, & phase_plasticity, &
@ -597,7 +599,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
case (PLASTICITY_isotropic_ID) plasticityType case (PLASTICITY_isotropic_ID) plasticityType
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(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 case default plasticityType
my_Li = 0.0_pReal my_Li = 0.0_pReal
my_dLi_dS = 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 : & use math, only : &
math_mul33x33, & math_mul33x33, &
math_mul3333xx33, & math_mul3333xx33, &
math_Mandel66to3333, & math_66toSym3333, &
math_I3 math_I3
use material, only: & use material, only: &
material_phase, & material_phase, &
@ -749,7 +751,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip
i, j i, j
ho = material_homogenizationAt(el) 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)) DegradationLoop: do d = 1_pInt, phase_NstiffnessDegradations(material_phase(ipc,ip,el))
degradationType: select case(phase_stiffnessDegradation(d,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 debug_levelBasic
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_Mandel6to33, & math_6toSym33, &
math_Mandel33to6, & math_sym33to6, &
math_mul33x33 math_mul33x33
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
@ -860,7 +862,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,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))) 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) call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType 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) subdt,subfracArray,ip,el)
end select plasticityType end select plasticityType
@ -920,7 +922,7 @@ end subroutine constitutive_collectDotState
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> 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: & use prec, only: &
pReal, & pReal, &
pLongInt pLongInt
@ -929,8 +931,7 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use math, only: & use math, only: &
math_Mandel6to33, & math_sym33to6, &
math_Mandel33to6, &
math_mul33x33 math_mul33x33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
@ -954,18 +955,17 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola Kirchhoff stress
Fe, & !< elastic deformation gradient Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Mp Mp
integer(pInt) :: & integer(pInt) :: &
s, & !< counter in source loop i, &
instance, of 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))) 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) call plastic_kinehardening_deltaState(Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType 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 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 case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & 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: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
math_Mandel6to33, & math_6toSym33, &
math_mul33x33 math_mul33x33
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
@ -1076,7 +1076,7 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
constitutive_postResults = 0.0_pReal 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) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)

515
src/crystallite.f90
View File

@ -1654,29 +1654,7 @@ subroutine integrateStateFPI()
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
call update_dependentState call update_dependentState
!$OMP PARALLEL call update_stress(1.0_pReal)
! --- 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_dotState(1.0_pReal) call update_dotState(1.0_pReal)
!$OMP PARALLEL !$OMP PARALLEL
! --- UPDATE STATE --- ! --- UPDATE STATE ---
@ -1874,17 +1852,6 @@ end subroutine integrateStateFPI
subroutine integrateStateEuler() subroutine integrateStateEuler()
use, intrinsic :: & use, intrinsic :: &
IEEE_arithmetic 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: & use mesh, only: &
mesh_element, & mesh_element, &
mesh_NcpElems mesh_NcpElems
@ -1896,7 +1863,6 @@ subroutine integrateStateEuler()
constitutive_microstructure constitutive_microstructure
implicit none implicit none
integer(pInt) :: & integer(pInt) :: &
e, & ! element index in element loop e, & ! element index in element loop
i, & ! integration point index in ip 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_dotState(1.0_pReal)
call update_State(1.0_pReal) call update_State(1.0_pReal)
call update_deltaState
!$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 call update_dependentState
call update_stress(1.0_pReal)
call setConvergenceFlag
!$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
! --- CHECK NON-LOCAL CONVERGENCE --- ! --- CHECK NON-LOCAL CONVERGENCE ---
@ -2109,74 +2022,13 @@ subroutine integrateStateAdaptiveEuler()
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP 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 !$OMP END PARALLEL
call update_deltaState
call update_dependentState call update_dependentState
call update_stress(1.0_pReal)
! --- 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_dotState(1.0_pReal)
!$OMP PARALLEL !$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) --- ! --- ERROR ESTIMATE FOR STATE (HEUN METHOD) ---
!$OMP SINGLE !$OMP SINGLE
@ -2365,46 +2217,9 @@ subroutine integrateStateRK4()
!$OMP END PARALLEL !$OMP END PARALLEL
call update_state(TIMESTEPFRACTION(n)) call update_state(TIMESTEPFRACTION(n))
call update_deltaState
!$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 call update_dependentState
!$OMP PARALLEL call update_stress(TIMESTEPFRACTION(n))
! --- 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
! --- dot state and RK dot state--- ! --- dot state and RK dot state---
@ -2414,14 +2229,7 @@ subroutine integrateStateRK4()
enddo enddo
call setConvergenceFlag
! --- 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
! --- CHECK NONLOCAL CONVERGENCE --- ! --- CHECK NONLOCAL CONVERGENCE ---
@ -2584,50 +2392,9 @@ subroutine integrateStateRKCK45()
!$OMP END PARALLEL !$OMP END PARALLEL
call update_state(1.0_pReal) !MD: 1.0 correct? call update_state(1.0_pReal) !MD: 1.0 correct?
call update_deltaState
!$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 call update_dependentState
!$OMP PARALLEL call update_stress(C(stage))
! --- 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_dotState(C(stage))
enddo enddo
@ -2735,56 +2502,12 @@ subroutine integrateStateRKCK45()
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
!$OMP END PARALLEL
call update_deltaState
! --- 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 call update_dependentState
call update_stress(1.0_pReal)
!$OMP PARALLEL call setConvergenceFlag
!--------------------------------------------------------------------------------------------------
! --- 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
! --- nonlocal convergence check --- ! --- nonlocal convergence check ---
@ -2798,6 +2521,67 @@ subroutine integrateStateRKCK45()
end 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 !> @brief tbd
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -2824,6 +2608,7 @@ subroutine update_dependentState()
end subroutine update_dependentState end subroutine update_dependentState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t) !> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -2886,7 +2671,6 @@ subroutine update_dotState(timeFraction)
constitutive_collectDotState constitutive_collectDotState
implicit none implicit none
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeFraction timeFraction
integer(pInt) :: & integer(pInt) :: &
@ -2897,15 +2681,17 @@ subroutine update_dotState(timeFraction)
c, & c, &
s s
logical :: & logical :: &
NaN NaN, &
nonlocalStop
!$OMP PARALLEL nonlocalStop = .false.
!$OMP DO PRIVATE (p,c,NaN)
!$OMP PARALLEL DO PRIVATE (p,c,NaN)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e)) do g = 1,homogenization_Ngrains(mesh_element(3,e))
!$OMP FLUSH(crystallite_todo) !$OMP FLUSH(nonlocalStop)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then 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), & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, & crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
@ -2918,20 +2704,94 @@ subroutine update_dotState(timeFraction)
enddo enddo
if (NaN) then if (NaN) then
crystallite_todo(g,i,e) = .false. ! this one done (and broken) crystallite_todo(g,i,e) = .false. ! this one done (and broken)
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a local... if (.not. crystallite_localPlasticity(g,i,e)) nonlocalStop = .True.
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
!$OMP END CRITICAL (checkTodo)
endif
endif endif
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP END PARALLEL DO
!$OMP END PARALLEL
if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity
end subroutine update_DotState 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 !> @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 !> 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 phaseAt, phasememberAt
use constitutive, only: & use constitutive, only: &
constitutive_collectDeltaState constitutive_collectDeltaState
use math, only: &
math_6toSym33
implicit none implicit none
integer(pInt), intent(in):: & integer(pInt), intent(in):: &
@ -2969,57 +2831,50 @@ logical function stateJump(ipc,ip,el)
c, & c, &
p, & p, &
mySource, & mySource, &
myOffsetPlasticDeltaState, & myOffset, &
myOffsetSourceDeltaState, & mySize
mySizePlasticDeltaState, &
mySizeSourceDeltaState
c = phasememberAt(ipc,ip,el) c = phasememberAt(ipc,ip,el)
p = phaseAt(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_Fe(1:3,1:3,ipc,ip,el), &
crystallite_Fi(1:3,1:3,ipc,ip,el), & crystallite_Fi(1:3,1:3,ipc,ip,el), &
ipc,ip,el) ipc,ip,el)
myOffsetPlasticDeltaState = plasticState(p)%offsetDeltaState myOffset = plasticState(p)%offsetDeltaState
mySizePlasticDeltaState = plasticState(p)%sizeDeltaState 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. stateJump = .false.
return return
endif endif
plasticState(p)%state(myOffsetPlasticDeltaState + 1_pInt : & plasticState(p)%state(myOffset + 1_pInt:myOffset + mySize,c) = &
myOffsetPlasticDeltaState + mySizePlasticDeltaState,c) = & plasticState(p)%state(myOffset + 1_pInt:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c)
plasticState(p)%state(myOffsetPlasticDeltaState + 1_pInt : &
myOffsetPlasticDeltaState + mySizePlasticDeltaState,c) + &
plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)
do mySource = 1_pInt, phase_Nsources(p) do mySource = 1_pInt, phase_Nsources(p)
myOffsetSourceDeltaState = sourceState(p)%p(mySource)%offsetDeltaState myOffset = sourceState(p)%p(mySource)%offsetDeltaState
mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState mySize = sourceState(p)%p(mySource)%sizeDeltaState
if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)))) then ! NaN occured in deltaState if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then ! NaN occured in deltaState
stateJump = .false. stateJump = .false.
return return
endif endif
sourceState(p)%p(mySource)%state(myOffsetSourceDeltaState + 1_pInt : & sourceState(p)%p(mySource)%state(myOffset + 1_pInt: myOffset + mySize,c) = &
myOffsetSourceDeltaState + mySizeSourceDeltaState,c) = & sourceState(p)%p(mySource)%state(myOffset + 1_pInt: myOffset + mySize,c) + &
sourceState(p)%p(mySource)%state(myOffsetSourceDeltaState + 1_pInt : & sourceState(p)%p(mySource)%deltaState(1:mySize,c)
myOffsetSourceDeltaState + mySizeSourceDeltaState,c) + &
sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)
enddo enddo
#ifdef DEBUG #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. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .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,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', & write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state', &
plasticState(p)%state(myOffsetPlasticDeltaState + 1_pInt : & plasticState(p)%state(myOffset + 1_pInt : &
myOffsetPlasticDeltaState + mySizePlasticDeltaState,c) myOffset + mySize,c)
endif endif
#endif #endif

2311
src/lattice.f90
View File

File diff suppressed because it is too large Load Diff

4
src/material.f90
View File

@ -918,7 +918,8 @@ end subroutine material_parseTexture
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates the plastic state of a phase !> @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) Nslip,Ntwin,Ntrans)
use numerics, only: & use numerics, only: &
numerics_integrator2 => numerics_integrator ! compatibility hack numerics_integrator2 => numerics_integrator ! compatibility hack
@ -939,6 +940,7 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,sizeState,sizeDotState
plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeState = sizeState
plasticState(phase)%sizeDotState = sizeDotState plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizeDeltaState = sizeDeltaState plasticState(phase)%sizeDeltaState = sizeDeltaState
plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
plasticState(phase)%Nslip = Nslip plasticState(phase)%Nslip = Nslip
plasticState(phase)%Ntwin = Ntwin plasticState(phase)%Ntwin = Ntwin
plasticState(phase)%Ntrans= Ntrans plasticState(phase)%Ntrans= Ntrans

34
src/mesh_FEM.f90
View File

@ -12,7 +12,7 @@ module mesh
#include <petsc/finclude/petscis.h> #include <petsc/finclude/petscis.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use prec, only: pReal, pInt use prec, only: pReal, pInt
use mesh_base
use PETScdmplex use PETScdmplex
use PETScdmda use PETScdmda
use PETScis use PETScis
@ -80,6 +80,17 @@ use PETScis
integer(pInt), dimension(1_pInt), parameter, public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type integer(pInt), dimension(1_pInt), parameter, public :: FE_NipNeighbors = & !< number of ip neighbors / cell faces in a specific cell type
int([6],pInt) int([6],pInt)
type, public, extends(tMesh) :: tMesh_FEM
contains
procedure :: init => tMesh_FEM_init
end type tMesh_FEM
type(tMesh_FEM), public, protected :: theMesh
public :: & public :: &
mesh_init, & mesh_init, &
@ -89,6 +100,23 @@ use PETScis
contains contains
subroutine tMesh_FEM_init(self,dimen,order)
implicit none
integer(pInt), intent(in) :: dimen,order
class(tMesh_FEM) :: self
if (dimen == 2_pInt) then
if (order == 1_pInt) call self%elem%init(1_pInt)
if (order == 2_pInt) call self%elem%init(2_pInt)
elseif(dimen == 3_pInt) then
if (order == 1_pInt) call self%elem%init(6_pInt)
if (order == 2_pInt) call self%elem%init(8_pInt)
endif
end subroutine tMesh_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief initializes the mesh by calling all necessary private routines the mesh module !> @brief initializes the mesh by calling all necessary private routines the mesh module
@ -213,6 +241,8 @@ subroutine mesh_init()
FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder) FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
mesh_maxNips = FE_Nips(1_pInt) mesh_maxNips = FE_Nips(1_pInt)
write(6,*) 'mesh_maxNips',mesh_maxNips
call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
call mesh_FEM_build_ipVolumes(dimPlex) call mesh_FEM_build_ipVolumes(dimPlex)
@ -243,7 +273,7 @@ subroutine mesh_init()
mesh_homogenizationAt = mesh_element(3,:) mesh_homogenizationAt = mesh_element(3,:)
mesh_microstructureAt = mesh_element(4,:) mesh_microstructureAt = mesh_element(4,:)
!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!
call theMesh%init(dimplex,integrationOrder)
end subroutine mesh_init end subroutine mesh_init

84
src/plastic_disloUCLA.f90
View File

@ -28,8 +28,7 @@ module plastic_disloUCLA
shearrate_ID, & shearrate_ID, &
accumulatedshear_ID, & accumulatedshear_ID, &
mfp_ID, & mfp_ID, &
thresholdstress_ID, & thresholdstress_ID
dipoledistance_ID
end enum end enum
type, private :: tParameters type, private :: tParameters
@ -73,7 +72,7 @@ module plastic_disloUCLA
integer(kind(undefined_ID)), allocatable, dimension(:) :: & integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID !< ID of each post result output outputID !< ID of each post result output
logical :: & logical :: &
dipoleformation dipoleFormation !< flag indicating consideration of dipole formation
end type !< container type for internal constitutive parameters end type !< container type for internal constitutive parameters
type, private :: tDisloUCLAState type, private :: tDisloUCLAState
@ -93,7 +92,7 @@ module plastic_disloUCLA
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! containers for parameters and state ! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param type(tParameters), allocatable, dimension(:), private :: param
type(tDisloUCLAState ), allocatable, dimension(:), private :: & type(tDisloUCLAState), allocatable, dimension(:), private :: &
dotState, & dotState, &
state state
type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState
@ -127,7 +126,6 @@ subroutine plastic_disloUCLA_init()
debug_constitutive,& debug_constitutive,&
debug_levelBasic debug_levelBasic
use math, only: & use math, only: &
math_mul3x3, &
math_expand math_expand
use IO, only: & use IO, only: &
IO_error, & IO_error, &
@ -148,8 +146,6 @@ subroutine plastic_disloUCLA_init()
implicit none implicit none
integer(pInt) :: & integer(pInt) :: &
index_myFamily, index_otherFamily, &
f,j,k,o, &
Ninstance, & Ninstance, &
p, i, & p, i, &
NipcMyPhase, & NipcMyPhase, &
@ -164,7 +160,6 @@ subroutine plastic_disloUCLA_init()
outputID outputID
character(len=pStringLen) :: & character(len=pStringLen) :: &
structure = '',&
extmsg = '' extmsg = ''
character(len=65536), dimension(:), allocatable :: & character(len=65536), dimension(:), allocatable :: &
outputs outputs
@ -197,8 +192,6 @@ subroutine plastic_disloUCLA_init()
dst => dependentState(phase_plasticityInstance(p)), & dst => dependentState(phase_plasticityInstance(p)), &
config => config_phase(p)) config => config_phase(p))
structure = config%getString('lattice_structure')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined ! optional parameters that need to be defined
prm%mu = lattice_mu(p) prm%mu = lattice_mu(p)
@ -213,9 +206,10 @@ subroutine plastic_disloUCLA_init()
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip) prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray) defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
@ -224,25 +218,29 @@ subroutine plastic_disloUCLA_init()
prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid prm%nonSchmid_neg = prm%Schmid
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), & config%getFloats('interaction_slipslip'), &
structure(1:3)) config%getString('lattice_structure'))
prm%rho0 = config%getFloats('rhoedge0', requiredShape=shape(prm%Nslip)) prm%forestProjectionEdge = lattice_forestProjection(prm%Nslip,config%getString('lattice_structure'),&
prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredShape=shape(prm%Nslip)) config%getFloat('c/a',defaultVal=0.0_pReal))
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))
prm%clambda = config%getFloats('clambdaslip', requiredShape=shape(prm%Nslip)) prm%rho0 = config%getFloats('rhoedge0', requiredSize=size(prm%Nslip))
prm%tau_Peierls = config%getFloats('tau_peierls', requiredShape=shape(prm%Nslip)) ! ToDo: Deprecated prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%Nslip))
prm%p = config%getFloats('p_slip', requiredShape=shape(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))]) 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))]) defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%kink_height = config%getFloats('kink_height', requiredShape=shape(prm%Nslip)) prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%Nslip))
prm%w = config%getFloats('kink_width', requiredShape=shape(prm%Nslip)) prm%w = config%getFloats('kink_width', requiredSize=size(prm%Nslip))
prm%omega = config%getFloats('omega', requiredShape=shape(prm%Nslip)) prm%omega = config%getFloats('omega', requiredSize=size(prm%Nslip))
prm%B = config%getFloats('friction_coeff', requiredShape=shape(prm%Nslip)) prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%Nslip))
prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! ToDo: Deprecated prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! ToDo: Deprecated
prm%grainSize = config%getFloat('grainsize') prm%grainSize = config%getFloat('grainsize')
@ -250,7 +248,7 @@ subroutine plastic_disloUCLA_init()
prm%Qsd = config%getFloat('qsd') prm%Qsd = config%getFloat('qsd')
prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers**3.0_pReal prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers**3.0_pReal
prm%minDipDistance = config%getFloat('cedgedipmindistance') * prm%burgers 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 ! expand: family => system
prm%rho0 = math_expand(prm%rho0, prm%Nslip) 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) outputID = merge(mfp_ID,undefined_ID,prm%totalNslip>0_pInt)
case ('threshold_stress','threshold_stress_slip') case ('threshold_stress','threshold_stress_slip')
outputID = merge(thresholdstress_ID,undefined_ID,prm%totalNslip>0_pInt) 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 end select
@ -336,24 +332,6 @@ subroutine plastic_disloUCLA_init()
prm%totalNslip,0_pInt,0_pInt) prm%totalNslip,0_pInt,0_pInt)
plasticState(p)%sizePostResults = sum(plastic_disloUCLA_sizePostResult(:,phase_plasticityInstance(p))) 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 ! locally defined state aliases and initialization of state0 and aTolState
startIndex = 1_pInt startIndex = 1_pInt
@ -374,7 +352,7 @@ subroutine plastic_disloUCLA_init()
endIndex = endIndex + prm%totalNslip endIndex = endIndex + prm%totalNslip
stt%accshear=>plasticState(p)%state(startIndex:endIndex,:) stt%accshear=>plasticState(p)%state(startIndex:endIndex,:)
dot%accshear=>plasticState(p)%dotState(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 ! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%accumulatedSlip => plasticState(p)%state(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) postResults(c+1_pInt:c+prm%totalNslip) = dst%mfp(1_pInt:prm%totalNslip, of)
case (thresholdstress_ID) case (thresholdstress_ID)
postResults(c+1_pInt:c+prm%totalNslip) = dst%threshold_stress(1_pInt:prm%totalNslip,of) 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 end select

1414
src/plastic_dislotwin.f90
View File

File diff suppressed because it is too large Load Diff

25
src/plastic_kinematichardening.f90
View File

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

20
src/plastic_phenopowerlaw.f90
View File

@ -153,7 +153,6 @@ subroutine plastic_phenopowerlaw_init
outputID outputID
character(len=pStringLen) :: & character(len=pStringLen) :: &
structure = '',&
extmsg = '' extmsg = ''
character(len=65536), dimension(:), allocatable :: & character(len=65536), dimension(:), allocatable :: &
outputs outputs
@ -181,8 +180,6 @@ subroutine plastic_phenopowerlaw_init
stt => state(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), &
config => config_phase(p)) config => config_phase(p))
structure = config%getString('lattice_structure')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined ! optional parameters that need to be defined
prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal) prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal)
@ -204,9 +201,10 @@ subroutine plastic_phenopowerlaw_init
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip) prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray) defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
@ -217,7 +215,7 @@ subroutine plastic_phenopowerlaw_init
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), & 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_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
@ -250,12 +248,12 @@ subroutine plastic_phenopowerlaw_init
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%totalNtwin = sum(prm%Ntwin) prm%totalNtwin = sum(prm%Ntwin)
twinActive: if (prm%totalNtwin > 0_pInt) then 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)) config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,& prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,&
config%getFloats('interaction_twintwin'), & config%getFloats('interaction_twintwin'), &
structure(1:3)) config%getString('lattice_structure'))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),& prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a')) config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) 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 slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then
prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,& prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,&
config%getFloats('interaction_sliptwin'), & config%getFloats('interaction_sliptwin'), &
structure(1:3)) config%getString('lattice_structure'))
prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,& prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,&
config%getFloats('interaction_twinslip'), & config%getFloats('interaction_twinslip'), &
structure(1:3)) config%getString('lattice_structure'))
else slipAndTwinActive else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0 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 allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0