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

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
This commit is contained in:
Christoph Kords 2009-08-24 08:16:01 +00:00
parent 387195e036
commit 1dbd0865db
3 changed files with 215 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

178
code/constitutive_nonlocal.f90
View File

@ -138,6 +138,7 @@ integer(pInt) section, &
s, & ! index of my slip system s, & ! index of my slip system
s1, & ! index of my slip system s1, & ! index of my slip system
s2, & ! index of my slip system s2, & ! index of my slip system
it, & ! index of my interaction type
output, & output, &
mySize mySize
character(len=64) tag character(len=64) tag
@ -250,7 +251,7 @@ do
case ('burgers') case ('burgers')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_burgersBySlipFamily(f,i) = IO_floatValue(line,positions,1+f) forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_burgersBySlipFamily(f,i) = IO_floatValue(line,positions,1+f)
case ('interaction_slipslip') 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 end select
endif endif
enddo enddo
@ -267,7 +268,6 @@ enddo
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)
do f = 1,lattice_maxNslipFamily do f = 1,lattice_maxNslipFamily
if (constitutive_nonlocal_Nslip(f,i) > 0) then if (constitutive_nonlocal_Nslip(f,i) > 0) then
if (constitutive_nonlocal_rhoEdgePos0(f,i) < 0.0_pReal) call IO_error(220) if (constitutive_nonlocal_rhoEdgePos0(f,i) < 0.0_pReal) call IO_error(220)
@ -279,6 +279,7 @@ enddo
if (constitutive_nonlocal_lambda0BySlipFamily(f,i) <= 0.0_pReal) call IO_error(-1) if (constitutive_nonlocal_lambda0BySlipFamily(f,i) <= 0.0_pReal) call IO_error(-1)
endif endif
enddo enddo
if (sum(constitutive_nonlocal_interactionSlipSlip(:,i)) <= 0) call IO_error(-1)
!*** determine total number of active slip systems !*** determine total number of active slip systems
@ -331,10 +332,15 @@ do i = 1,maxNinstance
do o = 1,maxval(phase_Noutput) do o = 1,maxval(phase_Noutput)
select case(constitutive_nonlocal_output(o,i)) select case(constitutive_nonlocal_output(o,i))
case( 'dislocationdensity', & case( 'rho', &
'shearrate_slip', & 'rho_edge', &
'resolvedstress_slip', & 'rho_screw', &
'resistance_slip') 'excess_rho_edge', &
'excess_rho_screw', &
'rho_forest', &
'shearrate', &
'resolvedstress', &
'resistance')
mySize = constitutive_nonlocal_totalNslip(i) mySize = constitutive_nonlocal_totalNslip(i)
case default case default
mySize = 0_pInt 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_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_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) constitutive_nonlocal_nu(i) = constitutive_nonlocal_C12(i) / constitutive_nonlocal_C11(i)
@ -559,7 +565,6 @@ use math, only: math_Plain3333to99, &
math_mul33x33, & math_mul33x33, &
math_mul3x3, & math_mul3x3, &
math_mul33x3, & math_mul33x3, &
math_transpose3x3, &
pi pi
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
mesh_maxNips, & 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 ! 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)) ) (mesh_ipCenterOfGravity(:,ip,el) - mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el)) )
neighboring_rhoEdgePos = state(1, neighboring_ip, neighboring_el)%p( 1: ns) 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 sigma(3,1) = 0.0_pReal
! coordinate transformation from the slip coordinate system to the lattice coordinate system ! 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
enddo enddo
@ -831,11 +836,10 @@ backStress_v = state(g,ip,el)%p(6*ns+1:6*ns+6)
!*** loop over slip systems !*** loop over slip systems
do s = 1,constitutive_nonlocal_totalNslip(myInstance) do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
!*** Calculation of Lp !*** Calculation of Lp
tauSlip(s) = math_mul6x6( Tstar_v + backStress_v, lattice_Sslip_v(:,sLattice,myStructure) ) 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) gdotSlip(s) = sum(rho(:,s)) * constitutive_nonlocal_burgersBySlipSystem(s,myInstance) * v(s)
Lp = Lp + gdotSlip(s) * lattice_Sslip(:,:,sLattice,myStructure) 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 !*** Calculation of the tangent of Lp
@ -863,9 +867,12 @@ enddo
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) 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,'(a23,3(i3,x),/,12(e10.3,x),/)') 'dislocation density at ',g,ip,el, rho
!if (debugger) write(6,'(a15,3(i3,x),/,12(f10.3,x),/)') 'tauSlip / MPa at ',g,ip,el, tauSlip/1e6 ! 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(e10.3,x),/)') 'gdotSlip at ',g,ip,el, gdotSlip ! 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 endsubroutine
@ -990,9 +997,11 @@ enddo
invLambda = sqrt(rhoForest) / constitutive_nonlocal_lambda0BySlipSystem(:,myInstance) invLambda = sqrt(rhoForest) / constitutive_nonlocal_lambda0BySlipSystem(:,myInstance)
forall (t = 1:4, s = 1:ns) & forall (t = 1:4) &
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) & 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(s,myInstance) / 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 * !* 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 prec, only: pReal, &
use mesh, only: mesh_NcpElems,mesh_maxNips pInt, &
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput p_vec
use lattice, only: lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & use math, only: math_mul6x6
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin 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 implicit none
!* Definition of variables !*** input variables
integer(pInt), intent(in) :: ipc,ip,el integer(pInt), intent(in) :: g, & ! current grain number
real(pReal), intent(in) :: dt,Temperature ip, & ! current integration point
real(pReal), dimension(6), intent(in) :: Tstar_v el ! current element number
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state real(pReal), intent(in) :: dt, & ! time increment
integer(pInt) matID,structID,ns,nt,f,o,i,c,j,index_myFamily Temperature ! temperature
real(pReal) sumf,tau real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & type(p_vec), dimension(homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems), intent(in) :: &
constitutive_nonlocal_postResults state ! microstructural state
!* Shortened notation !*** output variables
matID = phase_constitutionInstance(material_phase(ipc,ip,el)) real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
structID = constitutive_nonlocal_structure(matID) constitutive_nonlocal_postResults
ns = constitutive_nonlocal_totalNslip(matID)
!*** 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 c = 0_pInt
constitutive_nonlocal_postResults = 0.0_pReal constitutive_nonlocal_postResults = 0.0_pReal
do o = 1,phase_Noutput(material_phase(ipc,ip,el)) do o = 1,phase_Noutput(material_phase(g,ip,el))
select case(constitutive_nonlocal_output(o,matID)) select case(constitutive_nonlocal_output(o,myInstance))
case ('dislocationdensity') case ('rho')
constitutive_nonlocal_postResults(c+1:c+ns) = state(ipc,ip,el)%p(1:ns) 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 + 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 end select
enddo enddo

93
code/crystallite.f90
View File

@ -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_temperatureConverged: ', shape(crystallite_temperatureConverged)
write(6,'(a35,x,7(i5,x))') 'crystallite_nonfinished: ', shape(crystallite_nonfinished) write(6,'(a35,x,7(i5,x))') 'crystallite_nonfinished: ', shape(crystallite_nonfinished)
write(6,*) 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) call flush(6)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
@ -347,6 +347,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) 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 i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_converged(g,i,e)) then if (crystallite_converged(g,i,e)) then
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -424,7 +425,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) 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 i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains 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 enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
!$OMP 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 do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) 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 i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains 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)) & ! all undone crystallites if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites
crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e) crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e)
enddo enddo
enddo enddo
enddo enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -494,7 +496,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) 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 i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains 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 enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
!$OMP PARALLEL DO !$OMP PARALLEL DO
@ -502,6 +505,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) 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 i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites if (crystallite_nonfinished(g,i,e)) & ! all undone crystallites
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) 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_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state
crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature 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) 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 if (crystallite_converged(g,i,e)) then
!$OMP CRITICAL (distributionState) !$OMP CRITICAL (distributionState)
debug_CrystalliteStateLoopDistribution(NiterationState) = & debug_CrystalliteStateLoopDistribution(NiterationState) = &
@ -530,7 +533,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo enddo
!$OMPEND PARALLEL DO !$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 enddo ! crystallite convergence loop
@ -538,6 +542,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo ! cutback loop 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 ------ ! ------ check for non-converged crystallites ------
@ -614,7 +620,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
NiterationState = NiterationState + 1_pInt NiterationState = NiterationState + 1_pInt
onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe)
if (onTrack) then 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) crystallite_invFp(:,:,g,i,e), crystallite_Temperature(g,i,e), g, i, e)
stateConverged = crystallite_updateState(g,i,e) ! update state stateConverged = crystallite_updateState(g,i,e) ! update state
@ -645,8 +651,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_P(:,:,g,i,e) = myP crystallite_P(:,:,g,i,e) = myP
!$OMP CRITICAL (out) !$OMP CRITICAL (out)
debug_StiffnessStateLoopDistribution(NiterationState) = & debug_StiffnessStateLoopDistribution(NiterationState) = &
debug_StiffnessstateLoopDistribution(NiterationState) + 1 debug_StiffnessstateLoopDistribution(NiterationState) + 1
if (nState < NiterationState) write(6,*) 'ohh shit!! stiffenss state loop debugging exceeded',NiterationState
!$OMPEND CRITICAL (out) !$OMPEND CRITICAL (out)
enddo enddo
enddo enddo
@ -718,9 +723,13 @@ endsubroutine
! if NaN occured then return without changing the state ! if NaN occured then return without changing the state
if (any(residuum/=residuum)) then if (any(residuum/=residuum)) then
crystallite_updateState = .false. ! indicate state update failed crystallite_updateState = .false. ! indicate state update failed
!$OMP CRITICAL (write2out) crystallite_onTrack(g,i,e) = .false. ! no need to calculate any further
write(6,*) '::: updateState encountered NaN',e,i,g crystallite_onTrack = crystallite_onTrack .and. crystallite_localConstitution ! all nonlocal crystallites have to be redone
!$OMPEND CRITICAL (write2out) if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: updateState encountered NaN',g,i,e
!$OMPEND CRITICAL (write2out)
endif
return return
endif endif
@ -731,7 +740,19 @@ endsubroutine
crystallite_updateState = all(constitutive_state(g,i,e)%p(1:mySize) == 0.0_pReal .or. & 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))) abs(residuum) < rTol_crystalliteState*abs(constitutive_state(g,i,e)%p(1:mySize)))
if (debugger) then 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 endif
return return
@ -787,7 +808,7 @@ endsubroutine
if (residuum/=residuum) then if (residuum/=residuum) then
crystallite_updateTemperature = .false. ! indicate update failed crystallite_updateTemperature = .false. ! indicate update failed
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) '::: updateTemperature encountered NaN',e,i,g write(6,*) '::: updateTemperature encountered NaN',g,i,e
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
return return
endif endif
@ -795,7 +816,7 @@ endsubroutine
! update the microstructure ! update the microstructure
crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum 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. & crystallite_updateTemperature = crystallite_Temperature(g,i,e) == 0.0_pReal .or. &
abs(residuum) < rTol_crystalliteTemperature*crystallite_Temperature(g,i,e) abs(residuum) < rTol_crystalliteTemperature*crystallite_Temperature(g,i,e)
@ -910,9 +931,13 @@ endsubroutine
! inversion of Fp_current... ! inversion of Fp_current...
invFp_current = math_inv3x3(Fp_current) invFp_current = math_inv3x3(Fp_current)
if (all(invFp_current == 0.0_pReal)) then ! ... failed? if (all(invFp_current == 0.0_pReal)) then ! ... failed?
!$OMP CRITICAL (write2out) if (debugger) then
write(6,*) '::: integrateStress failed on invFp_current inversion',e,i,g !$OMP CRITICAL (write2out)
!$OMPEND 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 return
endif endif
@ -973,9 +998,11 @@ LpLoop: do
! NaN occured at regular speed? ! NaN occured at regular speed?
if (any(residuum/=residuum) .and. leapfrog == 1.0) then if (any(residuum/=residuum) .and. leapfrog == 1.0) then
!$OMP CRITICAL (write2out) if (debugger) then
write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',e,i,g !$OMP CRITICAL (write2out)
!$OMPEND CRITICAL (write2out) write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',g,i,e
!$OMPEND CRITICAL (write2out)
endif
return return
! something went wrong at accelerated speed? ! something went wrong at accelerated speed?
@ -1008,11 +1035,11 @@ LpLoop: do
if (error) then if (error) then
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress 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,*)
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,'(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 write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'Lpguess at ',g,i,e,Lpguess
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
return return
@ -1040,7 +1067,9 @@ LpLoop: do
if (error) then if (error) then
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$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) !$OMPEND CRITICAL (write2out)
endif endif
return return
@ -1068,12 +1097,12 @@ LpLoop: do
write(6,*) 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)/))') '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)/))') 'Lp',crystallite_Lp(:,:,g,i,e)
write(6,'(a,/,3(3(f12.7,x)/))') 'Fp',crystallite_Fp(:,:,g,i,e)
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
endif endif
!$OMP CRITICAL (distributionStress) !$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(NiterationStress) = debug_StressLoopDistribution(NiterationStress) + 1 debug_StressLoopDistribution(NiterationStress) = debug_StressLoopDistribution(NiterationStress) + 1
if (nStress < NiterationStress) write(6,*) 'ohh shit!! debug loop of stress exceeded',NiterationStress
!$OMPEND CRITICAL (distributionStress) !$OMPEND CRITICAL (distributionStress)
return return

21
code/homogenization.f90
View File

@ -213,7 +213,8 @@ subroutine materialpoint_stressAndItsTangent(&
crystallite_requested, & crystallite_requested, &
crystallite_converged, & crystallite_converged, &
crystallite_stressAndItsTangent crystallite_stressAndItsTangent
use debug, only: debug_MaterialpointLoopDistribution, & use debug, only: debugger, &
debug_MaterialpointLoopDistribution, &
debug_MaterialpointStateLoopDistribution debug_MaterialpointStateLoopDistribution
implicit none implicit none
@ -270,6 +271,8 @@ subroutine materialpoint_stressAndItsTangent(&
myNgrains = homogenization_Ngrains(mesh_element(3,e)) 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 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 our materialpoint converged then we are either finished or have to wind forward
if (materialpoint_converged(i,e)) then 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_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)) 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 ! still stepping needed
if (materialpoint_subStep(i,e) > subStepMin) then if (materialpoint_subStep(i,e) > subStepMin) then
@ -302,6 +313,13 @@ subroutine materialpoint_stressAndItsTangent(&
materialpoint_subStep(i,e) = 0.5_pReal * materialpoint_subStep(i,e) 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... ! restore...
crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures
crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads
@ -415,7 +433,6 @@ elementLoop: do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
write (6,*) 'Material Point finished' 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? ! how to deal with stiffness?
return return