extreme values of stress and jacobian now recorded in CPFEM_general. variable declaration and generation of output moved to debug module.

This commit is contained in:
Christoph Kords 2011-03-17 13:13:13 +00:00
parent 235266b169
commit 9d7ede7e03
3 changed files with 204 additions and 186 deletions

View File

@ -218,8 +218,15 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
debug_i, &
debug_e, &
debugger, &
selectiveDebugger, &
verboseDebugger
verboseDebugger, &
debug_stressMaxLocation, &
debug_stressMinLocation, &
debug_jacobianMaxLocation, &
debug_jacobianMinLocation, &
debug_stressMax, &
debug_stressMin, &
debug_jacobianMax, &
debug_jacobianMin
use FEsolving, only: parallelExecution, &
outdatedFFN1, &
terminallyIll, &
@ -236,7 +243,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
math_transpose3x3, &
math_I3, &
math_Mandel3333to66, &
math_Mandel33to6
math_Mandel66to3333, &
math_Mandel33to6, &
math_Mandel6to33
use mesh, only: mesh_FEasCP, &
mesh_NcpElems, &
mesh_maxNips, &
@ -291,15 +300,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
!*** output variables ***!
real(pReal), dimension(6), intent(out) :: cauchyStress ! stress vector in Mandel notation
real(pReal), dimension(6,6), intent(out) :: jacobian ! jacobian in Mandel notation
real(pReal), dimension (3,3), intent(out) :: pstress ! Piola-Kirchhoff stress in Matrix notation
real(pReal), dimension (3,3), intent(out) :: pstress ! Piola-Kirchhoff stress in Matrix notation
real(pReal), dimension (3,3,3,3), intent(out) :: dPdF !
!*** local variables ***!
real(pReal) J_inverse, & ! inverse of Jacobian
rnd
real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress in Matrix notation
real(pReal), dimension (3,3) :: Kirchhoff, & ! Piola-Kirchhoff stress in Matrix notation
cauchyStress33 ! stress vector in Matrix notation
real(pReal), dimension (3,3,3,3) :: H_sym, &
H
H, &
jacobian3333 ! jacobian in Matrix notation
integer(pInt) cp_en, & ! crystal plasticity element number
i, &
j, &
@ -553,6 +564,28 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
!$OMP END CRITICAL (write2out)
endif
! remember extreme values of stress and jacobian
if (mode < 3) then
cauchyStress33 = math_Mandel6to33(cauchyStress)
if (maxval(cauchyStress33) > debug_stressMax) then
debug_stressMaxLocation = (/cp_en, IP/)
debug_stressMax = maxval(cauchyStress33)
endif
if (minval(cauchyStress33) < debug_stressMin) then
debug_stressMinLocation = (/cp_en, IP/)
debug_stressMin = minval(cauchyStress33)
endif
jacobian3333 = math_Mandel66to3333(jacobian)
if (maxval(jacobian3333) > debug_jacobianMax) then
debug_jacobianMaxLocation = (/cp_en, IP/)
debug_jacobianMax = maxval(jacobian3333)
endif
if (minval(jacobian3333) < debug_jacobianMin) then
debug_jacobianMinLocation = (/cp_en, IP/)
debug_jacobianMin = minval(jacobian3333)
endif
endif
! return temperature
if (theTime > 0.0_pReal) Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition.

View File

@ -1,33 +1,41 @@
!* $Id$
!##############################################################
MODULE debug
MODULE debug
!##############################################################
use prec
use prec
implicit none
character(len=64), parameter :: debug_configFile = 'debug.config' ! name of configuration file
implicit none
character(len=64), parameter :: debug_configFile = 'debug.config' ! name of configuration file
integer(pInt), dimension(:,:), allocatable :: debug_StressLoopDistribution
integer(pInt), dimension(:,:), allocatable :: debug_LeapfrogBreakDistribution
integer(pInt), dimension(:,:), allocatable :: debug_StateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_MaterialpointStateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_MaterialpointLoopDistribution
integer(pLongInt) :: debug_cumLpTicks = 0_pInt
integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt
integer(pLongInt) :: debug_cumDotTemperatureTicks = 0_pInt
integer(pInt) :: debug_cumLpCalls = 0_pInt
integer(pInt) :: debug_cumDotStateCalls = 0_pInt
integer(pInt) :: debug_cumDotTemperatureCalls = 0_pInt
integer(pInt) :: debug_e = 1_pInt
integer(pInt) :: debug_i = 1_pInt
integer(pInt) :: debug_g = 1_pInt
logical :: selectiveDebugger = .true.
logical :: verboseDebugger = .false.
logical :: debugger = .true.
logical :: distribution_init = .false.
integer(pInt), dimension(:,:), allocatable :: debug_StressLoopDistribution
integer(pInt), dimension(:,:), allocatable :: debug_LeapfrogBreakDistribution
integer(pInt), dimension(:,:), allocatable :: debug_StateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_MaterialpointStateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_MaterialpointLoopDistribution
integer(pLongInt) :: debug_cumLpTicks = 0_pInt
integer(pLongInt) :: debug_cumDotStateTicks = 0_pInt
integer(pLongInt) :: debug_cumDotTemperatureTicks = 0_pInt
integer(pInt) :: debug_cumLpCalls = 0_pInt
integer(pInt) :: debug_cumDotStateCalls = 0_pInt
integer(pInt) :: debug_cumDotTemperatureCalls = 0_pInt
integer(pInt) :: debug_e = 1_pInt
integer(pInt) :: debug_i = 1_pInt
integer(pInt) :: debug_g = 1_pInt
integer(pInt), dimension(2) :: debug_stressMaxLocation = 0_pInt
integer(pInt), dimension(2) :: debug_stressMinLocation = 0_pInt
integer(pInt), dimension(2) :: debug_jacobianMaxLocation = 0_pInt
integer(pInt), dimension(2) :: debug_jacobianMinLocation = 0_pInt
real(pReal) :: debug_stressMax
real(pReal) :: debug_stressMin
real(pReal) :: debug_jacobianMax
real(pReal) :: debug_jacobianMin
logical :: selectiveDebugger = .true.
logical :: verboseDebugger = .false.
logical :: debugger = .true.
logical :: distribution_init = .false.
CONTAINS
CONTAINS
!********************************************************************
@ -159,125 +167,146 @@ subroutine debug_reset()
debug_cumLpCalls = 0_pInt
debug_cumDotStateCalls = 0_pInt
debug_cumDotTemperatureCalls = 0_pInt
debug_stressMaxLocation = 0_pInt
debug_stressMinLocation = 0_pInt
debug_jacobianMaxLocation = 0_pInt
debug_jacobianMinLocation = 0_pInt
debug_stressMax = - 1.0_pReal / 0.0_pReal
debug_stressMin = + 1.0_pReal / 0.0_pReal
debug_jacobianMax = - 1.0_pReal / 0.0_pReal
debug_jacobianMin = + 1.0_pReal / 0.0_pReal
endsubroutine
!********************************************************************
! write debug statements to standard out
!********************************************************************
subroutine debug_info()
subroutine debug_info()
use prec
use numerics, only: nStress, &
nState, &
nCryst, &
nMPstate, &
nHomog
implicit none
use prec
use numerics, only: nStress, &
nState, &
nCryst, &
nMPstate, &
nHomog
implicit none
integer(pInt) i,integral
integer(pLongInt) tickrate
integer(pInt) i,integral
integer(pLongInt) tickrate
call system_clock(count_rate=tickrate)
call system_clock(count_rate=tickrate)
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) 'DEBUG Info'
write(6,*)
write(6,'(a33,x,i12)') 'total calls to LpAndItsTangent :',debug_cumLpCalls
if (debug_cumLpCalls > 0_pInt) then
write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumLpTicks)/tickrate
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumLpTicks)*1.0e6_pReal/tickrate/debug_cumLpCalls
endif
write(6,*)
write(6,'(a33,x,i12)') 'total calls to collectDotState :',debug_cumDotStateCalls
if (debug_cumdotStateCalls > 0_pInt) then
write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumDotStateTicks)/tickrate
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumDotStateTicks)*1.0e6_pReal/tickrate/debug_cumDotStateCalls
endif
write(6,*)
write(6,'(a33,x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls
if (debug_cumdotTemperatureCalls > 0_pInt) then
write(6,'(a33,x,f12.3)') 'total CPU time/s :', dble(debug_cumDotTemperatureTicks)/tickrate
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumDotTemperatureTicks)*1.0e6_pReal/tickrate/debug_cumDotTemperatureCalls
endif
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) 'DEBUG Info (from previous cycle)'
write(6,*)
write(6,'(a33,x,i12)') 'total calls to LpAndItsTangent :',debug_cumLpCalls
if (debug_cumLpCalls > 0_pInt) then
write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumLpTicks)/tickrate
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumLpTicks)*1.0e6_pReal/tickrate/debug_cumLpCalls
endif
write(6,*)
write(6,'(a33,x,i12)') 'total calls to collectDotState :',debug_cumDotStateCalls
if (debug_cumdotStateCalls > 0_pInt) then
write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumDotStateTicks)/tickrate
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumDotStateTicks)*1.0e6_pReal/tickrate/debug_cumDotStateCalls
endif
write(6,*)
write(6,'(a33,x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls
if (debug_cumdotTemperatureCalls > 0_pInt) then
write(6,'(a33,x,f12.3)') 'total CPU time/s :', dble(debug_cumDotTemperatureTicks)/tickrate
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumDotTemperatureTicks)*1.0e6_pReal/tickrate/debug_cumDotTemperatureCalls
endif
integral = 0_pInt
write(6,*)
write(6,*)
write(6,*) 'distribution_StressLoop : stress frogbreak stiffness frogbreak'
do i=1,nStress
if (any(debug_StressLoopDistribution(i,:) /= 0_pInt ) .or. &
any(debug_LeapfrogBreakDistribution(i,:) /= 0_pInt ) ) then
integral = integral + i*debug_StressLoopDistribution(i,1) + i*debug_StressLoopDistribution(i,2)
write(6,'(i25,x,i10,x,i10,x,i10,x,i10)') i,debug_StressLoopDistribution(i,1),debug_LeapfrogBreakDistribution(i,1), &
debug_StressLoopDistribution(i,2),debug_LeapfrogBreakDistribution(i,2)
endif
enddo
write(6,'(a15,i10,x,i10,12x,i10)') ' total',integral,&
sum(debug_StressLoopDistribution(:,1)), &
sum(debug_StressLoopDistribution(:,2))
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteStateLoop :'
do i=1,nState
if (any(debug_StateLoopDistribution(i,:) /= 0)) then
integral = integral + i*debug_StateLoopDistribution(i,1) + i*debug_StateLoopDistribution(i,2)
write(6,'(i25,x,i10,12x,i10)') i,debug_StateLoopDistribution(i,1),debug_StateLoopDistribution(i,2)
endif
enddo
write(6,'(a15,i10,x,i10,12x,i10)') ' total',integral,&
sum(debug_StateLoopDistribution(:,1)), &
sum(debug_StateLoopDistribution(:,2))
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteCutbackLoop :'
do i=1,nCryst+1
if (debug_CrystalliteLoopDistribution(i) /= 0) then
integral = integral + i*debug_CrystalliteLoopDistribution(i)
if (i <= nCryst) then
write(6,'(i25,x,i10)') i,debug_CrystalliteLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1,'+',debug_CrystalliteLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_MaterialpointStateLoop :'
do i=1,nMPstate
if (debug_MaterialpointStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointStateLoopDistribution(i)
write(6,'(i25,x,i10)') i,debug_MaterialpointStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_MaterialpointCutbackLoop :'
do i=1,nHomog+1
if (debug_MaterialpointLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointLoopDistribution(i)
if (i <= nHomog) then
write(6,'(i25,x,i10)') i,debug_MaterialpointLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1,'+',debug_MaterialpointLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
write(6,*)
write(6,*)
write(6,*) 'Extreme values of returned stress and jacobian'
write(6,*)
write(6,'(a39)') ' value el ip'
write(6,'(a14,x,e12.3,x,i6,x,i4)') 'stress min :', debug_stressMin, debug_stressMinLocation
write(6,'(a14,x,e12.3,x,i6,x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation
write(6,'(a14,x,e12.3,x,i6,x,i4)') 'jacobian min :', debug_jacobianMin, debug_jacobianMinLocation
write(6,'(a14,x,e12.3,x,i6,x,i4)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_StressLoop : stress frogbreak stiffness frogbreak'
do i=1,nStress
if (any(debug_StressLoopDistribution(i,:) /= 0_pInt ) .or. &
any(debug_LeapfrogBreakDistribution(i,:) /= 0_pInt ) ) then
integral = integral + i*debug_StressLoopDistribution(i,1) + i*debug_StressLoopDistribution(i,2)
write(6,'(i25,x,i10,x,i10,x,i10,x,i10)') i,debug_StressLoopDistribution(i,1),debug_LeapfrogBreakDistribution(i,1), &
debug_StressLoopDistribution(i,2),debug_LeapfrogBreakDistribution(i,2)
endif
enddo
write(6,'(a15,i10,x,i10,12x,i10)') ' total',integral,&
sum(debug_StressLoopDistribution(:,1)), &
sum(debug_StressLoopDistribution(:,2))
write(6,*)
!$OMP END CRITICAL (write2out)
endsubroutine
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteStateLoop :'
do i=1,nState
if (any(debug_StateLoopDistribution(i,:) /= 0)) then
integral = integral + i*debug_StateLoopDistribution(i,1) + i*debug_StateLoopDistribution(i,2)
write(6,'(i25,x,i10,12x,i10)') i,debug_StateLoopDistribution(i,1),debug_StateLoopDistribution(i,2)
endif
enddo
write(6,'(a15,i10,x,i10,12x,i10)') ' total',integral,&
sum(debug_StateLoopDistribution(:,1)), &
sum(debug_StateLoopDistribution(:,2))
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_CrystalliteCutbackLoop :'
do i=1,nCryst+1
if (debug_CrystalliteLoopDistribution(i) /= 0) then
integral = integral + i*debug_CrystalliteLoopDistribution(i)
if (i <= nCryst) then
write(6,'(i25,x,i10)') i,debug_CrystalliteLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1,'+',debug_CrystalliteLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
integral = 0_pInt
write(6,*)
write(6,*)
write(6,*) 'distribution_MaterialpointStateLoop :'
do i=1,nMPstate
if (debug_MaterialpointStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointStateLoopDistribution(i)
write(6,'(i25,x,i10)') i,debug_MaterialpointStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_MaterialpointCutbackLoop :'
do i=1,nHomog+1
if (debug_MaterialpointLoopDistribution(i) /= 0) then
integral = integral + i*debug_MaterialpointLoopDistribution(i)
if (i <= nHomog) then
write(6,'(i25,x,i10)') i,debug_MaterialpointLoopDistribution(i)
else
write(6,'(i25,a1,i10)') i-1,'+',debug_MaterialpointLoopDistribution(i)
endif
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
write(6,*)
!$OMP END CRITICAL (write2out)
endsubroutine
END MODULE debug
END MODULE debug

View File

@ -240,11 +240,7 @@ subroutine hypela2(&
real(pReal), dimension (3,3) :: pstress ! dummy argument for call of cpfem_general (used by mpie_spectral)
real(pReal), dimension (3,3,3,3) :: dPdF ! dummy argument for call of cpfem_general (used by mpie_spectral)
integer(pInt) computationMode, i, cp_en, &
s_max_e, s_max_i, s_min_e, s_min_i, &
d_max_e, d_max_i, d_min_e, d_min_i
real(pReal) s_max, s_min, &
d_max, d_min
integer(pInt) computationMode, i, cp_en
! OpenMP variable
!$ integer(pInt) defaultNumThreadsInt ! default value set by Marc
@ -253,7 +249,6 @@ subroutine hypela2(&
!$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc
if (.not. CPFEM_init_done) call CPFEM_initAll(t(1),n(1),nn)
cp_en = mesh_FEasCP('elem',n(1))
!$ call omp_set_num_threads(mpieNumThreadsInt) ! set number of threads for parallel execution set by MPIE_NUM_THREADS
@ -264,6 +259,7 @@ subroutine hypela2(&
computationMode = 6 ! --> just return known tangent
endif
else ! stress requested (lovl == 6)
cp_en = mesh_FEasCP('elem',n(1))
if (cptim > theTime .or. inc /= theInc) then ! reached "convergence"
terminallyIll = .false.
cycleCounter = -1 ! first calc step increments this to cycle = 0
@ -308,10 +304,6 @@ subroutine hypela2(&
call debug_reset() ! resets debugging
outdatedFFN1 = .false.
cycleCounter = cycleCounter + 1_pInt
s_max = - 1.0_pReal / 0.0_pReal ! reset stored max/min values
s_min = + 1.0_pReal / 0.0_pReal
d_max = - 1.0_pReal / 0.0_pReal
d_min = + 1.0_pReal / 0.0_pReal
endif
if ( outdatedByNewInc ) then
outdatedByNewInc = .false. ! reset flag
@ -323,17 +315,6 @@ subroutine hypela2(&
if ( lastMode /= calcMode(nn,cp_en) .and. &
.not. terminallyIll ) then
call debug_info() ! first after ping pong reports (meaningful) debugging
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) 'EXTREME VALUES OF RETURNED VARIABLES (from previous cycle)'
write(6,*)
write(6,'(a39)') ' value el ip'
write(6,'(a14,x,e12.3,x,i6,x,i4)') 'stress min :', s_min, s_min_e, s_min_i
write(6,'(a14,x,e12.3,x,i6,x,i4)') ' max :', s_max, s_max_e, s_max_i
write(6,'(a14,x,e12.3,x,i6,x,i4)') 'jacobian min :', d_min, d_min_e, d_min_i
write(6,'(a14,x,e12.3,x,i6,x,i4)') ' max :', d_max, d_max_e, d_max_i
write(6,*)
!$OMP END CRITICAL (write2out)
endif
if ( lastIncConverged ) then
lastIncConverged = .false. ! reset flag
@ -359,33 +340,8 @@ subroutine hypela2(&
s(1:ngens) = stress(1:ngens)*invnrmMandel(1:ngens)
if(symmetricSolver) d(1:ngens,1:ngens) = 0.5_pReal*(d(1:ngens,1:ngens)+transpose(d(1:ngens,1:ngens)))
if (calcMode(nn,cp_en)) then
if (maxval(s(1:ngens)) > s_max) then ! remember extreme values of stress and jacobian
s_max_e = cp_en
s_max_i = nn
s_max = maxval(s(1:ngens))
endif
if (minval(s(1:ngens)) < s_min) then
s_min_e = cp_en
s_min_i = nn
s_min = minval(s(1:ngens))
endif
if (maxval(d(1:ngens,1:ngens)) > d_max) then
d_max_e = cp_en
d_max_i = nn
d_max = maxval(d(1:ngens,1:ngens))
endif
if (minval(d(1:ngens,1:ngens)) < d_min) then
d_min_e = cp_en
d_min_i = nn
d_min = minval(d(1:ngens,1:ngens))
endif
endif
!$ call omp_set_num_threads(defaultNumThreadsInt) ! reset number of threads to stored default value
return
end subroutine