From b259fbd9c6d9144ace84b0020ad4a06b908b6fd3 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 30 Sep 2010 07:31:53 +0000 Subject: [PATCH] 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 --- code/crystallite.f90 | 32 ++++++++++++++++++-------------- code/debug.f90 | 25 +++++++++++++++++-------- code/lattice.f90 | 27 ++++++++++++++++++++++++--- code/mesh.f90 | 27 ++++++++++++++------------- 4 files changed, 73 insertions(+), 38 deletions(-) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index d2a83b834..a729363d5 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -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 @@ -1448,7 +1451,9 @@ LpLoop: do ! restore old residuum and Lp 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 diff --git a/code/debug.f90 b/code/debug.f90 index 992c01ded..baf104261 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -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,*) diff --git a/code/lattice.f90 b/code/lattice.f90 index e446ac914..165d94954 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -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() !************************************** diff --git a/code/mesh.f90 b/code/mesh.f90 index 366c71749..efd1912ba 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -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