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

Merge branch 'nonlocal-standard-access' into 'development' 

Using consistent access pattern

See merge request damask/DAMASK!520
This commit is contained in:
Sharan Roongta 2022-02-13 17:33:35 +00:00
parent 5f7a18c32d b43695067d
commit 26979da585
10 changed files with 253 additions and 300 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/CPFEM.f90
View File

@ -192,11 +192,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
else validCalculation
if (debugCPFEM%extensive) print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP])
call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip,(elCP-1)*discretization_nIPs + ip)
if (.not. terminallyIll) &
call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP])
terminalIllness: if (terminallyIll) then
call random_number(rnd)

4
src/grid/spectral_utilities.f90
View File

@ -812,9 +812,9 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field
call homogenization_mechanical_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) ! calculate P field
call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3) ! calculate P field
if (.not. terminallyIll) &
call homogenization_thermal_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3)
if (.not. terminallyIll) &
call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3])

106
src/homogenization.f90
View File

@ -222,14 +222,13 @@ end subroutine homogenization_init
!--------------------------------------------------------------------------------------------------
!> @brief
!--------------------------------------------------------------------------------------------------
subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
subroutine homogenization_mechanical_response(Delta_t,cell_start,cell_end)
real(pReal), intent(in) :: Delta_t !< time increment
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
integer, intent(in) :: &
cell_start, cell_end
integer :: &
NiterationMPstate, &
ip, & !< integration point number
el, & !< element number
co, ce, ho, en
logical :: &
converged
@ -237,45 +236,42 @@ subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving
doneAndHappy
!$OMP PARALLEL DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
!$OMP PARALLEL DO PRIVATE(en,ho,co,NiterationMPstate,converged,doneAndHappy)
do ce = cell_start, cell_end
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip
en = material_homogenizationEntry(ce)
ho = material_homogenizationID(ce)
en = material_homogenizationEntry(ce)
ho = material_homogenizationID(ce)
call phase_restore(ce,.false.) ! wrong name (is more a forward function)
call phase_restore(ce,.false.) ! wrong name (is more a forward function)
if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en)
call damage_partition(ce)
if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en)
call damage_partition(ce)
doneAndHappy = [.false.,.true.]
doneAndHappy = [.false.,.true.]
NiterationMPstate = 0
convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) &
.and. NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1
NiterationMPstate = 0
convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) &
.and. NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1
call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
if (converged) then
doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce)
converged = all(doneAndHappy)
else
doneAndHappy = [.true.,.false.]
endif
enddo convergenceLooping
call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
converged = all([(phase_mechanical_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))])
if (converged) then
doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce)
converged = all(doneAndHappy)
else
doneAndHappy = [.true.,.false.]
end if
end do convergenceLooping
converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ce),co=1,homogenization_Nconstituents(ho))])
if (.not. converged) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
endif
enddo
enddo
if (.not. converged) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
end if
end do
!$OMP END PARALLEL DO
end subroutine homogenization_mechanical_response
@ -284,31 +280,27 @@ end subroutine homogenization_mechanical_response
!--------------------------------------------------------------------------------------------------
!> @brief
!--------------------------------------------------------------------------------------------------
subroutine homogenization_thermal_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end)
real(pReal), intent(in) :: Delta_t !< time increment
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
integer, intent(in) :: &
cell_start, cell_end
integer :: &
ip, & !< integration point number
el, & !< element number
co, ce, ho
!$OMP PARALLEL DO PRIVATE(ho,ce)
do el = FEsolving_execElem(1),FEsolving_execElem(2)
!$OMP PARALLEL DO PRIVATE(ho)
do ce = cell_start, cell_end
if (terminallyIll) continue
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip
ho = material_homogenizationID(ce)
call thermal_partition(ce)
do co = 1, homogenization_Nconstituents(ho)
if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
endif
enddo
enddo
enddo
ho = material_homogenizationID(ce)
call thermal_partition(ce)
do co = 1, homogenization_Nconstituents(ho)
if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then
if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
terminallyIll = .true.
end if
end do
end do
!$OMP END PARALLEL DO
end subroutine homogenization_thermal_response
@ -331,11 +323,11 @@ subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolvin
elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip
ho = material_homogenizationID(ce)
do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el)
ho = material_homogenizationID(ce)
do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el)
enddo
call mechanical_homogenize(Delta_t,ce)
call mechanical_homogenize(Delta_t,ce)
enddo IpLooping3
enddo elementLooping3
!$OMP END PARALLEL DO

81
src/material.f90
View File

@ -17,16 +17,17 @@ module material
implicit none
private
type :: tRotationContainer
type(tRotation), dimension(:), allocatable :: data
end type
type :: tTensorContainer
type, public :: tRotationContainer
type(tRotation), dimension(:), allocatable :: data
end type tRotationContainer
type, public :: tTensorContainer
real(pReal), dimension(:,:,:), allocatable :: data
end type
end type tTensorContainer
type(tRotationContainer), dimension(:), allocatable :: material_O_0
type(tTensorContainer), dimension(:), allocatable :: material_F_i_0
type(tRotationContainer), dimension(:), allocatable, public, protected :: material_O_0
type(tTensorContainer), dimension(:), allocatable, public, protected :: material_F_i_0
integer, dimension(:), allocatable, public, protected :: &
homogenization_Nconstituents !< number of grains in each homogenization
@ -37,20 +38,14 @@ module material
material_name_phase, & !< name of each phase
material_name_homogenization !< name of each homogenization
integer, dimension(:), allocatable, public, protected :: & ! (elem)
material_homogenizationID, & !< per cell TODO: material_ID_homogenization
material_homogenizationEntry !< per cell TODO: material_entry_homogenization
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt, & !< phase ID of each element TODO: remove
material_phaseID, & !< per (constituent,cell) TODO: material_ID_phase
material_phaseEntry !< per (constituent,cell) TODO: material_entry_phase
integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem)
material_phaseMemberAt !TODO: remove
integer, dimension(:), allocatable, public, protected :: & ! (cell)
material_homogenizationID, & ! TODO: rename to material_ID_homogenization
material_homogenizationEntry ! TODO: rename to material_entry_homogenization
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,cell)
material_phaseID, & ! TODO: rename to material_ID_phase
material_phaseEntry ! TODO: rename to material_entry_phase
public :: &
tTensorContainer, &
tRotationContainer, &
material_F_i_0, &
material_O_0, &
material_init
contains
@ -97,11 +92,12 @@ subroutine parse()
counterPhase, &
counterHomogenization
real(pReal) :: &
frac
real(pReal) :: v
integer :: &
el, ip, co, ma, &
h, ce
el, ip, &
ho, ph, &
co, ce, &
ma
materials => config_material%get('material')
phases => config_material%get('phase')
@ -118,51 +114,48 @@ subroutine parse()
#endif
allocate(homogenization_Nconstituents(homogenizations%length))
do h=1, homogenizations%length
homogenization => homogenizations%get(h)
homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents')
do ho=1, homogenizations%length
homogenization => homogenizations%get(ho)
homogenization_Nconstituents(ho) = homogenization%get_asInt('N_constituents')
end do
homogenization_maxNconstituents = maxval(homogenization_Nconstituents)
allocate(counterPhase(phases%length),source=0)
allocate(counterHomogenization(homogenizations%length),source=0)
allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0)
allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0)
allocate(material_homogenizationID(discretization_nIPs*discretization_Nelems),source=0)
allocate(material_homogenizationEntry(discretization_nIPs*discretization_Nelems),source=0)
allocate(material_phaseID(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
do el = 1, discretization_Nelems
material => materials%get(discretization_materialAt(el))
constituents => material%get('constituents')
material => materials%get(discretization_materialAt(el))
ho = homogenizations%getIndex(material%get_asString('homogenization'))
do ip = 1, discretization_nIPs
ce = (el-1)*discretization_nIPs + ip
material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization'))
counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1
material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce))
material_homogenizationID(ce) = ho
counterHomogenization(ho) = counterHomogenization(ho) + 1
material_homogenizationEntry(ce) = counterHomogenization(ho)
end do
frac = 0.0_pReal
v = 0.0_pReal
constituents => material%get('constituents')
do co = 1, constituents%length
constituent => constituents%get(co)
frac = frac + constituent%get_asFloat('v')
v = v + constituent%get_asFloat('v')
material_phaseAt(co,el) = phases%getIndex(constituent%get_asString('phase'))
ph = phases%getIndex(constituent%get_asString('phase'))
do ip = 1, discretization_nIPs
ce = (el-1)*discretization_nIPs + ip
counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1
material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el))
material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el))
material_phaseID(co,ce) = material_phaseAt(co,el)
material_phaseID(co,ce) = ph
counterPhase(ph) = counterPhase(ph) + 1
material_phaseEntry(co,ce) = counterPhase(ph)
end do
end do
if (dNeq(frac,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent')
if (dNeq(v,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent')
end do

2
src/mesh/FEM_utilities.f90
View File

@ -150,7 +150,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
print'(/,1x,a)', '... evaluating constitutive response ......................................'
call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) ! calculate P field
call homogenization_mechanical_response(timeinc,1,mesh_maxNips*mesh_NcpElems) ! calculate P field
if (.not. terminallyIll) &
call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])
cutBack = .false.

35
src/phase.f90
View File

@ -228,15 +228,15 @@ module phase
end function phase_thermal_constitutive
module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
module function phase_damage_constitutive(Delta_t,co,ce) result(converged_)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co, ip, el
integer, intent(in) :: co, ce
logical :: converged_
end function phase_damage_constitutive
module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co, ip, el
integer, intent(in) :: co, ce
logical :: converged_
end function phase_mechanical_constitutive
@ -264,19 +264,18 @@ module phase
real(pReal) :: f
end function phase_f_T
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el)
integer, intent(in) :: &
ph, &
i, &
e
ip, &
el
type(tRotationContainer), dimension(:), intent(in) :: orientation
end subroutine plastic_nonlocal_updateCompatibility
module subroutine plastic_dependentState(co,ip,el)
module subroutine plastic_dependentState(ph,en)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
ph, &
en
end subroutine plastic_dependentState
@ -503,8 +502,8 @@ subroutine crystallite_init()
ce, &
co, & !< counter in integration point component loop
ip, & !< counter in integration point loop
el !< counter in element loop
el, & !< counter in element loop
en, ph
class(tNode), pointer :: &
num_crystallite, &
phases
@ -543,13 +542,15 @@ subroutine crystallite_init()
phases => config_material%get('phase')
!$OMP PARALLEL DO PRIVATE(ce)
!$OMP PARALLEL DO PRIVATE(ce,ph,en)
do el = 1, discretization_Nelems
do ip = 1, discretization_nIPs
ce = (el-1)*discretization_nIPs + ip
do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
en = material_phaseEntry(co,ce)
ph = material_phaseID(co,ce)
call crystallite_orientations(co,ip,el)
call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
call plastic_dependentState(ph,en) ! update dependent state variables to be consistent with basic states
end do
end do
end do
@ -577,8 +578,8 @@ subroutine crystallite_orientations(co,ip,el)
call phase_O(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en))))
if (plasticState(material_phaseAt(1,el))%nonlocal) &
call plastic_nonlocal_updateCompatibility(phase_O,material_phaseAt(1,el),ip,el)
if (plasticState(material_phaseID(1,(el-1)*discretization_nIPs + ip))%nonlocal) &
call plastic_nonlocal_updateCompatibility(phase_O,material_phaseID(1,(el-1)*discretization_nIPs + ip),ip,el)
end subroutine crystallite_orientations

9
src/phase_damage.f90
View File

@ -127,21 +127,20 @@ end subroutine damage_init
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
module function phase_damage_constitutive(Delta_t,co,ce) result(converged_)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
co, &
ip, &
el
ce
logical :: converged_
integer :: &
ph, en
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
ph = material_phaseID(co,ce)
en = material_phaseEntry(co,ce)
converged_ = .not. integrateDamageState(Delta_t,ph,en)

132
src/phase_mechanical.f90
View File

@ -79,15 +79,12 @@ submodule(phase) mechanical
en
end subroutine plastic_isotropic_LiAndItsTangent
module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
module function plastic_dotState(subdt,ph,en) result(dotState)
integer, intent(in) :: &
co, & !< constituent
ip, & !< integration point
el, & !< element
ph, &
en
real(pReal), intent(in) :: &
subdt !< timestep
subdt !< timestep
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function plastic_dotState
@ -365,13 +362,11 @@ end subroutine mechanical_results
!> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction
!--------------------------------------------------------------------------------------------------
function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
real(pReal), dimension(3,3), intent(in) :: F,subFp0,subFi0
real(pReal), intent(in) :: Delta_t
integer, intent(in):: el, & ! element index
ip, & ! integration point index
co ! grain index
integer, intent(in) :: ph, en
real(pReal), dimension(3,3):: Fp_new, & ! plastic deformation gradient at end of timestep
invFp_new, & ! inverse of Fp_new
@ -419,19 +414,13 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
ierr, & ! error indicator for LAPACK
o, &
p, &
ph, &
en, &
jacoCounterLp, &
jacoCounterLi ! counters to check for Jacobian update
logical :: error,broken
broken = .true.
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
call plastic_dependentState(co,ip,el)
call plastic_dependentState(ph,en)
Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en) ! take as first guess
Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en) ! take as first guess
@ -579,37 +568,32 @@ end function integrateStress
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
ph, &
en
logical :: &
broken
integer :: &
NiterationState, & !< number of iterations in state loop
ph, &
en, &
sizeDotState
real(pReal) :: &
zeta
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
r, & ! state residuum
dotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState,2) :: &
dotState_last
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
@ -620,10 +604,10 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1)
dotState_last(1:sizeDotState,1) = dotState
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
if(broken) exit iteration
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) exit iteration
zeta = damper(dotState,dotState_last(1:sizeDotState,1),dotState_last(1:sizeDotState,2))
@ -672,31 +656,26 @@ end function integrateStateFPI
!--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
ph, &
en !< grain index in grain loop
logical :: &
broken
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
integer :: &
ph, &
en, &
sizeDotState
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
@ -706,7 +685,7 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
broken = plastic_deltaState(ph,en)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
end function integrateStateEuler
@ -714,32 +693,27 @@ end function integrateStateEuler
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
ph, &
en
logical :: &
broken
integer :: &
ph, &
en, &
sizeDotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
r, &
dotState
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
@ -751,10 +725,10 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
broken = plastic_deltaState(ph,en)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
if(broken) return
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return
broken = .not. converged(r + 0.5_pReal * dotState * Delta_t, &
@ -767,12 +741,12 @@ end function integrateStateAdaptiveEuler
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the classic Runge Kutta method
!---------------------------------------------------------------------------------------------------
function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co,ip,el
integer, intent(in) :: ph, en
logical :: broken
real(pReal), dimension(3,3), parameter :: &
@ -787,7 +761,7 @@ function integrateStateRK4(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal]
broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C)
broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C)
end function integrateStateRK4
@ -795,12 +769,12 @@ end function integrateStateRK4
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the Cash-Carp method
!---------------------------------------------------------------------------------------------------
function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) result(broken)
function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: co,ip,el
integer, intent(in) :: ph, en
logical :: broken
real(pReal), dimension(5,5), parameter :: &
@ -822,7 +796,7 @@ function integrateStateRKCK45(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) re
13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal]
broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB)
broken = integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB)
end function integrateStateRKCK45
@ -831,7 +805,7 @@ end function integrateStateRKCK45
!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
!! embedded explicit Runge-Kutta method
!--------------------------------------------------------------------------------------------------
function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,DB) result(broken)
function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,ph,en,A,B,C,DB) result(broken)
real(pReal), intent(in),dimension(3,3) :: F_0,F,subFp0,subFi0
real(pReal), intent(in),dimension(:) :: subState0
@ -840,28 +814,23 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
real(pReal), dimension(:), intent(in) :: B, C
real(pReal), dimension(:), intent(in), optional :: DB
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
ph, &
en
logical :: broken
integer :: &
stage, & ! stage index in integration stage loop
n, &
ph, &
en, &
sizeDotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,size(B)) :: &
real(pReal), dimension(plasticState(ph)%sizeDotState,size(B)) :: &
plastic_RKdotState
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
broken = .true.
dotState = plastic_dotState(Delta_t, co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t,ph,en)
if (any(IEEE_is_NaN(dotState))) return
sizeDotState = plasticState(ph)%sizeDotState
@ -879,10 +848,10 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
plasticState(ph)%state(1:sizeDotState,en) = subState0 &
+ dotState * Delta_t
broken = integrateStress(F_0 + (F-F_0) * Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage),co,ip,el)
broken = integrateStress(F_0+(F-F_0)*Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage), ph,en)
if(broken) exit
dotState = plastic_dotState(Delta_t*C(stage), co,ip,el,ph,en)
dotState = plastic_dotState(Delta_t*C(stage), ph,en)
if (any(IEEE_is_NaN(dotState))) exit
enddo
@ -904,7 +873,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
broken = plastic_deltaState(ph,en)
if(broken) return
broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el)
broken = integrateStress(F,subFp0,subFi0,Delta_t,ph,en)
end function integrateStateRK
@ -1006,13 +975,12 @@ end subroutine mechanical_forward
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
module function phase_mechanical_constitutive(Delta_t,co,ce) result(converged_)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
co, &
ip, &
el
ce
logical :: converged_
real(pReal) :: &
@ -1028,16 +996,15 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
subLi0, &
subF0, &
subF
real(pReal), dimension(:), allocatable :: subState0
real(pReal), dimension(plasticState(material_phaseID(co,ce))%sizeState) :: subState0
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
sizeDotState = plasticState(ph)%sizeDotState
ph = material_phaseID(co,ce)
en = material_phaseEntry(co,ce)
subState0 = plasticState(ph)%state0(:,en)
subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
allocate(subState0,source=plasticState(ph)%State0(:,en))
subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en)
@ -1070,7 +1037,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
subStep = num%subStepSizeCryst * subStep
phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0
phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0
phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en) ! why no subS0 ? is S0 of any use?
phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)
if (subStep < 1.0_pReal) then ! actual (not initial) cutback
phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0
phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0
@ -1082,9 +1049,10 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
!--------------------------------------------------------------------------------------------------
! prepare for integration
if (todo) then
sizeDotState = plasticState(ph)%sizeDotState
subF = subF0 &
+ subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el)
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,ph,en)
endif
enddo cutbackLooping

33
src/phase_mechanical_plastic.f90
View File

@ -160,16 +160,14 @@ submodule(phase:mechanical) plastic
dotState
end function dislotungsten_dotState
module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el)
module subroutine nonlocal_dotState(Mp,timestep,ph,en)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
en, &
ip, & !< current integration point
el !< current element number
en
end subroutine nonlocal_dotState
module subroutine dislotwin_dependentState(T,ph,en)
@ -186,12 +184,10 @@ submodule(phase:mechanical) plastic
en
end subroutine dislotungsten_dependentState
module subroutine nonlocal_dependentState(ph, en, ip, el)
module subroutine nonlocal_dependentState(ph,en)
integer, intent(in) :: &
ph, &
en, &
ip, & !< current integration point
el !< current element number
en
end subroutine nonlocal_dependentState
module subroutine plastic_kinehardening_deltaState(Mp,ph,en)
@ -301,12 +297,9 @@ end subroutine plastic_LpAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
module function plastic_dotState(subdt,ph,en) result(dotState)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
en
real(pReal), intent(in) :: &
@ -340,7 +333,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState)
dotState = dislotungsten_dotState(Mp,ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType
call nonlocal_dotState(Mp,subdt,ph,en,ip,el)
call nonlocal_dotState(Mp,subdt,ph,en)
dotState = plasticState(ph)%dotState(:,en)
end select plasticType
@ -352,21 +345,13 @@ end function plastic_dotState
!--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different plasticity constitutive models
!--------------------------------------------------------------------------------------------------
module subroutine plastic_dependentState(co, ip, el)
module subroutine plastic_dependentState(ph,en)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el !< element
integer :: &
ph, &
en
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
plasticType: select case (phase_plasticity(ph))
case (PLASTIC_DISLOTWIN_ID) plasticType
@ -376,7 +361,7 @@ module subroutine plastic_dependentState(co, ip, el)
call dislotungsten_dependentState(ph,en)
case (PLASTIC_NONLOCAL_ID) plasticType
call nonlocal_dependentState(ph,en,ip,el)
call nonlocal_dependentState(ph,en)
end select plasticType
@ -422,7 +407,7 @@ module function plastic_deltaState(ph, en) result(broken)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
if (.not. broken) then
mySize = plasticState(ph)%sizeDeltaState
mySize = plasticState(ph)%sizeDeltaState
plasticState(ph)%deltaState2(1:mySize,en) = plasticState(ph)%deltaState2(1:mySize,en) &
+ plasticState(ph)%deltaState(1:mySize,en)
end if

148
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -15,6 +15,9 @@ submodule(phase:plastic) nonlocal
type :: tGeometry
real(pReal), dimension(:), allocatable :: V_0
integer, dimension(:,:,:), allocatable :: IPneighborhood
real(pReal), dimension(:,:), allocatable :: IParea, IPcoordinates
real(pReal), dimension(:,:,:), allocatable :: IPareaNormal
end type tGeometry
type(tGeometry), dimension(:), allocatable :: geom
@ -123,8 +126,10 @@ submodule(phase:plastic) nonlocal
type :: tNonlocalDependentState
real(pReal), allocatable, dimension(:,:) :: &
tau_pass, &
tau_Back
tau_pass, &
tau_Back
real(pReal), allocatable, dimension(:,:,:,:,:) :: &
compatibility
end type tNonlocalDependentState
type :: tNonlocalState
@ -160,7 +165,6 @@ submodule(phase:plastic) nonlocal
state0
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters
type(tNonlocalDependentState), dimension(:), allocatable :: dependentState
contains
@ -406,6 +410,10 @@ module function plastic_nonlocal_init() result(myPlasticity)
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention
allocate(geom(ph)%V_0(Nmembers))
allocate(geom(ph)%IPneighborhood(3,nIPneighbors,Nmembers))
allocate(geom(ph)%IPareaNormal(3,nIPneighbors,Nmembers))
allocate(geom(ph)%IParea(nIPneighbors,Nmembers))
allocate(geom(ph)%IPcoordinates(3,Nmembers))
call storeGeometry(ph)
if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
@ -489,6 +497,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,Nmembers),source=0.0_pReal)
end associate
if (Nmembers > 0) call stateInit(ini,ph,Nmembers)
@ -543,13 +552,11 @@ end function plastic_nonlocal_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure
!--------------------------------------------------------------------------------------------------
module subroutine nonlocal_dependentState(ph, en, ip, el)
module subroutine nonlocal_dependentState(ph, en)
integer, intent(in) :: &
ph, &
en, &
ip, &
el
en
integer :: &
no, & !< neighbor offset
@ -650,12 +657,12 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
nRealNeighbors = 0.0_pReal
neighbor_rhoTotal = 0.0_pReal
do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el)
neighbor_ip = IPneighborhood(2,n,ip,el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
if (neighbor_el > 0 .and. neighbor_ip > 0) then
if (material_phaseAt(1,neighbor_el) == ph) then
neighbor_el = geom(ph)%IPneighborhood(1,n,en)
neighbor_ip = geom(ph)%IPneighborhood(2,n,en)
if (neighbor_el > 0 .and. neighbor_ip > 0) then
if (material_phaseID(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip) == ph) then
no = material_phaseEntry(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip)
nRealNeighbors = nRealNeighbors + 1.0_pReal
rho_neighbor0 = getRho0(ph,no)
@ -665,11 +672,11 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
neighbor_rhoTotal(1,:,n) = sum(abs(rho_neighbor0(:,edg)),2)
neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor0(:,scr)),2)
connection_latticeConf(1:3,n) = matmul(invFe, discretization_IPcoords(1:3,neighbor_el+neighbor_ip-1) &
- discretization_IPcoords(1:3,el+neighbor_ip-1))
normal_latticeConf = matmul(transpose(invFp), IPareaNormal(1:3,n,ip,el))
connection_latticeConf(1:3,n) = matmul(invFe, geom(ph)%IPcoordinates(1:3,no) &
- geom(ph)%IPcoordinates(1:3,en))
normal_latticeConf = matmul(transpose(invFp), geom(ph)%IPareaNormal(1:3,n,en))
if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image
connection_latticeConf(1:3,n) = normal_latticeConf * IPvolume(ip,el)/IParea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
connection_latticeConf(1:3,n) = normal_latticeConf * geom(ph)%V_0(en)/geom(ph)%IParea(n,en) ! instead take the surface normal scaled with the diameter of the cell
else
! local neighbor or different lattice structure or different constitution instance -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pReal
@ -939,7 +946,7 @@ end subroutine plastic_nonlocal_deltaState
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
module subroutine nonlocal_dotState(Mp,timestep, &
ph,en,ip,el)
ph,en)
real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress
@ -947,9 +954,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
en, &
ip, & !< current integration point
el !< current element number
en
integer :: &
c, & !< character of dislocation
@ -1099,7 +1104,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = rhoDotFlux(timestep, ph,en,ip,el) &
rhoDot = rhoDotFlux(timestep, ph,en) &
+ rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation &
@ -1110,7 +1115,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
.or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at ph ',ph,' en ',en
print'(a)', '<< CONST >> enforcing cutback !!!'
end if
#endif
@ -1129,20 +1134,17 @@ end subroutine nonlocal_dotState
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
#if __INTEL_COMPILER >= 2020
non_recursive function rhoDotFlux(timestep,ph,en,ip,el)
non_recursive function rhoDotFlux(timestep,ph,en)
#else
function rhoDotFlux(timestep,ph,en,ip,el)
function rhoDotFlux(timestep,ph,en)
#endif
real(pReal), intent(in) :: &
timestep !< substepped crystallite time increment
integer, intent(in) :: &
ph, &
en, &
ip, & !< current integration point
el !< current element number
en
integer :: &
neighbor_ph, & !< phase of my neighbor's plasticity
ns, & !< short notation for the total number of active slip systems
c, & !< character of dislocation
n, & !< index of my current neighbor
@ -1215,14 +1217,14 @@ function rhoDotFlux(timestep,ph,en,ip,el)
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux
if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ...
.and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
> geom(ph)%V_0(en)/ maxval(geom(ph)%IParea(:,en)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at ph ',ph,' en ',en
print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
maxval(abs(v0), abs(dot_gamma) > 0.0_pReal &
.and. prm%C_CFL * abs(v0) * timestep &
> IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
> geom(ph)%V_0(en) / maxval(geom(ph)%IParea(:,en))), &
' at a timestep of ',timestep
print*, '<< CONST >> enforcing cutback !!!'
end if
@ -1245,19 +1247,18 @@ function rhoDotFlux(timestep,ph,en,ip,el)
neighbors: do n = 1,nIPneighbors
neighbor_el = IPneighborhood(1,n,ip,el)
neighbor_ip = IPneighborhood(2,n,ip,el)
neighbor_n = IPneighborhood(3,n,ip,el)
np = material_phaseAt(1,neighbor_el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el)
neighbor_el = geom(ph)%IPneighborhood(1,n,en)
neighbor_ip = geom(ph)%IPneighborhood(2,n,en)
neighbor_n = geom(ph)%IPneighborhood(3,n,en)
np = material_phaseID(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip)
no = material_phaseEntry(1,(neighbor_el-1)*discretization_nIPs + neighbor_ip)
opposite_neighbor = n + mod(n,2) - mod(n+1,2)
opposite_el = IPneighborhood(1,opposite_neighbor,ip,el)
opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el)
opposite_n = IPneighborhood(3,opposite_neighbor,ip,el)
opposite_el = geom(ph)%IPneighborhood(1,opposite_neighbor,en)
opposite_ip = geom(ph)%IPneighborhood(2,opposite_neighbor,en)
opposite_n = geom(ph)%IPneighborhood(3,opposite_neighbor,en)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_ph = material_phaseAt(1,neighbor_el)
neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no)
neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no)))
Favg = 0.5_pReal * (my_F + neighbor_F)
@ -1276,12 +1277,12 @@ function rhoDotFlux(timestep,ph,en,ip,el)
!* The entering flux from my neighbor will be distributed on my slip systems according to the
!* compatibility
if (neighbor_n > 0) then
if (phase_plasticity(material_phaseAt(1,neighbor_el)) == PLASTIC_NONLOCAL_ID .and. &
any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) then
if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID .and. &
any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pReal)) then
forall (s = 1:ns, t = 1:4)
neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,neighbor_ph),no)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_ph),no),0.0_pReal)
neighbor_v0(s,t) = plasticState(np)%state0(iV (s,t,np),no)
neighbor_rhoSgl0(s,t) = max(plasticState(np)%state0(iRhoU(s,t,np),no),0.0_pReal)
endforall
where (neighbor_rhoSgl0 * IPvolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%rho_min &
@ -1301,12 +1302,12 @@ function rhoDotFlux(timestep,ph,en,ip,el)
.and. v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density
lineLength = neighbor_rhoSgl0(s,t) * neighbor_v0(s,t) &
* math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (compatibility(c,:,s,n,ip,el) > 0.0_pReal) &
where (dependentState(ph)%compatibility(c,:,s,n,en) > 0.0_pReal) &
rhoDotFlux(:,t) = rhoDotFlux(1:ns,t) &
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2 ! transferring to equally signed mobile dislocation type
where (compatibility(c,:,s,n,ip,el) < 0.0_pReal) &
+ lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to equally signed mobile dislocation type
where (dependentState(ph)%compatibility(c,:,s,n,en) < 0.0_pReal) &
rhoDotFlux(:,topp) = rhoDotFlux(:,topp) &
+ lineLength/IPvolume(ip,el)*compatibility(c,:,s,n,ip,el)**2 ! transferring to opposite signed mobile dislocation type
+ lineLength/geom(ph)%V_0(en)*dependentState(ph)%compatibility(c,:,s,n,en)**2 ! transferring to opposite signed mobile dislocation type
end if
end do
@ -1322,28 +1323,28 @@ function rhoDotFlux(timestep,ph,en,ip,el)
!* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
!* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
if (opposite_n > 0) then
if (phase_plasticity(material_phaseAt(1,opposite_el)) == PLASTIC_NONLOCAL_ID) then
if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID) then
normal_me2neighbor_defConf = math_det33(Favg) &
* matmul(math_inv33(transpose(Favg)),IPareaNormal(1:3,n,ip,el)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor)
* matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en)) ! normal of the interface in (average) deformed configuration (pointing en => neighbor)
normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) &
/ math_det33(my_Fe) ! interface normal in my lattice configuration
area = IParea(n,ip,el) * norm2(normal_me2neighbor)
area = geom(ph)%IParea(n,en) * norm2(normal_me2neighbor)
normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length
do s = 1,ns
do t = 1,4
c = (t + 1) / 2
if (v0(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from en to my neighbor == leaving flux for en (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
if (v0(s,t) * neighbor_v0(s,t) >= 0.0_pReal) then ! no sign change in flux density
transmissivity = sum(compatibility(c,:,s,n,ip,el)**2) ! overall transmissivity from this slip system to my neighbor
transmissivity = sum(dependentState(ph)%compatibility(c,:,s,n,en)**2) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pReal
end if
lineLength = my_rhoSgl0(s,t) * v0(s,t) &
* math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / IPvolume(ip,el) ! subtract dislocation flux from current type
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / geom(ph)%V_0(en) ! subtract dislocation flux from current type
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
+ lineLength / IPvolume(ip,el) * (1.0_pReal - transmissivity) &
+ lineLength / geom(ph)%V_0(en) * (1.0_pReal - transmissivity) &
* sign(1.0_pReal, v0(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
end if
end do
@ -1364,14 +1365,14 @@ end function rhoDotFlux
! plane normals and signed cosine of the angle between the slip directions. Only the largest values
! that sum up to a total of 1 are considered, all others are set to zero.
!--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el)
type(tRotationContainer), dimension(:), intent(in) :: &
orientation ! crystal orientation
integer, intent(in) :: &
ph, &
i, &
e
ip, &
el
integer :: &
n, & ! neighbor index
@ -1397,16 +1398,16 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
associate(prm => param(ph))
ns = prm%sum_N_sl
en = material_phaseMemberAt(1,i,e)
en = material_phaseEntry(1,(el-1)*discretization_nIPs + ip)
!*** start out fully compatible
my_compatibility = 0.0_pReal
forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
neighbors: do n = 1,nIPneighbors
neighbor_e = IPneighborhood(1,n,i,e)
neighbor_i = IPneighborhood(2,n,i,e)
neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
neighbor_phase = material_phaseAt(1,neighbor_e)
neighbor_e = IPneighborhood(1,n,ip,el)
neighbor_i = IPneighborhood(2,n,ip,el)
neighbor_me = material_phaseEntry(1,(neighbor_e-1)*discretization_nIPs + neighbor_i)
neighbor_phase = material_phaseID(1,(neighbor_e-1)*discretization_nIPs + neighbor_i)
if (neighbor_e <= 0 .or. neighbor_i <= 0) then
!* FREE SURFACE
@ -1465,7 +1466,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
end do neighbors
compatibility(:,:,:,:,i,e) = my_compatibility
dependentState(ph)%compatibility(:,:,:,:,material_phaseEntry(1,(el-1)*discretization_nIPs + ip)) = my_compatibility
end associate
@ -1756,14 +1757,29 @@ subroutine storeGeometry(ph)
integer, intent(in) :: ph
integer :: ce, co
integer :: ce, co, nCell
real(pReal), dimension(:), allocatable :: V
integer, dimension(:,:,:), allocatable :: neighborhood
real(pReal), dimension(:,:), allocatable :: area, coords
real(pReal), dimension(:,:,:), allocatable :: areaNormal
nCell = product(shape(IPvolume))
V = reshape(IPvolume,[nCell])
neighborhood = reshape(IPneighborhood,[3,nIPneighbors,nCell])
area = reshape(IParea,[nIPneighbors,nCell])
areaNormal = reshape(IPareaNormal,[3,nIPneighbors,nCell])
coords = reshape(discretization_IPcoords,[3,nCell])
V = reshape(IPvolume,[product(shape(IPvolume))])
do ce = 1, size(material_homogenizationEntry,1)
do co = 1, homogenization_maxNconstituents
if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
if (material_phaseID(co,ce) == ph) then
geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
geom(ph)%IPneighborhood(:,:,material_phaseEntry(co,ce)) = neighborhood(:,:,ce)
geom(ph)%IParea(:,material_phaseEntry(co,ce)) = area(:,ce)
geom(ph)%IPareaNormal(:,:,material_phaseEntry(co,ce)) = areaNormal(:,:,ce)
geom(ph)%IPcoordinates(:,material_phaseEntry(co,ce)) = coords(:,ce)
end if
end do
end do