From 2aae7b7574d86ec66672e2cb51fb478a332335cf Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Wed, 12 Aug 2009 11:22:02 +0000 Subject: [PATCH] added dislocation multiplication to dotState in constitutive_nonlocal.f90 cleaned up debugging output statements to *.out file --- code/constitutive_nonlocal.f90 | 151 +++++++++++++++++++-------------- code/crystallite.f90 | 95 +++++++++++++-------- code/math.f90 | 3 +- code/mesh.f90 | 89 ++++++++++--------- 4 files changed, 192 insertions(+), 146 deletions(-) diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 8bbda8251..5a520844f 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -57,8 +57,10 @@ real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_ constitutive_nonlocal_rhoEdgeNeg0, & ! initial edge_neg dislocation density per slip system for each family and instance constitutive_nonlocal_rhoScrewPos0, & ! initial screw_pos dislocation density per slip system for each family and instance constitutive_nonlocal_rhoScrewNeg0, & ! initial screw_neg dislocation density per slip system for each family and instance - constitutive_nonlocal_v0BySlipFamily, & ! initial dislocation velocity [m/s] for each family and instance - constitutive_nonlocal_v0BySlipSystem, & ! initial dislocation velocity [m/s] for each slip system and instance + constitutive_nonlocal_v0BySlipFamily, & ! dislocation velocity prefactor [m/s] for each family and instance + constitutive_nonlocal_v0BySlipSystem, & ! dislocation velocity prefactor [m/s] for each slip system and instance + constitutive_nonlocal_lambda0BySlipFamily, & ! mean free path prefactor for each family and instance + constitutive_nonlocal_lambda0BySlipSystem, & ! mean free path prefactor for each slip system and instance constitutive_nonlocal_burgersBySlipFamily, & ! absolute length of burgers vector [m] for each family and instance constitutive_nonlocal_burgersBySlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance @@ -184,6 +186,8 @@ allocate(constitutive_nonlocal_v0BySlipFamily(lattice_maxNslipFamily, maxNinstan constitutive_nonlocal_v0BySlipFamily = 0.0_pReal allocate(constitutive_nonlocal_burgersBySlipFamily(lattice_maxNslipFamily, maxNinstance)); constitutive_nonlocal_burgersBySlipFamily = 0.0_pReal +allocate(constitutive_nonlocal_Lambda0BySlipFamily(lattice_maxNslipFamily, maxNinstance)); + constitutive_nonlocal_lambda0BySlipFamily = 0.0_pReal allocate(constitutive_nonlocal_interactionSlipSlip(lattice_maxNinteraction, maxNinstance)) constitutive_nonlocal_interactionSlipSlip = 0.0_pReal @@ -241,6 +245,8 @@ do forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoScrewNeg0(f,i) = IO_floatValue(line,positions,1+f) case ('v0') forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_v0BySlipFamily(f,i) = IO_floatValue(line,positions,1+f) + case ('lambda0') + forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_lambda0BySlipFamily(f,i) = IO_floatValue(line,positions,1+f) case ('burgers') forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_burgersBySlipFamily(f,i) = IO_floatValue(line,positions,1+f) case ('interaction_slipslip') @@ -256,29 +262,23 @@ enddo lattice_initializeStructure(constitutive_nonlocal_structureName(i), constitutive_nonlocal_CoverA(i)) ! our lattice structure is defined in the material.config file by the structureName (and the c/a ratio) -!*** sanity checks +!*** sanity checks +!*** !!! not yet complete !!! - if (constitutive_nonlocal_structure(i) < 1 .or. constitutive_nonlocal_structure(i) > 3) call IO_error(205) + 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) - if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0) call IO_error(225) - - if (any( constitutive_nonlocal_rhoEdgePos0(:,i) < 0.0_pReal & - .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(220) - - if (any( constitutive_nonlocal_rhoEdgeNeg0(:,i) < 0.0_pReal & - .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(220) - - if (any( constitutive_nonlocal_rhoScrewPos0(:,i) < 0.0_pReal & - .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(220) - - if (any( constitutive_nonlocal_rhoScrewNeg0(:,i) < 0.0_pReal & - .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(220) - - if (any( constitutive_nonlocal_burgersBySlipFamily(:,i) <= 0.0_pReal & - .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(221) - - if (any( constitutive_nonlocal_v0BySlipFamily(:,i) <= 0.0_pReal & - .and. constitutive_nonlocal_Nslip(:,i) > 0)) call IO_error(-1) + 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) + if (constitutive_nonlocal_rhoEdgeNeg0(f,i) < 0.0_pReal) call IO_error(220) + if (constitutive_nonlocal_rhoScrewPos0(f,i) < 0.0_pReal) call IO_error(220) + if (constitutive_nonlocal_rhoScrewNeg0(f,i) < 0.0_pReal) call IO_error(220) + if (constitutive_nonlocal_burgersBySlipFamily(f,i) <= 0.0_pReal) call IO_error(221) + if (constitutive_nonlocal_v0BySlipFamily(f,i) <= 0.0_pReal) call IO_error(-1) + if (constitutive_nonlocal_lambda0BySlipFamily(f,i) <= 0.0_pReal) call IO_error(-1) + endif + enddo !*** determine total number of active slip systems @@ -297,6 +297,8 @@ allocate(constitutive_nonlocal_burgersBySlipSystem(maxTotalNslip, maxNinstance)) constitutive_nonlocal_burgersBySlipSystem = 0.0_pReal allocate(constitutive_nonlocal_v0BySlipSystem(maxTotalNslip, maxNinstance)) constitutive_nonlocal_v0BySlipSystem = 0.0_pReal +allocate(constitutive_nonlocal_lambda0BySlipSystem(maxTotalNslip, maxNinstance)) + constitutive_nonlocal_lambda0BySlipSystem = 0.0_pReal allocate(constitutive_nonlocal_forestProjectionEdge(maxTotalNslip, maxTotalNslip, maxNinstance)) constitutive_nonlocal_forestProjectionEdge = 0.0_pReal allocate(constitutive_nonlocal_forestProjectionScrew(maxTotalNslip, maxTotalNslip, maxNinstance)) @@ -375,12 +377,18 @@ do i = 1,maxNinstance constitutive_nonlocal_nu(i) = constitutive_nonlocal_C12(i) / constitutive_nonlocal_C11(i) -!*** burgers vector and initial dislocation velocity for each slip system +!*** burgers vector, dislocation velocity prefactor and mean free path prefactor for each slip system do s = 1,constitutive_nonlocal_totalNslip(i) + constitutive_nonlocal_burgersBySlipSystem(s,i) & = constitutive_nonlocal_burgersBySlipFamily( constitutive_nonlocal_slipFamily(s,i), i ) + constitutive_nonlocal_v0BySlipSystem(s,i) = constitutive_nonlocal_v0BySlipFamily(constitutive_nonlocal_slipFamily(s,i),i) + + constitutive_nonlocal_lambda0BySlipSystem(s,i) & + = constitutive_nonlocal_lambda0BySlipFamily( constitutive_nonlocal_slipFamily(s,i), i ) + enddo @@ -852,12 +860,13 @@ do s = 1,constitutive_nonlocal_totalNslip(myInstance) dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) + dgdot_dtauSlip(s) * lattice_Sslip(i,j, sLattice,myStructure) & * lattice_Sslip(k,l, sLattice,myStructure) 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 -dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) - endsubroutine @@ -882,7 +891,6 @@ use mesh, only: mesh_NcpElems, & FE_NipNeighbors, & mesh_ipNeighborhood, & mesh_ipVolume, & - mesh_ipCenterOfGravity, & mesh_ipArea, & mesh_ipAreaNormal use material, only: homogenization_maxNgrains, & @@ -922,30 +930,32 @@ integer(pInt) myInstance, & ! current neighboring_ip, & ! integration point of my neighbor n, & ! index of my current neighbor t, & ! type of dislocation - s, & ! index of my current slip system - sLattice ! index of my current slip system as specified by lattice + s ! index of my current slip system real(pReal), dimension(4,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & rho, & ! dislocation densities - gdot ! shear rates + gdot, & ! shear rates + lineLength ! dislocation line length leaving the current interface real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & rhoForest, & ! forest dislocation density tauSlipThreshold, & ! threshold shear stress tauSlip, & ! resolved shear stress - v ! dislocation velocity + v, & ! dislocation velocity + invLambda ! inverse of mean free path for dislocations real(pReal), dimension(3,4,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & m ! direction of dislocation motion real(pReal), dimension(6) :: backStress_v ! backstress resulting from the neighboring excess dislocation densities as 2nd Piola-Kirchhoff stress -real(pReal), dimension(4) :: leavingRho ! dislocation densities leaving the current interface real(pReal), dimension(3) :: surfaceNormal ! surface normal of the current interface real(pReal) norm_surfaceNormal, & ! euclidic norm of the surface normal area ! area of the current interface - myInstance = phase_constitutionInstance(material_phase(g,ip,el)) myStructure = constitutive_nonlocal_structure(myInstance) ns = constitutive_nonlocal_totalNslip(myInstance) +tauSlip = 0.0_pReal +v = 0.0_pReal +gdot = 0.0_pReal !*** shortcut to state variables @@ -955,65 +965,80 @@ tauSlipThreshold = state(g,ip,el)%p(5*ns+1:6*ns) backStress_v = state(g,ip,el)%p(6*ns+1:6*ns+6) +!**************************************************************************** +!*** Calculate shear rate + do s = 1,ns -!*** Calculate gdot for each type of dislocation - - sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) - tauSlip(s) = math_mul6x6( Tstar_v + backStress_v, lattice_Sslip_v(:,sLattice,myStructure) ) - - v(s) = constitutive_nonlocal_v0BySlipSystem(s,myInstance) & - * exp( - ( tauSlipThreshold(s) - abs(tauSlip(s)) ) & - * constitutive_nonlocal_burgersBySlipSystem(s,myInstance)**2.0_pReal & - / ( kB * Temperature * sqrt(rhoForest(s)) ) ) & - * sign(1.0_pReal,tauSlip(s)) - - forall (t = 1:4) gdot(t,s) = rho(t,s) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance) * v(s) + tauSlip(s) = math_mul6x6( Tstar_v + backStress_v, & + lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) ) - -!*** Direction of dislocation motion - - m(:,1,s) = lattice_sd(:, sLattice, myStructure) - m(:,2,s) = -lattice_sd(:, sLattice, myStructure) - m(:,3,s) = lattice_st(:, sLattice, myStructure) - m(:,4,s) = -lattice_st(:, sLattice, myStructure) + forall (s = 1:ns, rhoForest(s) > 0.0_pReal) & + v(s) = constitutive_nonlocal_v0BySlipSystem(s,myInstance) & + * exp( - ( tauSlipThreshold(s) - abs(tauSlip(s)) ) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance)**2.0_pReal & + / ( kB * Temperature * sqrt(rhoForest(s)) ) ) & + * sign(1.0_pReal,tauSlip(s)) + + forall (t = 1:4, 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) -!*** loop through my neighbors if existent and calculate the area and the surface normal of the interface +!**************************************************************************** +!*** calculate dislocation fluxes +! Direction of dislocation motion +m(:,1,:) = lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) +m(:,2,:) = -lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) +m(:,3,:) = lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) +m(:,4,:) = -lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) + +! loop through my neighbors do n = 1,FE_NipNeighbors(mesh_element(2,el)) neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) + ! calculate the area and the surface normal of the interface surfaceNormal = math_mul33x3(math_transpose3x3(invFp), mesh_ipAreaNormal(:,n,ip,el)) norm_surfaceNormal = math_norm3(surfaceNormal) surfaceNormal = surfaceNormal / norm_surfaceNormal area = mesh_ipArea(n,ip,el) / norm_surfaceNormal - -!*** loop through my interfaces -!*** subtract dislocation densities that leave through an interface from my dotState and add them to the neighboring dotState - + ! loop through my interfaces do s = 1,ns - leavingRho = 0.0_pReal + lineLength = 0.0_pReal + + ! loop through dislocation types do t = 1,4 if ( sign(1.0_pReal,math_mul3x3(m(:,t,s),surfaceNormal)) == sign(1.0_pReal,gdot(t,s)) ) then - leavingRho(t) = gdot(t,s) / constitutive_nonlocal_burgersBySlipSystem(s,myInstance) & - * math_mul3x3(m(:,t,s),surfaceNormal) * area + ! dislocation line length that leaves this interface per second + lineLength(t,s) = gdot(t,s) / constitutive_nonlocal_burgersBySlipSystem(s,myInstance) & + * math_mul3x3(m(:,t,s),surfaceNormal) * area - dotState(1,ip,el)%p((t-1)*ns+s) = dotState(1,ip,el)%p((t-1)*ns+s) - leavingRho(t) + ! subtract dislocation density rate (= line length over volume) that leaves through an interface from my dotState ... + dotState(1,ip,el)%p((t-1)*ns+s) = dotState(1,ip,el)%p((t-1)*ns+s) - lineLength(t,s) / mesh_ipVolume(ip,el) + ! ... and add them to the neighboring dotState (if neighbor exists) if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then !***************************************************************************************************** !*** OMP locking for this neighbor missing !***************************************************************************************************** - dotState(1,neighboring_ip,neighboring_el)%p((t-1)*ns+s) = & - dotState(1,neighboring_ip,neighboring_el)%p((t-1)*ns+s) + leavingRho(t) + dotState(1,neighboring_ip,neighboring_el)%p((t-1)*ns+s) = dotState(1,neighboring_ip,neighboring_el)%p((t-1)*ns+s) & + + lineLength(t,s) / mesh_ipVolume(neighboring_ip,neighboring_el) endif endif enddo diff --git a/code/crystallite.f90 b/code/crystallite.f90 index f69226f3e..d05120b45 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -410,6 +410,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) .and. crystallite_onTrack & .and. .not. crystallite_converged) + ! --+>> preguess for state <<+-- ! ! incrementing by crystallite_subdt @@ -449,6 +450,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo; enddo; enddo !$OMPEND PARALLEL DO + ! --+>> state loop <<+-- NiterationState = 0_pInt @@ -458,6 +460,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) NiterationState = NiterationState + 1_pInt + ! --+>> stress integration <<+-- ! ! incrementing by crystallite_subdt @@ -469,16 +472,15 @@ 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 - crystallite_nonfinished = crystallite_nonfinished .and. crystallite_onTrack ! --+>> state integration <<+-- ! @@ -488,6 +490,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! first loop for collection of state evolution based on old state ! second loop for updating to new state + crystallite_nonfinished = crystallite_nonfinished .and. crystallite_onTrack + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) @@ -501,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) @@ -511,12 +516,11 @@ 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) + ! debugger = (e == 1 .and. i == 1 .and. g == 1) if (crystallite_nonfinished(g,i,e)) then ! all undone crystallites 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) = & @@ -569,7 +573,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 = (g == 1 .and. i == 1 .and. e == 1) + ! debugger = (g == 1 .and. i == 1 .and. e == 1) if (crystallite_converged(g,i,e)) then ! grain converged in above iteration mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain mySizeDotState = constitutive_sizeDotState(g,i,e) ! number of dotStates for this grain @@ -639,8 +643,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 @@ -701,8 +704,6 @@ endsubroutine mySize = constitutive_sizeDotState(g,i,e) ! calculate the residuum - if (any(abs(constitutive_dotState(g,i,e)%p) > 1e10_pReal)) & - write(6,*) 'dotState: crap found at',g,i,e,constitutive_dotState(g,i,e)%p call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) - & crystallite_subdt(g,i,e) * constitutive_dotState(g,i,e)%p(1:mySize) @@ -714,20 +715,33 @@ 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) + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: updateState encountered NaN',g,i,e + write(6,*) + !$OMPEND CRITICAL (write2out) + endif return endif ! update the microstructure constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum - ! setting flag to true if state is below relative tolerance, otherwise set it to false <<>> + ! setting flag to true if state is below relative tolerance, otherwise set it to false 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 did not converge',g,i,e + write(6,*) + else + write(6,*) '::: updateState converged',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 @@ -783,7 +797,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 @@ -791,7 +805,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) @@ -906,9 +920,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 @@ -969,9 +987,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? @@ -1004,11 +1024,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 @@ -1036,7 +1056,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 @@ -1060,16 +1082,15 @@ LpLoop: do crystallite_integrateStress = .true. if (debugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress converged at iteration', NiterationStress - 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,*) '::: integrateStress converged at iteration', NiterationStress + 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) !$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/math.f90 b/code/math.f90 index fb4858808..a25f0de6d 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -707,8 +707,9 @@ function math_transpose3x3(A) real(pReal),dimension(3,3),intent(in) :: A real(pReal),dimension(3,3) :: math_transpose3x3 + integer(pInt) i,j - math_transpose3x3 = reshape( (/A(1,1), A(1,2), A(1,3), A(2,1), A(2,2), A(2,3), A(3,1), A(3,2), A(3,3) /),(/3,3/) ) + forall(i=1:3, j=1:3) math_transpose3x3(i,j) = A(j,i) return endfunction diff --git a/code/mesh.f90 b/code/mesh.f90 index e167b7213..b0d3569fc 100644 --- a/code/mesh.f90 +++ b/code/mesh.f90 @@ -1697,51 +1697,50 @@ subroutine mesh_get_nodeElemDimensions (unit) !$OMP CRITICAL (write2out) - write (6,*) - write (6,*) "Input Parser: IP NEIGHBORHOOD" - write (6,*) - write (6,"(a10,x,a10,x,a10,x,a3,x,a13,x,a13)") "elem","IP","neighbor","","elemNeighbor","ipNeighbor" - do e = 1,mesh_NcpElems ! loop over cpElems - t = mesh_element(2,e) ! get elemType - do i = 1,FE_Nips(t) ! loop over IPs of elem - do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP - write (6,"(i10,x,i10,x,i10,x,a3,x,i13,x,i13)") e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) - enddo - enddo - enddo - write (6,*) - write (6,*) "Input Parser: ELEMENT VOLUME" - write (6,*) - write (6,"(a13,x,e15.8)") "total volume", sum(mesh_ipVolume) - write (6,*) - write (6,"(a5,x,a5,x,a15,x,a5,x,a15,x,a16)") "elem","IP","volume","face","area","-- normal --" - do e = 1,mesh_NcpElems - do i = 1,FE_Nips(mesh_element(2,e)) - write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e) - do f = 1,FE_NipNeighbors(mesh_element(2,e)) - write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) - enddo - enddo - enddo - !write (6,*) - !write (6,*) "Input Parser: SUBNODE COORDINATES" - !write (6,*) - !write(6,'(a5,x,a5,x,a15,x,a15,x,a20,3(x,a8))') 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z' - !do e = 1,mesh_NcpElems ! loop over cpElems - ! t = mesh_element(2,e) ! get elemType - ! do i = 1,FE_Nips(t) ! loop over IPs of elem - ! do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP - ! do n = 1,FE_NipFaceNodes ! loop over nodes on interface - ! write(6,'(i5,x,i5,x,i15,x,i15,x,i20,3(x,f8.3))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),& - ! mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),& - ! mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),& - ! mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e) - ! enddo - ! enddo - ! enddo - !enddo - write (6,*) - write (6,*) + ! write (6,*) + ! write (6,*) "Input Parser: IP NEIGHBORHOOD" + ! write (6,*) + ! write (6,"(a10,x,a10,x,a10,x,a3,x,a13,x,a13)") "elem","IP","neighbor","","elemNeighbor","ipNeighbor" + ! do e = 1,mesh_NcpElems ! loop over cpElems + ! t = mesh_element(2,e) ! get elemType + ! do i = 1,FE_Nips(t) ! loop over IPs of elem + ! do n = 1,FE_NipNeighbors(t) ! loop over neighbors of IP + ! write (6,"(i10,x,i10,x,i10,x,a3,x,i13,x,i13)") e,i,n,'-->',mesh_ipNeighborhood(1,n,i,e),mesh_ipNeighborhood(2,n,i,e) + ! enddo + ! enddo + ! enddo + ! write (6,*) + ! write (6,*) "Input Parser: ELEMENT VOLUME" + ! write (6,*) + ! write (6,"(a13,x,e15.8)") "total volume", sum(mesh_ipVolume) + ! write (6,*) + ! write (6,"(a5,x,a5,x,a15,x,a5,x,a15,x,a16)") "elem","IP","volume","face","area","-- normal --" + ! do e = 1,mesh_NcpElems + ! do i = 1,FE_Nips(mesh_element(2,e)) + ! write (6,"(i5,x,i5,x,e15.8)") e,i,mesh_IPvolume(i,e) + ! do f = 1,FE_NipNeighbors(mesh_element(2,e)) + ! write (6,"(i33,x,e15.8,x,3(f6.3,x))") f,mesh_ipArea(f,i,e),mesh_ipAreaNormal(:,f,i,e) + ! enddo + ! enddo + ! enddo + ! write (6,*) + ! write (6,*) "Input Parser: SUBNODE COORDINATES" + ! write (6,*) + ! write(6,'(a5,x,a5,x,a15,x,a15,x,a20,3(x,a8))') 'elem','IP','IP neighbor','IPFaceNodes','subNodeOnIPFace','x','y','z' + ! do e = 1,mesh_NcpElems ! loop over cpElems + ! t = mesh_element(2,e) ! get elemType + ! do i = 1,FE_Nips(t) ! loop over IPs of elem + ! do f = 1,FE_NipNeighbors(t) ! loop over interfaces of IP + ! do n = 1,FE_NipFaceNodes ! loop over nodes on interface + ! write(6,'(i5,x,i5,x,i15,x,i15,x,i20,3(x,f8.3))') e,i,f,n,FE_subNodeOnIPFace(n,f,i,t),& + ! mesh_subNodeCoord(1,FE_subNodeOnIPFace(n,f,i,t),e),& + ! mesh_subNodeCoord(2,FE_subNodeOnIPFace(n,f,i,t),e),& + ! mesh_subNodeCoord(3,FE_subNodeOnIPFace(n,f,i,t),e) + ! enddo + ! enddo + ! enddo + ! enddo + ! write (6,*) write (6,*) "Input Parser: STATISTICS" write (6,*) write (6,*) mesh_Nelems, " : total number of elements in mesh"