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

new subroutine crystallite_orientations: calculates the orientation matrix and the euler angles and - in case of single grain ips - the misorientation for each neighboring ip. 

The calculation of the misorientation is now done once in crystallite init and at the end of every FE increment. This saves a lot of time compared to doing it in dotState for every crystallite subinc.
This commit is contained in:
Christoph Kords 2009-12-18 15:46:33 +00:00
parent 3ab5cdc770
commit 1f7aebfa4d
4 changed files with 221 additions and 118 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
code/constitutive.f90
View File

@ -363,7 +363,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, misorientation, 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 *
@ -380,7 +380,8 @@ use prec, only: pReal, pInt
use debug, only: debug_cumDotStateCalls, & use debug, only: debug_cumDotStateCalls, &
debug_cumDotStateTicks debug_cumDotStateTicks
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
mesh_maxNips mesh_maxNips, &
mesh_maxNipNeighbors
use material, only: phase_constitution, & use material, only: phase_constitution, &
material_phase, & material_phase, &
homogenization_maxNgrains homogenization_maxNgrains
@ -395,6 +396,8 @@ 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) :: &
misorientation
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
@ -421,7 +424,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, subdt, & call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, misorientation, subdt, &
constitutive_state, constitutive_subState0, subdt, ipc, ip, el) constitutive_state, constitutive_subState0, subdt, ipc, ip, el)
end select end select
@ -481,7 +484,7 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
endfunction endfunction
function constitutive_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, dt, subdt, ipc, ip, el) function constitutive_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, misorientation, dt, subdt, ipc, ip, el)
!********************************************************************* !*********************************************************************
!* return array of constitutive results * !* return array of constitutive results *
!* INPUT: * !* INPUT: *
@ -493,7 +496,8 @@ function constitutive_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, dt,
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
mesh_maxNips mesh_maxNips, &
mesh_maxNipNeighbors
use material, only: phase_constitution, & use material, only: phase_constitution, &
material_phase, & material_phase, &
homogenization_maxNgrains homogenization_maxNgrains
@ -507,6 +511,7 @@ function constitutive_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, dt,
integer(pInt), intent(in) :: ipc,ip,el integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: dt, Temperature, subdt real(pReal), intent(in) :: dt, Temperature, subdt
real(pReal), dimension(6), intent(in) :: Tstar_v, subTstar0_v real(pReal), dimension(6), intent(in) :: Tstar_v, subTstar0_v
real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: misorientation
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(constitutive_sizePostResults(ipc,ip,el)) :: constitutive_postResults real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: constitutive_postResults
@ -523,8 +528,8 @@ function constitutive_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, dt,
constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, dt, subdt, & constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, subTstar0_v, Fe, Fp, Temperature, misorientation, &
constitutive_state, constitutive_subState0, & dt, subdt, constitutive_state, constitutive_subState0, &
constitutive_dotstate, ipc, ip, el) constitutive_dotstate, ipc, ip, el)
end select end select

90
code/constitutive_nonlocal.f90
View File

@ -1071,7 +1071,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, misorientation, dt_previous, &
state, previousState, timestep, g,ip,el) state, previousState, timestep, g,ip,el)
use prec, only: pReal, & use prec, only: pReal, &
@ -1089,6 +1089,7 @@ use math, only: math_norm3, &
pi pi
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
mesh_maxNips, & mesh_maxNips, &
mesh_maxNipNeighbors, &
mesh_element, & mesh_element, &
FE_NipNeighbors, & FE_NipNeighbors, &
mesh_ipNeighborhood, & mesh_ipNeighborhood, &
@ -1116,6 +1117,8 @@ 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) :: &
misorientation ! crystal misorientation between me and my neighbor (axis, angle pair)
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
@ -1331,7 +1334,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
lineLength(s,t) = gdot(s,t) / constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) & lineLength(s,t) = gdot(s,t) / constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) &
* math_mul3x3(m(:,s,t), surfaceNormal) * area & * math_mul3x3(m(:,s,t), surfaceNormal) * area &
* constitutive_nonlocal_transmissivity(Fe, ip, el, neighboring_ip, neighboring_el) ! dislocation line length that leaves this interface per second; weighted by a transmissivity factor * constitutive_nonlocal_transmissivity(misorientation(4,n), misorientation(1:3,n)) ! dislocation line length that leaves this interface per second; weighted by a transmissivity factor
thisRhoDot(s,t) = thisRhoDot(s,t) - lineLength(s,t) / mesh_ipVolume(ip,el) ! subtract dislocation density rate (= line length over volume) that leaves through an interface from my dotState ... thisRhoDot(s,t) = thisRhoDot(s,t) - lineLength(s,t) / mesh_ipVolume(ip,el) ! subtract dislocation density rate (= line length over volume) that leaves through an interface from my dotState ...
@ -1463,92 +1466,32 @@ endsubroutine
!********************************************************************* !*********************************************************************
!* transmissivity of IP interface * !* transmissivity of IP interface *
!********************************************************************* !*********************************************************************
function constitutive_nonlocal_transmissivity(Fe, ip,el, neighboring_ip,neighboring_el) function constitutive_nonlocal_transmissivity(misorientationAngle, misorientationAxis)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt
p_vec
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
material_phase, &
phase_constitutionInstance
use math, only: math_inv3x3, &
math_mul33x33, &
math_pDecomposition, &
math_misorientation
use IO, only: IO_warning
implicit none implicit none
!* input variables !* input variables
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & real(pReal), dimension(3), intent(in) :: misorientationAxis ! misorientation axis
Fe ! elastic deformation gradient real(pReal), intent(in) :: misorientationAngle ! misorientation angle
integer(pInt), intent(in) :: ip, & ! my IP
el, & ! my element number
neighboring_ip, & ! neighboring IP
neighboring_el ! neighboring elment number
!* output variables !* output variables
real(pReal) constitutive_nonlocal_transmissivity ! factor determining the transmissivity of an IP interface for dislocations real(pReal) constitutive_nonlocal_transmissivity ! transmissivity of an IP interface for dislocations
!* local variables !* local variables
real(pReal), dimension(3,3) :: U, R, & ! polar decomposition of my Fe
neighboring_U, neighboring_R, & ! polar decomposition of my neighbor's Fe
dR ! net misorientation
real(pReal), dimension(3) :: axis ! rotation axis
real(pReal) angle ! misorientation angle
integer(pInt) myInstance, & ! current instance of this constitution
symmetryType
logical error
! initialize with perfect transmissivity
constitutive_nonlocal_transmissivity = 1.0_pReal
! only change transmissivity when we have a valid neighbor
if ((neighboring_el > 0) .and. (neighboring_ip > 0)) then
! calculate orientations from elastic deformation gradients
! write(6,'(a,/,3(3(f12.8),/))') 'my Fe', Fe(:,:,1,ip,el)
! write(6,'(a,/,3(3(f12.8),/))') 'neighboring Fe', Fe(:,:,1,neighboring_ip,neighboring_el)
call math_pDecomposition(Fe(:,:,1,ip,el), U, R, error) ! polar decomposition of Fe
if (.not. error) &
call math_pDecomposition(Fe(:,:,1,neighboring_ip,neighboring_el), neighboring_U, neighboring_R, error) ! polar decomposition of my neighbors Fe
if (error) then ! polar decomposition failed?
call IO_warning(650, el, ip, 1)
else
! choose symmetry type of my crystal structure
myInstance = phase_constitutionInstance(material_phase(1,ip,el))
select case (constitutive_nonlocal_structureName(myInstance))
case ('fcc','bcc')
symmetryType = 1_pInt
case ('hex')
symmetryType = 2_pInt
case default
symmetryType = 0_pInt
end select
! calculate misorientation
R = transpose(R) ! transpose defines orientation (rotation from current into lattice conf)
neighboring_R = transpose(neighboring_R)
call math_misorientation(axis, angle, dR, R, neighboring_R, symmetryType)
! write(6,'(a,2(x,i3),/,3(3(f12.8),/))') 'my R at',ip,el, R
! write(6,'(a,2(x,i3),/,3(3(f12.8),/))') 'neighboring R at',neighboring_ip, neighboring_el, neighboring_R
! write(6,'(a,x,i3,x,i3,x,f10.3,x,3(f8.5,x),/)') 'constitutive_nonlocal_transmissivity', ip,el, angle, axis
! transmissivity depends on misorientation angle ! transmissivity depends on misorientation angle
if (angle < 3.0_pReal) then if (misorientationAngle < 3.0_pReal) then
constitutive_nonlocal_transmissivity = 1.0_pReal constitutive_nonlocal_transmissivity = 1.0_pReal
elseif (angle < 10.0_pReal) then elseif (misorientationAngle < 10.0_pReal) then
constitutive_nonlocal_transmissivity = 0.5_pReal constitutive_nonlocal_transmissivity = 0.5_pReal
else else
constitutive_nonlocal_transmissivity = 0.1_pReal constitutive_nonlocal_transmissivity = 0.1_pReal
endif endif
endif
endif
endfunction endfunction
@ -1589,7 +1532,7 @@ endfunction
!********************************************************************* !*********************************************************************
!* return array of constitutive results * !* return array of constitutive results *
!********************************************************************* !*********************************************************************
function constitutive_nonlocal_postResults(Tstar_v, previousTstar_v, Fe, Fp, Temperature, dt, dt_previous, & function constitutive_nonlocal_postResults(Tstar_v, previousTstar_v, Fe, Fp, Temperature, misorientation, dt, dt_previous, &
state, previousState, dotState, g,ip,el) state, previousState, dotState, g,ip,el)
use prec, only: pReal, & use prec, only: pReal, &
@ -1606,6 +1549,7 @@ use math, only: math_norm3, &
pi pi
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
mesh_maxNips, & mesh_maxNips, &
mesh_maxNipNeighbors, &
mesh_element, & mesh_element, &
FE_NipNeighbors, & FE_NipNeighbors, &
mesh_ipNeighborhood, & mesh_ipNeighborhood, &
@ -1634,6 +1578,8 @@ 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) :: &
misorientation ! crystal misorientation between me and my neighbor (axis, angle pair)
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
@ -1789,7 +1735,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
if ( sign(1.0_pReal,math_mul3x3(m(:,s,t),surfaceNormal)) == sign(1.0_pReal,gdot(s,t)) ) & if ( sign(1.0_pReal,math_mul3x3(m(:,s,t),surfaceNormal)) == sign(1.0_pReal,gdot(s,t)) ) &
fluxes(s,n,t) = gdot(s,t) / constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) & fluxes(s,n,t) = gdot(s,t) / constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) &
* math_mul3x3(m(:,s,t), surfaceNormal) * area & * math_mul3x3(m(:,s,t), surfaceNormal) * area &
* constitutive_nonlocal_transmissivity(Fe, ip, el, neighboring_ip, neighboring_el) & * constitutive_nonlocal_transmissivity(misorientation(4,n), misorientation(1:3,n)) &
/ mesh_ipVolume(ip,el) / mesh_ipVolume(ip,el)
enddo enddo
enddo enddo

200
code/crystallite.f90
View File

@ -31,7 +31,8 @@ real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, &
real(pReal), dimension (:,:,:,:), allocatable :: crystallite_Tstar_v, & ! current 2nd Piola-Kirchhoff stress vector (end of converged time step) real(pReal), dimension (:,:,:,:), allocatable :: crystallite_Tstar_v, & ! current 2nd Piola-Kirchhoff stress vector (end of converged time step)
crystallite_Tstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_Tstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_partionedTstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of homog inc crystallite_partionedTstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of homog inc
crystallite_subTstar0_v ! 2nd Piola-Kirchhoff stress vector at start of crystallite inc crystallite_subTstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of crystallite inc
crystallite_eulerangles ! euler angles phi1 Phi phi2
real(pReal), dimension (:,:,:,:,:), allocatable :: crystallite_Fe, & ! current "elastic" def grad (end of converged time step) real(pReal), dimension (:,:,:,:,:), allocatable :: crystallite_Fe, & ! current "elastic" def grad (end of converged time step)
crystallite_Fp, & ! current plastic def grad (end of converged time step) crystallite_Fp, & ! current plastic def grad (end of converged time step)
crystallite_invFp, & ! inverse of current plastic def grad (end of converged time step) crystallite_invFp, & ! inverse of current plastic def grad (end of converged time step)
@ -47,7 +48,9 @@ real(pReal), dimension (:,:,:,:,:), allocatable :: crystallite_Fe, &
crystallite_Lp0, & ! plastic velocitiy grad at start of FE inc crystallite_Lp0, & ! plastic velocitiy grad at start of FE inc
crystallite_partionedLp0,& ! plastic velocity grad at start of homog inc crystallite_partionedLp0,& ! plastic velocity grad at start of homog inc
crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc
crystallite_P ! 1st Piola-Kirchhoff stress per grain crystallite_P, & ! 1st Piola-Kirchhoff stress per grain
crystallite_R, & ! crystal orientation (rotation matrix current -> lattice conf)
crystallite_misorientation ! misorientation between two neighboring ips (only calculated for single grain IPs)
real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: crystallite_dPdF, & ! individual dPdF per grain real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: crystallite_dPdF, & ! individual dPdF per grain
crystallite_fallbackdPdF ! dPdF fallback for non-converged grains (elastic prediction) crystallite_fallbackdPdF ! dPdF fallback for non-converged grains (elastic prediction)
real(pReal) crystallite_statedamper ! damping for state update real(pReal) crystallite_statedamper ! damping for state update
@ -78,7 +81,8 @@ subroutine crystallite_init(Temperature)
FEsolving_execIP FEsolving_execIP
use mesh, only: mesh_element, & use mesh, only: mesh_element, &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips, &
mesh_maxNipNeighbors
use material, only: homogenization_Ngrains, & use material, only: homogenization_Ngrains, &
homogenization_maxNgrains, & homogenization_maxNgrains, &
material_EulerAngles, & material_EulerAngles, &
@ -98,11 +102,13 @@ subroutine crystallite_init(Temperature)
gMax, & ! maximum number of grains gMax, & ! maximum number of grains
iMax, & ! maximum number of integration points iMax, & ! maximum number of integration points
eMax, & ! maximum number of elements eMax, & ! maximum number of elements
nMax, & ! maximum number of ip neighbors
myNgrains myNgrains
gMax = homogenization_maxNgrains gMax = homogenization_maxNgrains
iMax = mesh_maxNips iMax = mesh_maxNips
eMax = mesh_NcpElems eMax = mesh_NcpElems
nMax = mesh_maxNipNeighbors
allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = Temperature allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = Temperature
allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal
@ -126,6 +132,9 @@ subroutine crystallite_init(Temperature)
allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal
allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal
allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal
allocate(crystallite_R(3,3,gMax,iMax,eMax)); crystallite_R = 0.0_pReal
allocate(crystallite_eulerangles(3,gMax,iMax,eMax)); crystallite_eulerangles = 0.0_pReal
allocate(crystallite_misorientation(4,nMax,gMax,iMax,eMax)); crystallite_misorientation = 0.0_pReal
allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal
allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 0.0_pReal allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 0.0_pReal
allocate(crystallite_fallbackdPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_fallbackdPdF = 0.0_pReal allocate(crystallite_fallbackdPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_fallbackdPdF = 0.0_pReal
@ -148,6 +157,7 @@ subroutine crystallite_init(Temperature)
do g = 1,myNgrains do g = 1,myNgrains
crystallite_partionedTemperature0(g,i,e) = Temperature ! isothermal assumption crystallite_partionedTemperature0(g,i,e) = Temperature ! isothermal assumption
crystallite_Fp0(:,:,g,i,e) = math_EulerToR(material_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation crystallite_Fp0(:,:,g,i,e) = math_EulerToR(material_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation
crystallite_Fe(:,:,g,i,e) = transpose(crystallite_Fp0(:,:,g,i,e))
crystallite_F0(:,:,g,i,e) = math_I3 crystallite_F0(:,:,g,i,e) = math_I3
crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e) crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e)
crystallite_partionedF0(:,:,g,i,e) = crystallite_F0(:,:,g,i,e) crystallite_partionedF0(:,:,g,i,e) = crystallite_F0(:,:,g,i,e)
@ -159,6 +169,7 @@ subroutine crystallite_init(Temperature)
enddo enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
call crystallite_orientations()
call crystallite_stressAndItsTangent(.true.) ! request elastic answers call crystallite_stressAndItsTangent(.true.) ! request elastic answers
crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback
@ -194,6 +205,9 @@ subroutine crystallite_init(Temperature)
write(6,'(a35,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v) write(6,'(a35,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v)
write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF)
write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF)
write(6,'(a35,x,7(i5,x))') 'crystallite_R: ', shape(crystallite_R)
write(6,'(a35,x,7(i5,x))') 'crystallite_eulerangles: ', shape(crystallite_eulerangles)
write(6,'(a35,x,7(i5,x))') 'crystallite_misorientation: ', shape(crystallite_misorientation)
write(6,'(a35,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt) write(6,'(a35,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt)
write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt)
write(6,'(a35,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac) write(6,'(a35,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac)
@ -316,7 +330,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
storedTemperature storedTemperature
real(pReal), dimension(constitutive_maxSizeState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(constitutive_maxSizeState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
storedState, & storedState
real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
storedDotState storedDotState
logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
storedConvergenceFlag storedConvergenceFlag
@ -371,7 +386,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
! debugger = (e == 1 .and. i == 1 .and. g == 1) debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_converged(g,i,e)) then if (crystallite_converged(g,i,e)) then
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -464,7 +479,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
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_subdt(g,i,e), g, i, e) crystallite_misorientation(:,:,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
@ -557,7 +572,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
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_subdt(g,i,e), g, i, e) crystallite_misorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e)
delta_dotState1 = constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p delta_dotState1 = constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p
delta_dotState2 = constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p delta_dotState2 = constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p
dot_prod12 = dot_product(delta_dotState1, delta_dotState2) dot_prod12 = dot_product(delta_dotState1, delta_dotState2)
@ -595,6 +610,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) count(crystallite_converged(:,:,:)),'grains converged after state integration no.', NiterationState write(6,*) count(crystallite_converged(:,:,:)),'grains converged after state integration no.', NiterationState
write(6,*)
!write(6,'(8(L,x))') crystallite_converged(:,:,:)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged? if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged?
@ -705,7 +722,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
constitutive_dotState(g,i,e)%p = 0.0_pReal constitutive_dotState(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_Fp, crystallite_Temperature(g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), &
crystallite_subdt(g,i,e), g, i, e) crystallite_misorientation(:,:,g,i,e), crystallite_subdt(g,i,e), &
g, i, e)
stateConverged = crystallite_updateState(g,i,e) ! update state stateConverged = crystallite_updateState(g,i,e) ! update state
temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature
converged = stateConverged .and. temperatureConverged converged = stateConverged .and. temperatureConverged
@ -768,7 +786,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
constitutive_dotState(g,i,e)%p = 0.0_pReal constitutive_dotState(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_Fp, crystallite_Temperature(g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), &
crystallite_subdt(g,i,e), g, i, e) crystallite_misorientation(:,:,g,i,e), crystallite_subdt(g,i,e), &
g, i, e)
stateConverged = crystallite_updateState(g,i,e) ! update state stateConverged = crystallite_updateState(g,i,e) ! update state
temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature
converged = stateConverged .and. temperatureConverged converged = stateConverged .and. temperatureConverged
@ -873,7 +892,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if (crystallite_todo(gg,ii,ee)) & if (crystallite_todo(gg,ii,ee)) &
call constitutive_collectDotState(crystallite_Tstar_v(:,gg,ii,ee), crystallite_subTstar0_v(:,gg,ii,ee), & call constitutive_collectDotState(crystallite_Tstar_v(:,gg,ii,ee), crystallite_subTstar0_v(:,gg,ii,ee), &
crystallite_Fe, crystallite_Fp, crystallite_Temperature(gg,ii,ee), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(gg,ii,ee), &
crystallite_subdt(gg,ii,ee), gg, ii, ee) ! collect dot state crystallite_misorientation(:,:,g,i,e), crystallite_subdt(gg,ii,ee), &
gg, ii, ee) ! collect dot state
enddo; enddo; enddo enddo; enddo; enddo
do ee = FEsolving_execElem(1),FEsolving_execElem(2) do ee = FEsolving_execElem(1),FEsolving_execElem(2)
@ -928,10 +948,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo; enddo; enddo ! element,ip,grain loop (e,i,g) enddo; enddo; enddo ! element,ip,grain loop (e,i,g)
crystallite_converged = storedConvergenceFlag crystallite_converged = storedConvergenceFlag
endif endif
endif ! jacobian calculation endif ! jacobian calculation
@ -1396,6 +1412,147 @@ LpLoop: do
!********************************************************************
! calculates orientations and misorientations (in case of single grain ips)
!********************************************************************
subroutine crystallite_orientations()
!*** variables and functions from other modules ***!
use prec, only: pInt, &
pReal
use math, only: math_pDecomposition, &
math_RtoEuler, &
math_misorientation, &
inDeg
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP
use IO, only: IO_warning
use material, only: material_phase, &
homogenization_Ngrains, &
phase_constitution, &
phase_constitutionInstance
use mesh, only: mesh_element, &
mesh_ipNeighborhood, &
FE_NipNeighbors
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, &
constitutive_phenopowerlaw_structure
use constitutive_dislotwin, only: constitutive_dislotwin_label, &
constitutive_dislotwin_structure
use constitutive_nonlocal, only: constitutive_nonlocal_label, &
constitutive_nonlocal_structure
implicit none
!*** input variables ***!
!*** output variables ***!
!*** local variables ***!
integer(pInt) e, & ! element index
i, & ! integration point index
g, & ! grain index
n, & ! neighbor index
myPhase, & ! phase
myStructure, & ! lattice structure
neighboring_e, & ! element index of my neighbor
neighboring_i, & ! integration point index of my neighbor
neighboringPhase, & ! phase of my neighbor
neighboringStructure, & ! lattice structure of my neighbor
symmetryType ! type of crystal symmetry
real(pReal), dimension(3,3) :: U, R, & ! polar decomposition of Fe
netRotation ! net rotation between two orientations
logical error
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,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
if (error) then
call IO_warning(650, e, i, g)
crystallite_R(:,:,g,i,e) = 0.0_pReal
crystallite_eulerangles(:,g,i,e) = (/400.0, 400.0, 400.0/) ! fake orientation
else
crystallite_R(:,:,g,i,e) = transpose(R)
crystallite_eulerangles(:,g,i,e) = math_RtoEuler(crystallite_R(:,:,g,i,e)) * inDeg
endif
enddo
enddo
enddo
!$OMPEND PARALLEL DO
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
if (homogenization_Ngrains(mesh_element(3,e)) == 1_pInt) then ! if single grain ip
myPhase = material_phase(1,i,e) ! get my crystal structure
select case (phase_constitution(myPhase))
case (constitutive_phenopowerlaw_label)
myStructure = constitutive_phenopowerlaw_structure(phase_constitutionInstance(myPhase))
case (constitutive_dislotwin_label)
myStructure = constitutive_dislotwin_structure(phase_constitutionInstance(myPhase))
case (constitutive_nonlocal_label)
myStructure = constitutive_nonlocal_structure(phase_constitutionInstance(myPhase))
case default
myStructure = ''
end select
do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors
neighboring_e = mesh_ipNeighborhood(1,n,i,e)
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists
neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's crystal structure
select case (phase_constitution(neighboringPhase))
case (constitutive_phenopowerlaw_label)
neighboringStructure = constitutive_phenopowerlaw_structure(phase_constitutionInstance(neighboringPhase))
case (constitutive_dislotwin_label)
neighboringStructure = constitutive_dislotwin_structure(phase_constitutionInstance(neighboringPhase))
case (constitutive_nonlocal_label)
neighboringStructure = constitutive_nonlocal_structure(phase_constitutionInstance(neighboringPhase))
case default
neighboringStructure = ''
end select
if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure as me
select case (myStructure) ! get type of symmetry
case (1_pInt, 2_pInt) ! fcc and bcc:
symmetryType = 1_pInt ! -> cubic symmetry
case (3_pInt) ! hex:
symmetryType = 2_pInt ! -> hexagonal symmetry
case default
symmetryType = 0_pInt
end select
call math_misorientation( crystallite_misorientation(1:3,n,1,i,e), &
crystallite_misorientation(4,n,1,i,e), &
netRotation, &
crystallite_R(:,:,1,i,e), &
crystallite_R(:,:,1,neighboring_i,neighboring_e), &
symmetryType) ! calculate misorientation
else
call IO_warning(-1, e, i, g, 'couldnt calculate misorientation because of two different lattice structures')
endif
endif
enddo
endif
enddo
enddo
!$OMPEND PARALLEL DO
endsubroutine
!******************************************************************** !********************************************************************
! return results of particular grain ! return results of particular grain
@ -1410,10 +1567,6 @@ function crystallite_postResults(&
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pInt, & use prec, only: pInt, &
pReal pReal
use math, only: math_pDecomposition, &
math_RtoEuler, &
inDeg
use IO, only: IO_warning
use material, only: material_phase, & use material, only: material_phase, &
material_volume material_volume
use constitutive, only: constitutive_sizePostResults, & use constitutive, only: constitutive_sizePostResults, &
@ -1443,13 +1596,7 @@ function crystallite_postResults(&
c = c+2_pInt c = c+2_pInt
endif endif
if (crystallite_Nresults >= 5) then if (crystallite_Nresults >= 5) then
call math_pDecomposition(crystallite_Fe(:,:,g,i,e),U,R,error) ! polar decomposition of Fe crystallite_postResults(c+1:c+3) = crystallite_eulerangles(:,i,e,g) ! fake orientation
if (error) then
call IO_warning(650,e,i,g)
crystallite_postResults(c+1:c+3) = (/400.0,400.0,400.0/) ! fake orientation
else
crystallite_postResults(c+1:c+3) = math_RtoEuler(transpose(R))*inDeg ! orientation
endif
c = c+3_pInt c = c+3_pInt
endif endif
if (crystallite_Nresults >= 14) then ! deformation gradient if (crystallite_Nresults >= 14) then ! deformation gradient
@ -1460,7 +1607,8 @@ function crystallite_postResults(&
crystallite_postResults(c+1) = constitutive_sizePostResults(g,i,e); c = c+1_pInt ! size of constitutive results crystallite_postResults(c+1) = constitutive_sizePostResults(g,i,e); c = c+1_pInt ! size of constitutive results
crystallite_postResults(c+1:c+constitutive_sizePostResults(g,i,e)) = & crystallite_postResults(c+1:c+constitutive_sizePostResults(g,i,e)) = &
constitutive_postResults(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & constitutive_postResults(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(g,i,e), dt, crystallite_subdt(g,i,e), g, i, e) crystallite_Temperature(g,i,e), crystallite_misorientation(:,:,g,i,e), dt, &
crystallite_subdt(g,i,e), g, i, e)
c = c + constitutive_sizePostResults(g,i,e) c = c + constitutive_sizePostResults(g,i,e)
return return

6
code/homogenization.f90
View File

@ -217,7 +217,8 @@ subroutine materialpoint_stressAndItsTangent(&
crystallite_dt, & crystallite_dt, &
crystallite_requested, & crystallite_requested, &
crystallite_converged, & crystallite_converged, &
crystallite_stressAndItsTangent crystallite_stressAndItsTangent, &
crystallite_orientations
use debug, only: debugger, & use debug, only: debugger, &
debug_MaterialpointLoopDistribution, & debug_MaterialpointLoopDistribution, &
debug_MaterialpointStateLoopDistribution debug_MaterialpointStateLoopDistribution
@ -429,6 +430,9 @@ subroutine materialpoint_stressAndItsTangent(&
enddo ! cutback loop enddo ! cutback loop
! calculate crystal orientations
call crystallite_orientations()
! check for non-performer: any(.not. converged) ! check for non-performer: any(.not. converged)
! replace everybody with odd response ? ! replace everybody with odd response ?