in preguess of FPintegrator: resetting of previous dotstates can be done in same loop as collectdotstate

polishing: indentation level and capitalization
This commit is contained in:
Christoph Kords 2013-11-21 10:58:41 +00:00
parent 0337b4f319
commit 883669bd77
1 changed files with 638 additions and 536 deletions

View File

@ -411,7 +411,7 @@ subroutine crystallite_init(temperature)
!--------------------------------------------------------------------------------------------------
! debug output
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Temperature: ', shape(crystallite_temperature)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_temperature: ', shape(crystallite_temperature)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_heat: ', shape(crystallite_heat)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp: ', shape(crystallite_Fp)
@ -1047,6 +1047,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+--
elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
@ -1084,6 +1085,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
computeJacobian: if(updateJaco) then
jacobianMethod: if (analyticJaco) then
! --- ANALYTIC JACOBIAN ---
!$OMP PARALLEL DO PRIVATE(dFedF,dSdF,dSdFe,myNgrains)
elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
@ -1103,9 +1107,16 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
enddo; enddo
enddo elementLooping6
!$OMP END PARALLEL DO
else jacobianMethod
! --- STANDARD (PERTURBATION METHOD) FOR JACOBIAN ---
numerics_integrationMode = 2_pInt
! --- BACKUP ---
elementLooping7: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
@ -1124,8 +1135,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
P_backup = crystallite_P
convergenceFlag_backup = crystallite_converged
! --- CALCULATE STATE AND STRESS FOR PERTURBATION ---
dPdF_perturbation1 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
dPdF_perturbation2 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
do perturbation = 1,2 ! forward and backward perturbation
@ -1139,6 +1150,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
endif
! --- INITIALIZE UNPERTURBED STATE ---
select case(numerics_integrator(numerics_integrationMode))
case(1_pInt) ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state
do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -1172,7 +1184,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
crystallite_Tstar_v = crystallite_subTstar0_v
end select
! --- PERTURB EITHER FORWARD OR BACKWARD ---
crystallite_subF = F_backup
@ -1214,8 +1225,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
endif
enddo ! perturbation direction
! --- STIFFNESS ACCORDING TO PERTURBATION METHOD AND CONVERGENCE ---
elementLooping9: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
select case(pert_method)
@ -1239,6 +1250,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
enddo elementLooping9
! --- RESTORE ---
elementLooping10: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
@ -1258,6 +1270,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
crystallite_converged = convergenceFlag_backup
endif jacobianMethod
rateSensitivity: if (rate_sensitivity) then
!$OMP PARALLEL DO PRIVATE(dFedFdot,dSdFdot,dSdFe,Fpinv_rate,FDot_inv,counter,dFp_invdFdot,myNgrains)
elementLooping11: do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -1380,8 +1393,9 @@ 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
constitutive_RK4dotState(g,i,e)%p = 0.0_pReal ! initialize Runge-Kutta dotState
if (crystallite_todo(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_temperature(i,e), &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -1416,10 +1430,11 @@ subroutine crystallite_integrateStateRK4()
if (crystallite_todo(g,i,e)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
if (n < 4) then
constitutive_RK4dotState(g,i,e)%p = constitutive_RK4dotState(g,i,e)%p + weight(n)*constitutive_dotState(g,i,e)%p
constitutive_RK4dotState(g,i,e)%p = constitutive_RK4dotState(g,i,e)%p &
+ weight(n)*constitutive_dotState(g,i,e)%p
elseif (n == 4) then
constitutive_dotState(g,i,e)%p = (constitutive_RK4dotState(g,i,e)%p + &
weight(n)*constitutive_dotState(g,i,e)%p) / 6.0_pReal ! use weighted RKdotState for final integration
constitutive_dotState(g,i,e)%p = (constitutive_RK4dotState(g,i,e)%p &
+ weight(n)*constitutive_dotState(g,i,e)%p) / 6.0_pReal ! use weighted RKdotState for final integration
endif
endif
enddo; enddo; enddo
@ -1429,7 +1444,8 @@ subroutine crystallite_integrateStateRK4()
if (crystallite_todo(g,i,e)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) * timeStepFraction(n)
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
* crystallite_subdt(g,i,e) * timeStepFraction(n)
if (n == 4) then ! final integration step
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
@ -1470,7 +1486,7 @@ subroutine crystallite_integrateStateRK4()
!$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
if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
call constitutive_microstructure(crystallite_temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
endif
enddo; enddo; enddo
@ -1501,8 +1517,9 @@ subroutine crystallite_integrateStateRK4()
!$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
if (crystallite_todo(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(i,e), timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_temperature(i,e), &
timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
crystallite_subFrac, g,i,e)
endif
enddo; enddo; enddo
@ -1560,7 +1577,8 @@ end subroutine crystallite_integrateStateRK4
!> adaptive step size (use 5th order solution to advance = "local extrapolation")
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRKCK45()
use debug, only: debug_level, &
use debug, only: &
debug_level, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
@ -1569,16 +1587,21 @@ subroutine crystallite_integrateStateRKCK45()
debug_i, &
debug_g, &
debug_StateLoopDistribution
use numerics, only: rTol_crystalliteState, &
use numerics, only: &
rTol_crystalliteState, &
numerics_integrationMode
use FEsolving, only: FEsolving_execElem, &
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
use mesh, only: mesh_element, &
use mesh, only: &
mesh_element, &
mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_Ngrains, &
use material, only: &
homogenization_Ngrains, &
homogenization_maxNgrains
use constitutive, only: constitutive_sizeDotState, &
use constitutive, only: &
constitutive_sizeDotState, &
constitutive_maxSizeDotState, &
constitutive_state, &
constitutive_aTolState, &
@ -1591,17 +1614,9 @@ subroutine crystallite_integrateStateRKCK45()
constitutive_microstructure
implicit none
integer(pInt) :: &
e, & ! element index in element loop
i, & ! integration point index in ip loop
g, & ! grain index in grain loop
n, & ! stage index in integration stage loop
mySizeDotState, & ! size of dot State
s ! state index
integer(pInt), dimension(2) :: eIter ! bounds for element iteration
integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
gIter ! bounds for grain iteration
real(pReal), dimension(5,5), parameter :: A = reshape([&
real(pReal), dimension(5,5), parameter :: &
A = reshape([&
.2_pReal, .075_pReal, .3_pReal, -11.0_pReal/54.0_pReal, 1631.0_pReal/55296.0_pReal, &
.0_pReal, .225_pReal, -.9_pReal, 2.5_pReal, 175.0_pReal/512.0_pReal, &
.0_pReal, .0_pReal, 1.2_pReal, -70.0_pReal/27.0_pReal, 575.0_pReal/13824.0_pReal, &
@ -1609,20 +1624,34 @@ subroutine crystallite_integrateStateRKCK45()
.0_pReal, .0_pReal, .0_pReal, .0_pReal, 253.0_pReal/4096.0_pReal], &
[5,5], order=[2,1]) !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
real(pReal), dimension(6), parameter :: B = &
real(pReal), dimension(6), parameter :: &
B = &
[37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, &
125.0_pReal/594.0_pReal, .0_pReal, 512.0_pReal/1771.0_pReal], & !< coefficients in Butcher tableau (used for final integration and error estimate)
DB = B - &
[2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,&
13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 0.25_pReal] !< coefficients in Butcher tableau (used for final integration and error estimate)
real(pReal), dimension(5), parameter :: &
C = [0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
real(pReal), dimension(5), parameter :: C = &
[0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
integer(pInt) :: &
e, & ! element index in element loop
i, & ! integration point index in ip loop
g, & ! grain index in grain loop
n, & ! stage index in integration stage loop
mySizeDotState, & ! size of dot State
s ! state index
integer(pInt), dimension(2) :: &
eIter ! bounds for element iteration
integer(pInt), dimension(2,mesh_NcpElems) :: &
iIter, & ! bounds for ip iteration
gIter ! bounds for grain iteration
real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
stateResiduum, & ! residuum from evolution in micrstructure
relStateResiduum ! relative residuum from evolution in microstructure
logical :: singleRun ! flag indicating computation for single (g,i,e) triple
logical :: &
singleRun ! flag indicating computation for single (g,i,e) triple
! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP ---
eIter = FEsolving_execElem(1:2)
@ -1646,8 +1675,9 @@ subroutine crystallite_integrateStateRKCK45()
!$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
if (crystallite_todo(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_temperature(i,e), &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -1689,25 +1719,25 @@ 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
if (crystallite_todo(g,i,e)) then
if (n == 1) then ! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION (CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE)
constitutive_dotState(g,i,e)%p = a(1,1) * constitutive_RKCK45dotState(1,g,i,e)%p
constitutive_dotState(g,i,e)%p = A(1,1) * constitutive_RKCK45dotState(1,g,i,e)%p
elseif (n == 2) then
constitutive_dotState(g,i,e)%p = a(1,2) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ a(2,2) * constitutive_RKCK45dotState(2,g,i,e)%p
constitutive_dotState(g,i,e)%p = A(1,2) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ A(2,2) * constitutive_RKCK45dotState(2,g,i,e)%p
elseif (n == 3) then
constitutive_dotState(g,i,e)%p = a(1,3) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ a(2,3) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ a(3,3) * constitutive_RKCK45dotState(3,g,i,e)%p
constitutive_dotState(g,i,e)%p = A(1,3) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ A(2,3) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ A(3,3) * constitutive_RKCK45dotState(3,g,i,e)%p
elseif (n == 4) then
constitutive_dotState(g,i,e)%p = a(1,4) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ a(2,4) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ a(3,4) * constitutive_RKCK45dotState(3,g,i,e)%p &
+ a(4,4) * constitutive_RKCK45dotState(4,g,i,e)%p
constitutive_dotState(g,i,e)%p = A(1,4) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ A(2,4) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ A(3,4) * constitutive_RKCK45dotState(3,g,i,e)%p &
+ A(4,4) * constitutive_RKCK45dotState(4,g,i,e)%p
elseif (n == 5) then
constitutive_dotState(g,i,e)%p = a(1,5) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ a(2,5) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ a(3,5) * constitutive_RKCK45dotState(3,g,i,e)%p &
+ a(4,5) * constitutive_RKCK45dotState(4,g,i,e)%p &
+ a(5,5) * constitutive_RKCK45dotState(5,g,i,e)%p
constitutive_dotState(g,i,e)%p = A(1,5) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ A(2,5) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ A(3,5) * constitutive_RKCK45dotState(3,g,i,e)%p &
+ A(4,5) * constitutive_RKCK45dotState(4,g,i,e)%p &
+ A(5,5) * constitutive_RKCK45dotState(5,g,i,e)%p
endif
endif
enddo; enddo; enddo
@ -1717,7 +1747,8 @@ subroutine crystallite_integrateStateRKCK45()
if (crystallite_todo(g,i,e)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
* crystallite_subdt(g,i,e)
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -1746,7 +1777,7 @@ subroutine crystallite_integrateStateRKCK45()
!$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
if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
call constitutive_microstructure(crystallite_temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
endif
enddo; enddo; enddo
@ -1780,8 +1811,9 @@ subroutine crystallite_integrateStateRKCK45()
!$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
if (crystallite_todo(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(i,e), c(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_temperature(i,e), &
C(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
crystallite_subFrac, g,i,e)
endif
enddo; enddo; enddo
@ -1831,22 +1863,22 @@ subroutine crystallite_integrateStateRKCK45()
! CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE
stateResiduum(1:mySizeDotState,g,i,e) = &
( db(1) * constitutive_RKCK45dotState(1,g,i,e)%p(1:mySizeDotState) &
+ db(2) * constitutive_RKCK45dotState(2,g,i,e)%p(1:mySizeDotState) &
+ db(3) * constitutive_RKCK45dotState(3,g,i,e)%p(1:mySizeDotState) &
+ db(4) * constitutive_RKCK45dotState(4,g,i,e)%p(1:mySizeDotState) &
+ db(5) * constitutive_RKCK45dotState(5,g,i,e)%p(1:mySizeDotState) &
+ db(6) * constitutive_RKCK45dotState(6,g,i,e)%p(1:mySizeDotState)) &
( DB(1) * constitutive_RKCK45dotState(1,g,i,e)%p(1:mySizeDotState) &
+ DB(2) * constitutive_RKCK45dotState(2,g,i,e)%p(1:mySizeDotState) &
+ DB(3) * constitutive_RKCK45dotState(3,g,i,e)%p(1:mySizeDotState) &
+ DB(4) * constitutive_RKCK45dotState(4,g,i,e)%p(1:mySizeDotState) &
+ DB(5) * constitutive_RKCK45dotState(5,g,i,e)%p(1:mySizeDotState) &
+ DB(6) * constitutive_RKCK45dotState(6,g,i,e)%p(1:mySizeDotState)) &
* crystallite_subdt(g,i,e)
! --- dot state ---
constitutive_dotState(g,i,e)%p = b(1) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ b(2) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ b(3) * constitutive_RKCK45dotState(3,g,i,e)%p &
+ b(4) * constitutive_RKCK45dotState(4,g,i,e)%p &
+ b(5) * constitutive_RKCK45dotState(5,g,i,e)%p &
+ b(6) * constitutive_RKCK45dotState(6,g,i,e)%p
constitutive_dotState(g,i,e)%p = B(1) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ B(2) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ B(3) * constitutive_RKCK45dotState(3,g,i,e)%p &
+ B(4) * constitutive_RKCK45dotState(4,g,i,e)%p &
+ B(5) * constitutive_RKCK45dotState(5,g,i,e)%p &
+ B(6) * constitutive_RKCK45dotState(6,g,i,e)%p
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -1858,7 +1890,8 @@ subroutine crystallite_integrateStateRKCK45()
if (crystallite_todo(g,i,e)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
* crystallite_subdt(g,i,e)
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -1872,7 +1905,6 @@ subroutine crystallite_integrateStateRKCK45()
forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) &
relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s)
!$OMP FLUSH(relStateResiduum)
crystallite_todo(g,i,e) = &
( all( abs(relStateResiduum(:,g,i,e)) < rTol_crystalliteState &
.or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) ))
@ -1918,7 +1950,7 @@ subroutine crystallite_integrateStateRKCK45()
!$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
if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
call constitutive_microstructure(crystallite_temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
endif
enddo; enddo; enddo
@ -1981,7 +2013,8 @@ end subroutine crystallite_integrateStateRKCK45
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateAdaptiveEuler()
use debug, only: debug_level, &
use debug, only: &
debug_level, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
@ -1990,16 +2023,21 @@ subroutine crystallite_integrateStateAdaptiveEuler()
debug_i, &
debug_g, &
debug_StateLoopDistribution
use numerics, only: rTol_crystalliteState, &
use numerics, only: &
rTol_crystalliteState, &
numerics_integrationMode
use FEsolving, only: FEsolving_execElem, &
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
use mesh, only: mesh_element, &
use mesh, only: &
mesh_element, &
mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_Ngrains, &
use material, only: &
homogenization_Ngrains, &
homogenization_maxNgrains
use constitutive, only: constitutive_sizeDotState, &
use constitutive, only: &
constitutive_sizeDotState, &
constitutive_maxSizeDotState, &
constitutive_state, &
constitutive_aTolState, &
@ -2009,19 +2047,23 @@ subroutine crystallite_integrateStateAdaptiveEuler()
constitutive_microstructure
implicit none
!*** local variables ***!
integer(pInt) e, & ! element index in element loop
integer(pInt) :: &
e, & ! element index in element loop
i, & ! integration point index in ip loop
g, & ! grain index in grain loop
mySizeDotState, & ! size of dot State
s ! state index
integer(pInt), dimension(2) :: eIter ! bounds for element iteration
integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
integer(pInt), dimension(2) :: &
eIter ! bounds for element iteration
integer(pInt), dimension(2,mesh_NcpElems) :: &
iIter, & ! bounds for ip iteration
gIter ! bounds for grain iteration
real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
stateResiduum, & ! residuum from evolution in micrstructure
relStateResiduum ! relative residuum from evolution in microstructure
logical :: singleRun ! flag indicating computation for single (g,i,e) triple
logical :: &
singleRun ! flag indicating computation for single (g,i,e) triple
! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP ---
@ -2045,8 +2087,9 @@ subroutine crystallite_integrateStateAdaptiveEuler()
!$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
if (crystallite_todo(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_temperature(i,e), &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -2074,9 +2117,11 @@ 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
if (crystallite_todo(g,i,e)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
stateResiduum(1:mySizeDotState,g,i,e) = - 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state
stateResiduum(1:mySizeDotState,g,i,e) = - 0.5_pReal * constitutive_dotState(g,i,e)%p &
* crystallite_subdt(g,i,e) ! contribution to absolute residuum in state
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
* crystallite_subdt(g,i,e)
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -2105,7 +2150,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
!$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
if (crystallite_todo(g,i,e)) &
call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
call constitutive_microstructure(crystallite_temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
enddo; enddo; enddo
!$OMP ENDDO
@ -2137,8 +2182,9 @@ subroutine crystallite_integrateStateAdaptiveEuler()
!$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
if (crystallite_todo(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_temperature(i,e), &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
enddo; enddo; enddo
!$OMP ENDDO
!$OMP DO
@ -2173,7 +2219,8 @@ subroutine crystallite_integrateStateAdaptiveEuler()
! --- contribution of heun step to absolute residui ---
stateResiduum(1:mySizeDotState,g,i,e) = stateResiduum(1:mySizeDotState,g,i,e) &
+ 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state
+ 0.5_pReal * constitutive_dotState(g,i,e)%p &
* crystallite_subdt(g,i,e) ! contribution to absolute residuum in state
!$OMP FLUSH(stateResiduum)
! --- relative residui ---
@ -2259,9 +2306,8 @@ end subroutine crystallite_integrateStateAdaptiveEuler
!> @brief integrate stress, and state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateEuler()
use numerics, only: numerics_integrationMode, &
numerics_timeSyncing
use debug, only: debug_level, &
use debug, only: &
debug_level, &
debug_crystallite, &
debug_levelBasic, &
debug_levelExtensive, &
@ -2270,12 +2316,19 @@ subroutine crystallite_integrateStateEuler()
debug_i, &
debug_g, &
debug_StateLoopDistribution
use FEsolving, only: FEsolving_execElem, &
use numerics, only: &
numerics_integrationMode, &
numerics_timeSyncing
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
use mesh, only: mesh_element, &
use mesh, only: &
mesh_element, &
mesh_NcpElems
use material, only: homogenization_Ngrains
use constitutive, only: constitutive_sizeDotState, &
use material, only: &
homogenization_Ngrains
use constitutive, only: &
constitutive_sizeDotState, &
constitutive_state, &
constitutive_subState0, &
constitutive_dotState, &
@ -2283,15 +2336,20 @@ subroutine crystallite_integrateStateEuler()
constitutive_microstructure
implicit none
!*** local variables ***!
integer(pInt) e, & ! element index in element loop
integer(pInt) :: &
e, & ! element index in element loop
i, & ! integration point index in ip loop
g, & ! grain index in grain loop
mySizeDotState
integer(pInt), dimension(2) :: eIter ! bounds for element iteration
integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
mySizeDotState ! size of dot State
integer(pInt), dimension(2) :: &
eIter ! bounds for element iteration
integer(pInt), dimension(2,mesh_NcpElems) :: &
iIter, & ! bounds for ip iteration
gIter ! bounds for grain iteration
logical singleRun ! flag indicating computation for single (g,i,e) triple
logical :: &
singleRun ! flag indicating computation for single (g,i,e) triple
eIter = FEsolving_execElem(1:2)
do e = eIter(1),eIter(2)
@ -2311,8 +2369,9 @@ eIter = FEsolving_execElem(1:2)
!$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
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_temperature(i,e), &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
enddo; enddo; enddo
!$OMP ENDDO
!$OMP DO
@ -2340,7 +2399,8 @@ eIter = FEsolving_execElem(1:2)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
* crystallite_subdt(g,i,e)
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
@ -2379,7 +2439,7 @@ eIter = FEsolving_execElem(1:2)
!$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
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
call constitutive_microstructure(crystallite_temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
enddo; enddo; enddo
!$OMP ENDDO
@ -2441,7 +2501,8 @@ end subroutine crystallite_integrateStateEuler
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateFPI()
use debug, only: debug_e, &
use debug, only: &
debug_e, &
debug_i, &
debug_g, &
debug_level,&
@ -2450,15 +2511,20 @@ subroutine crystallite_integrateStateFPI()
debug_levelExtensive, &
debug_levelSelective, &
debug_StateLoopDistribution
use numerics, only: nState, &
use numerics, only: &
nState, &
numerics_integrationMode, &
rTol_crystalliteState
use FEsolving, only: FEsolving_execElem, &
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
use mesh, only: mesh_element, &
use mesh, only: &
mesh_element, &
mesh_NcpElems
use material, only: homogenization_Ngrains
use constitutive, only: constitutive_subState0, &
use material, only: &
homogenization_Ngrains
use constitutive, only: &
constitutive_subState0, &
constitutive_state, &
constitutive_sizeDotState, &
constitutive_maxSizeDotState, &
@ -2477,16 +2543,20 @@ subroutine crystallite_integrateStateFPI()
i, & !< integration point index in ip loop
g, & !< grain index in grain loop
mySizeDotState
integer(pInt), dimension(2) :: eIter ! bounds for element iteration
integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
integer(pInt), dimension(2) :: &
eIter ! bounds for element iteration
integer(pInt), dimension(2,mesh_NcpElems) :: &
iIter, & ! bounds for ip iteration
gIter ! bounds for grain iteration
real(pReal) dot_prod12, &
real(pReal) :: &
dot_prod12, &
dot_prod22, &
stateDamper ! damper for integration of state
real(pReal), dimension(constitutive_maxSizeDotState) :: &
stateResiduum, &
tempState
logical :: singleRun ! flag indicating computation for single (g,i,e) triple
logical :: &
singleRun ! flag indicating computation for single (g,i,e) triple
eIter = FEsolving_execElem(1:2)
do e = eIter(1),eIter(2)
@ -2499,23 +2569,19 @@ subroutine crystallite_integrateStateFPI()
! --+>> PREGUESS FOR STATE <<+--
! --- DOT STATES ---
!$OMP PARALLEL
!$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
constitutive_previousDotState(g,i,e)%p = 0.0_pReal
constitutive_previousDotState2(g,i,e)%p = 0.0_pReal
enddo; enddo; enddo
!$OMP ENDDO
! --- DOT STATES ---
!$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
if (crystallite_todo(g,i,e)) & ! ToDo: Put in loop above?
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
if (crystallite_todo(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_temperature(i,e), &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
endif
enddo; enddo; enddo
!$OMP ENDDO
!$OMP DO
@ -2525,7 +2591,7 @@ subroutine crystallite_integrateStateFPI()
if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) ) then ! NaN occured in dotState
if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity !< @ToDo: no (g,i,e) index here? ...all non-locals done (and broken)
crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
!$OMP END CRITICAL (checkTodo)
else ! broken one was local...
crystallite_todo(g,i,e) = .false. ! ... done (and broken)
@ -2543,7 +2609,8 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
+ constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e)
+ constitutive_dotState(g,i,e)%p &
* crystallite_subdt(g,i,e)
endif
enddo; enddo; enddo
!$OMP ENDDO
@ -2565,7 +2632,7 @@ subroutine crystallite_integrateStateFPI()
!$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
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_microstructure(crystallite_Temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
call constitutive_microstructure(crystallite_temperature(i,e), crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! remember previous dotState
constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p ! remember current dotState
@ -2604,8 +2671,9 @@ subroutine crystallite_integrateStateFPI()
!$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
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_temperature(i,e), &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
enddo; enddo; enddo
!$OMP ENDDO
!$OMP DO
@ -2651,10 +2719,12 @@ subroutine crystallite_integrateStateFPI()
stateResiduum(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
- constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
- (constitutive_dotState(g,i,e)%p * statedamper &
+ constitutive_previousDotState(g,i,e)%p * (1.0_pReal - statedamper)) * crystallite_subdt(g,i,e)
+ constitutive_previousDotState(g,i,e)%p &
* (1.0_pReal - statedamper)) * crystallite_subdt(g,i,e)
! --- correct state with residuum ---
tempState(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) - stateResiduum(1:mySizeDotState) ! need to copy to local variable, since we cant flush a pointer in openmp
tempState(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
- stateResiduum(1:mySizeDotState) ! need to copy to local variable, since we cant flush a pointer in openmp
#ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt &
@ -2670,7 +2740,8 @@ subroutine crystallite_integrateStateFPI()
! --- store corrected dotState --- (cannot do this before state update, because not sure how to flush pointers in openmp)
constitutive_dotState(g,i,e)%p = constitutive_dotState(g,i,e)%p * statedamper &
+ constitutive_previousDotState(g,i,e)%p * (1.0_pReal - statedamper)
+ constitutive_previousDotState(g,i,e)%p &
* (1.0_pReal - statedamper)
! --- converged ? ---
@ -2756,23 +2827,28 @@ logical function crystallite_stateJump(g,i,e)
debug_e, &
debug_i, &
debug_g
use FEsolving, only: FEsolving_execElem, &
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
use mesh, only: mesh_element, &
use mesh, only: &
mesh_element, &
mesh_NcpElems
use material, only: homogenization_Ngrains
use constitutive, only: constitutive_sizeDotState, &
use material, only: &
homogenization_Ngrains
use constitutive, only: &
constitutive_sizeDotState, &
constitutive_state, &
constitutive_deltaState, &
constitutive_collectDeltaState, &
constitutive_microstructure
implicit none
integer(pInt), intent(in):: e, & ! element index
integer(pInt), intent(in):: &
e, & ! element index
i, & ! integration point index
g ! grain index
integer(pInt) :: mySizeDotState
integer(pInt) :: &
mySizeDotState
crystallite_stateJump = .false.
@ -2780,8 +2856,10 @@ logical function crystallite_stateJump(g,i,e)
call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,g,i,e), g,i,e)
mySizeDotState = constitutive_sizeDotState(g,i,e)
if (any(constitutive_deltaState(g,i,e)%p(1:mySizeDotState) /= constitutive_deltaState(g,i,e)%p(1:mySizeDotState))) &
if (any(constitutive_deltaState(g,i,e)%p(1:mySizeDotState) &
/= constitutive_deltaState(g,i,e)%p(1:mySizeDotState))) then
return
endif
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
+ constitutive_deltaState(g,i,e)%p(1:mySizeDotState)
@ -2923,7 +3001,8 @@ logical function crystallite_integrateStress(&
if (present(timeFraction)) then
dt = crystallite_subdt(g,i,e) * timeFraction
Fg_new = crystallite_subF0(1:3,1:3,g,i,e) + (crystallite_subF(1:3,1:3,g,i,e) - crystallite_subF0(1:3,1:3,g,i,e)) * timeFraction
Fg_new = crystallite_subF0(1:3,1:3,g,i,e) &
+ (crystallite_subF(1:3,1:3,g,i,e) - crystallite_subF0(1:3,1:3,g,i,e)) * timeFraction
else
dt = crystallite_subdt(g,i,e)
Fg_new = crystallite_subF(1:3,1:3,g,i,e)
@ -2984,10 +3063,12 @@ logical function crystallite_integrateStress(&
!* calculate plastic velocity gradient and its tangent from constitutive law
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
endif
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, Tstar_v, crystallite_Temperature(i,e), g, i, e)
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, Tstar_v, &
crystallite_temperature(i,e), g, i, e)
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
@ -3110,7 +3191,8 @@ logical function crystallite_integrateStress(&
!* add volumetric component to 2nd Piola-Kirchhoff stress and calculate 1st Piola-Kirchhoff stress
forall (n=1_pInt:3_pInt) Tstar_v(n) = Tstar_v(n) + p_hydro
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_new, math_mul33x33(math_Mandel6to33(Tstar_v), math_transpose33(invFp_new)))
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_new, math_mul33x33(math_Mandel6to33(Tstar_v), &
math_transpose33(invFp_new)))
!* store local values in global variables
@ -3152,27 +3234,34 @@ end function crystallite_integrateStress
!> @brief calculates orientations and disorientations (in case of single grain ips)
!--------------------------------------------------------------------------------------------------
subroutine crystallite_orientations
use math, only: math_pDecomposition, &
use math, only: &
math_pDecomposition, &
math_RtoQ, &
math_qDisorientation, &
math_qConj
use FEsolving, only: FEsolving_execElem, &
use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP
use IO, only: IO_warning
use material, only: material_phase, &
use IO, only: &
IO_warning
use material, only: &
material_phase, &
homogenization_Ngrains, &
phase_localPlasticity, &
phase_plasticityInstance
use mesh, only: mesh_element, &
use mesh, only: &
mesh_element, &
mesh_ipNeighborhood, &
FE_NipNeighbors, &
FE_geomtype, &
FE_celltype
use constitutive_nonlocal, only: constitutive_nonlocal_structure, &
use constitutive_nonlocal, only: &
constitutive_nonlocal_structure, &
constitutive_nonlocal_updateCompatibility
implicit none
integer(pInt) e, & ! element index
integer(pInt) &
e, & ! element index
i, & ! integration point index
g, & ! grain index
n, & ! neighbor index
@ -3184,9 +3273,13 @@ subroutine crystallite_orientations
neighboringInstance, &
myStructure, & ! lattice structure
neighboringStructure
real(pReal), dimension(3,3) :: U, R
real(pReal), dimension(4) :: orientation
logical error
real(pReal), dimension(3,3) :: &
U, &
R
real(pReal), dimension(4) :: &
orientation
logical &
error
! --- CALCULATE ORIENTATION AND LATTICE ROTATION ---
@ -3303,17 +3396,26 @@ function crystallite_postResults(ipc, ip, el)
constitutive_homogenizedC
implicit none
integer(pInt), intent(in):: el, & !< element index
integer(pInt), intent(in):: &
el, & !< element index
ip, & !< integration point index
ipc !< grain index
real(pReal), dimension(1+crystallite_sizePostResults(microstructure_crystallite(mesh_element(4,el)))+ &
1+constitutive_sizePostResults(ipc,ip,el)) :: crystallite_postResults
real(pReal), dimension(3,3) :: Ee
real(pReal), dimension(4) :: rotation
real(pReal) :: detF
integer(pInt) :: o,c,crystID,mySize,n
1+constitutive_sizePostResults(ipc,ip,el)) :: &
crystallite_postResults
real(pReal), dimension(3,3) :: &
Ee
real(pReal), dimension(4) :: &
rotation
real(pReal) :: &
detF
integer(pInt) :: &
o, &
c, &
crystID, &
mySize, &
n
crystID = microstructure_crystallite(mesh_element(4,el))
@ -3334,8 +3436,8 @@ function crystallite_postResults(ipc, ip, el)
case ('volume')
mySize = 1_pInt
detF = math_det33(crystallite_partionedF(1:3,1:3,ipc,ip,el)) ! V_current = det(F) * V_reference
crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) / &
homogenization_Ngrains(mesh_element(3,el)) ! grain volume (not fraction but absolute)
crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) &
/ homogenization_Ngrains(mesh_element(3,el)) ! grain volume (not fraction but absolute)
case ('heat')
mySize = 1_pInt
crystallite_postResults(c+1) = crystallite_heat(ipc,ip,el) ! heat production
@ -3344,8 +3446,8 @@ function crystallite_postResults(ipc, ip, el)
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,ipc,ip,el) ! grain orientation as quaternion
case ('eulerangles')
mySize = 3_pInt
crystallite_postResults(c+1:c+mySize) = inDeg * &
math_qToEuler(crystallite_orientation(1:4,ipc,ip,el)) ! grain orientation as Euler angles in degree
crystallite_postResults(c+1:c+mySize) = inDeg &
* math_qToEuler(crystallite_orientation(1:4,ipc,ip,el)) ! grain orientation as Euler angles in degree
case ('grainrotation')
mySize = 4_pInt
crystallite_postResults(c+1:c+mySize) = &
@ -3423,7 +3525,7 @@ function crystallite_postResults(ipc, ip, el)
if (constitutive_sizePostResults(ipc,ip,el) > 0_pInt) &
crystallite_postResults(c+1:c+constitutive_sizePostResults(ipc,ip,el)) = &
constitutive_postResults(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe, &
crystallite_Temperature(ip,el), ipc, ip, el)
crystallite_temperature(ip,el), ipc, ip, el)
c = c + constitutive_sizePostResults(ipc,ip,el)
end function crystallite_postResults