1) added distribution of leapfrog breaks

2) lattice_symmetryType is now a function (former lookup array was buggy)
3) stricter check of state var values (>0!) and memory deallocation done
This commit is contained in:
Philip Eisenlohr 2010-09-30 07:31:53 +00:00
parent e562df35a9
commit b259fbd9c6
4 changed files with 73 additions and 38 deletions

View File

@ -95,7 +95,7 @@ subroutine crystallite_init(Temperature)
mesh_maxNipNeighbors
use IO
use material
use lattice, only: lattice_symmetryTypes
use lattice, only: lattice_symmetryType
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, &
constitutive_phenopowerlaw_structure
use constitutive_titanmod, only: constitutive_titanmod_label, &
@ -302,7 +302,7 @@ subroutine crystallite_init(Temperature)
myStructure = -1_pInt ! does this happen for j2 material?
end select
if (myStructure>0_pInt) then
crystallite_symmetryID(g,i,e)=lattice_symmetryTypes(myStructure) ! structure = 1(fcc) or 2(bcc) => 1; 3(hex)=>2
crystallite_symmetryID(g,i,e) = lattice_symmetryType(myStructure) ! structure = 1(fcc) or 2(bcc) => 1; 3+(hex)=>2
endif
enddo
enddo
@ -1257,7 +1257,8 @@ endsubroutine
verboseDebugger, &
debug_cumLpCalls, &
debug_cumLpTicks, &
debug_StressLoopDistribution
debug_StressLoopDistribution, &
debug_LeapfrogBreakDistribution
use constitutive, only: constitutive_homogenizedC, &
constitutive_LpAndItsTangent
use math, only: math_mul33x33, &
@ -1421,9 +1422,11 @@ LpLoop: do
! NaN occured at regular speed?
if (any(residuum/=residuum) .and. leapfrog == 1.0) then
if (verboseDebugger) then
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress encountered NaN at ',g,i,e,' ; iteration ', NiterationStress
write(6,'(a,i3,x,i2,x,i5,x,a,i3,x,a)') '::: integrateStress encountered NaN at ',g,i,e,&
'; iteration ', NiterationStress, &
'>> returning..!'
!$OMPEND CRITICAL (write2out)
endif
return
@ -1437,8 +1440,8 @@ LpLoop: do
) then
if (verboseDebugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress encountered high-speed crash at ',g,i,e,' ; iteration ', &
NiterationStress
write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress encountered high-speed crash at ',g,i,e,&
'; iteration ', NiterationStress
!$OMPEND CRITICAL (write2out)
endif
maxleap = 0.5_pReal * leapfrog ! limit next acceleration
@ -1449,6 +1452,8 @@ LpLoop: do
Lpguess = Lpguess_old
residuum = residuum_old
debug_LeapfrogBreakDistribution(NiterationStress,mode) = debug_LeapfrogBreakDistribution(NiterationStress,mode) + 1
! residuum got better
else
! calculate Jacobian for correction term
@ -1465,8 +1470,8 @@ LpLoop: do
if (error) then
if (verboseDebugger .and. selectiveDebugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on dR/dLp inversion at ',g,i,e, &
' ; iteration ', NiterationStress
write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress failed on dR/dLp inversion at ',g,i,e, &
'; iteration ', NiterationStress
write(6,*)
write(6,'(a,/,9(9(e15.3,x)/))') 'dRdLp',dRdLp
write(6,'(a,/,9(9(e15.3,x)/))') 'dLpdT_constitutive',dLpdT_constitutive
@ -1476,10 +1481,10 @@ LpLoop: do
endif
return
else
if (.false. .and. verboseDebugger .and. selectiveDebugger) then
if (verboseDebugger .and. selectiveDebugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress did dR/dLp inversion at ',g,i,e, &
' ; iteration ', NiterationStress
write(6,'(a,i3,x,i2,x,i5,x,a,i3)') '::: integrateStress did dR/dLp inversion at ',g,i,e, &
'; iteration ', NiterationStress
write(6,*)
write(6,'(a,/,9(9(e15.3,x)/))') 'dRdLp',dRdLp
write(6,'(a,/,9(9(e15.3,x)/))') 'dLpdT_constitutive',dLpdT_constitutive
@ -1499,7 +1504,6 @@ LpLoop: do
! leapfrog to updated Lp
do k=1,3; do l=1,3; do m=1,3; do n=1,3
! forall (k=1:3,l=1:3,m=1:3,n=1:3) &
Lpguess(k,l) = Lpguess(k,l) - leapfrog*invdRdLp(3*(k-1)+l,3*(m-1)+n)*residuum(m,n)
enddo; enddo; enddo; enddo
enddo LpLoop

View File

@ -8,6 +8,7 @@
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_CrystalliteStateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_StiffnessStateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution
@ -22,8 +23,8 @@
integer(pInt) :: debug_e = 1_pInt
integer(pInt) :: debug_i = 1_pInt
integer(pInt) :: debug_g = 1_pInt
logical :: selectiveDebugger = .false.
logical :: verboseDebugger = .true.
logical :: selectiveDebugger = .true.
logical :: verboseDebugger = .false.
logical :: debugger = .true.
logical :: distribution_init = .false.
@ -68,6 +69,7 @@ subroutine debug_init()
write(6,*)
allocate(debug_StressLoopDistribution(nStress,2)) ; debug_StressLoopDistribution = 0_pInt
allocate(debug_LeapfrogBreakDistribution(nStress,2)) ; debug_LeapfrogBreakDistribution = 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
@ -120,6 +122,10 @@ subroutine debug_init()
write(6,'(a24,x,i8)') ' element: ',debug_e
write(6,'(a24,x,i8)') ' ip: ',debug_i
write(6,'(a24,x,i8)') ' grain: ',debug_g
else
debug_e = 0_pInt ! switch off selective debugging
debug_i = 0_pInt
debug_g = 0_pInt
endif
@ -134,6 +140,7 @@ subroutine debug_reset()
implicit none
debug_StressLoopDistribution = 0_pInt ! initialize debugging data
debug_LeapfrogBreakDistribution = 0_pInt
debug_CrystalliteStateLoopDistribution = 0_pInt
debug_StiffnessStateLoopDistribution = 0_pInt
debug_CrystalliteLoopDistribution = 0_pInt
@ -192,16 +199,18 @@ endsubroutine
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_StressLoop :'
write(6,*) 'distribution_StressLoop : stress frogbreak stiffness frogbreak'
do i=1,nStress
if (debug_StressLoopDistribution(i,1) /= 0 .or. debug_StressLoopDistribution(i,2) /= 0) then
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)') i,debug_StressLoopDistribution(i,1),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,x,i10)') ' total',integral,&
sum(debug_StressLoopDistribution(:,1)), &
sum(debug_StressLoopDistribution(:,2))
write(6,'(a15,i10,x,i10,12x,i10)') ' total',integral,&
sum(debug_StressLoopDistribution(:,1)), &
sum(debug_StressLoopDistribution(:,2))
integral = 0_pInt
write(6,*)

View File

@ -26,9 +26,6 @@ integer(pInt), parameter :: lattice_maxNslip = 48 ! max # of sli
integer(pInt), parameter :: lattice_maxNtwin = 24 ! max # of twin systems over lattice structures
integer(pInt), parameter :: lattice_maxNinteraction = 20 ! max # of interaction types (in hardening matrix part)
integer(pInt), parameter, dimension(3) :: lattice_symmetryTypes =(/1, 1, 2/) ! maps crystal structures to symmetry tpyes
integer(pInt), pointer, dimension(:,:) :: interactionSlipSlip, &
interactionSlipTwin, &
interactionTwinSlip, &
@ -661,6 +658,30 @@ CONTAINS
!* - lattice_initializeStructure
!****************************************
pure function lattice_symmetryType(structID)
!**************************************
!* maps structure to symmetry type *
!* fcc(1) and bcc(2) are cubic(1) *
!* hex(3+) is hexagonal(2) *
!**************************************
implicit none
integer(pInt), intent(in) :: structID
integer(pInt) lattice_symmetryType
select case(structID)
case (1,2)
lattice_symmetryType = 1_pInt
case (3:)
lattice_symmetryType = 2_pInt
case default
lattice_symmetryType = 0_pInt
end select
return
end function
subroutine lattice_init()
!**************************************

View File

@ -2522,10 +2522,10 @@ subroutine mesh_marc_count_cpSizes (unit)
pos = IO_stringPos(line,maxNchunks) ! limit to 64 nodes max (plus ID, type)
e = mesh_FEasCP('elem',IO_intValue(line,pos,1))
if (e /= 0) then ! disregard non CP elems
mesh_element (1,e) = IO_IntValue (line,pos,1) ! FE id
mesh_element (2,e) = FE_mapElemtype(IO_StringValue(line,pos,2)) ! elem type
mesh_element(1,e) = IO_IntValue (line,pos,1) ! FE id
mesh_element(2,e) = FE_mapElemtype(IO_StringValue(line,pos,2)) ! elem type
forall (j = 1:FE_Nnodes(mesh_element(2,e))) &
mesh_element(j+4,e) = IO_IntValue(line,pos,j+2) ! copy FE ids of nodes
mesh_element(j+4,e) = IO_IntValue(line,pos,j+2) ! copy FE ids of nodes
call IO_skipChunks(unit,FE_NoriginalNodes(mesh_element(2,e))-(pos(1)-2)) ! read on if FE_Nnodes exceeds node count present on current line
endif
enddo
@ -2978,18 +2978,19 @@ subroutine mesh_marc_count_cpSizes (unit)
integer(pInt) i,e,n,f,t
if (mesh_maxValStateVar(1) == 0) call IO_error(110) ! no materials specified
if (mesh_maxValStateVar(2) == 0) call IO_error(120) ! no textures specified
if (mesh_maxValStateVar(1) < 1_pInt) call IO_error(110) ! no homogenization specified
if (mesh_maxValStateVar(2) < 1_pInt) call IO_error(120) ! no microstructure specified
allocate (mesh_HomogMicro(mesh_maxValStateVar(1),mesh_maxValStateVar(2))); mesh_HomogMicro = 0_pInt
do i=1,mesh_NcpElems
mesh_HomogMicro(mesh_element(3,i),mesh_element(4,i)) = &
mesh_HomogMicro(mesh_element(3,i),mesh_element(4,i)) + 1 ! count combinations of homogenization and microstructure
do e = 1,mesh_NcpElems
if (mesh_element(3,e) < 1_pInt) call IO_error(110,e) ! no homogenization specified
if (mesh_element(4,e) < 1_pInt) call IO_error(120,e) ! no homogenization specified
mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) = &
mesh_HomogMicro(mesh_element(3,e),mesh_element(4,e)) + 1 ! count combinations of homogenization and microstructure
enddo
if (verboseDebugger) then
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) 'Input Parser: IP COORDINATES'
write(6,'(a5,x,a5,3(x,a12))') 'elem','IP','x','y','z'
@ -3041,12 +3042,11 @@ subroutine mesh_marc_count_cpSizes (unit)
enddo
enddo
enddo
!$OMP END CRITICAL (write2out)
endif
!$OMP END CRITICAL (write2out)
endif
!$OMP CRITICAL (write2out)
write (6,*)
write (6,*) "Input Parser: STATISTICS"
write (6,*)
@ -3072,8 +3072,9 @@ endif
enddo
write (6,*)
call flush(6)
!$OMP END CRITICAL (write2out)
!$OMP END CRITICAL (write2out)
deallocate(mesh_HomogMicro)
return
endsubroutine