From 1dbd0865dbedfbee92bb7eba0ddc6a2d4da243d8 Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Mon, 24 Aug 2009 08:16:01 +0000 Subject: [PATCH] constitutive_nonlocal.f90 - completed postResults output function - connecting vector of neighboring material points is mapped to intermediate configuration of my neighbor crystallite.f90 - zero out dotState only when crystallite is non-finished - set nonfinished flag to false if crystallite is not on Track after state update - in updateState: set onTrack flag to false if encounter NaN - removed some old debugging outputs and added others homogenization.f90 - in debugging mode now telling when a cutback happens --- code/constitutive_nonlocal.f90 | 186 +++++++++++++++++++++++---------- code/crystallite.f90 | 99 +++++++++++------- code/homogenization.f90 | 23 +++- 3 files changed, 215 insertions(+), 93 deletions(-) diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 5a520844f..806402475 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -138,6 +138,7 @@ integer(pInt) section, & s, & ! index of my slip system s1, & ! index of my slip system s2, & ! index of my slip system + it, & ! index of my interaction type output, & mySize character(len=64) tag @@ -250,7 +251,7 @@ do case ('burgers') forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_burgersBySlipFamily(f,i) = IO_floatValue(line,positions,1+f) case ('interaction_slipslip') - forall (f = 1:lattice_maxNinteraction) constitutive_nonlocal_interactionSlipSlip(f,i) = IO_floatValue(line,positions,1+f) + forall (it = 1:lattice_maxNinteraction) constitutive_nonlocal_interactionSlipSlip(it,i) = IO_floatValue(line,positions,1+it) end select endif enddo @@ -267,7 +268,6 @@ enddo if (constitutive_nonlocal_structure(i) < 1 .or. constitutive_nonlocal_structure(i) > 3) call IO_error(205) if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0) call IO_error(225) - do f = 1,lattice_maxNslipFamily if (constitutive_nonlocal_Nslip(f,i) > 0) then if (constitutive_nonlocal_rhoEdgePos0(f,i) < 0.0_pReal) call IO_error(220) @@ -279,7 +279,8 @@ enddo if (constitutive_nonlocal_lambda0BySlipFamily(f,i) <= 0.0_pReal) call IO_error(-1) endif enddo - + if (sum(constitutive_nonlocal_interactionSlipSlip(:,i)) <= 0) call IO_error(-1) + !*** determine total number of active slip systems @@ -331,10 +332,15 @@ do i = 1,maxNinstance do o = 1,maxval(phase_Noutput) select case(constitutive_nonlocal_output(o,i)) - case( 'dislocationdensity', & - 'shearrate_slip', & - 'resolvedstress_slip', & - 'resistance_slip') + case( 'rho', & + 'rho_edge', & + 'rho_screw', & + 'excess_rho_edge', & + 'excess_rho_screw', & + 'rho_forest', & + 'shearrate', & + 'resolvedstress', & + 'resistance') mySize = constitutive_nonlocal_totalNslip(i) case default mySize = 0_pInt @@ -373,7 +379,7 @@ do i = 1,maxNinstance constitutive_nonlocal_Cslip_66(:,:,i) = math_Mandel3333to66(math_Voigt66to3333(constitutive_nonlocal_Cslip_66(:,:,i))) constitutive_nonlocal_Cslip_3333(:,:,:,:,i) = math_Voigt66to3333(constitutive_nonlocal_Cslip_66(:,:,i)) - constitutive_nonlocal_Gmod(i) = constitutive_nonlocal_C44(i) ! shear modulus is given by elastic constant C44 + constitutive_nonlocal_Gmod(i) = constitutive_nonlocal_C44(i) constitutive_nonlocal_nu(i) = constitutive_nonlocal_C12(i) / constitutive_nonlocal_C11(i) @@ -559,7 +565,6 @@ use math, only: math_Plain3333to99, & math_mul33x33, & math_mul3x3, & math_mul33x3, & - math_transpose3x3, & pi use mesh, only: mesh_NcpElems, & mesh_maxNips, & @@ -680,7 +685,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! calculate the connecting vector between me and my neighbor and his excess dislocation density - connectingVector = math_mul33x3( Fp(:,:,g,ip,el), & + connectingVector = math_mul33x3( Fp(:,:,g,neighboring_ip,neighboring_el), & (mesh_ipCenterOfGravity(:,ip,el) - mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el)) ) neighboring_rhoEdgePos = state(1, neighboring_ip, neighboring_el)%p( 1: ns) @@ -732,7 +737,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) sigma(3,1) = 0.0_pReal ! coordinate transformation from the slip coordinate system to the lattice coordinate system - backStress_v = backStress_v + math_Mandel33to6( math_mul33x33(math_transpose3x3(transform), math_mul33x33(sigma, transform) ) ) + backStress_v = backStress_v + math_Mandel33to6( math_mul33x33(transpose(transform), math_mul33x33(sigma, transform) ) ) enddo enddo @@ -830,12 +835,11 @@ backStress_v = state(g,ip,el)%p(6*ns+1:6*ns+6) ! if (debugger) write(6,'(a15,3(i3,x),/,3(3(f12.3,x)/))') 'Tstar / MPa at ',g,ip,el, math_Mandel6to33(Tstar_v/1e6) !*** loop over slip systems - -do s = 1,constitutive_nonlocal_totalNslip(myInstance) + +do s = 1,ns sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) - !*** Calculation of Lp tauSlip(s) = math_mul6x6( Tstar_v + backStress_v, lattice_Sslip_v(:,sLattice,myStructure) ) @@ -849,7 +853,7 @@ do s = 1,constitutive_nonlocal_totalNslip(myInstance) gdotSlip(s) = sum(rho(:,s)) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance) * v(s) Lp = Lp + gdotSlip(s) * lattice_Sslip(:,:,sLattice,myStructure) - + ! if (debugger) write(6,'(a4,i2,a3,/,3(3(f15.7)/))') 'dLp(',s,'): ',gdotSlip(s) * lattice_Sslip(:,:,sLattice,myStructure) !*** Calculation of the tangent of Lp @@ -863,9 +867,12 @@ enddo dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) -!if (debugger) write(6,'(a26,3(i3,x),/,12(f10.3,x),/)') 'tauSlipThreshold / MPa at ',g,ip,el, tauSlipThreshold/1e6 -!if (debugger) write(6,'(a15,3(i3,x),/,12(f10.3,x),/)') 'tauSlip / MPa at ',g,ip,el, tauSlip/1e6 -!if (debugger) write(6,'(a15,3(i3,x),/,12(e10.3,x),/)') 'gdotSlip at ',g,ip,el, gdotSlip +! if (debugger) write(6,'(a23,3(i3,x),/,12(e10.3,x),/)') 'dislocation density at ',g,ip,el, rho +! if (debugger) write(6,'(a26,3(i3,x),/,12(f10.5,x),/)') 'tauSlipThreshold / MPa at ',g,ip,el, tauSlipThreshold/1e6 +! if (debugger) write(6,'(a15,3(i3,x),/,12(f10.5,x),/)') 'tauSlip / MPa at ',g,ip,el, tauSlip/1e6 +! if (debugger) write(6,'(a5,3(i3,x),/,12(e10.3,x),/)') 'v at ',g,ip,el, v +! if (debugger) write(6,'(a15,3(i3,x),/,12(e10.3,x),/)') 'gdotSlip at ',g,ip,el, gdotSlip +! if (debugger) write(6,'(a6,3(i3,x),/,3(3(f15.7)/))') 'Lp at ',g,ip,el, Lp endsubroutine @@ -983,16 +990,18 @@ do s = 1,ns gdot(t,s) = rho(t,s) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance) * v(s) enddo - + !**************************************************************************** !*** calculate dislocation multiplication invLambda = sqrt(rhoForest) / constitutive_nonlocal_lambda0BySlipSystem(:,myInstance) -forall (t = 1:4, s = 1:ns) & - dotState(1,ip,el)%p((t-1)*ns+s) = dotState(1,ip,el)%p((t-1)*ns+s) + 0.25_pReal * sum(gdot(:,s)) * invLambda(s) & - / constitutive_nonlocal_burgersBySlipSystem(s,myInstance) +forall (t = 1:4) & + dotState(1,ip,el)%p((t-1)*ns+1:t*ns) = dotState(1,ip,el)%p((t-1)*ns+1:t*ns) + 0.25_pReal * sum(abs(gdot),1) * invLambda & + / constitutive_nonlocal_burgersBySlipSystem(:,myInstance) +! if (debugger) write(6,'(a30,3(i3,x),/,12(e10.3,x),/)') 'dislocation multiplication at ',g,ip,el, & + ! 0.25_pReal * sum(abs(gdot),1) * invLambda / constitutive_nonlocal_burgersBySlipSystem(:,myInstance) !**************************************************************************** @@ -1086,49 +1095,116 @@ endfunction !********************************************************************* !* return array of constitutive results * -!* INPUT: * -!* - Temperature : temperature * -!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * -!* - dt : current time increment * -!* - ipc : component-ID at current integration point * -!* - ip : current integration point * -!* - el : current element * !********************************************************************* -pure function constitutive_nonlocal_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) +pure function constitutive_nonlocal_postResults(Tstar_v, Temperature, dt, state, g, ip, el) -use prec, only: pReal,pInt,p_vec -use mesh, only: mesh_NcpElems,mesh_maxNips -use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput -use lattice, only: lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin +use prec, only: pReal, & + pInt, & + p_vec +use math, only: math_mul6x6 +use mesh, only: mesh_NcpElems, & + mesh_maxNips +use material, only: homogenization_maxNgrains, & + material_phase, & + phase_constitutionInstance, & + phase_Noutput +use lattice, only: lattice_Sslip_v, & + lattice_NslipSystem implicit none -!* Definition of variables -integer(pInt), intent(in) :: ipc,ip,el -real(pReal), intent(in) :: dt,Temperature -real(pReal), dimension(6), intent(in) :: Tstar_v -type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state -integer(pInt) matID,structID,ns,nt,f,o,i,c,j,index_myFamily -real(pReal) sumf,tau -real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & - constitutive_nonlocal_postResults +!*** input variables +integer(pInt), intent(in) :: g, & ! current grain number + ip, & ! current integration point + el ! current element number +real(pReal), intent(in) :: dt, & ! time increment + Temperature ! temperature +real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation +type(p_vec), dimension(homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems), intent(in) :: & + state ! microstructural state -!* Shortened notation -matID = phase_constitutionInstance(material_phase(ipc,ip,el)) -structID = constitutive_nonlocal_structure(matID) -ns = constitutive_nonlocal_totalNslip(matID) +!*** output variables +real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_constitutionInstance(material_phase(g,ip,el)))) :: & + constitutive_nonlocal_postResults + +!*** local variables +integer(pInt) myInstance, & ! current instance of this constitution + myStructure, & ! current lattice structure + ns, & ! short notation for the total number of active slip systems + o, & ! index of current output + s, & ! index of current slip system + sLattice, & ! index of current slip system as specified by lattice + c +real(pReal) tau, & ! resolved shear stress on current slip system + v ! dislocation velocity on current slip system + + +myInstance = phase_constitutionInstance(material_phase(g,ip,el)) +myStructure = constitutive_nonlocal_structure(myInstance) +ns = constitutive_nonlocal_totalNslip(myInstance) -!* Required output c = 0_pInt constitutive_nonlocal_postResults = 0.0_pReal -do o = 1,phase_Noutput(material_phase(ipc,ip,el)) - select case(constitutive_nonlocal_output(o,matID)) - - case ('dislocationdensity') - constitutive_nonlocal_postResults(c+1:c+ns) = state(ipc,ip,el)%p(1:ns) - c = c + ns +do o = 1,phase_Noutput(material_phase(g,ip,el)) + select case(constitutive_nonlocal_output(o,myInstance)) + + case ('rho') + constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns) + state(g,ip,el)%p(ns+1:2*ns) & + + state(g,ip,el)%p(2*ns+1:3*ns) + state(g,ip,el)%p(3*ns+1:4*ns) + c = c + ns + + case ('rho_edge') + constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns) + state(g,ip,el)%p(ns+1:2*ns) + c = c + ns + + case ('rho_screw') + constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(2*ns+1:3*ns) + state(g,ip,el)%p(3*ns+1:4*ns) + c = c + ns + + case ('excess_rho_edge') + constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns) - state(g,ip,el)%p(ns+1:2*ns) + c = c + ns + + case ('excess_rho_screw') + constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(2*ns+1:3*ns) - state(g,ip,el)%p(3*ns+1:4*ns) + c = c + ns + + case ('rho_forest') + constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(4*ns+1:5*ns) + c = c + ns + + case ('shearrate') + do s = 1,ns + sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) + tau = math_mul6x6( Tstar_v + state(g,ip,el)%p(6*ns+1:6*ns+6), lattice_Sslip_v(:,sLattice,myStructure) ) + + if (state(g,ip,el)%p(4*ns+s) > 0.0_pReal) then + v = constitutive_nonlocal_v0BySlipSystem(s,myInstance) & + * exp( - ( state(g,ip,el)%p(5*ns+s) - abs(tau) ) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance)**2.0_pReal & + / ( kB * Temperature * sqrt(state(g,ip,el)%p(4*ns+s)) ) ) & + * sign(1.0_pReal,tau) + else + v = 0.0_pReal + endif + + constitutive_nonlocal_postResults(c+s) = ( state(g,ip,el)%p(s) + state(g,ip,el)%p(ns+s) & + + state(g,ip,el)%p(2*ns+s) + state(g,ip,el)%p(3*ns+s) ) & + * constitutive_nonlocal_burgersBySlipSystem(s,myInstance) * v + enddo + c = c + ns + + case ('resolvedstress') + do s = 1,ns + sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) + constitutive_nonlocal_postResults(c+s) = math_mul6x6( Tstar_v + state(g,ip,el)%p(6*ns+1:6*ns+6), & + lattice_Sslip_v(:,sLattice,myStructure) ) + enddo + c = c + ns + + case ('resistance') + constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(5*ns+1:6*ns) + c = c + ns end select enddo diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 28eda5eb6..1a05340e8 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -204,7 +204,7 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_temperatureConverged: ', shape(crystallite_temperatureConverged) write(6,'(a35,x,7(i5,x))') 'crystallite_nonfinished: ', shape(crystallite_nonfinished) write(6,*) - write(6,*) 'Number of non-local grains: ',count(.not. crystallite_localConstitution) + write(6,*) 'Number of nonlocal grains: ',count(.not. crystallite_localConstitution) call flush(6) !$OMPEND CRITICAL (write2out) @@ -347,6 +347,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains + debugger = (e == 1 .and. i == 1 .and. g == 1) if (crystallite_converged(g,i,e)) then if (debugger) then !$OMP CRITICAL (write2out) @@ -424,7 +425,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState + if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites + constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState enddo; enddo; enddo !$OMPEND PARALLEL DO !$OMP PARALLEL DO @@ -470,11 +472,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains - debugger = (e == 1 .and. i == 1 .and. g == 1) - if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites - crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e) - enddo + do g = 1,myNgrains + debugger = (e == 1 .and. i == 1 .and. g == 1) + if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites + crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e) + enddo enddo enddo !$OMPEND PARALLEL DO @@ -494,7 +496,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains - constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState + if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites + constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState enddo; enddo; enddo !$OMPEND PARALLEL DO !$OMP PARALLEL DO @@ -502,6 +505,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains + debugger = (e == 1 .and. i == 1 .and. g == 1) if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), & crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e) @@ -517,7 +521,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) .and. crystallite_temperatureConverged(g,i,e) - if (debugger) write (6,*) g,i,e,'converged after updState',crystallite_converged(g,i,e) if (crystallite_converged(g,i,e)) then !$OMP CRITICAL (distributionState) debug_CrystalliteStateLoopDistribution(NiterationState) = & @@ -529,15 +532,18 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMPEND PARALLEL DO - - crystallite_nonfinished = crystallite_nonfinished .and. .not. crystallite_converged + + write(6,*) 'NiterationState: ',NiterationState + crystallite_nonfinished = crystallite_nonfinished .and. crystallite_onTrack .and. .not. crystallite_converged enddo ! crystallite convergence loop NiterationCrystallite = NiterationCrystallite + 1 - + enddo ! cutback loop - + + ! write (6,'(a,/,8(L,x))') 'crystallite_nonfinished',crystallite_nonfinished + ! write (6,'(a,/,8(L,x))') 'crystallite_converged',crystallite_converged ! ------ check for non-converged crystallites ------ @@ -614,7 +620,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) NiterationState = NiterationState + 1_pInt onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) if (onTrack) then - call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), & + call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_Fp(:,:,g,i,e), & crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e) stateConverged = crystallite_updateState(g,i,e) ! update state @@ -645,8 +651,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_P(:,:,g,i,e) = myP !$OMP CRITICAL (out) debug_StiffnessStateLoopDistribution(NiterationState) = & - debug_StiffnessstateLoopDistribution(NiterationState) + 1 - if (nState < NiterationState) write(6,*) 'ohh shit!! stiffenss state loop debugging exceeded',NiterationState + debug_StiffnessstateLoopDistribution(NiterationState) + 1 !$OMPEND CRITICAL (out) enddo enddo @@ -718,9 +723,13 @@ endsubroutine ! if NaN occured then return without changing the state if (any(residuum/=residuum)) then crystallite_updateState = .false. ! indicate state update failed - !$OMP CRITICAL (write2out) - write(6,*) '::: updateState encountered NaN',e,i,g - !$OMPEND CRITICAL (write2out) + crystallite_onTrack(g,i,e) = .false. ! no need to calculate any further + crystallite_onTrack = crystallite_onTrack .and. crystallite_localConstitution ! all nonlocal crystallites have to be redone + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: updateState encountered NaN',g,i,e + !$OMPEND CRITICAL (write2out) + endif return endif @@ -731,7 +740,19 @@ endsubroutine crystallite_updateState = all(constitutive_state(g,i,e)%p(1:mySize) == 0.0_pReal .or. & abs(residuum) < rTol_crystalliteState*abs(constitutive_state(g,i,e)%p(1:mySize))) if (debugger) then - write(6,'(a,/,12(f10.5,x))') 'resid tolerance',abs(residuum/rTol_crystalliteState/constitutive_state(g,i,e)%p(1:mySize)) + !$OMP CRITICAL (write2out) + if (crystallite_updateState) then + write(6,*) '::: updateState converged',g,i,e + write(6,*) + write(6,'(a10,/,12(e12.3,x))') 'new state ',constitutive_state(g,i,e)%p(1:mySize) + write(6,*) + else + write(6,*) '::: updateState did not converge',g,i,e + write(6,*) + endif + write(6,'(a,/,12(f10.5,x))') 'resid tolerance',abs(residuum/rTol_crystalliteState/constitutive_state(g,i,e)%p(1:mySize)) + write(6,*) + !$OMPEND CRITICAL (write2out) endif return @@ -787,7 +808,7 @@ endsubroutine if (residuum/=residuum) then crystallite_updateTemperature = .false. ! indicate update failed !$OMP CRITICAL (write2out) - write(6,*) '::: updateTemperature encountered NaN',e,i,g + write(6,*) '::: updateTemperature encountered NaN',g,i,e !$OMPEND CRITICAL (write2out) return endif @@ -795,7 +816,7 @@ endsubroutine ! update the microstructure crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum - ! setting flag to true if residuum is below relative tolerance (or zero Kelvin), otherwise set it to false <<>> + ! setting flag to true if residuum is below relative tolerance (or zero Kelvin), otherwise set it to false crystallite_updateTemperature = crystallite_Temperature(g,i,e) == 0.0_pReal .or. & abs(residuum) < rTol_crystalliteTemperature*crystallite_Temperature(g,i,e) @@ -910,9 +931,13 @@ endsubroutine ! inversion of Fp_current... invFp_current = math_inv3x3(Fp_current) if (all(invFp_current == 0.0_pReal)) then ! ... failed? - !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress failed on invFp_current inversion',e,i,g - !$OMPEND CRITICAL (write2out) + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: integrateStress failed on invFp_current inversion',g,i,e + write(6,*) + write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new + !$OMPEND CRITICAL (write2out) + endif return endif @@ -973,9 +998,11 @@ LpLoop: do ! NaN occured at regular speed? if (any(residuum/=residuum) .and. leapfrog == 1.0) then - !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',e,i,g - !$OMPEND CRITICAL (write2out) + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',g,i,e + !$OMPEND CRITICAL (write2out) + endif return ! something went wrong at accelerated speed? @@ -1008,11 +1035,11 @@ LpLoop: do if (error) then if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress - write(6,'(a9,3(i3,x),/,9(9(e12.2,x)/))') 'dTdLp at ',g,i,e,dTdLp - write(6,'(a20,3(i3,x),/,9(9(e12.2,x)/))') 'dLp_constitutive at ',g,i,e,dLp_constitutive - write(6,'(a9,3(i3,x),/,9(9(f12.7,x)/))') 'dRdLp at ',g,i,e,dRdLp - write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'Lpguess at ',g,i,e,Lpguess + write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress + write(6,*) + write(6,'(a9,3(i3,x),/,9(9(f12.7,x)/))') 'dRdLp at ',g,i,e,dRdLp + write(6,'(a20,3(i3,x),/,9(9(e12.2,x)/))') 'dLp_constitutive at ',g,i,e,dLp_constitutive + write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'Lpguess at ',g,i,e,Lpguess !$OMPEND CRITICAL (write2out) endif return @@ -1040,7 +1067,9 @@ LpLoop: do if (error) then if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress + write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress + write(6,*) + write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new !$OMPEND CRITICAL (write2out) endif return @@ -1068,12 +1097,12 @@ LpLoop: do write(6,*) write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',crystallite_P(:,:,g,i,e)/1e6 write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',crystallite_Lp(:,:,g,i,e) + write(6,'(a,/,3(3(f12.7,x)/))') 'Fp',crystallite_Fp(:,:,g,i,e) !$OMP CRITICAL (write2out) endif !$OMP CRITICAL (distributionStress) debug_StressLoopDistribution(NiterationStress) = debug_StressLoopDistribution(NiterationStress) + 1 - if (nStress < NiterationStress) write(6,*) 'ohh shit!! debug loop of stress exceeded',NiterationStress !$OMPEND CRITICAL (distributionStress) return diff --git a/code/homogenization.f90 b/code/homogenization.f90 index e9af941c8..1ccfb98e3 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -213,7 +213,8 @@ subroutine materialpoint_stressAndItsTangent(& crystallite_requested, & crystallite_converged, & crystallite_stressAndItsTangent - use debug, only: debug_MaterialpointLoopDistribution, & + use debug, only: debugger, & + debug_MaterialpointLoopDistribution, & debug_MaterialpointStateLoopDistribution implicit none @@ -270,6 +271,8 @@ subroutine materialpoint_stressAndItsTangent(& myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + ! debugger = (e == 1 .and. i == 1) + ! if our materialpoint converged then we are either finished or have to wind forward if (materialpoint_converged(i,e)) then @@ -277,6 +280,14 @@ subroutine materialpoint_stressAndItsTangent(& materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e) materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), 2.0_pReal * materialpoint_subStep(i,e)) + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,'(a21,f10.8,a34,f10.8,a37,/)') 'winding forward from ', & + materialpoint_subFrac(i,e) - materialpoint_subStep(i,e),' to current materialpoint_subFrac ', & + materialpoint_subFrac(i,e),' in materialpoint_stressAndItsTangent' + !$OMPEND CRITICAL (write2out) + endif + ! still stepping needed if (materialpoint_subStep(i,e) > subStepMin) then @@ -301,6 +312,13 @@ subroutine materialpoint_stressAndItsTangent(& else materialpoint_subStep(i,e) = 0.5_pReal * materialpoint_subStep(i,e) + + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,'(a82,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 @@ -415,8 +433,7 @@ elementLoop: do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate !$OMP END PARALLEL DO write (6,*) 'Material Point finished' - write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 1 1',crystallite_Lp(1:3,:,1,1,1) - + ! how to deal with stiffness? return