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

homogenization_RGC.f90 >> just switching-off one of the debugging statement. 

crystallite.f90 >> optimizing the perturbation algorithm for local tangent.
This commit is contained in:
Denny Tjahjanto 2009-12-22 12:28:02 +00:00
parent 1f7aebfa4d
commit 043356e8a9
2 changed files with 94 additions and 154 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

154
code/crystallite.f90
View File

@ -19,7 +19,7 @@ implicit none
! **************************************************************** ! ****************************************************************
! *** General variables for the crystallite calculation *** ! *** General variables for the crystallite calculation ***
! **************************************************************** ! ****************************************************************
integer(pInt), parameter :: crystallite_Nresults = 5_pInt ! phaseID, volume, Euler angles, def gradient integer(pInt), parameter :: crystallite_Nresults = 14_pInt ! phaseID, volume, Euler angles, def gradient
real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, & ! requested time increment of each grain real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, & ! requested time increment of each grain
crystallite_subdt, & ! substepped time increment of each grain crystallite_subdt, & ! substepped time increment of each grain
@ -292,7 +292,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!*** output variables ***! !*** output variables ***!
!*** local variables ***! !*** local variables ***!
real(pReal) myTemperature ! local copy of the temperature real(pReal) myTemperature, & ! local copy of the temperature
myPert ! perturbation with correct sign
real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient
Fe_guess, & ! guess for elastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient
Tstar ! 2nd Piola-Kirchhoff stress tensor Tstar ! 2nd Piola-Kirchhoff stress tensor
@ -303,6 +304,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
g, gg, & ! grain index g, gg, & ! grain index
k, & k, &
l, & l, &
perturbation , & ! loop counter for forward,backward perturbation mode
comp, & comp, &
myNgrains, & myNgrains, &
mySizeState, & mySizeState, &
@ -313,7 +315,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
stateConverged, & ! flag indicating if state converged stateConverged, & ! flag indicating if state converged
converged ! flag indicating if iteration converged converged ! flag indicating if iteration converged
real(pReal), dimension(9,9) :: dPdF99 real(pReal), dimension(9,9) :: dPdF99
real(pReal), dimension(3,3,3,3) :: dPdF_pos,dPdF_neg real(pReal), dimension(3,3,3,3,2) :: dPdF_perturbation
real(pReal), dimension(constitutive_maxSizeDotState) :: delta_dotState1, & ! difference between current and previous dotstate real(pReal), dimension(constitutive_maxSizeDotState) :: delta_dotState1, & ! difference between current and previous dotstate
delta_dotState2 ! difference between previousDotState and previousDotState2 delta_dotState2 ! difference between previousDotState and previousDotState2
real(pReal) dot_prod12, & real(pReal) dot_prod12, &
@ -386,7 +388,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) 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 do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1) ! debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_converged(g,i,e)) then if (crystallite_converged(g,i,e)) then
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -447,8 +449,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
.and. crystallite_onTrack & .and. crystallite_onTrack &
.and. .not. crystallite_converged) .and. .not. crystallite_converged)
crystallite_statedamper = 1.0_pReal
! --+>> preguess for state <<+-- ! --+>> preguess for state <<+--
! !
! incrementing by crystallite_subdt ! incrementing by crystallite_subdt
@ -483,7 +483,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
crystallite_statedamper = 1.0_pReal crystallite_statedamper = 1.0_pReal
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
@ -506,7 +508,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
.and. NiterationState < nState) ! convergence loop for crystallite .and. NiterationState < nState) ! convergence loop for crystallite
NiterationState = NiterationState + 1_pInt NiterationState = NiterationState + 1_pInt
! --+>> stress integration <<+-- ! --+>> stress integration <<+--
! !
! incrementing by crystallite_subdt ! incrementing by crystallite_subdt
@ -522,10 +523,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! debugger = (e == 1 .and. i == 1 .and. g == 1) ! debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_todo(g,i,e)) & ! all undone crystallites if (crystallite_todo(g,i,e)) & ! all undone crystallites
crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e) crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e)
enddo enddo; enddo; enddo
enddo
enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) count(crystallite_onTrack(:,:,:)),'grains onTrack after stress integration' write(6,*) count(crystallite_onTrack(:,:,:)),'grains onTrack after stress integration'
@ -533,6 +533,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif endif
crystallite_todo = crystallite_todo .and. crystallite_onTrack ! continue with non-broken grains crystallite_todo = crystallite_todo .and. crystallite_onTrack ! continue with non-broken grains
if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & ! any non-local is broken? if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & ! any non-local is broken?
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped
@ -562,7 +563,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
crystallite_statedamper = 1.0_pReal crystallite_statedamper = 1.0_pReal
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
@ -585,6 +588,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
@ -628,13 +632,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif endif
enddo ! crystallite convergence loop enddo ! crystallite convergence loop
NiterationCrystallite = NiterationCrystallite + 1 NiterationCrystallite = NiterationCrystallite + 1
enddo ! cutback loop enddo ! cutback loop
! ------ check for non-converged crystallites ------ ! ------ check for non-converged crystallites ------
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
@ -689,6 +691,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) 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 do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
! debugger = (e==1 .and. i==1 .and. g==1)
if (crystallite_requested(g,i,e)) then ! first check whether is requested at all! if (crystallite_requested(g,i,e)) then ! first check whether is requested at all!
if (crystallite_converged(g,i,e)) then ! grain converged in above iteration if (crystallite_converged(g,i,e)) then ! grain converged in above iteration
if (debugger) then if (debugger) then
@ -702,10 +705,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
if (pert_method == 1_pInt .or. pert_method == 3_pInt) then ! <<< when forward or central difference is desired >>> do perturbation = 1,2
if (iand(pert_method,perturbation) > 0) then
myPert = -pert_Fg * (-1.0_pReal)**perturbation ! forward or backward perturbation
do k = 1,3 ! perturbation... do k = 1,3 ! perturbation...
do l = 1,3 ! ...components to the positive direction do l = 1,3 ! ...components to the positive direction
crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + myPert ! perturb single component (either forward or backward)
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,'(i1,x,i1)') k,l write (6,'(i1,x,i1)') k,l
@ -723,7 +728,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), &
crystallite_misorientation(:,:,g,i,e), crystallite_subdt(g,i,e), & crystallite_misorientation(:,:,g,i,e), crystallite_subdt(g,i,e), &
g, i, e) g,i,e)
stateConverged = crystallite_updateState(g,i,e) ! update state stateConverged = crystallite_updateState(g,i,e) ! update state
temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature
converged = stateConverged .and. temperatureConverged converged = stateConverged .and. temperatureConverged
@ -732,106 +737,41 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write (6,*) '-------------' write (6,*) '-------------'
write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged
write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 write (6,'(a,/,3(3(f12.4,x)/))') 'pertP/MPa of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6
write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-storedP(1:3,:,g,i,e))/1e6 write (6,'(a,/,3(3(f12.4,x)/))') 'DP/MPa of 1 1 1',(crystallite_P(1:3,:,g,i,e)-storedP(1:3,:,g,i,e))/1e6
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
enddo enddo
if (converged) then ! converged state warrants stiffness update if (converged) & ! converged state warrants stiffness update
dPdF_pos(:,:,k,l) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/pert_Fg ! tangent dP_ij/dFg_kl dPdF_perturbation(:,:,k,l,perturbation) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/myPert ! tangent dP_ij/dFg_kl
if (pert_method == 1_pInt) crystallite_dPdF(:,:,k,l,g,i,e) = dPdF_pos(:,:,k,l)
endif mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain
do ee = FEsolving_execElem(1),FEsolving_execElem(2) mySizeDotState = constitutive_sizeDotState(g,i,e) ! number of dotStates for this grain
myNgrains = homogenization_Ngrains(mesh_element(3,ee)) constitutive_state(g,i,e)%p = storedState(1:mySizeState,g,i,e)
do ii = FEsolving_execIP(1,ee),FEsolving_execIP(2,ee) constitutive_dotState(g,i,e)%p = storedDotState(1:mySizeDotState,g,i,e)
do gg = 1,myNgrains crystallite_Temperature(g,i,e) = storedTemperature(g,i,e)
mySizeState = constitutive_sizeState(gg,ii,ee) ! number of state variables for this grain crystallite_subF(:,:,g,i,e) = storedF(:,:,g,i,e)
mySizeDotState = constitutive_sizeDotState(gg,ii,ee) ! number of dotStates for this grain crystallite_Fp(:,:,g,i,e) = storedFp(:,:,g,i,e)
constitutive_state(gg,ii,ee)%p = storedState(1:mySizeState,gg,ii,ee) crystallite_invFp(:,:,g,i,e) = storedInvFp(:,:,g,i,e)
constitutive_dotState(gg,ii,ee)%p = storedDotState(1:mySizeDotState,gg,ii,ee) crystallite_Fe(:,:,g,i,e) = storedFe(:,:,g,i,e)
crystallite_Temperature(gg,ii,ee) = storedTemperature(gg,ii,ee) crystallite_Lp(:,:,g,i,e) = storedLp(:,:,g,i,e)
crystallite_subF(:,:,gg,ii,ee) = storedF(:,:,gg,ii,ee) crystallite_Tstar_v(:,g,i,e) = storedTstar_v(:,g,i,e)
crystallite_Fp(:,:,gg,ii,ee) = storedFp(:,:,gg,ii,ee) crystallite_P(:,:,g,i,e) = storedP(:,:,g,i,e)
crystallite_invFp(:,:,gg,ii,ee) = storedInvFp(:,:,gg,ii,ee)
crystallite_Fe(:,:,gg,ii,ee) = storedFe(:,:,gg,ii,ee)
crystallite_Lp(:,:,gg,ii,ee) = storedLp(:,:,gg,ii,ee)
crystallite_Tstar_v(:,gg,ii,ee) = storedTstar_v(:,gg,ii,ee)
crystallite_P(:,:,gg,ii,ee) = storedP(:,:,gg,ii,ee)
enddo; enddo; enddo
!$OMP CRITICAL (out) !$OMP CRITICAL (out)
debug_StiffnessStateLoopDistribution(NiterationState) = & debug_StiffnessStateLoopDistribution(NiterationState) = &
debug_StiffnessstateLoopDistribution(NiterationState) + 1 debug_StiffnessstateLoopDistribution(NiterationState) + 1
!$OMPEND CRITICAL (out) !$OMPEND CRITICAL (out)
enddo enddo; enddo
enddo
endif endif
enddo ! perturbation direction
if (pert_method == 2_pInt .or. pert_method == 3_pInt) then ! <<< when backward or central difference is desired >>> select case(pert_method)
do k = 1,3 ! perturbation... case (1)
do l = 1,3 ! ...components to the negative direction crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,1)
crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) - pert_Fg ! perturb single component case (2)
if (debugger) then crystallite_dPdF(:,:,:,:,g,i,e) = dPdF_perturbation(:,:,:,:,2)
!$OMP CRITICAL (write2out) case (3)
write (6,'(i1,x,i1)') k,l crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_perturbation(:,:,:,:,1)+dPdF_perturbation(:,:,:,:,2))
write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) end select
!$OMPEND CRITICAL (write2out)
endif
onTrack = .true.
converged = .false.
NiterationState = 0_pInt
do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged)
NiterationState = NiterationState + 1_pInt
onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe)
if (onTrack) then
constitutive_dotState(g,i,e)%p = 0.0_pReal
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), &
crystallite_misorientation(:,:,g,i,e), crystallite_subdt(g,i,e), &
g, i, e)
stateConverged = crystallite_updateState(g,i,e) ! update state
temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature
converged = stateConverged .and. temperatureConverged
endif
if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) '-------------'
write (6,'(a,x,l,x,l)') 'ontrack + converged:',onTrack,converged
write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6
write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-storedP(1:3,:,g,i,e))/1e6
!$OMPEND CRITICAL (write2out)
endif
enddo
if (converged) then ! converged state warrants stiffness update
dPdF_neg(:,:,k,l) = (storedP(:,:,g,i,e) - crystallite_P(:,:,g,i,e))/pert_Fg ! tangent dP_ij/dFg_kl
if (pert_method == 2_pInt) crystallite_dPdF(:,:,k,l,g,i,e) = dPdF_neg(:,:,k,l)
endif
do ee = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,ee))
do ii = FEsolving_execIP(1,ee),FEsolving_execIP(2,ee)
do gg = 1,myNgrains
mySizeState = constitutive_sizeState(gg,ii,ee) ! number of state variables for this grain
mySizeDotState = constitutive_sizeDotState(gg,ii,ee) ! number of dotStates for this grain
constitutive_state(gg,ii,ee)%p = storedState(1:mySizeState,gg,ii,ee)
constitutive_dotState(gg,ii,ee)%p = storedDotState(1:mySizeDotState,gg,ii,ee)
crystallite_Temperature(gg,ii,ee) = storedTemperature(gg,ii,ee)
crystallite_subF(:,:,gg,ii,ee) = storedF(:,:,gg,ii,ee)
crystallite_Fp(:,:,gg,ii,ee) = storedFp(:,:,gg,ii,ee)
crystallite_invFp(:,:,gg,ii,ee) = storedInvFp(:,:,gg,ii,ee)
crystallite_Fe(:,:,gg,ii,ee) = storedFe(:,:,gg,ii,ee)
crystallite_Lp(:,:,gg,ii,ee) = storedLp(:,:,gg,ii,ee)
crystallite_Tstar_v(:,gg,ii,ee) = storedTstar_v(:,gg,ii,ee)
crystallite_P(:,:,gg,ii,ee) = storedP(:,:,gg,ii,ee)
enddo; enddo; enddo
!$OMP CRITICAL (out)
debug_StiffnessStateLoopDistribution(NiterationState) = &
debug_StiffnessstateLoopDistribution(NiterationState) + 1
!$OMPEND CRITICAL (out)
enddo
enddo
endif
if (pert_method == 3_pInt) crystallite_dPdF(:,:,:,:,g,i,e) = 0.5_pReal*(dPdF_neg + dPdF_pos)
else ! grain did not converge else ! grain did not converge
crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use (elastic) fallback
endif ! grain convergence endif ! grain convergence

2
code/homogenization_RGC.f90
View File

@ -334,7 +334,7 @@ function homogenization_RGC_updateState(&
call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID) call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID)
!* Debugging the mismatch, stress and penalties of grains !* Debugging the mismatch, stress and penalties of grains
if (el == 1 .and. ip == 1) then if (RGCdebug) then
do iGrain = 1,nGrain do iGrain = 1,nGrain
write(6,'(x,a30,x,i3,x,a4,x,e14.8)')'Mismatch magnitude of grain(',iGrain,') :',NN(iGrain) write(6,'(x,a30,x,i3,x,a4,x,e14.8)')'Mismatch magnitude of grain(',iGrain,') :',NN(iGrain)
write(6,*)' ' write(6,*)' '