223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

added dislocation multiplication to dotState in constitutive_nonlocal.f90 

cleaned up debugging output statements to *.out file
This commit is contained in:
Christoph Kords 2009-08-12 11:22:02 +00:00
parent 6171361c7e
commit 2aae7b7574
4 changed files with 192 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

151
code/constitutive_nonlocal.f90
View File

@ -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

95
code/crystallite.f90
View File

@ -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 <<<updated 31.07.2009>>>
! 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 <<<updated 31.07.2009>>>
! 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

3
code/math.f90
View File

@ -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

89
code/mesh.f90
View File

@ -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"