diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index b13e5c6b9..4d8b1e829 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -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. diff --git a/code/debug.f90 b/code/debug.f90 index e7efbdf98..71e1a392c 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -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 diff --git a/code/mpie_cpfem_marc.f90 b/code/mpie_cpfem_marc.f90 index 27257e9cf..f376ec28e 100644 --- a/code/mpie_cpfem_marc.f90 +++ b/code/mpie_cpfem_marc.f90 @@ -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