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

polishing

This commit is contained in:
Philip Eisenlohr 2012-11-22 09:58:36 +00:00
parent e0e5683386
commit 7d196fbb25
2 changed files with 110 additions and 110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

218
code/crystallite.f90
View File

@ -583,97 +583,97 @@ NiterationCrystallite = 0_pInt
numerics_integrationMode = 1_pInt numerics_integrationMode = 1_pInt
do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2)))) ! cutback loop for crystallites do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2)))) ! cutback loop for crystallites
!$OMP PARALLEL DO PRIVATE(myNgrains,formerSubStep) !$OMP PARALLEL DO PRIVATE(myNgrains,formerSubStep)
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))
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
! --- wind forward --- ! --- wind forward ---
if (crystallite_converged(g,i,e)) then if (crystallite_converged(g,i,e)) then
formerSubStep = crystallite_subStep(g,i,e) formerSubStep = crystallite_subStep(g,i,e)
crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e) crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e)
!$OMP FLUSH(crystallite_subFrac) !$OMP FLUSH(crystallite_subFrac)
crystallite_subStep(g,i,e) = min( 1.0_pReal - crystallite_subFrac(g,i,e), & crystallite_subStep(g,i,e) = min( 1.0_pReal - crystallite_subFrac(g,i,e), &
stepIncreaseCryst * crystallite_subStep(g,i,e) ) stepIncreaseCryst * crystallite_subStep(g,i,e) )
!$OMP FLUSH(crystallite_subStep) !$OMP FLUSH(crystallite_subStep)
if (crystallite_subStep(g,i,e) > 0.0_pReal) then if (crystallite_subStep(g,i,e) > 0.0_pReal) then
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
!$OMP FLUSH(crystallite_subF0) !$OMP FLUSH(crystallite_subF0)
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_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
crystallite_todo(g,i,e) = .true. crystallite_todo(g,i,e) = .true.
!$OMP FLUSH(crystallite_todo) !$OMP FLUSH(crystallite_todo)
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt & if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,f12.8,a,f12.8,a)') '<< CRYST >> winding forward from ', &
crystallite_subFrac(g,i,e)-formerSubStep,' to current crystallite_subfrac ', &
crystallite_subFrac(g,i,e),' in crystallite_stressAndItsTangent'
write(6,*)
endif
#endif
elseif (formerSubStep > 0.0_pReal) then ! this crystallite just converged
crystallite_todo(g,i,e) = .false. ! so done here
!$OMP FLUSH(crystallite_todo)
if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionCrystallite)
debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) = &
debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) + 1_pInt
!$OMP END CRITICAL (distributionCrystallite)
endif
endif
! --- cutback ---
else
crystallite_subStep(g,i,e) = subStepSizeCryst * crystallite_subStep(g,i,e) ! cut step in half and restore...
!$OMP FLUSH(crystallite_subStep)
crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature
crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) ! ...plastic def grad
!$OMP FLUSH(crystallite_Fp)
crystallite_invFp(1:3,1:3,g,i,e) = math_inv33(crystallite_Fp(1:3,1:3,g,i,e))
!$OMP FLUSH(crystallite_invFp)
crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad
constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure
crystallite_Tstar_v(1:6,g,i,e) = crystallite_subTstar0_v(1:6,g,i,e) ! ...2nd PK stress
! cant restore dotState here, since not yet calculated in first cutback after initialization
crystallite_todo(g,i,e) = crystallite_subStep(g,i,e) > subStepMinCryst ! still on track or already done (beyond repair)
!$OMP FLUSH(crystallite_todo)
#ifndef _OPENMP
if (crystallite_todo(g,i,e) &
.and. iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,f12.8)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& write(6,'(a,f12.8,a,f12.8,a)') '<< CRYST >> winding forward from ', &
crystallite_subStep(g,i,e) crystallite_subFrac(g,i,e)-formerSubStep,' to current crystallite_subfrac ', &
crystallite_subFrac(g,i,e),' in crystallite_stressAndItsTangent'
write(6,*) write(6,*)
endif endif
#endif #endif
elseif (formerSubStep > 0.0_pReal) then ! this crystallite just converged
crystallite_todo(g,i,e) = .false. ! so done here
!$OMP FLUSH(crystallite_todo)
if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (distributionCrystallite)
debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) = &
debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) + 1_pInt
!$OMP END CRITICAL (distributionCrystallite)
endif
endif endif
! --- cutback ---
else
crystallite_subStep(g,i,e) = subStepSizeCryst * crystallite_subStep(g,i,e) ! cut step in half and restore...
!$OMP FLUSH(crystallite_subStep)
crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature
crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) ! ...plastic def grad
!$OMP FLUSH(crystallite_Fp)
crystallite_invFp(1:3,1:3,g,i,e) = math_inv33(crystallite_Fp(1:3,1:3,g,i,e))
!$OMP FLUSH(crystallite_invFp)
crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad
constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure
crystallite_Tstar_v(1:6,g,i,e) = crystallite_subTstar0_v(1:6,g,i,e) ! ...2nd PK stress
! cant restore dotState here, since not yet calculated in first cutback after initialization
crystallite_todo(g,i,e) = crystallite_subStep(g,i,e) > subStepMinCryst ! still on track or already done (beyond repair)
!$OMP FLUSH(crystallite_todo)
#ifndef _OPENMP
if (crystallite_todo(g,i,e) &
.and. iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,f12.8)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',&
crystallite_subStep(g,i,e)
write(6,*)
endif
#endif
endif
! --- prepare for integration --- ! --- prepare for integration ---
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
crystallite_subF(1:3,1:3,g,i,e) = 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) &
+ crystallite_subStep(g,i,e) & + crystallite_subStep(g,i,e) &
* (crystallite_partionedF(1:3,1:3,g,i,e) - crystallite_partionedF0(1:3,1:3,g,i,e)) * (crystallite_partionedF(1:3,1:3,g,i,e) - crystallite_partionedF0(1:3,1:3,g,i,e))
!$OMP FLUSH(crystallite_subF) !$OMP FLUSH(crystallite_subF)
crystallite_Fe(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)) crystallite_Fe(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))
crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e)
crystallite_converged(g,i,e) = .false. ! start out non-converged crystallite_converged(g,i,e) = .false. ! start out non-converged
endif endif
enddo enddo ! grains
enddo enddo ! IPs
enddo enddo ! elements
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
! --- integrate --- requires fully defined state array (basic + dependent state) ! --- integrate --- requires fully defined state array (basic + dependent state)
@ -689,15 +689,15 @@ do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2)))
call crystallite_integrateStateRK4() call crystallite_integrateStateRK4()
case(5_pInt) case(5_pInt)
call crystallite_integrateStateRKCK45() call crystallite_integrateStateRKCK45()
endselect end select
endif endif
where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) & where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) & ! do not try non-converged & fully cutbacked any further
crystallite_todo = .true. crystallite_todo = .true.
NiterationCrystallite = NiterationCrystallite + 1_pInt NiterationCrystallite = NiterationCrystallite + 1_pInt
enddo ! cutback loop enddo ! cutback loop
! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+-- ! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+--
@ -709,7 +709,7 @@ do e = FEsolving_execElem(1),FEsolving_execElem(2)
if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway)
if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> no convergence : respond fully elastic at el ip g ',e,i,g write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> no convergence: respond fully elastic at el ip g ',e,i,g
write(6,*) write(6,*)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -743,7 +743,6 @@ if(updateJaco) then
numerics_integrationMode = 2_pInt numerics_integrationMode = 2_pInt
! --- BACKUP --- ! --- BACKUP ---
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
@ -755,14 +754,14 @@ if(updateJaco) then
constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) ! ... dotStates, ... constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) ! ... dotStates, ...
endforall endforall
enddo enddo
Temperature_backup = crystallite_Temperature ! ... Temperature, ... Temperature_backup = crystallite_Temperature ! ... Temperature, ...
F_backup = crystallite_subF ! ... and kinematics F_backup = crystallite_subF ! ... and kinematics
Fp_backup = crystallite_Fp Fp_backup = crystallite_Fp
InvFp_backup = crystallite_invFp InvFp_backup = crystallite_invFp
Fe_backup = crystallite_Fe Fe_backup = crystallite_Fe
Lp_backup = crystallite_Lp Lp_backup = crystallite_Lp
Tstar_v_backup = crystallite_Tstar_v Tstar_v_backup = crystallite_Tstar_v
P_backup = crystallite_P P_backup = crystallite_P
convergenceFlag_backup = crystallite_converged convergenceFlag_backup = crystallite_converged
@ -771,7 +770,7 @@ if(updateJaco) then
dPdF_perturbation1 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment 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 dPdF_perturbation2 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
do perturbation = 1,2 ! forward and backward perturbation do perturbation = 1,2 ! forward and backward perturbation
if (iand(pert_method,perturbation) > 0) then ! mask for desired direction if (iand(pert_method,perturbation) > 0_pInt) 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
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
@ -796,11 +795,11 @@ if(updateJaco) then
endforall endforall
enddo enddo
crystallite_Temperature = Temperature_backup crystallite_Temperature = Temperature_backup
crystallite_Fp = Fp_backup crystallite_Fp = Fp_backup
crystallite_invFp = InvFp_backup crystallite_invFp = InvFp_backup
crystallite_Fe = Fe_backup crystallite_Fe = Fe_backup
crystallite_Lp = Lp_backup crystallite_Lp = Lp_backup
crystallite_Tstar_v = Tstar_v_backup crystallite_Tstar_v = Tstar_v_backup
case(2_pInt,3_pInt) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step case(2_pInt,3_pInt) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step
case(4_pInt,5_pInt) ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress case(4_pInt,5_pInt) ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
@ -813,9 +812,9 @@ if(updateJaco) then
endforall endforall
enddo enddo
crystallite_Temperature = crystallite_subTemperature0 crystallite_Temperature = crystallite_subTemperature0
crystallite_Fp = crystallite_subFp0 crystallite_Fp = crystallite_subFp0
crystallite_Fe = crystallite_subFe0 crystallite_Fe = crystallite_subFe0
crystallite_Tstar_v = crystallite_subTstar0_v crystallite_Tstar_v = crystallite_subTstar0_v
end select end select
@ -855,7 +854,8 @@ if(updateJaco) then
end select end select
enddo enddo
enddo; enddo ! k,l loop enddo; enddo ! k,l component perturbation loop
endif endif
enddo ! perturbation direction enddo ! perturbation direction
@ -898,14 +898,14 @@ if(updateJaco) then
endforall endforall
enddo enddo
crystallite_Temperature = Temperature_backup crystallite_Temperature = Temperature_backup
crystallite_subF = F_backup crystallite_subF = F_backup
crystallite_Fp = Fp_backup crystallite_Fp = Fp_backup
crystallite_invFp = InvFp_backup crystallite_invFp = InvFp_backup
crystallite_Fe = Fe_backup crystallite_Fe = Fe_backup
crystallite_Lp = Lp_backup crystallite_Lp = Lp_backup
crystallite_Tstar_v = Tstar_v_backup crystallite_Tstar_v = Tstar_v_backup
crystallite_P = P_backup crystallite_P = P_backup
crystallite_converged = convergenceFlag_backup crystallite_converged = convergenceFlag_backup
else ! Calculate Jacobian using analytical expression else ! Calculate Jacobian using analytical expression
@ -943,14 +943,14 @@ if(updateJaco) then
FDot_inv = crystallite_subF(1:3,1:3,g,i,e) - crystallite_F0(1:3,1:3,g,i,e) FDot_inv = crystallite_subF(1:3,1:3,g,i,e) - crystallite_F0(1:3,1:3,g,i,e)
counter = 0.0_pReal counter = 0.0_pReal
do p=1_pInt,3_pInt; do o=1_pInt,3_pInt do p=1_pInt,3_pInt; do o=1_pInt,3_pInt
if (abs(FDot_inv(o,p)) .lt. relevantStrain) then if (abs(FDot_inv(o,p)) < relevantStrain) then
FDot_inv(o,p) = 0.0_pReal FDot_inv(o,p) = 0.0_pReal
else else
counter = counter + 1.0_pReal counter = counter + 1.0_pReal
FDot_inv(o,p) = crystallite_dt(g,i,e)/FDot_inv(o,p) FDot_inv(o,p) = crystallite_dt(g,i,e)/FDot_inv(o,p)
endif endif
enddo; enddo enddo; enddo
if (counter .gt. 0.0_pReal) FDot_inv = FDot_inv/counter if (counter > 0.0_pReal) FDot_inv = FDot_inv/counter
do p=1_pInt,3_pInt; do o=1_pInt,3_pInt do p=1_pInt,3_pInt; do o=1_pInt,3_pInt
dFp_invdFdot(o,p,1:3,1:3) = Fpinv_rate(o,p)*FDot_inv dFp_invdFdot(o,p,1:3,1:3) = Fpinv_rate(o,p)*FDot_inv
enddo; enddo enddo; enddo

2
code/math.f90
View File

@ -239,7 +239,7 @@ subroutine math_init
enddo enddo
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'size of random seed: ', randSize write(6,*) 'size of random seed: ', randSize
do i =1, randSize do i =1, randSize
write(6,*) 'value of random seed: ', i, randInit(i) write(6,*) 'value of random seed: ', i, randInit(i)
enddo enddo