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

1) fixed terminallyIll bug occurring with Ngrains=1 homogenization 

2) brushed up the output to be more easily readable/understandable
This commit is contained in:
Philip Eisenlohr 2010-09-01 21:04:02 +00:00
parent 4f29e8c2fe
commit c19524f264
4 changed files with 130 additions and 115 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
code/CPFEM.f90
View File

@ -217,13 +217,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
cp_en = mesh_FEasCP('elem',element)
selectiveDebugger = (cp_en == debug_e .and. IP == debug_i)
if (debugger .and. selectiveDebugger) then
if (selectiveDebugger) then
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '#####################################'
write(6,'(a10,x,f8.4,x,a10,x,f8.4,x,a10,x,i6,x,a10,x,i3,x,a16,x,i2,x,a16,x,i2)') &
'theTime',theTime,'theDelta',theDelta,'theInc',theInc,'cycleCounter',cycleCounter,'computationMode',mode
write(6,*) '#####################################'
write(6,'(a)') '#######################################################'
write(6,'(a32,x,i5,x,i2)') 'reporting for element, ip:',cp_en,IP
write(6,'(a32,x,f15.7)') 'theTime',theTime
write(6,'(a32,x,f15.7)') 'theDelta',theDelta
write(6,'(a32,x,i8)') 'theInc',theInc
write(6,'(a32,x,i8)') 'cycleCounter',cycleCounter
write(6,'(a32,x,i8)') 'computationMode',mode
write(6,'(a)') '#######################################################'
call flush (6)
!$OMP END CRITICAL (write2out)
endif
@ -243,9 +247,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
j = 1:mesh_maxNips, &
k = 1:mesh_NcpElems ) &
constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites
if (debugger .and. selectiveDebugger) then
if (selectiveDebugger) then
!$OMP CRITICAL (write2out)
write(6,'(a16,i2,x,i5,/,4(3(e20.8,x),/))') 'aged state at 1 ',IP,cp_en, constitutive_state(1,IP,cp_en)%p
write(6,'(a32,x,i8,x,i2,/,4(3(e20.8,x),/))') '°°° AGED state of grain 1, element ip',cp_en,IP, &
constitutive_state(1,IP,cp_en)%p
!$OMP END CRITICAL (write2out)
endif
do k = 1,mesh_NcpElems
@ -265,14 +270,15 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
! deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one
if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,IP,cp_en)) > defgradTolerance)) then
if (.not. terminallyIll .and. .not. outdatedFFN1) then
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a11,x,i2,x,i5,x,a12,/,3(3(f10.6,x),/))') 'outdated at',IP,cp_en,'; FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3)
write(6,'(a32,x,i5,x,i2)') '°°° OUTDATED at element ip',cp_en,IP
write(6,'(a32,/,3(3(f10.6,x),/))') ' FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3)
!$OMP END CRITICAL (write2out)
endif
outdatedFFN1 = .true.
endif
CPFEM_cs(:,IP,cp_en) = CPFEM_odd_stress
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = 1.0_pReal - rnd
CPFEM_cs(:,IP,cp_en) = rnd*CPFEM_odd_stress
CPFEM_dcsde(:,:,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6)
! deformation gradient is not outdated
@ -307,8 +313,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
CPFEM_calc_done = .true.
endif
if (terminallyIll) then
CPFEM_cs(:,IP,cp_en) = CPFEM_odd_stress
if ( terminallyIll ) then
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = 1.0_pReal - rnd
CPFEM_cs(:,IP,cp_en) = rnd*CPFEM_odd_stress
CPFEM_dcsde(:,:,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6)
else
! translate from P to CS
@ -334,7 +342,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
endif
endif
! --+>> COLLECTION OF FEM INPUT WITH RETURNING OF ODD STRESS AND JACOBIAN <<+--
! --+>> COLLECTION OF FEM INPUT WITH RETURNING OF RANDOMIZED ODD STRESS AND JACOBIAN <<+--
case (3,4,5)
if (mode == 4) then
CPFEM_dcsde_knownGood = CPFEM_dcsde ! --+>> BACKUP JACOBIAN FROM FORMER CONVERGED INC
@ -353,8 +361,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
! --+>> RECYCLING OF FORMER RESULTS (MARC SPECIALTY) <<+--
case (6)
! do nothing
! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC
case (7)
CPFEM_dcsde = CPFEM_dcsde_knownGood ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC
CPFEM_dcsde = CPFEM_dcsde_knownGood
end select
@ -365,10 +374,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
! copy P and dPdF to the output variables
pstress(:,:) = materialpoint_P(:,:,IP,cp_en)
dPdF(:,:,:,:) = materialpoint_dPdF(:,:,:,:,IP,cp_en)
if (debugger .and. selectiveDebugger) then
if ((debugger .and. selectiveDebugger) .and. &
mode < 6) then
!$OMP CRITICAL (write2out)
write(6,'(a16,x,i2,x,a2,x,i4,/,6(f10.3,x)/)') 'stress/MPa at ip', IP, 'el', cp_en, cauchyStress/1e6
write(6,'(a18,x,i2,x,a2,x,i4,/,6(6(f10.3,x)/))') 'jacobian/GPa at ip', IP, 'el', cp_en, jacobian/1e9
write(6,'(a,x,i2,x,a,x,i4,/,6(f10.3,x)/)') 'stress/MPa at ip', IP, 'el', cp_en, cauchyStress/1e6
write(6,'(a,x,i2,x,a,x,i4,/,6(6(f10.3,x)/))') 'jacobian/GPa at ip', IP, 'el', cp_en, jacobian/1e9
call flush(6)
!$OMP END CRITICAL (write2out)
endif

31
code/debug.f90
View File

@ -47,7 +47,7 @@ subroutine debug_init()
write(6,*)
allocate(debug_StressLoopDistribution(nStress)) ; debug_StressLoopDistribution = 0_pInt
allocate(debug_CrystalliteStateLoopDistribution(nState)); debug_CrystalliteStateLoopDistribution = 0_pInt
allocate(debug_CrystalliteStateLoopDistribution(nState)) ; debug_CrystalliteStateLoopDistribution = 0_pInt
allocate(debug_StiffnessStateLoopDistribution(nState)) ; debug_StiffnessStateLoopDistribution = 0_pInt
allocate(debug_CrystalliteLoopDistribution(nCryst+1)) ; debug_CrystalliteLoopDistribution = 0_pInt
allocate(debug_MaterialpointStateLoopDistribution(nMPstate)) ; debug_MaterialpointStateLoopDistribution = 0_pInt
@ -143,18 +143,7 @@ endsubroutine
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_StiffnessStateLoop :'
do i=1,nState
if (debug_StiffnessStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_StiffnessStateLoopDistribution(i)
write(6,'(i25,x,i10)') i,debug_StiffnessStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_StiffnessStateLoopDistribution)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteLoop :'
write(6,*) 'distribution_CrystalliteCutbackLoop :'
do i=1,nCryst+1
if (debug_CrystalliteLoopDistribution(i) /= 0) then
integral = integral + i*debug_CrystalliteLoopDistribution(i)
@ -166,11 +155,22 @@ endsubroutine
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_StiffnessStateLoop :'
do i=1,nState
if (debug_StiffnessStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_StiffnessStateLoopDistribution(i)
write(6,'(i25,x,i10)') i,debug_StiffnessStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_StiffnessStateLoopDistribution)
!* Material point loop counter <<<updated 31.07.2009>>>
integral = 0_pInt
write(6,*)
write(6,*)
write(6,*) 'distribution_MaterialpointStateLoop :'
do i=1,nMPstate
if (debug_MaterialpointStateLoopDistribution(i) /= 0) then
@ -179,11 +179,10 @@ endsubroutine
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution)
write(6,*)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_MaterialpointLoop :'
write(6,*) 'distribution_MaterialpointCutbackLoop :'
do i=1,nHomog+1
if (debug_MaterialpointLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointLoopDistribution(i)
@ -195,8 +194,8 @@ endsubroutine
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
write(6,*)
write(6,*)
endsubroutine

66
code/homogenization.f90
View File

@ -278,7 +278,7 @@ subroutine materialpoint_stressAndItsTangent(&
if (debugger) then
write (6,*)
write (6,*) 'Material Point start'
write (6,'(a,/,(f12.7,x))') 'Temp0 of 1 1' ,materialpoint_Temperature(1,1)
write (6,'(a,/,(f12.7,x))') 'Temp0 of 1 1',materialpoint_Temperature(1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'F0 of 1 1',materialpoint_F0(1:3,:,1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'F of 1 1',materialpoint_F(1:3,:,1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp0 of 1 1 1',crystallite_Fp0(1:3,:,1,1,1)
@ -316,9 +316,9 @@ subroutine materialpoint_stressAndItsTangent(&
! ------ cutback loop ------
do while (any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinHomog)) ! cutback loop for material points
do while (.not. terminallyIll .and. &
any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMinHomog)) ! cutback loop for material points
! write(6,'(a,/,125(8(f8.5,x),/))') 'mp_subSteps',materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2))
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
@ -326,13 +326,12 @@ subroutine materialpoint_stressAndItsTangent(&
selectiveDebugger = (e == debug_e .and. i == debug_i)
! if our materialpoint converged or consists of only one single grain then we are either finished or have to wind forward
if ( materialpoint_converged(i,e) .or. (myNgrains == 1_pInt .and. materialpoint_subStep(i,e) <= 1.0_pReal) ) then
if ( materialpoint_converged(i,e) ) then
if (verboseDebugger .and. selectiveDebugger) then
!$OMP CRITICAL (write2out)
write(6,'(a21,f10.8,a34,f10.8,a37,/)') 'winding forward from ', &
materialpoint_subFrac(i,e), ' to current materialpoint_subFrac ', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),' in materialpoint_stressAndItsTangent'
write(6,'(a,x,f10.8,x,a,x,f10.8,x,a,/)') '°°° winding forward from', &
materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent'
!$OMPEND CRITICAL (write2out)
endif
@ -363,26 +362,31 @@ subroutine materialpoint_stressAndItsTangent(&
! materialpoint didn't converge, so we need a cutback here
else
if ( (myNgrains == 1_pInt .and. materialpoint_subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
subStepSizeHomog * materialpoint_subStep(i,e) <= subStepMinHomog ) then ! would require too small subStep
! cutback makes no sense and...
terminallyIll = .true. ! ...one kills all
else ! cutback makes sense
materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
! <<modified to add more flexibility in cutback>>
if (verboseDebugger .and. selectiveDebugger) then
!$OMP CRITICAL (write2out)
write(6,'(a82,f10.8,/)') 'cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep: ',&
write(6,'(a,x,f10.8,/)') '°°° cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',&
materialpoint_subStep(i,e)
!$OMPEND CRITICAL (write2out)
endif
! restore...
crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures
! ...initial def grad unchanged
crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads
crystallite_Lp(:,:,1:myNgrains,i,e) = crystallite_partionedLp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_Tstar_v(:,1:myNgrains,i,e) = crystallite_partionedTstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress
forall (g = 1:myNgrains) constitutive_state(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructures
if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_state(i,e)%p = homogenization_subState0(i,e)%p ! ...internal state of homog scheme
endif
endif
materialpoint_requested(i,e) = materialpoint_subStep(i,e) > subStepMinHomog
@ -392,8 +396,8 @@ subroutine materialpoint_stressAndItsTangent(&
materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt
materialpoint_doneAndHappy(:,i,e) = (/.false.,.true./)
endif
enddo
enddo
enddo ! loop IPs
enddo ! loop elements
!$OMP END PARALLEL DO
!* Checks for cutback/substepping loops: added <<<updated 31.07.2009>>>
@ -406,9 +410,11 @@ subroutine materialpoint_stressAndItsTangent(&
NiterationMPstate = 0_pInt
do while (any( materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) &
do while (.not. terminallyIll .and. &
any( materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) &
.and. .not. materialpoint_doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
) .and. NiterationMPstate < nMPstate) ! convergence loop for materialpoint
) .and. &
NiterationMPstate < nMPstate) ! convergence loop for materialpoint
NiterationMPstate = NiterationMPstate + 1
! write(6,'(a,/,125(8(l,x),/))') 'material point request and not done', &
@ -476,33 +482,29 @@ subroutine materialpoint_stressAndItsTangent(&
enddo ! cutback loop
! calculate crystal orientations
call crystallite_orientations()
! check for non-performer: any(.not. converged)
! replace everybody with odd response ?
if (.not. terminallyIll ) then
!$OMP PARALLEL DO
elementLoop: do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
if (materialpoint_converged(i,e)) then
call homogenization_averageStressAndItsTangent(i,e)
call homogenization_averageTemperature(i,e)
else
terminallyIll = .true.
write(6,'(a48,i4,i4,/)') 'homogenization terminally-ill ',i,e
exit elementLoop
endif
enddo
enddo elementLoop
!$OMP END PARALLEL DO
enddo; enddo
!$OMP END PARALLEL DO
if (debugger) then
write (6,*)
write (6,*) 'Material Point end'
write (6,'(a)') '°°° Material Point end'
write (6,*)
endif
else
write (6,*)
write (6,'(a)') '°°° Material Point terminally ill'
write (6,*)
endif
return
endsubroutine

12
code/mpie_cpfem_marc.f90
View File

@ -51,6 +51,8 @@ subroutine mpie_interface_init()
write(6,*) '<<<+- mpie_cpfem_marc init -+>>>'
write(6,*) '$Id$'
write(6,*)
write(6,*)
write(6,*)
return
end subroutine
@ -247,7 +249,7 @@ subroutine hypela2(&
terminallyIll = .false.
cycleCounter = 0
!$OMP CRITICAL (write2out)
write (6,'(i6,x,i2,x,a)') n(1),nn,'<< hypela2 >> lastIncConverged + outdated'; call flush(6)
write (6,'(i6,x,i2,x,a)') n(1),nn,'<< hypela2 >> former increment converged..!'; call flush(6)
!$OMP END CRITICAL (write2out)
else if ( timinc < theDelta ) then ! cutBack
@ -261,10 +263,11 @@ subroutine hypela2(&
calcMode(nn,cp_en) = .not. calcMode(nn,cp_en) ! ping pong (calc <--> collect)
if ( calcMode(nn,cp_en) ) then ! now calc
if ( calcMode(nn,cp_en) ) then ! now --- CALC ---
if ( lastMode .neqv. calcMode(nn,cp_en) ) then ! first after ping pong
call debug_reset() ! resets debugging
outdatedFFN1 = .false.
terminallyIll = .false.
cycleCounter = cycleCounter + 1
endif
if ( outdatedByNewInc ) then
@ -273,8 +276,9 @@ subroutine hypela2(&
else
computationMode = 2 ! plain calc
endif
else ! now collect
if ( lastMode .neqv. calcMode(nn,cp_en) ) call debug_info() ! first after ping pong reports debugging
else ! now --- COLLECT ---
if ( lastMode .neqv. calcMode(nn,cp_en) .and. &
.not. terminallyIll ) call debug_info() ! first after ping pong reports (meaningful) debugging
if ( lastIncConverged ) then
lastIncConverged = .false.
computationMode = 4 ! collect and backup Jacobian after convergence