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

* had to correct nonlocal slip system compatibility and flux calculation. last update gave wrong results.

This commit is contained in:
Christoph Kords 2010-10-15 13:19:26 +00:00
parent 61f8a5fcbe
commit fffe731447
3 changed files with 178 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
code/constitutive.f90
View File

@ -508,7 +508,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ip
endsubroutine endsubroutine
subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, ipc, ip, el) subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, orientation, ipc, ip, el)
!********************************************************************* !*********************************************************************
!* This subroutine contains the constitutive equation for * !* This subroutine contains the constitutive equation for *
!* calculating the rate of change of microstructure * !* calculating the rate of change of microstructure *
@ -550,6 +550,8 @@ real(pReal), intent(in) :: Temperature, &
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe, & Fe, &
Fp Fp
real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
orientation
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v, & Tstar_v, &
subTstar0_v subTstar0_v
@ -577,7 +579,8 @@ select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, & call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, &
constitutive_state, constitutive_subState0, constitutive_relevantState, subdt, ipc, ip, el) constitutive_state, constitutive_subState0, constitutive_relevantState, subdt, &
orientation, ipc, ip, el)
end select end select
@ -709,7 +712,7 @@ function constitutive_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, mis
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, misorientation, & constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, misorientation, &
dt, subdt, constitutive_state, constitutive_subState0, & dt, subdt, constitutive_state, constitutive_subState0, &
constitutive_dotstate, ipc, ip, el) constitutive_relevantState, constitutive_dotstate, ipc, ip, el)
end select end select
return return

203
code/constitutive_nonlocal.f90
View File

@ -81,8 +81,7 @@ real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_
constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance
real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_v, & ! dislocation velocity real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_v, & ! dislocation velocity
constitutive_nonlocal_rhoDotFlux ! dislocation convection term constitutive_nonlocal_rhoDotFlux ! dislocation convection term
real(pReal), dimension(:,:,:,:,:,:), allocatable :: constitutive_nonlocal_compatibility, & ! slip system compatibility between me and my neighbors real(pReal), dimension(:,:,:,:,:,:), allocatable :: constitutive_nonlocal_compatibility ! slip system compatibility between me and my neighbors
constitutive_nonlocal_transmissivity ! transmissivity between me and my neighbors
real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance
constitutive_nonlocal_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance constitutive_nonlocal_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance
constitutive_nonlocal_interactionMatrixSlipSlip ! interaction matrix of the different slip systems for each instance constitutive_nonlocal_interactionMatrixSlipSlip ! interaction matrix of the different slip systems for each instance
@ -425,12 +424,9 @@ constitutive_nonlocal_v = 0.0_pReal
allocate(constitutive_nonlocal_rhoDotFlux(maxTotalNslip, 8, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems)) allocate(constitutive_nonlocal_rhoDotFlux(maxTotalNslip, 8, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
constitutive_nonlocal_rhoDotFlux = 0.0_pReal constitutive_nonlocal_rhoDotFlux = 0.0_pReal
allocate(constitutive_nonlocal_compatibility(maxTotalNslip, 2, maxTotalNslip, FE_maxNipNeighbors, mesh_maxNips, mesh_NcpElems)) allocate(constitutive_nonlocal_compatibility(2,maxTotalNslip, maxTotalNslip, FE_maxNipNeighbors, mesh_maxNips, mesh_NcpElems))
constitutive_nonlocal_compatibility = 0.0_pReal constitutive_nonlocal_compatibility = 0.0_pReal
allocate(constitutive_nonlocal_transmissivity(maxTotalNslip, 2, maxTotalNslip, FE_maxNipNeighbors, mesh_maxNips, mesh_NcpElems))
constitutive_nonlocal_transmissivity = 0.0_pReal
do i = 1,maxNinstance do i = 1,maxNinstance
myStructure = constitutive_nonlocal_structure(i) ! lattice structure of this instance myStructure = constitutive_nonlocal_structure(i) ! lattice structure of this instance
@ -1224,7 +1220,7 @@ endsubroutine
!* rate of change of microstructure * !* rate of change of microstructure *
!********************************************************************* !*********************************************************************
subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, previousTstar_v, Fe, Fp, Temperature, dt_previous, & subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, previousTstar_v, Fe, Fp, Temperature, dt_previous, &
state, previousState, relevantState, timestep, g,ip,el) state, previousState, relevantState, timestep, orientation, g,ip,el)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
@ -1241,6 +1237,8 @@ use math, only: math_norm3, &
math_inv3x3, & math_inv3x3, &
math_det3x3, & math_det3x3, &
math_Mandel6to33, & math_Mandel6to33, &
math_QuaternionDisorientation, &
math_qRot, &
pi, & pi, &
NaN NaN
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
@ -1279,6 +1277,8 @@ real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe, & ! elastic deformation gradient Fe, & ! elastic deformation gradient
Fp ! plastic deformation gradient Fp ! plastic deformation gradient
real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
orientation ! crystal lattice orientation
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state, & ! current microstructural state state, & ! current microstructural state
previousState, & ! previous microstructural state previousState, & ! previous microstructural state
@ -1340,18 +1340,20 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan
previousDUpper, & ! previous maximum stable dipole distance for edges and screws previousDUpper, & ! previous maximum stable dipole distance for edges and screws
dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
m ! direction of dislocation motion m, & ! direction of dislocation motion
neighboring_m ! direction of dislocation motion for my neighbor in central MP's lattice configuration
real(pReal), dimension(3,3) :: F, & ! total deformation gradient real(pReal), dimension(3,3) :: F, & ! total deformation gradient
neighboring_F, & ! total deformation gradient of my neighbor neighboring_F, & ! total deformation gradient of my neighbor
Favg ! average total deformation gradient of me and my neighbor Favg ! average total deformation gradient of me and my neighbor
real(pReal), dimension(6) :: Tdislocation_v, & ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress real(pReal), dimension(6) :: Tdislocation_v, & ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
previousTdislocation_v ! previous dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress previousTdislocation_v ! previous dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor
real(pReal), dimension(3) :: surfaceNormal, & ! surface normal in lattice configuration real(pReal), dimension(3) :: surfaceNormal, & ! surface normal in lattice configuration
surfaceNormal_currentconf ! surface normal in current configuration surfaceNormal_currentconf ! surface normal in current configuration
real(pReal) area, & ! area of the current interface real(pReal) area, & ! area of the current interface
detFe, & ! determinant of elastic defornmation gradient detFe, & ! determinant of elastic defornmation gradient
transmissivity, & ! transmissivity of dislocation flux for different slip systems in neighboring material points for a specific charcter of dislocations transmissivity, & ! transmissivity of dislocation flux for different slip systems in neighboring material points
compatibility, & ! compatibility of different slip systems in neighboring material points for a specific charcter of dislocations compatibility, & ! compatibility of different slip systems in neighboring material points for a specific character of dislocations
lineLength, & ! dislocation line length leaving the current interface lineLength, & ! dislocation line length leaving the current interface
D ! self diffusion D ! self diffusion
logical, dimension(3) :: periodicSurfaceFlux ! flag indicating periodic fluxes at surfaces when surface normal points mainly in x, y and z direction respectively (in reference configuration) logical, dimension(3) :: periodicSurfaceFlux ! flag indicating periodic fluxes at surfaces when surface normal points mainly in x, y and z direction respectively (in reference configuration)
@ -1482,7 +1484,7 @@ rhoDotMultiplication(:,5:10) = 0.0_pReal ! used dislocation densities and d
!*** calculate dislocation fluxes !*** calculate dislocation fluxes
rhoDotFlux = 0.0_pReal rhoDotFlux = 0.0_pReal
periodicSurfaceFlux = (/.false.,.false.,.false./) periodicSurfaceFlux = (/.false.,.true.,.false./)
m(:,:,1) = lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) m(:,:,1) = lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure)
m(:,:,2) = -lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) m(:,:,2) = -lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure)
@ -1494,6 +1496,7 @@ detFe = math_det3x3(Fe(:,:,g,ip,el))
fluxdensity = rhoSgl(:,1:4) * constitutive_nonlocal_v(:,:,g,ip,el) fluxdensity = rhoSgl(:,1:4) * constitutive_nonlocal_v(:,:,g,ip,el)
!if (selectiveDebugger) write(6,*) '--> dislocation flux <---'
do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors
opposite_n = n - 1_pInt + 2_pInt*mod(n,2_pInt) opposite_n = n - 1_pInt + 2_pInt*mod(n,2_pInt)
@ -1505,8 +1508,8 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
if ( neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt ) then ! if neighbor exists ... if ( neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt ) then ! if neighbor exists ...
do neighboring_n = 1,FE_NipNeighbors(mesh_element(2,neighboring_el)) ! find neighboring index that points from my neighbor to myself do neighboring_n = 1,FE_NipNeighbors(mesh_element(2,neighboring_el)) ! find neighboring index that points from my neighbor to myself
if ( el == mesh_ipNeighborhood(1,opposite_n,ip,neighboring_el) & ! special case if no neighbor at all... if ( el == mesh_ipNeighborhood(1,neighboring_n,neighboring_ip,neighboring_el) & ! special case if no neighbor at all...
.and. ip == mesh_ipNeighborhood(2,opposite_n,ip,neighboring_el) ) & .and. ip == mesh_ipNeighborhood(2,neighboring_n,neighboring_ip,neighboring_el) ) &
exit ! ...exit without any flux calculation exit ! ...exit without any flux calculation
enddo enddo
endif endif
@ -1528,38 +1531,69 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
forall (t = 1:4) & ! ... then calculate neighboring flux density forall (t = 1:4) & ! ... then calculate neighboring flux density
neighboring_fluxdensity(:,t) = state(g,neighboring_ip,neighboring_el)%p((t-1)*ns+1:t*ns) & neighboring_fluxdensity(:,t) = state(g,neighboring_ip,neighboring_el)%p((t-1)*ns+1:t*ns) &
* constitutive_nonlocal_v(:,t,g,neighboring_ip,neighboring_el) * constitutive_nonlocal_v(:,t,g,neighboring_ip,neighboring_el)
absoluteMisorientation = math_QuaternionDisorientation( orientation(:,1,ip,el), &
orientation(:,1,neighboring_ip,neighboring_el), 0_pInt)
else ! ... and is of local constitution... else ! ... and is of local constitution...
neighboring_fluxdensity = fluxdensity ! ... then copy flux density to neighbor to ensure zero gradient in fluxdensity neighboring_fluxdensity = fluxdensity ! ... then copy flux density to neighbor to ensure zero gradient in fluxdensity
absoluteMisorientation = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/)
endif endif
else ! if no neighbor existent... else ! if no neighbor existent...
if ( all(periodicSurfaceFlux(maxloc(abs(mesh_ipAreaNormal(:,n,ip,el))))) ) then ! ... and we want periodic fluxes at surface... if ( all(periodicSurfaceFlux(maxloc(abs(mesh_ipAreaNormal(:,n,ip,el))))) ) then ! ... and we want periodic fluxes at surface...
forall (t = 1:4) & ! ... then mirror fluxes forall (t = 1:4) & ! ... then mirror fluxes
neighboring_fluxdensity(:,t) = fluxdensity(:,t-1_pInt+2_pInt*mod(t,2_pInt)) neighboring_fluxdensity(:,t) = fluxdensity(:,t-1+2*mod(t,2))
else ! ... and we have a free surface... else ! ... and we have a free surface...
neighboring_fluxdensity = 0.0_pReal ! ... assume zero density neighboring_fluxdensity = 0.0_pReal ! ... assume zero density
endif endif
absoluteMisorientation = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/)
endif endif
do t = 1,4
do s = 1,ns
neighboring_m(:,s,t) = math_qRot(absoluteMisorientation, m(:,s,t)) ! calculate neighboring dislocation movement directions in my lattice frame (we simply assume same crystal structure for both me and my neighbor)
enddo
enddo
! if (selectiveDebugger) then
! write(6,'(a,x,i1)') 'neighbor',n
! write(6,'(a,x,i2)') 'neighboring_ip',neighboring_ip
! write(6,'(a,x,i1)') 'neighboring_n',neighboring_n
! write(6,'(a,12(x,f5.3))') 'compatibility of neighbor with me', &
! constitutive_nonlocal_compatibility(:,:,:,neighboring_n,neighboring_ip,neighboring_el)
! endif
do s = 1,ns do s = 1,ns
! if (selectiveDebugger) write(6,'(a,x,i2)') ' system',s
do t = 1,4 do t = 1,4
c = (t + 1) / 2 ! dislocation character ! if (selectiveDebugger) write(6,'(a,x,i1)') ' type',t
if ( fluxdensity(s,t) * math_mul3x3(m(:,s,t), surfaceNormal) > 0.0_pReal ) then ! outgoing flux if ( fluxdensity(s,t) * math_mul3x3(m(:,s,t),surfaceNormal) > 0.0_pReal ) then ! outgoing flux
transmissivity = sum(constitutive_nonlocal_transmissivity(:,c,s,n,ip,el)) ! overall transmissivity between my system s and all neighboring systems s2 for this dislocation character transmissivity = sum(abs(constitutive_nonlocal_compatibility(1,:,s,n,ip,el))) ! ..overall transmissivity between my system s and all neighboring systems s2 for this dislocation character
lineLength = fluxdensity(s,t) * math_mul3x3(m(:,s,t), surfaceNormal) * area ! line length that wants to leave this interface lineLength = fluxdensity(s,t) * math_mul3x3(m(:,s,t),surfaceNormal) * area ! ..line length that wants to leave thrugh this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract positive dislocation flux that leaves the material point rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! ..subtract positive dislocation flux that leaves the material point
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) & rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) &
* sign(1.0_pReal, fluxdensity(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point * sign(1.0_pReal, fluxdensity(s,t)) ! ..dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
else ! incoming flux ! if (selectiveDebugger) write(6,'(a,x,e12.5,x,a,x,f10.5)') ' outgoing flux:', lineLength / mesh_ipVolume(ip,el),&
do s2 = 1,ns ! ' transmissivity:',transmissivity
compatibility = constitutive_nonlocal_compatibility(s,c,s2,neighboring_n,neighboring_ip,neighboring_el) ! compatibility of system s2 of my neighbor to system s in my material point else ! perhaps we get something from our neighbor?
transmissivity = constitutive_nonlocal_transmissivity(s,c,s2,neighboring_n,neighboring_ip,neighboring_el) c = (t + 1) / 2
lineLength = neighboring_fluxdensity(s2,t) * math_mul3x3(m(:,s,t), surfaceNormal) * area ! line length that wants to leave this interface topp = t + mod(t,2) - mod(t+1,2)
if (compatibility > 0.0_pReal) then ! compatible with equally signed dislocation density ! if (selectiveDebugger) write(6,'(a,12(x,f5.3))') ' compatibility', &
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) * transmissivity * compatibility ! subtract negative dislocation flux that enters the material point ! constitutive_nonlocal_compatibility(c,s,:,neighboring_n,neighboring_ip,neighboring_el)
elseif (compatibility < 0.0_pReal) then ! compatible with inversely signed dislocation density do s2 = 1,ns ! assuming same crystal structure at neighbor
topp = t + mod(t,2) - mod(t+1,2) compatibility = constitutive_nonlocal_compatibility(c,s,s2,neighboring_n,neighboring_ip,neighboring_el) ! ..compatibility of system s2 of my neighbor to system s in my material point
rhoDotFlux(s,topp) = rhoDotFlux(s,topp) - lineLength / mesh_ipVolume(ip,el) * transmissivity * abs(compatibility) ! subtract negative dislocation flux that enters the material point transmissivity = abs(constitutive_nonlocal_compatibility(1,s,s2,neighboring_n,neighboring_ip,neighboring_el))
if ( compatibility > 0.0_pReal ) then ! ..dislocation signs have same sense on neighboring system
if (neighboring_fluxdensity(s2,t) * math_mul3x3(neighboring_m(:,s2,t),surfaceNormal) < 0.0_pReal) then ! ....incoming flux
lineLength = neighboring_fluxdensity(s2,t) * math_mul3x3(neighboring_m(:,s2,t), surfaceNormal) * area ! ......line length that enters through this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) * transmissivity * abs(compatibility) ! ......subtract negative dislocation flux that enters the material point
endif
elseif ( compatibility < 0.0_pReal ) then ! ..dislocation signs have opposite sense on neighboring system, so consider opposite dislocation type
if (neighboring_fluxdensity(s2,topp) * math_mul3x3(neighboring_m(:,s2,topp),surfaceNormal) < 0.0_pReal) then ! ....incoming flux
lineLength = neighboring_fluxdensity(s2,topp) * math_mul3x3(neighboring_m(:,s2,topp), surfaceNormal) * area ! ......line length that enters through this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) * transmissivity * abs(compatibility) ! ......subtract negative dislocation flux that enters the material point
endif
endif endif
! if (selectiveDebugger) write(6,'(a,x,e12.5,x,a,x,f10.5)') ' incoming flux:', &
! - lineLength / mesh_ipVolume(ip,el) * transmissivity * abs(compatibility), &
! ' compatibility:',compatibility
enddo enddo
endif endif
enddo enddo
@ -1680,7 +1714,17 @@ do i = 1,10*ns
if (verboseDebugger) then if (verboseDebugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'negative dislocation density: enforced cutback at ',g,ip,el,i write(6,*) 'negative dislocation density: enforced cutback at ',g,ip,el,i
write(6,'(a12,x,e15.8)') 'dotState was',dotState(g,ip,el)%p(i) write(6,'(a,/,8(12(e12.5,x),/))') 'dislocation remobilization', rhoDotRemobilization(:,1:8) * timestep
write(6,'(a,/,4(12(e12.5,x),/))') 'dislocation multiplication', rhoDotMultiplication(:,1:4) * timestep
write(6,'(a,/,8(12(e12.5,x),/))') 'dislocation flux', rhoDotFlux(:,1:8) * timestep
write(6,'(a,/,10(12(e12.5,x),/))') 'dipole formation by glide', rhoDotSingle2DipoleGlide * timestep
write(6,'(a,/,2(12(e12.5,x),/))') 'athermal dipole annihilation', rhoDotAthermalAnnihilation(:,1:2) * timestep
write(6,'(a,/,2(12(e12.5,x),/))') 'thermally activated dipole annihilation', rhoDotThermalAnnihilation(:,9:10) * timestep
! write(6,'(a,/,10(12(e12.5,x),/))') 'dipole dissociation by stress increase', rhoDotDipole2SingleStressChange * timestep
! write(6,'(a,/,10(12(e12.5,x),/))') 'dipole formation by stress decrease', rhoDotSingle2DipoleStressChange * timestep
write(6,'(a,/,10(12(e12.5,x),/))') 'total density change', rhoDot * timestep
write(6,'(a,/,10(12(f12.7,x),/))') 'relative density change', rhoDot(:,1:8) * timestep / (abs(rhoSgl)+1.0_pReal), &
rhoDot(:,9:10) * timestep / (rhoDip+1.0_pReal)
write(6,*) write(6,*)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
@ -1700,6 +1744,8 @@ if (verboseDebugger .and. selectiveDebugger) then
! write(6,'(a,/,10(12(e12.5,x),/))') 'dipole dissociation by stress increase', rhoDotDipole2SingleStressChange * timestep ! write(6,'(a,/,10(12(e12.5,x),/))') 'dipole dissociation by stress increase', rhoDotDipole2SingleStressChange * timestep
! write(6,'(a,/,10(12(e12.5,x),/))') 'dipole formation by stress decrease', rhoDotSingle2DipoleStressChange * timestep ! write(6,'(a,/,10(12(e12.5,x),/))') 'dipole formation by stress decrease', rhoDotSingle2DipoleStressChange * timestep
write(6,'(a,/,10(12(e12.5,x),/))') 'total density change', rhoDot * timestep write(6,'(a,/,10(12(e12.5,x),/))') 'total density change', rhoDot * timestep
write(6,'(a,/,10(12(f12.7,x),/))') 'relative density change', rhoDot(:,1:8) * timestep / (abs(rhoSgl)+1.0_pReal), &
rhoDot(:,9:10) * timestep / (rhoDip+1.0_pReal)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
@ -1708,12 +1754,12 @@ endsubroutine
!********************************************************************* !*********************************************************************
!* TRANSMISSIVITY AND COMPATIBILITY UPDATE * !* COMPATIBILITY UPDATE *
!* Transmissivity is defined as absolute minimum of the cosine of * !* Compatibility is defined as normalized product of signed cosine *
!* the angle between the slip directions and the cosine of the angle * !* of the angle between the slip plane normals and signed cosine of *
!* between the slip plane normals. Compatibility is defined as * !* the angle between the slip directions. Only the largest values *
!* product of cosine of the angle between the slip plane normals and * !* that sum up to a total of 1 are considered, all others are set to *
!* signed cosine of the angle between the slip directions * !* zero. *
!********************************************************************* !*********************************************************************
subroutine constitutive_nonlocal_updateCompatibility(orientation,i,e) subroutine constitutive_nonlocal_updateCompatibility(orientation,i,e)
@ -1734,8 +1780,7 @@ use mesh, only: mesh_element, &
mesh_NcpElems mesh_NcpElems
use lattice, only: lattice_sn, & use lattice, only: lattice_sn, &
lattice_sd, & lattice_sd, &
lattice_st, & lattice_st
lattice_maxNslip
use debug, only: debugger, & use debug, only: debugger, &
debug_e, debug_i, debug_g, & debug_e, debug_i, debug_g, &
verboseDebugger, & verboseDebugger, &
@ -1763,31 +1808,36 @@ integer(pInt) n, & !
neighboringStructure, & neighboringStructure, &
myNSlipSystems, & ! number of active slip systems myNSlipSystems, & ! number of active slip systems
neighboringNSlipSystems, & neighboringNSlipSystems, &
c, & ! dislocation character index (1=edge, 2=screw)
s1, & ! slip system index (me) s1, & ! slip system index (me)
s2 ! slip system index (my neighbor) s2 ! slip system index (my neighbor)
integer(pInt), dimension(lattice_maxNslip) :: mySlipSystems, & ! slip system numbering according to lattice integer(pInt), dimension(maxval(constitutive_nonlocal_totalNslip)) :: &
mySlipSystems, & ! slip system numbering according to lattice
neighboringSlipSystems neighboringSlipSystems
real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor
real(pReal), dimension(3,lattice_maxNslip) :: myNormals, & ! slip plane normals real(pReal), dimension(3,maxval(constitutive_nonlocal_totalNslip)) :: &
neighboringNormals myNormals, & ! slip plane normals
real(pReal), dimension(3,lattice_maxNslip,2) :: mySlipDirections, & ! slip directions for positive edge and screws neighboringNormals, &
mySlipDirections, & ! slip directions
neighboringSlipDirections neighboringSlipDirections
real(pReal) compatibilitySum, &
compatibilityMax, &
compatibilityMaxCount
logical, dimension(maxval(constitutive_nonlocal_totalNslip)) :: &
compatibilityMask
selectiveDebugger = (debug_i==i .and. debug_e==e)
myPhase = material_phase(1,i,e) myPhase = material_phase(1,i,e)
myInstance = phase_constitutionInstance(myPhase) myInstance = phase_constitutionInstance(myPhase)
myStructure = constitutive_nonlocal_structure(myInstance) myStructure = constitutive_nonlocal_structure(myInstance)
myNSlipSystems = constitutive_nonlocal_totalNslip(myInstance) myNSlipSystems = constitutive_nonlocal_totalNslip(myInstance)
mySlipSystems(1:myNSlipSystems) = constitutive_nonlocal_slipSystemLattice(1:myNSlipSystems,myInstance) mySlipSystems(1:myNSlipSystems) = constitutive_nonlocal_slipSystemLattice(1:myNSlipSystems,myInstance)
myNormals = lattice_sn(:, mySlipSystems, myStructure) myNormals = lattice_sn(:, mySlipSystems, myStructure)
mySlipDirections(:,:,1) = lattice_sd(:, mySlipSystems, myStructure) ! direction of positive edges mySlipDirections = lattice_sd(:, mySlipSystems, myStructure)
mySlipDirections(:,:,2) = lattice_st(:, mySlipSystems, myStructure) ! direction of positive screws
do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors
neighboring_e = mesh_ipNeighborhood(1,n,i,e) neighboring_e = mesh_ipNeighborhood(1,n,i,e)
neighboring_i = mesh_ipNeighborhood(2,n,i,e) neighboring_i = mesh_ipNeighborhood(2,n,i,e)
if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists
neighboringPhase = material_phase(1,neighboring_i,neighboring_e) neighboringPhase = material_phase(1,neighboring_i,neighboring_e)
@ -1798,42 +1848,48 @@ do n = 1,FE_NipNeighbors(mesh_element(2,e))
neighboringSlipSystems(1:neighboringNSlipSystems) = constitutive_nonlocal_slipSystemLattice(1:neighboringNSlipSystems,& neighboringSlipSystems(1:neighboringNSlipSystems) = constitutive_nonlocal_slipSystemLattice(1:neighboringNSlipSystems,&
neighboringInstance) neighboringInstance)
neighboringNormals = lattice_sn(:, neighboringSlipSystems, neighboringStructure) neighboringNormals = lattice_sn(:, neighboringSlipSystems, neighboringStructure)
neighboringSlipDirections(:,:,1) = lattice_sd(:, neighboringSlipSystems, neighboringStructure) ! direction of positive edges neighboringSlipDirections = lattice_sd(:, neighboringSlipSystems, neighboringStructure)
neighboringSlipDirections(:,:,2) = lattice_st(:, neighboringSlipSystems, neighboringStructure) ! direction of positive screws
if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me
absoluteMisorientation = math_QuaternionDisorientation( orientation(:,1,i,e), & absoluteMisorientation = math_QuaternionDisorientation( orientation(:,1,i,e), &
orientation(:,1,neighboring_i,neighboring_e), 0_pInt) orientation(:,1,neighboring_i,neighboring_e), 0_pInt)
do s1 = 1,myNSlipSystems ! loop through my slip systems do s1 = 1,myNSlipSystems ! loop through my slip systems
do c = 1,2 ! loop through edge and screw character do s2 = 1,neighboringNSlipSystems ! loop through my neighbors' slip systems
do s2 = 1,neighboringNSlipSystems ! loop through my neighbors' slip systems constitutive_nonlocal_compatibility(1,s2,s1,n,i,e) = &
constitutive_nonlocal_transmissivity(s2,c,s1,n,i,e) = & math_mul3x3(myNormals(:,s1), math_qRot(absoluteMisorientation, neighboringNormals(:,s2))) &
min(abs(math_mul3x3(mySlipDirections(:,s1,c), & * abs(math_mul3x3(mySlipDirections(:,s1), math_qRot(absoluteMisorientation, neighboringSlipDirections(:,s2))))
math_qRot(absoluteMisorientation, neighboringSlipDirections(:,s2,c)))), & constitutive_nonlocal_compatibility(2,s2,s1,n,i,e) = &
abs(math_mul3x3(myNormals(:,s1), math_qRot(absoluteMisorientation, neighboringNormals(:,s2)))) ) math_mul3x3(myNormals(:,s1), math_qRot(absoluteMisorientation, neighboringNormals(:,s2))) &
constitutive_nonlocal_compatibility(s2,c,s1,n,i,e) = & * math_mul3x3(mySlipDirections(:,s1), math_qRot(absoluteMisorientation, neighboringSlipDirections(:,s2)))
abs(math_mul3x3(myNormals(:,s1), math_qRot(absoluteMisorientation, neighboringNormals(:,s2)))) &
* math_mul3x3(mySlipDirections(:,s1,c), math_qRot(absoluteMisorientation, neighboringSlipDirections(:,s2,c)))
enddo
if (any(abs(constitutive_nonlocal_compatibility(:,c,s1,n,i,e)) > 0.0_pReal)) then
constitutive_nonlocal_compatibility(:,c,s1,n,i,e) = constitutive_nonlocal_compatibility(:,c,s1,n,i,e) &
/ sum(abs(constitutive_nonlocal_compatibility(:,c,s1,n,i,e))) ! normalize to a total of one
endif
enddo enddo
compatibilitySum = 0.0_pReal
compatibilityMask = .true.
do while ( (1.0_pReal - compatibilitySum > 0.0_pReal) .and. any(compatibilityMask) ) ! only those largest values that sum up to 1 are considered (round off of the smallest considered values to ensure sum to be exactly 1)
compatibilityMax = maxval(abs(constitutive_nonlocal_compatibility(1,:,s1,n,i,e)), compatibilityMask)
compatibilityMaxCount = dble(count(abs(constitutive_nonlocal_compatibility(1,:,s1,n,i,e)) == compatibilityMax))
where (abs(constitutive_nonlocal_compatibility(1,:,s1,n,i,e)) >= compatibilityMax) compatibilityMask = .false.
if (compatibilitySum + compatibilityMax * compatibilityMaxCount > 1.0_pReal) &
where (abs(constitutive_nonlocal_compatibility(:,:,s1,n,i,e)) == compatibilityMax) &
constitutive_nonlocal_compatibility(:,:,s1,n,i,e) = sign((1.0_pReal - compatibilitySum) / compatibilityMaxCount, &
constitutive_nonlocal_compatibility(:,:,s1,n,i,e))
compatibilitySum = compatibilitySum + compatibilityMaxCount * compatibilityMax
enddo
where (compatibilityMask) constitutive_nonlocal_compatibility(1,:,s1,n,i,e) = 0.0_pReal
where (compatibilityMask) constitutive_nonlocal_compatibility(2,:,s1,n,i,e) = 0.0_pReal
enddo enddo
else ! neighbor has different crystal structure else ! neighbor has different crystal structure
constitutive_nonlocal_transmissivity(:,:,:,n,i,e) = 0.0_pReal ! no transmissivity... constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 0.0_pReal ! no compatibility
constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 0.0_pReal ! ...and compatibility
endif endif
else ! neighbor has local constitution else ! neighbor has local constitution
constitutive_nonlocal_transmissivity(:,:,:,n,i,e) = 1.0_pReal ! assume perfectly transmissive... constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 0.0_pReal
constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 1.0_pReal ! ...and compatible forall(s1 = 1:maxval(constitutive_nonlocal_totalNslip)) &
constitutive_nonlocal_compatibility(:,s1,s1,n,i,e) = 1.0_pReal ! assume perfect compatibility for equal slip system index
endif endif
else ! no neighbor present else ! no neighbor present
constitutive_nonlocal_transmissivity(:,:,:,n,i,e) = 1.0_pReal ! perfectly transmissive... constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 0.0_pReal
constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 1.0_pReal ! ...and compatible forall(s1 = 1:maxval(constitutive_nonlocal_totalNslip)) &
constitutive_nonlocal_compatibility(:,s1,s1,n,i,e) = 1.0_pReal ! perfect compatibility for equal slip system index
endif endif
enddo enddo
@ -1879,7 +1935,7 @@ endfunction
!* return array of constitutive results * !* return array of constitutive results *
!********************************************************************* !*********************************************************************
function constitutive_nonlocal_postResults(Tstar_v, previousTstar_v, Fe, Fp, Temperature, disorientation, dt, dt_previous, & function constitutive_nonlocal_postResults(Tstar_v, previousTstar_v, Fe, Fp, Temperature, disorientation, dt, dt_previous, &
state, previousState, dotState, g,ip,el) state, previousState, relevantState, dotState, g,ip,el)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
@ -1932,6 +1988,7 @@ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state, & ! current microstructural state state, & ! current microstructural state
previousState, & ! previous microstructural state previousState, & ! previous microstructural state
relevantState, &
dotState ! evolution rate of microstructural state dotState ! evolution rate of microstructural state
!*** output variables !*** output variables
@ -2023,8 +2080,8 @@ do t = 1,4
enddo enddo
enddo enddo
forall (t = 1:4) & forall (s = 1:ns, t = 1:4, rhoSgl(s,t) > relevantState(g,ip,el)%p((t-1)*ns+s)) & ! no shear rate contribution for densities below relevant state
gdot(:,t) = rhoSgl(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * constitutive_nonlocal_v(:,t,g,ip,el) gdot(s,t) = rhoSgl(s,t) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * constitutive_nonlocal_v(s,t,g,ip,el)
!* calculate limits for stable dipole height and its rate of change !* calculate limits for stable dipole height and its rate of change

63
code/crystallite.f90
View File

@ -905,6 +905,7 @@ endif
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, &
crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
@ -918,9 +919,10 @@ RK4dotTemperature = 0.0_pReal
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e) crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), &
crystallite_orientation, g,i,e)
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), &
crystallite_Temperature(g,i,e),g,i,e) crystallite_Temperature(g,i,e),g,i,e)
if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
@ -966,6 +968,7 @@ do n = 1,4
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), &
crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
@ -978,7 +981,7 @@ do n = 1,4
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
if (crystallite_integrateStress(mode,g,i,e,timeStepFraction(n))) then ! fraction of original times step if (crystallite_integrateStress(mode,g,i,e,timeStepFraction(n))) then ! fraction of original times step
if (n == 4) then ! final integration step if (n == 4) then ! final integration step
if (mode==1 .and. verboseDebugger .and. e == debug_e .and. i == debug_i .and. g == debug_g) then if (mode==1 .and. verboseDebugger .and. e == debug_e .and. i == debug_i .and. g == debug_g) then
@ -986,9 +989,9 @@ do n = 1,4
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) '::: updateState',g,i,e write(6,*) '::: updateState',g,i,e
write(6,*) write(6,*)
write(6,'(a,/,12(e14.8,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) write(6,'(a,/,12(e12.5,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,*) write(6,*)
write(6,'(a,/,12(e14.8,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:mySizeDotState) write(6,'(a,/,12(e12.5,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*) write(6,*)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
@ -1018,10 +1021,11 @@ do n = 1,4
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), &
timeStepFraction(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
crystallite_orientation, g,i,e)
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), &
crystallite_Temperature(g,i,e),g,i,e) crystallite_Temperature(g,i,e),g,i,e)
if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
@ -1212,7 +1216,8 @@ endif
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1)
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e) crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), &
crystallite_orientation, g,i,e)
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), &
crystallite_Temperature(g,i,e),g,i,e) crystallite_Temperature(g,i,e),g,i,e)
if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
@ -1299,7 +1304,8 @@ do n = 1,5
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), &
c(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep c(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
crystallite_orientation, g,i,e)
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), &
crystallite_Temperature(g,i,e),g,i,e) crystallite_Temperature(g,i,e),g,i,e)
if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
@ -1379,9 +1385,9 @@ relTemperatureResiduum = 0.0_pReal
! write(6,'(12(e14.8,x))') constitutive_RKCK45dotState(j,g,i,e)%p(1:sizeDotState) ! write(6,'(12(e14.8,x))') constitutive_RKCK45dotState(j,g,i,e)%p(1:sizeDotState)
! write(6,*) ! write(6,*)
! enddo ! enddo
write(6,'(a,/,12(e14.8,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:sizeDotState) write(6,'(a,/,12(e12.5,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:sizeDotState)
write(6,*) write(6,*)
write(6,'(a,/,12(e14.8,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:sizeDotState) write(6,'(a,/,12(e12.5,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:sizeDotState)
write(6,*) write(6,*)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
@ -1555,7 +1561,8 @@ endif
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1)
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e) crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), &
crystallite_orientation, g,i,e)
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), &
crystallite_Temperature(g,i,e),g,i,e) crystallite_Temperature(g,i,e),g,i,e)
@ -1631,7 +1638,8 @@ endif
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1)
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e) crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), &
crystallite_orientation, g,i,e)
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), &
crystallite_Temperature(g,i,e),g,i,e) crystallite_Temperature(g,i,e),g,i,e)
if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
@ -1684,10 +1692,10 @@ relTemperatureResiduum = 0.0_pReal
write(6,*) write(6,*)
write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance', relStateResiduum(1:sizeDotState,g,i,e) write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance', relStateResiduum(1:sizeDotState,g,i,e)
write(6,*) write(6,*)
write(6,'(a,/,12(e14.8,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:sizeDotState) & write(6,'(a,/,12(e12.5,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:sizeDotState) &
- 2.0_pReal * stateResiduum(1:sizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum - 2.0_pReal * stateResiduum(1:sizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum
write(6,*) write(6,*)
write(6,'(a,/,12(e14.8,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:sizeDotState) write(6,'(a,/,12(e12.5,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:sizeDotState)
write(6,*) write(6,*)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
@ -1794,6 +1802,7 @@ endif
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), &
crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
@ -1806,9 +1815,10 @@ endif
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e) crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), &
crystallite_orientation, g,i,e)
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), & crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e), &
crystallite_Temperature(g,i,e),g,i,e) crystallite_Temperature(g,i,e),g,i,e)
if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
@ -1831,7 +1841,7 @@ endif
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
mySizeDotState = constitutive_sizeDotState(g,i,e) mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) & constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
+ constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
@ -1842,9 +1852,9 @@ endif
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) '::: updateState',g,i,e write(6,*) '::: updateState',g,i,e
write(6,*) write(6,*)
write(6,'(a,/,12(e14.8,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) write(6,'(a,/,12(e12.5,x))') 'updateState: dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
write(6,*) write(6,*)
write(6,'(a,/,12(e14.8,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:mySizeDotState) write(6,'(a,/,12(e12.5,x))') 'updateState: new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*) write(6,*)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
@ -1858,6 +1868,7 @@ endif
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), &
crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
@ -1869,7 +1880,7 @@ endif
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e=eIter(1),eIter(2); do i=iIter(1,e),iIter(2,e); do g=gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode == 1) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
if (crystallite_integrateStress(mode,g,i,e)) then if (crystallite_integrateStress(mode,g,i,e)) then
crystallite_converged(g,i,e) = .true. crystallite_converged(g,i,e) = .true.
@ -1993,7 +2004,8 @@ endif
constitutive_previousDotState2(g,i,e)%p = 0.0_pReal constitutive_previousDotState2(g,i,e)%p = 0.0_pReal
constitutive_previousDotState(g,i,e)%p = 0.0_pReal constitutive_previousDotState(g,i,e)%p = 0.0_pReal
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g, i, e) crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), &
crystallite_orientation, g, i, e)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -2059,7 +2071,8 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! wind forward dotStates constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! wind forward dotStates
constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g, i, e) crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), &
crystallite_orientation, g, i, e)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -2216,9 +2229,9 @@ if (verboseDebugger .and. selectiveDebugger) then
write(6,*) write(6,*)
write(6,'(a,f6.1)') 'updateState: crystallite_statedamper',crystallite_statedamper(g,i,e) write(6,'(a,f6.1)') 'updateState: crystallite_statedamper',crystallite_statedamper(g,i,e)
write(6,*) write(6,*)
write(6,'(a,/,12(e14.8,x))') 'updateState: dotState',constitutive_dotState(g,i,e)%p(1:mySize) write(6,'(a,/,12(e12.5,x))') 'updateState: dotState',constitutive_dotState(g,i,e)%p(1:mySize)
write(6,*) write(6,*)
write(6,'(a,/,12(e14.8,x))') 'updateState: new state',constitutive_state(g,i,e)%p(1:mySize) write(6,'(a,/,12(e12.5,x))') 'updateState: new state',constitutive_state(g,i,e)%p(1:mySize)
write(6,*) write(6,*)
write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance',abs(residuum / rTol_crystalliteState & write(6,'(a,/,12(f12.1,x))') 'updateState: resid tolerance',abs(residuum / rTol_crystalliteState &
/ constitutive_state(g,i,e)%p(1:mySize)) / constitutive_state(g,i,e)%p(1:mySize))