Depending on which state integrator one uses for the stiffness calculation, the initial state has to be chosen accordingly: e.g. for FPI choose last converged state, but for explicit RK choose converged state from start of increment (in case of explicit euler no state integration at all, but only stress integration). For this purpose we also need to remember "Fe" which now follows the cutbacking procedure as it is used for "Fp".

This commit is contained in:
Christoph Kords 2011-11-04 12:44:50 +00:00
parent 4966231e5b
commit ca3d21a3b6
1 changed files with 209 additions and 147 deletions

View File

@ -66,6 +66,7 @@ real(pReal), dimension (:,:,:,:), allocatable :: &
crystallite_rotation ! grain rotation away from initial orientation as axis-angle (in degrees) crystallite_rotation ! grain rotation away from initial orientation as axis-angle (in degrees)
real(pReal), dimension (:,:,:,:,:), allocatable :: & real(pReal), dimension (:,:,:,:,:), allocatable :: &
crystallite_Fe, & ! current "elastic" def grad (end of converged time step) crystallite_Fe, & ! current "elastic" def grad (end of converged time step)
crystallite_subFe0,& ! "elastic" def grad at start of crystallite inc
crystallite_Fp, & ! current plastic def grad (end of converged time step) crystallite_Fp, & ! current plastic def grad (end of converged time step)
crystallite_invFp, & ! inverse of current plastic def grad (end of converged time step) crystallite_invFp, & ! inverse of current plastic def grad (end of converged time step)
crystallite_Fp0, & ! plastic def grad at start of FE inc crystallite_Fp0, & ! plastic def grad at start of FE inc
@ -202,6 +203,7 @@ allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crys
allocate(crystallite_Fp(3,3,gMax,iMax,eMax)); crystallite_Fp = 0.0_pReal allocate(crystallite_Fp(3,3,gMax,iMax,eMax)); crystallite_Fp = 0.0_pReal
allocate(crystallite_invFp(3,3,gMax,iMax,eMax)); crystallite_invFp = 0.0_pReal allocate(crystallite_invFp(3,3,gMax,iMax,eMax)); crystallite_invFp = 0.0_pReal
allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal
allocate(crystallite_subFe0(3,3,gMax,iMax,eMax)); crystallite_subFe0 = 0.0_pReal
allocate(crystallite_Lp0(3,3,gMax,iMax,eMax)); crystallite_Lp0 = 0.0_pReal allocate(crystallite_Lp0(3,3,gMax,iMax,eMax)); crystallite_Lp0 = 0.0_pReal
allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax)); crystallite_partionedLp0 = 0.0_pReal allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax)); crystallite_partionedLp0 = 0.0_pReal
allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal
@ -400,6 +402,7 @@ if (debug_verbosity > 0) then
write(6,'(a35,x,7(i8,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0) write(6,'(a35,x,7(i8,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0)
write(6,'(a35,x,7(i8,x))') 'crystallite_symmetryID: ', shape(crystallite_symmetryID) write(6,'(a35,x,7(i8,x))') 'crystallite_symmetryID: ', shape(crystallite_symmetryID)
write(6,'(a35,x,7(i8,x))') 'crystallite_subF0: ', shape(crystallite_subF0) write(6,'(a35,x,7(i8,x))') 'crystallite_subF0: ', shape(crystallite_subF0)
write(6,'(a35,x,7(i8,x))') 'crystallite_subFe0: ', shape(crystallite_subFe0)
write(6,'(a35,x,7(i8,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) write(6,'(a35,x,7(i8,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0)
write(6,'(a35,x,7(i8,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0) write(6,'(a35,x,7(i8,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0)
write(6,'(a35,x,7(i8,x))') 'crystallite_P: ', shape(crystallite_P) write(6,'(a35,x,7(i8,x))') 'crystallite_P: ', shape(crystallite_P)
@ -561,6 +564,8 @@ crystallite_subStep = 0.0_pReal
crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness
crystallite_subF0(1:3,1:3,g,i,e) = crystallite_partionedF0(1:3,1:3,g,i,e) ! ...def grad crystallite_subF0(1:3,1:3,g,i,e) = crystallite_partionedF0(1:3,1:3,g,i,e) ! ...def grad
crystallite_subTstar0_v(1:6,g,i,e) = crystallite_partionedTstar0_v(1:6,g,i,e) !...2nd PK stress crystallite_subTstar0_v(1:6,g,i,e) = crystallite_partionedTstar0_v(1:6,g,i,e) !...2nd PK stress
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF0(1:3,1:3,g,i,e), &
math_inv3x3(crystallite_subFp0(1:3,1:3,g,i,e))) ! only needed later on for stiffness calculation
crystallite_subFrac(g,i,e) = 0.0_pReal crystallite_subFrac(g,i,e) = 0.0_pReal
crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst
@ -607,6 +612,7 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward...
crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad
crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation
crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient
constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure
crystallite_subTstar0_v(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) ! ...2nd PK stress crystallite_subTstar0_v(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) ! ...2nd PK stress
@ -719,6 +725,7 @@ enddo
if(updateJaco) then ! Jacobian required if(updateJaco) then ! Jacobian required
numerics_integrationMode = 2_pInt numerics_integrationMode = 2_pInt
! --- BACKUP --- ! --- BACKUP ---
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState) !$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState)
@ -753,14 +760,63 @@ if(updateJaco) then
if (iand(pert_method,perturbation) > 0) then ! mask for desired direction if (iand(pert_method,perturbation) > 0) then ! mask for desired direction
myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step
do k = 1,3; do l = 1,3 ! ...alter individual components do k = 1,3; do l = 1,3 ! ...alter individual components
#ifndef _OPENMP
if (debug_verbosity> 5) then if (debug_verbosity> 5) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,2(x,i1),x,a)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]' write(6,'(a,2(x,i1),x,a)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]'
write(6,*) write(6,*)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
crystallite_subF(k,l,:,:,:) = crystallite_subF(k,l,:,:,:) + myPert ! perturb either forward or backward #endif
! --- INITIALIZE UNPERTURBED STATE ---
select case(numerics_integrator(numerics_integrationMode))
case(1) ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState)
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)
do g = 1,myNgrains
mySizeState = constitutive_sizeState(g,i,e)
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_state_backup(g,i,e)%p(1:mySizeState)
constitutive_dotState(g,i,e)%p(1:mySizeDotState) = constitutive_dotState_backup(g,i,e)%p(1:mySizeDotState)
enddo; enddo; enddo
!OMP END PARALLEL DO
crystallite_Temperature = Temperature_backup
crystallite_Fp = Fp_backup
crystallite_invFp = InvFp_backup
crystallite_Fe = Fe_backup
crystallite_Lp = Lp_backup
crystallite_Tstar_v = Tstar_v_backup
case(2,3) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step
case(4,5) ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState)
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)
do g = 1,myNgrains
mySizeState = constitutive_sizeState(g,i,e)
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_subState0(g,i,e)%p(1:mySizeState)
constitutive_dotState(g,i,e)%p(1:mySizeDotState) = constitutive_dotState_backup(g,i,e)%p(1:mySizeDotState)
enddo; enddo; enddo
!OMP END PARALLEL DO
crystallite_Temperature = crystallite_subTemperature0
crystallite_Fp = crystallite_subFp0
crystallite_Fe = crystallite_subFe0
crystallite_Tstar_v = crystallite_subTstar0_v
end select
! --- PERTURB EITHER FORWARD OR BACKWARD ---
crystallite_subF = F_backup
crystallite_subF(k,l,:,:,:) = crystallite_subF(k,l,:,:,:) + myPert
crystallite_converged = convergenceFlag_backup
crystallite_todo = crystallite_requested .and. crystallite_converged crystallite_todo = crystallite_requested .and. crystallite_converged
where (crystallite_todo) crystallite_converged = .false. ! start out non-converged where (crystallite_todo) crystallite_converged = .false. ! start out non-converged
@ -793,32 +849,6 @@ if(updateJaco) then
enddo; enddo; enddo enddo; enddo; enddo
!OMP END PARALLEL DO !OMP END PARALLEL DO
! --- RESTORE ---
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState)
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)
do g = 1,myNgrains
mySizeState = constitutive_sizeState(g,i,e)
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_state_backup(g,i,e)%p(1:mySizeState)
constitutive_dotState(g,i,e)%p(1:mySizeDotState) = constitutive_dotState_backup(g,i,e)%p(1:mySizeDotState)
enddo
enddo
enddo
!OMP END PARALLEL DO
crystallite_Temperature = Temperature_backup
crystallite_subF = F_backup
crystallite_Fp = Fp_backup
crystallite_invFp = InvFp_backup
crystallite_Fe = Fe_backup
crystallite_Lp = Lp_backup
crystallite_Tstar_v = Tstar_v_backup
crystallite_P = P_backup
crystallite_converged = convergenceFlag_backup
enddo; enddo ! k,l loop enddo; enddo ! k,l loop
endif endif
enddo ! perturbation direction enddo ! perturbation direction
@ -831,7 +861,7 @@ if(updateJaco) then
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,myNgrains do g = 1,myNgrains
if (crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) then ! central solution converged if (crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) then ! central solution converged
select case(pert_method) select case(pert_method)
case(1) case(1)
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e) crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e)
@ -841,7 +871,7 @@ if(updateJaco) then
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.5_pReal* ( dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e) & crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.5_pReal* ( dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e) &
+ dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e)) + dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e))
end select end select
elseif (crystallite_requested(g,i,e) .and. .not. crystallite_converged(g,i,e)) then ! central solution did not converge elseif (crystallite_requested(g,i,e) .and. .not. convergenceFlag_backup(g,i,e)) then ! central solution did not converge
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! use (elastic) fallback crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! use (elastic) fallback
endif endif
enddo enddo
@ -849,6 +879,30 @@ if(updateJaco) then
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
! --- RESTORE ---
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState)
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)
do g = 1,myNgrains
mySizeState = constitutive_sizeState(g,i,e)
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_state_backup(g,i,e)%p(1:mySizeState)
constitutive_dotState(g,i,e)%p(1:mySizeDotState) = constitutive_dotState_backup(g,i,e)%p(1:mySizeDotState)
enddo; enddo; enddo
!OMP END PARALLEL DO
crystallite_Temperature = Temperature_backup
crystallite_subF = F_backup
crystallite_Fp = Fp_backup
crystallite_invFp = InvFp_backup
crystallite_Fe = Fe_backup
crystallite_Lp = Lp_backup
crystallite_Tstar_v = Tstar_v_backup
crystallite_P = P_backup
crystallite_converged = convergenceFlag_backup
endif ! jacobian calculation endif ! jacobian calculation
endsubroutine endsubroutine
@ -1671,75 +1725,79 @@ else
endif endif
! --- RESET DOTSTATE ---
!$OMP PARALLEL PRIVATE(mySizeDotState) !$OMP PARALLEL PRIVATE(mySizeDotState)
!$OMP DO if (numerics_integrationMode < 2) then
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_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero ! --- RESET DOTSTATE ---
enddo; enddo; enddo
!$OMP ENDDO !$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_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero
enddo; enddo; enddo
!$OMP ENDDO
! --- DOT STATE AND TEMPERATURE (EULER INTEGRATION) --- ! --- DOT STATE AND TEMPERATURE (EULER INTEGRATION) ---
stateResiduum = 0.0_pReal stateResiduum = 0.0_pReal
!$OMP DO !$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 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 (crystallite_todo(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e)
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Temperature(g,i,e),g,i,e) crystallite_Temperature(g,i,e),g,i,e)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
!$OMP DO !$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 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 (crystallite_todo(g,i,e)) then
if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
.or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo) !$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo) !$OMP END CRITICAL (checkTodo)
else ! if broken local... else ! if broken local...
crystallite_todo(g,i,e) = .false. ! ... skip this one next time crystallite_todo(g,i,e) = .false. ! ... skip this one next time
endif
endif endif
endif endif
endif enddo; enddo; enddo
enddo; enddo; enddo !$OMP ENDDO
!$OMP ENDDO
! --- STATE UPDATE (EULER INTEGRATION) --- ! --- STATE UPDATE (EULER INTEGRATION) ---
!$OMP DO !$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 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 (crystallite_todo(g,i,e)) then
mySizeDotState = constitutive_sizeDotState(g,i,e) 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 and temperature 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 and temperature
temperatureResiduum(g,i,e) = - 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) temperatureResiduum(g,i,e) = - 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & 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)
crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
+ crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
! --- UPDATE DEPENDENT STATES (EULER INTEGRATION) --- ! --- UPDATE DEPENDENT STATES (EULER INTEGRATION) ---
!$OMP DO !$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 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 (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
endif
! --- STRESS INTEGRATION (EULER INTEGRATION) --- ! --- STRESS INTEGRATION (EULER INTEGRATION) ---
@ -1938,82 +1996,86 @@ else
endif endif
! --- RESET DOTSTATE ---
!$OMP PARALLEL PRIVATE(mySizeDotState) !$OMP PARALLEL PRIVATE(mySizeDotState)
!$OMP DO if (numerics_integrationMode < 2) then
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_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero ! --- RESET DOTSTATE ---
enddo; enddo; enddo
!$OMP ENDDO !$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_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero
enddo; enddo; enddo
!$OMP ENDDO
! --- DOT STATE AND TEMPERATURE --- ! --- DOT STATE AND TEMPERATURE ---
!$OMP DO !$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 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 (crystallite_todo(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e)
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Temperature(g,i,e),g,i,e) crystallite_Temperature(g,i,e),g,i,e)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
!$OMP DO !$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 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 (crystallite_todo(g,i,e)) then
if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
.or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local...
!$OMP CRITICAL (checkTodo) !$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo) !$OMP END CRITICAL (checkTodo)
else ! if broken local... else ! if broken local...
crystallite_todo(g,i,e) = .false. ! ... skip this one next time crystallite_todo(g,i,e) = .false. ! ... skip this one next time
endif
endif endif
endif endif
endif enddo; enddo; enddo
enddo; enddo; enddo !$OMP ENDDO
!$OMP ENDDO
! --- UPDATE STATE AND TEMPERATURE --- ! --- UPDATE STATE AND TEMPERATURE ---
!$OMP DO !$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 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 (crystallite_todo(g,i,e)) then
mySizeDotState = constitutive_sizeDotState(g,i,e) mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & 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)
crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) & crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
+ crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
#ifndef _OPENMP #ifndef _OPENMP
if (debug_verbosity > 5 & if (debug_verbosity > 5 &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*) write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,*) write(6,*)
write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState) write(6,'(a,/,(12(x),12(e12.5,x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*) write(6,*)
endif endif
#endif #endif
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
! --- UPDATE DEPENDENT STATES --- ! --- UPDATE DEPENDENT STATES ---
!$OMP DO !$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 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 (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
endif
! --- STRESS INTEGRATION --- ! --- STRESS INTEGRATION ---