Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

This commit is contained in:
Martin Diehl 2018-09-20 07:16:03 +02:00
commit fc016bbc6e
4 changed files with 122 additions and 179 deletions

View File

@ -1 +1 @@
v2.0.2-560-g20f0bee4
v2.0.2-565-g59043d58

View File

@ -5,15 +5,10 @@
(output) orientation # quaternion
(output) eulerangles # orientation as Bunge triple in degree
(output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) in crystal reference coordinates
(output) grainrotationx # deviation from initial orientation as angle in degrees around sample reference x axis
(output) grainrotationy # deviation from initial orientation as angle in degrees around sample reference y axis
(output) grainrotationz # deviation from initial orientation as angle in degrees around sample reference z axis
(output) f # deformation gradient tensor; synonyms: "defgrad"
(output) f # deformation gradient tensor
(output) fe # elastic deformation gradient tensor
(output) fp # plastic deformation gradient tensor
(output) e # total strain as Green-Lagrange tensor
(output) ee # elastic strain as Green-Lagrange tensor
(output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola"
(output) s # second Piola-Kichhoff stress tensor; synonyms: "tstar", "secondpiola", "2ndpiola"
(output) p # first Piola-Kichhoff stress tensor
(output) s # second Piola-Kichhoff stress tensor
(output) lp # plastic velocity gradient tensor
(output) elasmatrix # elastic stiffness matrix

View File

@ -18,8 +18,6 @@ mech none
(output) f # deformation gradient tensor; synonyms: "defgrad"
(output) fe # elastic deformation gradient tensor
(output) fp # plastic deformation gradient tensor
(output) e # total strain as Green-Lagrange tensor
(output) ee # elastic strain as Green-Lagrange tensor
(output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola"
(output) lp # plastic velocity gradient tensor

View File

@ -90,9 +90,6 @@ module crystallite
phase_ID, &
texture_ID, &
volume_ID, &
grainrotationx_ID, &
grainrotationy_ID, &
grainrotationz_ID, &
orientation_ID, &
grainrotation_ID, &
eulerangles_ID, &
@ -102,8 +99,6 @@ module crystallite
fi_ID, &
lp_ID, &
li_ID, &
e_ID, &
ee_ID, &
p_ID, &
s_ID, &
elasmatrix_ID, &
@ -122,13 +117,13 @@ module crystallite
crystallite_postResults
private :: &
integrateState, &
crystallite_integrateStateFPI, &
crystallite_integrateStateEuler, &
crystallite_integrateStateAdaptiveEuler, &
crystallite_integrateStateRK4, &
crystallite_integrateStateRKCK45, &
crystallite_integrateStress, &
crystallite_stateJump
integrateStateFPI, &
integrateStateEuler, &
integrateStateAdaptiveEuler, &
integrateStateRK4, &
integrateStateRKCK45, &
integrateStress, &
stateJump
contains
@ -142,12 +137,14 @@ subroutine crystallite_init
compiler_version, &
compiler_options
#endif
#ifdef DEBUG
use debug, only: &
debug_info, &
debug_reset, &
debug_level, &
debug_crystallite, &
debug_levelBasic
#endif
use numerics, only: &
numerics_integrator, &
worldrank, &
@ -175,8 +172,7 @@ subroutine crystallite_init
use config, only: &
config_deallocate, &
config_crystallite, &
crystallite_name, &
material_Nphase
crystallite_name
use constitutive, only: &
constitutive_initialFi, &
constitutive_microstructure ! derived (shortcut) quantities of given state
@ -190,7 +186,6 @@ subroutine crystallite_init
e, & !< counter in element loop
o = 0_pInt, & !< counter in output loop
r, &
ph, & !< counter in crystallite loop
cMax, & !< maximum number of integration point components
iMax, & !< maximum number of integration points
eMax, & !< maximum number of elements
@ -199,8 +194,6 @@ subroutine crystallite_init
mySize
character(len=65536), dimension(:), allocatable :: str
character(len=65536) :: &
tag = ''
write(6,'(/,a)') ' <<<+- crystallite init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -272,15 +265,15 @@ subroutine crystallite_init
select case(numerics_integrator(1))
case(1_pInt)
integrateState => crystallite_integrateStateFPI
integrateState => integrateStateFPI
case(2_pInt)
integrateState => crystallite_integrateStateEuler
integrateState => integrateStateEuler
case(3_pInt)
integrateState => crystallite_integrateStateAdaptiveEuler
integrateState => integrateStateAdaptiveEuler
case(4_pInt)
integrateState => crystallite_integrateStateRK4
integrateState => integrateStateRK4
case(5_pInt)
integrateState => crystallite_integrateStateRKCK45
integrateState => integrateStateRKCK45
end select
@ -302,12 +295,6 @@ subroutine crystallite_init
crystallite_outputID(o,c) = texture_ID
case ('volume') outputName
crystallite_outputID(o,c) = volume_ID
case ('grainrotationx') outputName
crystallite_outputID(o,c) = grainrotationx_ID
case ('grainrotationy') outputName
crystallite_outputID(o,c) = grainrotationy_ID
case ('grainrotationz') outputName
crystallite_outputID(o,c) = grainrotationx_ID
case ('orientation') outputName
crystallite_outputID(o,c) = orientation_ID
case ('grainrotation') outputName
@ -326,10 +313,6 @@ subroutine crystallite_init
crystallite_outputID(o,c) = lp_ID
case ('li') outputName
crystallite_outputID(o,c) = li_ID
case ('e') outputName
crystallite_outputID(o,c) = e_ID
case ('ee') outputName
crystallite_outputID(o,c) = ee_ID
case ('p','firstpiola','1stpiola') outputName
crystallite_outputID(o,c) = p_ID
case ('s','tstar','secondpiola','2ndpiola') outputName
@ -341,7 +324,7 @@ subroutine crystallite_init
case ('neighboringelement') outputName
crystallite_outputID(o,c) = neighboringelement_ID
case default outputName
call IO_error(105_pInt,ext_msg=tag//' (Crystallite)')
call IO_error(105_pInt,ext_msg=trim(str(o))//' (Crystallite)')
end select outputName
enddo
enddo
@ -350,13 +333,13 @@ subroutine crystallite_init
do r = 1_pInt,size(config_crystallite)
do o = 1_pInt,crystallite_Noutput(r)
select case(crystallite_outputID(o,r))
case(phase_ID,texture_ID,volume_ID,grainrotationx_ID,grainrotationy_ID,grainrotationz_ID)
case(phase_ID,texture_ID,volume_ID)
mySize = 1_pInt
case(orientation_ID,grainrotation_ID)
mySize = 4_pInt
case(eulerangles_ID)
mySize = 3_pInt
case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,e_ID,ee_ID,p_ID,s_ID)
case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,p_ID,s_ID)
mySize = 9_pInt
case(elasmatrix_ID)
mySize = 36_pInt
@ -437,59 +420,19 @@ subroutine crystallite_init
call crystallite_stressAndItsTangent(.true.) ! request elastic answers
!--------------------------------------------------------------------------------------------------
! debug output
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp: ', shape(crystallite_Fp)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fi: ', shape(crystallite_Fi)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp: ', shape(crystallite_Lp)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Li: ', shape(crystallite_Li)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_F0: ', shape(crystallite_F0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp0: ', shape(crystallite_Fp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fi0: ', shape(crystallite_Fi0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp0: ', shape(crystallite_Lp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Li0: ', shape(crystallite_Li0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF: ', shape(crystallite_partionedF)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFi0: ', shape(crystallite_partionedFi0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLi0: ', shape(crystallite_partionedLi0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF: ', shape(crystallite_subF)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF0: ', shape(crystallite_subF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFp0: ', shape(crystallite_subFp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFi0: ', shape(crystallite_subFi0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLp0: ', shape(crystallite_subLp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLi0: ', shape(crystallite_subLi0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_P: ', shape(crystallite_P)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedTstar0_v: ', shape(crystallite_partionedTstar0_v)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_dPdF: ', shape(crystallite_dPdF)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_dPdF0: ', shape(crystallite_dPdF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partioneddPdF0: ', shape(crystallite_partioneddPdF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_orientation: ', shape(crystallite_orientation)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_orientation0: ', shape(crystallite_orientation0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_rotation: ', shape(crystallite_rotation)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_disorientation: ', shape(crystallite_disorientation)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_dt: ', shape(crystallite_dt)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subdt: ', shape(crystallite_subdt)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFrac: ', shape(crystallite_subFrac)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subStep: ', shape(crystallite_subStep)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_localPlasticity: ', shape(crystallite_localPlasticity)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_requested: ', shape(crystallite_requested)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_todo: ', shape(crystallite_todo)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_converged: ', shape(crystallite_converged)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_sizePostResults: ', shape(crystallite_sizePostResults)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_sizePostResult: ', shape(crystallite_sizePostResult)
write(6,'(/,a35,1x,i10)') 'Number of nonlocal grains: ',count(.not. crystallite_localPlasticity)
write(6,'(a42,1x,i10)') ' # of elements: ', eMax
write(6,'(a42,1x,i10)') 'max # of integration points/element: ', iMax
write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax
write(6,'(a42,1x,i10)') 'max # of neigbours/integration point: ', nMax
write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity)
flush(6)
endif
call debug_info
call debug_reset
#endif
end subroutine crystallite_init
@ -506,6 +449,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
subStepSizeCryst, &
stepIncreaseCryst, &
numerics_timeSyncing
#ifdef DEBUG
use debug, only: &
debug_level, &
debug_crystallite, &
@ -515,6 +459,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
debug_e, &
debug_i, &
debug_g
#endif
use IO, only: &
IO_warning, &
IO_error
@ -591,7 +536,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
real(pReal), dimension(9,9):: temp_99
logical :: error
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt &
.and. FEsolving_execElem(1) <= debug_e &
.and. debug_e <= FEsolving_execElem(2)) then
@ -610,6 +555,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Li0', &
transpose(crystallite_partionedLi0(1:3,1:3,debug_g,debug_i,debug_e))
endif
#endif
!--------------------------------------------------------------------------------------------------
! initialize to starting condition
@ -653,10 +599,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
NiterationCrystallite = 0_pInt
cutbackLooping: do while (any(crystallite_todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2))))
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i6)') '<< CRYST >> crystallite iteration ',NiterationCrystallite
#endif
timeSyncing1: if (any(.not. crystallite_localPlasticity) .and. numerics_timeSyncing) then
! Time synchronization can only be used for nonlocal calculations, and only there it makes sense.
@ -675,6 +621,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! and its not clear how to fix this, so all nonlocals become terminally ill.
if (any(crystallite_syncSubFrac .and. .not. crystallite_converged(1,:,:))) then
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
@ -683,6 +630,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
enddo
endif
#endif
crystallite_syncSubFrac = .false.
where(.not. crystallite_localPlasticity)
crystallite_substep = 0.0_pReal
@ -697,8 +645,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
enddo
!$OMP END PARALLEL DO
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i6)') '<< CRYST >> time synchronization: wind forward'
#endif
endif
elseif (any(crystallite_syncSubFracCompleted)) then
@ -714,8 +664,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_clearToCutback(i,e) = crystallite_localPlasticity(1,i,e) .or. .not. crystallite_converged(1,i,e)
enddo
enddo
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i6)') '<< CRYST >> time synchronization: done, proceed with cutback'
#endif
else
! Normal calculation.
@ -741,8 +693,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if (all(crystallite_localPlasticity .or. crystallite_converged)) then
if (all(crystallite_localPlasticity .or. crystallite_subStep + crystallite_subFrac >= 1.0_pReal)) then
crystallite_clearToWindForward = .true. ! final wind forward
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i6)') '<< CRYST >> final wind forward'
#endif
else
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -751,8 +705,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
enddo
!$OMP END PARALLEL DO
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i6)') '<< CRYST >> wind forward'
#endif
endif
else
subFracIntermediate = maxval(crystallite_subFrac, mask=.not.crystallite_localPlasticity)
@ -838,8 +794,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
enddo
!$OMP END PARALLEL DO
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i6)') '<< CRYST >> time synchronization: cutback'
#endif
else
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -848,8 +806,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
enddo
!$OMP END PARALLEL DO
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i6)') '<< CRYST >> cutback'
#endif
endif
endif
endif
@ -979,6 +939,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
timeSyncing2: if(numerics_timeSyncing) then
if (any(.not. crystallite_localPlasticity .and. .not. crystallite_todo .and. .not. crystallite_converged &
.and. crystallite_subStep <= subStepMinCryst)) then ! no way of rescuing a nonlocal ip that violated the lower time step limit, ...
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then
elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
@ -990,6 +951,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
enddo elementLooping4
endif
#endif
where(.not. crystallite_localPlasticity)
crystallite_todo = .false. ! ... so let all nonlocal ips die peacefully
crystallite_subStep = 0.0_pReal
@ -997,6 +959,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif
endif timeSyncing2
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then
write(6,'(/,a,f8.5)') '<< CRYST >> min(subStep) ',minval(crystallite_subStep)
write(6,'(a,f8.5)') '<< CRYST >> max(subStep) ',maxval(crystallite_subStep)
@ -1009,9 +972,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
flush(6)
endif
endif
#endif
! --- integrate --- requires fully defined state array (basic + dependent state)
if (any(crystallite_todo)) call integrateState()
where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) & ! do not try non-converged & fully cutbacked any further
crystallite_todo = .true.
@ -1027,9 +990,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do c = 1,homogenization_Ngrains(mesh_element(3,e))
if (.not. crystallite_converged(c,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway)
#ifdef DEBUG
if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> no convergence: respond fully elastic at el (elFE) ip ipc ', &
e,'(',mesh_element(1,e),')',i,c
#endif
invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,c,i,e))
Fe_guess = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,c,i,e), invFp), &
math_inv33(crystallite_partionedFi0(1:3,1:3,c,i,e)))
@ -1037,6 +1002,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_P(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,c,i,e), invFp), &
math_mul33x33(Tstar,transpose(invFp)))
endif
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. c == debug_g) &
.or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) then
@ -1053,6 +1019,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
transpose(crystallite_Li(1:3,1:3,c,i,e))
flush(6)
endif
#endif
enddo
enddo
enddo elementLooping5
@ -1177,20 +1144,20 @@ end subroutine crystallite_stressAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 4th order explicit Runge Kutta method
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRK4()
subroutine integrateStateRK4()
use, intrinsic :: &
IEEE_arithmetic
use debug, only: &
#ifdef DEBUG
use debug, only: &
debug_e, &
debug_i, &
debug_g, &
#endif
debug_level, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective
#endif
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
@ -1358,7 +1325,7 @@ subroutine crystallite_integrateStateRK4()
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) = crystallite_stateJump(g,i,e)
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)
@ -1389,7 +1356,7 @@ subroutine crystallite_integrateStateRK4()
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) = crystallite_integrateStress(g,i,e,timeStepFraction(n)) ! fraction of original times step
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)
@ -1460,27 +1427,27 @@ subroutine crystallite_integrateStateRK4()
endif
endif
end subroutine crystallite_integrateStateRK4
end subroutine integrateStateRK4
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with
!> adaptive step size (use 5th order solution to advance = "local extrapolation")
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRKCK45()
subroutine integrateStateRKCK45()
use, intrinsic :: &
IEEE_arithmetic
use debug, only: &
#ifdef DEBUG
use debug, only: &
debug_e, &
debug_i, &
debug_g, &
#endif
debug_level, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective
#endif
use numerics, only: &
rTol_crystalliteState
use FEsolving, only: &
@ -1556,8 +1523,10 @@ subroutine crystallite_integrateStateRKCK45()
singleRun ! flag indicating computation for single (g,i,e) triple
eIter = FEsolving_execElem(1:2)
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,1x,i1)') '<< CRYST >> Runge--Kutta step',1
#endif
! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP ---
do e = eIter(1),eIter(2)
@ -1677,7 +1646,7 @@ subroutine crystallite_integrateStateRKCK45()
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) = crystallite_stateJump(g,i,e)
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)
@ -1708,7 +1677,7 @@ subroutine crystallite_integrateStateRKCK45()
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) = crystallite_integrateStress(g,i,e,C(stage)) ! fraction of original time step
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)
@ -1894,7 +1863,7 @@ subroutine crystallite_integrateStateRKCK45()
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) = crystallite_stateJump(g,i,e)
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)
@ -1925,7 +1894,7 @@ subroutine crystallite_integrateStateRKCK45()
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) = crystallite_integrateStress(g,i,e)
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)
@ -1949,32 +1918,33 @@ subroutine crystallite_integrateStateRKCK45()
! --- nonlocal convergence check ---
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged' ! if not requesting Integration of just a single IP
#endif
if ((.not. singleRun) .and. any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
end subroutine crystallite_integrateStateRKCK45
end subroutine integrateStateRKCK45
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateAdaptiveEuler()
subroutine integrateStateAdaptiveEuler()
use, intrinsic :: &
IEEE_arithmetic
use debug, only: &
#ifdef DEBUG
use debug, only: &
debug_e, &
debug_i, &
debug_g, &
#endif
debug_level, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective
#endif
use numerics, only: &
rTol_crystalliteState
use FEsolving, only: &
@ -2120,7 +2090,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
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) = crystallite_stateJump(g,i,e)
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)
@ -2152,7 +2122,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
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) = crystallite_integrateStress(g,i,e)
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)
@ -2273,33 +2243,34 @@ subroutine crystallite_integrateStateAdaptiveEuler()
! --- NONLOCAL CONVERGENCE CHECK ---
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged'
#endif
if ((.not. singleRun) .and. any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
end subroutine crystallite_integrateStateAdaptiveEuler
end subroutine integrateStateAdaptiveEuler
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, and state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateEuler()
subroutine integrateStateEuler()
use, intrinsic :: &
IEEE_arithmetic
use debug, only: &
#ifdef DEBUG
use debug, only: &
debug_e, &
debug_i, &
debug_g, &
#endif
debug_level, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective
#endif
use numerics, only: &
numerics_timeSyncing
use FEsolving, only: &
@ -2427,7 +2398,7 @@ eIter = FEsolving_execElem(1:2)
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) = crystallite_stateJump(g,i,e)
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) & ! if broken non-local...
.and. .not. numerics_timeSyncing) then
@ -2461,7 +2432,7 @@ eIter = FEsolving_execElem(1:2)
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) = crystallite_integrateStress(g,i,e)
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) & ! if broken non-local...
.and. .not. numerics_timeSyncing) then
@ -2493,27 +2464,27 @@ eIter = FEsolving_execElem(1:2)
crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
endif
end subroutine crystallite_integrateStateEuler
end subroutine integrateStateEuler
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateFPI()
subroutine integrateStateFPI()
use, intrinsic :: &
IEEE_arithmetic
use debug, only: &
#ifdef DEBUG
use debug, only: &
debug_e, &
debug_i, &
debug_g, &
#endif
debug_level,&
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective
#endif
use numerics, only: &
nState, &
rTol_crystalliteState
@ -2577,8 +2548,10 @@ subroutine crystallite_integrateStateFPI()
singleRun = (eIter(1) == eIter(2) .and. iIter(1,eIter(1)) == iIter(2,eIter(2)))
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo at start of state integration'
#endif
!--------------------------------------------------------------------------------------------------
! initialize dotState
@ -2633,8 +2606,10 @@ subroutine crystallite_integrateStateFPI()
NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
enddo
if (NaN) then ! NaN occured in any dotState
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,*) '<< CRYST >> dotstate ',plasticState(p)%dotState(:,c)
#endif
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
@ -2648,9 +2623,10 @@ subroutine crystallite_integrateStateFPI()
!$OMP ENDDO
! --- UPDATE STATE ---
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after preguess of state'
#endif
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c)
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
@ -2706,14 +2682,15 @@ subroutine crystallite_integrateStateFPI()
! --- 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) = crystallite_integrateStress(g,i,e)
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)
@ -2724,13 +2701,10 @@ subroutine crystallite_integrateStateFPI()
enddo; enddo; enddo
!$OMP ENDDO
!$OMP SINGLE
!$OMP CRITICAL (write2out)
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration'
!$OMP END CRITICAL (write2out)
!$OMP END SINGLE
#endif
! --- DOT STATE ---
@ -2904,7 +2878,7 @@ subroutine crystallite_integrateStateFPI()
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. crystallite_converged(g,i,e)) then ! converged and still alive...
crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
crystallite_todo(g,i,e) = stateJump(g,i,e)
!$OMP FLUSH(crystallite_todo)
if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken
crystallite_converged(g,i,e) = .false.
@ -2919,10 +2893,11 @@ subroutine crystallite_integrateStateFPI()
!$OMP ENDDO
!$OMP END PARALLEL
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
' grains converged after state integration #', NiterationState
#endif
! --- NON-LOCAL CONVERGENCE CHECK ---
@ -2931,12 +2906,15 @@ subroutine crystallite_integrateStateFPI()
crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
endif
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
' grains converged after non-local check'
write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_todo(:,:,:)), &
' grains todo after state integration #', NiterationState
endif
#endif
! --- CHECK IF DONE WITH INTEGRATION ---
doneWithIntegration = .true.
@ -2950,14 +2928,14 @@ subroutine crystallite_integrateStateFPI()
enddo elemLoop
enddo crystalliteLooping
end subroutine crystallite_integrateStateFPI
end subroutine integrateStateFPI
!--------------------------------------------------------------------------------------------------
!> @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
!--------------------------------------------------------------------------------------------------
logical function crystallite_stateJump(ipc,ip,el)
logical function stateJump(ipc,ip,el)
use, intrinsic :: &
IEEE_arithmetic
use prec, only: &
@ -3007,7 +2985,7 @@ logical function crystallite_stateJump(ipc,ip,el)
mySizePlasticDeltaState = plasticState(p)%sizeDeltaState
if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)))) then ! NaN occured in deltaState
crystallite_stateJump = .false.
stateJump = .false.
return
endif
@ -3021,7 +2999,7 @@ logical function crystallite_stateJump(ipc,ip,el)
myOffsetSourceDeltaState = sourceState(p)%p(mySource)%offsetDeltaState
mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState
if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)))) then ! NaN occured in deltaState
crystallite_stateJump = .false.
stateJump = .false.
return
endif
sourceState(p)%p(mySource)%state(myOffsetSourceDeltaState + 1_pInt : &
@ -3044,9 +3022,9 @@ logical function crystallite_stateJump(ipc,ip,el)
endif
#endif
crystallite_stateJump = .true.
stateJump = .true.
end function crystallite_stateJump
end function stateJump
!--------------------------------------------------------------------------------------------------
@ -3080,7 +3058,7 @@ end function crystallite_push33ToRef
!> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction
!--------------------------------------------------------------------------------------------------
logical function crystallite_integrateStress(&
logical function integrateStress(&
ipc,& ! grain number
ip,& ! integration point number
el,& ! element number
@ -3097,16 +3075,16 @@ logical function crystallite_integrateStress(&
iJacoLpresiduum, &
subStepSizeLp, &
subStepSizeLi
use debug, only: debug_level, &
#ifdef DEBUG
use debug, only: debug_level, &
debug_e, &
debug_i, &
debug_g, &
#endif
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
debug_levelSelective
#endif
use constitutive, only: constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangents, &
@ -3198,7 +3176,7 @@ logical function crystallite_integrateStress(&
dgesv
!* be pessimistic
crystallite_integrateStress = .false.
integrateStress = .false.
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
@ -3537,7 +3515,7 @@ logical function crystallite_integrateStress(&
!* set return flag to true
crystallite_integrateStress = .true.
integrateStress = .true.
#ifdef DEBUG
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
@ -3552,7 +3530,7 @@ logical function crystallite_integrateStress(&
endif
#endif
end function crystallite_integrateStress
end function integrateStress
!--------------------------------------------------------------------------------------------------
@ -3668,9 +3646,7 @@ function crystallite_postResults(ipc, ip, el)
math_det33, &
math_I3, &
inDeg, &
math_Mandel6to33, &
math_qMul, &
math_qConj
math_Mandel6to33
use mesh, only: &
mesh_element, &
mesh_ipVolume, &
@ -3701,10 +3677,6 @@ function crystallite_postResults(ipc, ip, el)
1+plasticState(material_phase(ipc,ip,el))%sizePostResults + &
sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: &
crystallite_postResults
real(pReal), dimension(3,3) :: &
Ee
real(pReal), dimension(4) :: &
rotation
real(pReal) :: &
detF
integer(pInt) :: &
@ -3748,18 +3720,6 @@ function crystallite_postResults(ipc, ip, el)
crystallite_postResults(c+1:c+mySize) = &
math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates
crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree
case (grainrotationx_ID)
mySize = 1_pInt
rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates
crystallite_postResults(c+1) = inDeg * rotation(1) * rotation(4) ! angle in degree
case (grainrotationy_ID)
mySize = 1_pInt
rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates
crystallite_postResults(c+1) = inDeg * rotation(2) * rotation(4) ! angle in degree
case (grainrotationz_ID)
mySize = 1_pInt
rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates
crystallite_postResults(c+1) = inDeg * rotation(3) * rotation(4) ! angle in degree
! remark: tensor output is of the form 11,12,13, 21,22,23, 31,32,33
! thus row index i is slow, while column index j is fast. reminder: "row is slow"
@ -3768,20 +3728,10 @@ function crystallite_postResults(ipc, ip, el)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_partionedF(1:3,1:3,ipc,ip,el)),[mySize])
case (e_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( &
transpose(crystallite_partionedF(1:3,1:3,ipc,ip,el)), &
crystallite_partionedF(1:3,1:3,ipc,ip,el)) - math_I3),[mySize])
case (fe_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &
reshape(transpose(crystallite_Fe(1:3,1:3,ipc,ip,el)),[mySize])
case (ee_ID)
Ee = 0.5_pReal *(math_mul33x33(transpose(crystallite_Fe(1:3,1:3,ipc,ip,el)), &
crystallite_Fe(1:3,1:3,ipc,ip,el)) - math_I3)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(Ee,[mySize])
case (fp_ID)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = &