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

* in nonlocal model: dislocation flux now with valid description also for large angle grain boundaries; transfer of dislocations from one slip system to the other according to new variable "constitutive_nonlocal_compatibility" which depends on the angle between slip plane normals and slip directions and is updated once per cycle in function "crystallite_orientations" 

* reactivated debugging functionality for "non-standard" integration methods
This commit is contained in:
Christoph Kords 2010-10-12 13:08:54 +00:00
parent efd92d9b51
commit 724960686f
3 changed files with 307 additions and 163 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
code/constitutive.f90
View File

@ -403,7 +403,7 @@ return
endfunction endfunction
subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,disorientation,ipc,ip,el) subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el)
!********************************************************************* !*********************************************************************
!* This function calculates from state needed variables * !* This function calculates from state needed variables *
!* INPUT: * !* INPUT: *
@ -432,7 +432,6 @@ integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature real(pReal), intent(in) :: Temperature
real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp
real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: disorientation
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
@ -449,7 +448,7 @@ real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: disorientation
call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el) call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Tstar_v, Fe, Fp, disorientation, ipc, ip, el) call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Tstar_v, Fe, Fp, ipc, ip, el)
end select end select
@ -509,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, disorientation, subdt, ipc, ip, el) subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, 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 *
@ -548,8 +547,6 @@ implicit none
integer(pInt), intent(in) :: ipc, ip, el integer(pInt), intent(in) :: ipc, ip, el
real(pReal), intent(in) :: Temperature, & real(pReal), intent(in) :: Temperature, &
subdt subdt
real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: &
disorientation
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
@ -579,7 +576,7 @@ select case (phase_constitution(material_phase(ipc,ip,el)))
constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, disorientation, 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, ipc, ip, el)
end select end select

231
code/constitutive_nonlocal.f90
View File

@ -81,6 +81,8 @@ 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
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
@ -116,7 +118,8 @@ use IO, only: IO_lc, &
IO_intValue, & IO_intValue, &
IO_error IO_error
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
mesh_maxNips mesh_maxNips, &
FE_maxNipNeighbors
use material, only: homogenization_maxNgrains, & use material, only: homogenization_maxNgrains, &
phase_constitution, & phase_constitution, &
phase_constitutionInstance, & phase_constitutionInstance, &
@ -422,6 +425,12 @@ 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))
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
@ -742,7 +751,7 @@ endfunction
!********************************************************************* !*********************************************************************
!* calculates quantities characterizing the microstructure * !* calculates quantities characterizing the microstructure *
!********************************************************************* !*********************************************************************
subroutine constitutive_nonlocal_microstructure(state, Temperature, Tstar_v, Fe, Fp, disorientation, g, ip, el) subroutine constitutive_nonlocal_microstructure(state, Temperature, Tstar_v, Fe, Fp, g, ip, el)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
@ -769,6 +778,7 @@ use mesh, only: mesh_NcpElems, &
mesh_ipCenterOfGravity mesh_ipCenterOfGravity
use material, only: homogenization_maxNgrains, & use material, only: homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_localConstitution, &
phase_constitutionInstance phase_constitutionInstance
use lattice, only: lattice_Sslip, & use lattice, only: lattice_Sslip, &
lattice_Sslip_v, & lattice_Sslip_v, &
@ -791,8 +801,6 @@ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)
Fp ! plastic deformation gradient Fp ! plastic deformation gradient
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: &
disorientation ! crystal disorientation between me and my neighbor as quaternion
!*** input/output variables !*** input/output variables
@ -825,6 +833,8 @@ real(pReal), dimension(3,3) :: sigma, & ! dislocation stre
invPositionDifference ! inverse of a 3x3 matrix containing finite differences of pairs of position vectors invPositionDifference ! inverse of a 3x3 matrix containing finite differences of pairs of position vectors
real(pReal), dimension(6) :: transmissivity, & ! transmissivity factor for each interface real(pReal), dimension(6) :: transmissivity, & ! transmissivity factor for each interface
Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress in Mandel notation Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
rhoExcess ! central excess density
real(pReal), dimension(6,2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(6,2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
neighboring_rhoExcess ! excess density for each neighbor, dislo character and slip system neighboring_rhoExcess ! excess density for each neighbor, dislo character and slip system
real(pReal), dimension(6,3) :: neighboring_position ! position vector of each neighbor when seen from the centreal material point's lattice frame real(pReal), dimension(6,3) :: neighboring_position ! position vector of each neighbor when seen from the centreal material point's lattice frame
@ -879,27 +889,48 @@ F = math_mul33x33(Fe(:,:,g,ip,el), Fp(:,:,g,ip,el))
invFe = math_inv3x3(Fe(:,:,g,ip,el)) invFe = math_inv3x3(Fe(:,:,g,ip,el))
nu = constitutive_nonlocal_nu(myInstance) nu = constitutive_nonlocal_nu(myInstance)
forall (s = 1:ns, c = 1:2) &
rhoExcess(c,s) = state(g,ip,el)%p((2*c-2)*ns+s) + abs(state(g,ip,el)%p((2*c+2)*ns+s)) &
- state(g,ip,el)%p((2*c-1)*ns+s) - abs(state(g,ip,el)%p((2*c+3)*ns+s))
do n = 1,6 do n = 1,6
neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
if ( neighboring_ip == 0 .or. constitutive_nonlocal_transmissivity(disorientation(:,n)) < 1.0_pReal ) then if ( neighboring_ip == 0 .or. neighboring_el == 0 ) then ! at free surfaces ...
neighboring_el = el neighboring_el = el ! ... use central values instead of neighboring values
neighboring_ip = ip neighboring_ip = ip
neighboring_F = F neighboring_position(n,:) = 0.0_pReal
neighboring_rhoExcess(n,:,:) = rhoExcess
elseif (.not. phase_localConstitution(material_phase(1,neighboring_ip,neighboring_el))) then ! for neighbors with local constitution
neighboring_el = el ! ... use central values instead of neighboring values
neighboring_ip = ip
neighboring_position(n,:) = 0.0_pReal
neighboring_rhoExcess(n,:,:) = rhoExcess
elseif (myStructure /= &
constitutive_nonlocal_structure(phase_constitutionInstance(material_phase(1,neighboring_ip,neighboring_el)))) then ! for neighbors with different crystal structure
neighboring_el = el ! ... use central values instead of neighboring values
neighboring_ip = ip
neighboring_position(n,:) = 0.0_pReal
neighboring_rhoExcess(n,:,:) = rhoExcess
else else
neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el))
endif
neighboring_position(n,:) = &
0.5_pReal * math_mul33x3( math_mul33x33(invFe,neighboring_F) + Fp(:,:,g,ip,el), &
mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el) )
forall (s = 1:ns, c = 1:2) & forall (s = 1:ns, c = 1:2) &
neighboring_rhoExcess(n,c,s) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*ns+s) & neighboring_rhoExcess(n,c,s) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*ns+s) &
+ abs(state(g,neighboring_ip,neighboring_el)%p((2*c+2)*ns+s)) & + abs(state(g,neighboring_ip,neighboring_el)%p((2*c+2)*ns+s)) &
- state(g,neighboring_ip,neighboring_el)%p((2*c-1)*ns+s) & - state(g,neighboring_ip,neighboring_el)%p((2*c-1)*ns+s) &
- abs(state(g,neighboring_ip,neighboring_el)%p((2*c+3)*ns+s)) - abs(state(g,neighboring_ip,neighboring_el)%p((2*c+3)*ns+s))
if ( any( neighboring_rhoExcess(n,:,:)*rhoExcess < 0.0_pReal &
.and. abs(neighboring_rhoExcess(n,:,:)) > 1.0_pReal ) ) then ! at grain boundary (=significant change of sign in any excess density) ...
neighboring_el = el ! ... use central values instead of neighboring values
neighboring_ip = ip
neighboring_position(n,:) = 0.0_pReal
neighboring_rhoExcess(n,:,:) = rhoExcess
else
neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el))
neighboring_position(n,:) = &
0.5_pReal * math_mul33x3( math_mul33x33(invFe,neighboring_F) + Fp(:,:,g,ip,el), &
mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el) )
endif
endif
enddo enddo
@ -1192,7 +1223,7 @@ endsubroutine
!********************************************************************* !*********************************************************************
!* rate of change of microstructure * !* rate of change of microstructure *
!********************************************************************* !*********************************************************************
subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, previousTstar_v, Fe, Fp, Temperature, disorientation, 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, g,ip,el)
use prec, only: pReal, & use prec, only: pReal, &
@ -1245,8 +1276,6 @@ real(pReal), intent(in) :: Temperature, & ! temperat
dt_previous ! time increment between previous and current state dt_previous ! time increment between previous and current state
real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation
previousTstar_v ! previous 2nd Piola-Kirchhoff stress in Mandel notation previousTstar_v ! previous 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: &
disorientation ! crystal disorientation between me and my neighbor as quaternion
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
@ -1265,6 +1294,7 @@ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), in
integer(pInt) myInstance, & ! current instance of this constitution integer(pInt) myInstance, & ! current instance of this constitution
myStructure, & ! current lattice structure myStructure, & ! current lattice structure
ns, & ! short notation for the total number of active slip systems ns, & ! short notation for the total number of active slip systems
neighboring_n, & ! neighbor index of myself when seen from my neighbor
neighboring_el, & ! element number of my neighbor neighboring_el, & ! element number of my neighbor
neighboring_ip, & ! integration point of my neighbor neighboring_ip, & ! integration point of my neighbor
c, & ! character of dislocation c, & ! character of dislocation
@ -1273,7 +1303,9 @@ integer(pInt) myInstance, & ! current
opposite_ip, & ! ip of my opposite neighbor opposite_ip, & ! ip of my opposite neighbor
opposite_el, & ! element index of my opposite neighbor opposite_el, & ! element index of my opposite neighbor
t, & ! type of dislocation t, & ! type of dislocation
topp, & ! type of dislocation with opposite sign to t
s, & ! index of my current slip system s, & ! index of my current slip system
s2, &
sLattice, & ! index of my current slip system according to lattice order sLattice, & ! index of my current slip system according to lattice order
i i
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),10) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),10) :: &
@ -1318,7 +1350,8 @@ real(pReal), dimension(3) :: surfaceNormal, & ! surface
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 interfaces for dislocation flux transmissivity, & ! transmissivity of dislocation flux for 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 charcter 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)
@ -1365,7 +1398,7 @@ previousTdislocation_v = previousState(g,ip,el)%p(12*ns+1:12*ns+6)
!*** sanity check for timestep !*** sanity check for timestep
if (timestep <= 0.0_pReal) then ! if illegal timestep... if (timestep <= 0.0_pReal) then ! if illegal timestep...
dotState(1,ip,el)%p(1:10*ns) = 0.0_pReal ! ...return without doing anything (-> zero dotState) dotState(g,ip,el)%p = 0.0_pReal ! ...return without doing anything (-> zero dotState)
return return
endif endif
@ -1470,6 +1503,14 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el) opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el)
opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el) opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el)
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
if ( el == mesh_ipNeighborhood(1,opposite_n,ip,neighboring_el) & ! special case if no neighbor at all...
.and. ip == mesh_ipNeighborhood(2,opposite_n,ip,neighboring_el) ) &
exit ! ...exit without any flux calculation
enddo
endif
if ( neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt ) then ! if neighbor exists, average deformation gradient if ( neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt ) then ! if neighbor exists, average deformation gradient
neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el)) neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el))
Favg = 0.5_pReal * (F + neighboring_F) Favg = 0.5_pReal * (F + neighboring_F)
@ -1482,8 +1523,6 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal) area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal)
surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) ! normalize the surface normal to unit length surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) ! normalize the surface normal to unit length
transmissivity = constitutive_nonlocal_transmissivity(disorientation(:,n))
if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then ! if neighbor exists... if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then ! if neighbor exists...
if ( .not. phase_localConstitution(material_phase(1,neighboring_ip,neighboring_el))) then ! ... and is of nonlocal constitution... if ( .not. phase_localConstitution(material_phase(1,neighboring_ip,neighboring_el))) then ! ... and is of nonlocal constitution...
forall (t = 1:4) & ! ... then calculate neighboring flux density forall (t = 1:4) & ! ... then calculate neighboring flux density
@ -1503,20 +1542,25 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
do s = 1,ns do s = 1,ns
do t = 1,4 do t = 1,4
c = (t + 1) / 2 ! dislocation character
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
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 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
! if (selectiveDebugger .and. s==1) &
! write(6,'(a22,i2,a15,i2,a3,e20.10)') 'outgoing flux of type ',t,' to neighbor ',n,' : ', &
! -lineLength / mesh_ipVolume(ip,el)
else ! incoming flux else ! incoming flux
lineLength = neighboring_fluxdensity(s,t) * math_mul3x3(m(:,s,t), surfaceNormal) * area ! line length that wants to leave this interface do s2 = 1,ns
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) * transmissivity ! subtract negative dislocation flux that enters the material point 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
! if (selectiveDebugger .and. s==1) & transmissivity = constitutive_nonlocal_transmissivity(s,c,s2,neighboring_n,neighboring_ip,neighboring_el)
! write(6,'(a22,i2,a15,i2,a3,e20.10)') 'incoming flux of type ',t,' from neighbor ',n,' : ', & lineLength = neighboring_fluxdensity(s2,t) * math_mul3x3(m(:,s,t), surfaceNormal) * area ! line length that wants to leave this interface
! -lineLength / mesh_ipVolume(ip,el) * transmissivity if (compatibility > 0.0_pReal) then ! compatible with equally signed dislocation density
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) * transmissivity * compatibility ! subtract negative dislocation flux that enters the material point
elseif (compatibility < 0.0_pReal) then ! compatible with inversely signed dislocation density
topp = t + mod(t,2) - mod(t+1,2)
rhoDotFlux(s,topp) = rhoDotFlux(s,topp) - lineLength / mesh_ipVolume(ip,el) * transmissivity * abs(compatibility) ! subtract negative dislocation flux that enters the material point
endif
enddo
endif endif
enddo enddo
enddo enddo
@ -1664,40 +1708,137 @@ endsubroutine
!********************************************************************* !*********************************************************************
!* transmissivity of IP interface * !* TRANSMISSIVITY AND COMPATIBILITY UPDATE *
!* Transmissivity is defined as absolute minimum of the cosine of *
!* the angle between the slip directions and the cosine of the angle *
!* between the slip plane normals. Compatibility is defined as *
!* product of cosine of the angle between the slip plane normals and *
!* signed cosine of the angle between the slip directions *
!********************************************************************* !*********************************************************************
function constitutive_nonlocal_transmissivity(disorientation) subroutine constitutive_nonlocal_updateCompatibility(orientation,i,e)
use prec, only: pReal, & use prec, only: pReal, &
pInt pInt
use math, only: math_QuaternionToAxisAngle use math, only: math_QuaternionDisorientation, &
math_mul3x3, &
math_qRot
use material, only: material_phase, &
phase_constitution, &
phase_localConstitution, &
phase_constitutionInstance, &
homogenization_maxNgrains
use mesh, only: mesh_element, &
mesh_ipNeighborhood, &
FE_NipNeighbors, &
mesh_maxNips, &
mesh_NcpElems
use lattice, only: lattice_sn, &
lattice_sd, &
lattice_st, &
lattice_maxNslip
use debug, only: debugger, &
debug_e, debug_i, debug_g, &
verboseDebugger, &
selectiveDebugger
implicit none implicit none
!* input variables !* input variables
real(pReal), dimension(4), intent(in) :: disorientation ! disorientation as quaternion integer(pInt), intent(in) :: i, & ! ip index
e ! element index
real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
orientation ! crystal orientation in quaternions
!* output variables !* output variables
real(pReal) constitutive_nonlocal_transmissivity ! transmissivity of an IP interface for dislocations
!* local variables !* local variables
real(pReal) disorientationAngle integer(pInt) n, & ! neighbor index
real(pReal), dimension(3) :: disorientationAxis neighboring_e, & ! element index of my neighbor
real(pReal), dimension(4) :: disorientationAxisAngle neighboring_i, & ! integration point index of my neighbor
myPhase, & ! phase
neighboringPhase, &
myInstance, & ! instance of constitution
neighboringInstance, &
myStructure, & ! lattice structure
neighboringStructure, &
myNSlipSystems, & ! number of active slip systems
neighboringNSlipSystems, &
c, & ! dislocation character index (1=edge, 2=screw)
s1, & ! slip system index (me)
s2 ! slip system index (my neighbor)
integer(pInt), dimension(lattice_maxNslip) :: mySlipSystems, & ! slip system numbering according to lattice
neighboringSlipSystems
real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor
real(pReal), dimension(3,lattice_maxNslip) :: myNormals, & ! slip plane normals
neighboringNormals
real(pReal), dimension(3,lattice_maxNslip,2) :: mySlipDirections, & ! slip directions for positive edge and screws
neighboringSlipDirections
disorientationAxisAngle = math_QuaternionToAxisAngle(disorientation) myPhase = material_phase(1,i,e)
myInstance = phase_constitutionInstance(myPhase)
myStructure = constitutive_nonlocal_structure(myInstance)
myNSlipSystems = constitutive_nonlocal_totalNslip(myInstance)
mySlipSystems(1:myNSlipSystems) = constitutive_nonlocal_slipSystemLattice(1:myNSlipSystems,myInstance)
myNormals = lattice_sn(:, mySlipSystems, myStructure)
mySlipDirections(:,:,1) = lattice_sd(:, mySlipSystems, myStructure) ! direction of positive edges
mySlipDirections(:,:,2) = lattice_st(:, mySlipSystems, myStructure) ! direction of positive screws
disorientationAxis = disorientationAxisAngle(1:3) do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors
disorientationAngle = disorientationAxisAngle(4) neighboring_e = mesh_ipNeighborhood(1,n,i,e)
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
if (disorientationAngle < 5.0_pReal) then if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists
constitutive_nonlocal_transmissivity = 1.0_pReal neighboringPhase = material_phase(1,neighboring_i,neighboring_e)
else
constitutive_nonlocal_transmissivity = 0.5_pReal if (.not. phase_localConstitution(neighboringPhase)) then ! neighbor got also nonlocal constitution
neighboringInstance = phase_constitutionInstance(neighboringPhase)
neighboringStructure = constitutive_nonlocal_structure(neighboringInstance)
neighboringNSlipSystems = constitutive_nonlocal_totalNslip(neighboringInstance)
neighboringSlipSystems(1:neighboringNSlipSystems) = constitutive_nonlocal_slipSystemLattice(1:neighboringNSlipSystems,&
neighboringInstance)
neighboringNormals = lattice_sn(:, neighboringSlipSystems, neighboringStructure)
neighboringSlipDirections(:,:,1) = lattice_sd(:, neighboringSlipSystems, neighboringStructure) ! direction of positive edges
neighboringSlipDirections(:,:,2) = lattice_st(:, neighboringSlipSystems, neighboringStructure) ! direction of positive screws
if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me
absoluteMisorientation = math_QuaternionDisorientation( orientation(:,1,i,e), &
orientation(:,1,neighboring_i,neighboring_e), 0_pInt)
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
constitutive_nonlocal_transmissivity(s2,c,s1,n,i,e) = &
min(abs(math_mul3x3(mySlipDirections(:,s1,c), &
math_qRot(absoluteMisorientation, neighboringSlipDirections(:,s2,c)))), &
abs(math_mul3x3(myNormals(:,s1), math_qRot(absoluteMisorientation, neighboringNormals(:,s2)))) )
constitutive_nonlocal_compatibility(s2,c,s1,n,i,e) = &
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
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 ! ...and compatibility
endif
else ! neighbor has local constitution
constitutive_nonlocal_transmissivity(:,:,:,n,i,e) = 1.0_pReal ! assume perfectly transmissive...
constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 1.0_pReal ! ...and compatible
endif
else ! no neighbor present
constitutive_nonlocal_transmissivity(:,:,:,n,i,e) = 1.0_pReal ! perfectly transmissive...
constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 1.0_pReal ! ...and compatible
endif endif
endfunction enddo
endsubroutine

176
code/crystallite.f90
View File

@ -184,7 +184,7 @@ subroutine crystallite_init(Temperature)
allocate(crystallite_statedamper(gMax,iMax,eMax)); crystallite_statedamper = 1.0_pReal allocate(crystallite_statedamper(gMax,iMax,eMax)); crystallite_statedamper = 1.0_pReal
allocate(crystallite_symmetryID(gMax,iMax,eMax)); crystallite_symmetryID = 0.0_pReal !NEW allocate(crystallite_symmetryID(gMax,iMax,eMax)); crystallite_symmetryID = 0.0_pReal !NEW
allocate(crystallite_orientation(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal allocate(crystallite_orientation(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal
allocate(crystallite_orientation0(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal allocate(crystallite_orientation0(4,gMax,iMax,eMax)); crystallite_orientation0 = 0.0_pReal
allocate(crystallite_rotation(4,gMax,iMax,eMax)); crystallite_rotation = 0.0_pReal allocate(crystallite_rotation(4,gMax,iMax,eMax)); crystallite_rotation = 0.0_pReal
allocate(crystallite_disorientation(4,nMax,gMax,iMax,eMax)); crystallite_disorientation = 0.0_pReal allocate(crystallite_disorientation(4,nMax,gMax,iMax,eMax)); crystallite_disorientation = 0.0_pReal
allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true. allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true.
@ -342,7 +342,7 @@ subroutine crystallite_init(Temperature)
write(6,'(a35,x,7(i5,x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0) write(6,'(a35,x,7(i5,x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0)
write(6,'(a35,x,7(i5,x))') 'crystallite_subF: ', shape(crystallite_subF) write(6,'(a35,x,7(i5,x))') 'crystallite_subF: ', shape(crystallite_subF)
write(6,'(a35,x,7(i5,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0) write(6,'(a35,x,7(i5,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0)
write(6,'(a35,x,7(i5,x))') 'crystallite_symmetryID: ', shape(crystallite_symmetryID) !NEW write(6,'(a35,x,7(i5,x))') 'crystallite_symmetryID: ', shape(crystallite_symmetryID)
write(6,'(a35,x,7(i5,x))') 'crystallite_subF0: ', shape(crystallite_subF0) write(6,'(a35,x,7(i5,x))') 'crystallite_subF0: ', shape(crystallite_subF0)
write(6,'(a35,x,7(i5,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) write(6,'(a35,x,7(i5,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0)
write(6,'(a35,x,7(i5,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0) write(6,'(a35,x,7(i5,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0)
@ -483,7 +483,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
logical forceLocalStiffnessCalculation ! flag indicating that stiffness calculation is always done locally logical forceLocalStiffnessCalculation ! flag indicating that stiffness calculation is always done locally
forceLocalStiffnessCalculation = .true. forceLocalStiffnessCalculation = .true.
! --+>> INITIALIZE TO STARTING CONDITION <<+-- ! --+>> INITIALIZE TO STARTING CONDITION <<+--
crystallite_subStep = 0.0_pReal crystallite_subStep = 0.0_pReal
@ -564,7 +563,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad
constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure
crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress
! canÕt restore dotState here, since not yet calculated in first cutback after initialization ! can<EFBFBD>t restore dotState here, since not yet calculated in first cutback after initialization
if (debugger .and. selectiveDebugger) then if (debugger .and. selectiveDebugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a78,f10.8)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& write(6,'(a78,f10.8)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',&
@ -907,7 +906,7 @@ endif
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
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, crystallite_disorientation(:,:,g,i,e), 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
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -919,10 +918,9 @@ 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 = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) 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), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e)
crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), 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
@ -968,8 +966,8 @@ 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
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_Fp, crystallite_disorientation(:,:,g,i,e), 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
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -980,7 +978,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 = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g .and. mode==1)
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
@ -1020,10 +1018,9 @@ 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 = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) 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), & 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), &
crystallite_disorientation(:,:,g,i,e), &
timeStepFraction(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep timeStepFraction(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep
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)
@ -1201,8 +1198,8 @@ 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
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_Fp, crystallite_disorientation(:,:,g,i,e), 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
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -1213,10 +1210,9 @@ 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 = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) 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), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e)
crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), 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
@ -1268,8 +1264,8 @@ do n = 1,5
!$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
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_Fp, crystallite_disorientation(:,:,g,i,e), 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
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -1279,7 +1275,7 @@ do n = 1,5
!$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 = .false. !(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
if (.not. crystallite_integrateStress(mode,g,i,e,c(n))) then ! fraction of original time step if (.not. crystallite_integrateStress(mode,g,i,e,c(n))) then ! fraction of original time step
if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local...
@ -1299,11 +1295,11 @@ do n = 1,5
!$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 = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) 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), & 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), &
crystallite_disorientation(:,:,g,i,e), c(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep c(n)*crystallite_subdt(g,i,e), g,i,e) ! fraction of original timestep
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
@ -1402,8 +1398,8 @@ relTemperatureResiduum = 0.0_pReal
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
sizeDotState = constitutive_sizeDotState(g,i,e) sizeDotState = constitutive_sizeDotState(g,i,e)
if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then
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_Fp, crystallite_disorientation(:,:,g,i,e), 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
endif endif
enddo; enddo; enddo enddo; enddo; enddo
@ -1415,7 +1411,7 @@ relTemperatureResiduum = 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 = .false. !(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)
sizeDotState = constitutive_sizeDotState(g,i,e) sizeDotState = constitutive_sizeDotState(g,i,e)
if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then if ( all(relStateResiduum(1:sizeDotState,g,i,e) < 1.0_pReal) .and. relTemperatureResiduum(g,i,e) < 1.0_pReal ) then
@ -1544,8 +1540,8 @@ 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
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_Fp, crystallite_disorientation(:,:,g,i,e), 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
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -1556,11 +1552,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 = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) 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), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e)
crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), 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)
@ -1602,8 +1597,8 @@ 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
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_Fp, crystallite_disorientation(:,:,g,i,e), 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
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -1613,7 +1608,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 = .false. !(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
if (.not. crystallite_integrateStress(mode,g,i,e)) then if (.not. crystallite_integrateStress(mode,g,i,e)) then
if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local... if (.not. crystallite_localConstitution(g,i,e)) then ! if broken non-local...
@ -1633,11 +1628,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
selectiveDebugger = .false. !(e == debug_e .and. i == debug_i .and. g == debug_g) 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), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e)
crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), 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
@ -1800,8 +1794,8 @@ 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
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_Fp, crystallite_disorientation(:,:,g,i,e), 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
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -1812,10 +1806,9 @@ 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 = .false. !(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), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g,i,e)
crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), 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
@ -1865,8 +1858,8 @@ 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
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_Fp, crystallite_disorientation(:,:,g,i,e), 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
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -1876,7 +1869,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 = .false. !(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
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.
@ -1984,8 +1977,8 @@ 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
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_Fp, crystallite_disorientation(:,:,g,i,e), 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
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -1995,13 +1988,12 @@ 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 = .false. !(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
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), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g, i, e)
crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -2039,7 +2031,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$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 = .false. !(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
crystallite_todo(g,i,e) = crystallite_integrateStress(mode,g,i,e) crystallite_todo(g,i,e) = crystallite_integrateStress(mode,g,i,e)
if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if broken non-local... if ( .not. crystallite_localConstitution(g,i,e) .and. .not. crystallite_todo(g,i,e)) then ! if broken non-local...
@ -2062,13 +2054,12 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$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 = .false. !(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
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), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & crystallite_Fp, crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g, i, e)
crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -2117,8 +2108,8 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$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
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_Fp, crystallite_disorientation(:,:,g,i,e), 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
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -2656,7 +2647,9 @@ use FEsolving, only: FEsolving_execElem, &
use IO, only: IO_warning use IO, only: IO_warning
use material, only: material_phase, & use material, only: material_phase, &
homogenization_Ngrains, & homogenization_Ngrains, &
phase_constitution phase_constitution, &
phase_localConstitution, &
phase_constitutionInstance
use mesh, only: mesh_element, & use mesh, only: mesh_element, &
mesh_ipNeighborhood, & mesh_ipNeighborhood, &
FE_NipNeighbors FE_NipNeighbors
@ -2664,7 +2657,8 @@ use debug, only: debugger, &
debug_e, debug_i, debug_g, & debug_e, debug_i, debug_g, &
verboseDebugger, & verboseDebugger, &
selectiveDebugger selectiveDebugger
use constitutive_nonlocal, only: constitutive_nonlocal_label use constitutive_nonlocal, only: constitutive_nonlocal_structure, &
constitutive_nonlocal_updateCompatibility
implicit none implicit none
@ -2677,21 +2671,24 @@ integer(pInt) e, & ! element index
i, & ! integration point index i, & ! integration point index
g, & ! grain index g, & ! grain index
n, & ! neighbor index n, & ! neighbor index
myPhase, & ! phase
neighboring_e, & ! element index of my neighbor neighboring_e, & ! element index of my neighbor
neighboring_i, & ! integration point index of my neighbor neighboring_i, & ! integration point index of my neighbor
neighboringPhase, & ! phase of my neighbor myPhase, & ! phase
neighboringStructure ! lattice structure of my neighbor neighboringPhase, &
myInstance, & ! instance of constitution
neighboringInstance, &
myStructure, & ! lattice structure
neighboringStructure
real(pReal), dimension(3,3) :: U, R real(pReal), dimension(3,3) :: U, R
logical error logical error
! --- CALCULATE ORIENTATION AND LATTICE ROTATION ---
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e)) do g = 1,homogenization_Ngrains(mesh_element(3,e))
! calculate orientation in terms of rotation matrix and euler angles
call math_pDecomposition(crystallite_Fe(:,:,g,i,e), U, R, error) ! polar decomposition of Fe call math_pDecomposition(crystallite_Fe(:,:,g,i,e), U, R, error) ! polar decomposition of Fe
if (error) then if (error) then
call IO_warning(650, e, i, g) call IO_warning(650, e, i, g)
@ -2710,42 +2707,51 @@ logical error
enddo enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL ---
! --- we use crystallite_orientation from above, so need a seperate loop
!$OMP PARALLEL DO !$OMP PARALLEL DO
! Another loop for nonlocal material which uses the orientations from the first one.
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
selectiveDebugger = (e == debug_e .and. i == debug_i) selectiveDebugger = (e == debug_e .and. i == debug_i)
myPhase = material_phase(1,i,e) ! get my crystal structure myPhase = material_phase(1,i,e) ! get my phase
if (phase_constitution(myPhase) == constitutive_nonlocal_label) then ! if nonlocal model if (.not. phase_localConstitution(myPhase)) then ! if nonlocal model
myInstance = phase_constitutionInstance(myPhase)
myStructure = constitutive_nonlocal_structure(myInstance) ! get my crystal structure
! --- calculate disorientation between me and my neighbor ---
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) ! get my neighbor's phase
neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's crystal structure if (.not. phase_localConstitution(neighboringPhase)) then ! neighbor got also nonlocal constitution
if (myPhase == neighboringPhase) then ! if my neighbor has same phase like me neighboringInstance = phase_constitutionInstance(neighboringPhase)
neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure
if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me
crystallite_disorientation(:,n,1,i,e) = & crystallite_disorientation(:,n,1,i,e) = &
math_QuaternionDisorientation( crystallite_orientation(:,1,i,e), & math_QuaternionDisorientation( crystallite_orientation(:,1,i,e), &
crystallite_orientation(:,1,neighboring_i,neighboring_e), & crystallite_orientation(:,1,neighboring_i,neighboring_e), &
crystallite_symmetryID(1,i,e)) ! calculate disorientation crystallite_symmetryID(1,i,e)) ! calculate disorientation
else ! for neighbor with different phase else ! for neighbor with different phase
crystallite_disorientation(:,n,1,i,e) = (/0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal/) ! 180 degree rotation about 100 axis crystallite_disorientation(:,n,1,i,e) = (/0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal/) ! 180 degree rotation about 100 axis
endif
else ! for neighbor with local constitution
crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity
endif endif
else ! no existing neighbor else ! no existing neighbor
crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity
endif endif
if (verboseDebugger .and. selectiveDebugger) then
!$OMP CRITICAL (write2out)
write(6,'(a27,i2,a3,4(f12.5,x))') 'disorientation to neighbor ',n,' : ',crystallite_disorientation(:,n,1,i,e)
!$OMP END CRITICAL (write2out)
endif
enddo enddo
! --- calculate compatibility and transmissivity between me and my neighbor ---
call constitutive_nonlocal_updateCompatibility(crystallite_orientation,i,e)
endif endif
enddo enddo
enddo enddo