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

named some loops and made BUTCHER TABLEAUs parameters in crystallite_integrateStateRKCK45

This commit is contained in:
Martin Diehl 2013-04-29 11:17:30 +00:00
parent 9afaa9f4d0
commit 533fcec96f
1 changed files with 829 additions and 842 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

403
code/crystallite.f90
View File

@ -26,7 +26,9 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module crystallite module crystallite
use prec, only: pReal, pInt use prec, only: &
pReal, &
pInt
implicit none implicit none
private private
@ -41,26 +43,26 @@ module crystallite
integer(pInt), dimension(:,:,:), allocatable, private :: & integer(pInt), dimension(:,:,:), allocatable, private :: &
crystallite_symmetryID !< crystallographic symmetry 1=cubic 2=hexagonal, needed in all orientation calcs crystallite_symmetryID !< crystallographic symmetry 1=cubic 2=hexagonal, needed in all orientation calcs
real(pReal), dimension (:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:), allocatable, public :: &
crystallite_dt, & !< requested time increment of each grain crystallite_dt, & !< requested time increment of each grain
crystallite_Temperature, & !< Temp of each grain crystallite_Temperature, & !< Temp of each grain
crystallite_partionedTemperature0 !< Temp of each grain at start of homog inc crystallite_partionedTemperature0 !< Temp of each grain at start of homog inc
real(pReal), dimension (:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:), allocatable, private :: &
crystallite_subdt, & !< substepped time increment of each grain crystallite_subdt, & !< substepped time increment of each grain
crystallite_subFrac, & !< already calculated fraction of increment crystallite_subFrac, & !< already calculated fraction of increment
crystallite_subStep, & !< size of next integration step crystallite_subStep, & !< size of next integration step
crystallite_subTemperature0, & !< Temp of each grain at start of crystallite inc crystallite_subTemperature0, & !< Temp of each grain at start of crystallite inc
crystallite_dotTemperature !< evolution of Temperature of each grain crystallite_dotTemperature !< evolution of Temperature of each grain
real(pReal), dimension (:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:), allocatable, public :: &
crystallite_Tstar_v, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) 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
real(pReal), dimension (:,:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:,:), allocatable, private :: &
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_orientation, & !< orientation as quaternion crystallite_orientation, & !< orientation as quaternion
crystallite_orientation0, & !< initial orientation as quaternion crystallite_orientation0, & !< initial orientation as quaternion
crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees) in crystal reference frame crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees) in crystal reference frame
real(pReal), dimension (:,:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
crystallite_Fp, & !< current plastic def grad (end of converged time step) crystallite_Fp, & !< current plastic def grad (end of converged time step)
crystallite_Fp0, & !< plastic def grad at start of FE inc crystallite_Fp0, & !< plastic def grad at start of FE inc
crystallite_partionedFp0,& !< plastic def grad at start of homog inc crystallite_partionedFp0,& !< plastic def grad at start of homog inc
@ -71,7 +73,7 @@ module crystallite
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_P !< 1st Piola-Kirchhoff stress per grain crystallite_P !< 1st Piola-Kirchhoff stress per grain
real(pReal), dimension (:,:,:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:,:,:), allocatable, private :: &
crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc
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)
@ -80,20 +82,20 @@ module crystallite
crystallite_subF0, & !< def grad at start of crystallite inc crystallite_subF0, & !< def grad at start of crystallite inc
crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc
crystallite_disorientation !< disorientation between two neighboring ips (only calculated for single grain IPs) crystallite_disorientation !< disorientation between two neighboring ips (only calculated for single grain IPs)
real(pReal), dimension (:,:,:,:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:,:,:,:), allocatable, public :: &
crystallite_dPdF, & !< current individual dPdF per grain (end of converged time step) crystallite_dPdF, & !< current individual dPdF per grain (end of converged time step)
crystallite_dPdF0, & !< individual dPdF per grain at start of FE inc crystallite_dPdF0, & !< individual dPdF per grain at start of FE inc
crystallite_partioneddPdF0 !< individual dPdF per grain at start of homog inc crystallite_partioneddPdF0 !< individual dPdF per grain at start of homog inc
real(pReal), dimension (:,:,:,:,:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:,:,:,:,:), allocatable, private :: &
crystallite_fallbackdPdF !< dPdF fallback for non-converged grains (elastic prediction) crystallite_fallbackdPdF !< dPdF fallback for non-converged grains (elastic prediction)
logical, dimension (:,:,:), allocatable, public :: & logical, dimension(:,:,:), allocatable, public :: &
crystallite_requested !< flag to request crystallite calculation crystallite_requested !< flag to request crystallite calculation
logical, dimension (:,:,:), allocatable, public, protected :: & logical, dimension(:,:,:), allocatable, public, protected :: &
crystallite_converged, & !< convergence flag crystallite_converged, & !< convergence flag
crystallite_localPlasticity !< indicates this grain to have purely local constitutive law crystallite_localPlasticity !< indicates this grain to have purely local constitutive law
logical, dimension (:,:,:), allocatable, private :: & logical, dimension(:,:,:), allocatable, private :: &
crystallite_todo !< flag to indicate need for further computation crystallite_todo !< flag to indicate need for further computation
logical, dimension (:,:), allocatable, private :: & logical, dimension(:,:), allocatable, private :: &
crystallite_clearToWindForward, & !< description now available crystallite_clearToWindForward, & !< description now available
crystallite_clearToCutback, & !< description now available crystallite_clearToCutback, & !< description now available
crystallite_syncSubFrac, & !< description now available crystallite_syncSubFrac, & !< description now available
@ -122,37 +124,56 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_init(Temperature) subroutine crystallite_init(Temperature)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: debug_info, & use debug, only: &
debug_info, &
debug_reset, & debug_reset, &
debug_level, & debug_level, &
debug_crystallite, & debug_crystallite, &
debug_levelBasic debug_levelBasic
use numerics, only: & use numerics, only: &
usePingPong usePingPong
use math, only: math_I3, & use math, only: &
math_I3, &
math_EulerToR, & math_EulerToR, &
math_inv33, & math_inv33, &
math_transpose33, & math_transpose33, &
math_mul33xx33, & math_mul33xx33, &
math_mul33x33 math_mul33x33
use FEsolving, only: FEsolving_execElem, & use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP FEsolving_execIP
use mesh, only: mesh_element, & use mesh, only: &
mesh_element, &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips, & mesh_maxNips, &
mesh_maxNipNeighbors mesh_maxNipNeighbors
use IO use IO, only: &
IO_timeStamp, &
IO_open_jobFile_stat, &
IO_open_file, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_write_jobFile, &
IO_error
use material use material
use lattice, only: lattice_symmetryType use lattice, only: &
lattice_symmetryType
use constitutive, only: constitutive_microstructure use constitutive, only: &
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, & constitutive_microstructure
use constitutive_phenopowerlaw, only: &
constitutive_phenopowerlaw_label, &
constitutive_phenopowerlaw_structureName constitutive_phenopowerlaw_structureName
use constitutive_titanmod, only: constitutive_titanmod_label, & use constitutive_titanmod, only: &
constitutive_titanmod_label, &
constitutive_titanmod_structureName constitutive_titanmod_structureName
use constitutive_dislotwin, only: constitutive_dislotwin_label, & use constitutive_dislotwin, only: &
constitutive_dislotwin_label, &
constitutive_dislotwin_structureName constitutive_dislotwin_structureName
use constitutive_nonlocal, only: constitutive_nonlocal_label, & use constitutive_nonlocal, only: &
constitutive_nonlocal_label, &
constitutive_nonlocal_structureName constitutive_nonlocal_structureName
implicit none implicit none
@ -160,16 +181,16 @@ subroutine crystallite_init(Temperature)
integer(pInt), parameter :: myFile = 200_pInt, & integer(pInt), parameter :: myFile = 200_pInt, &
maxNchunks = 2_pInt maxNchunks = 2_pInt
!*** local variables ***!
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) g, & ! grain number integer(pInt) :: &
i, & ! integration point number g, & !< grain number
e, & ! element number i, & !< integration point number
gMax, & ! maximum number of grains e, & !< element number
iMax, & ! maximum number of integration points gMax, & !< maximum number of grains
eMax, & ! maximum number of elements iMax, & !< maximum number of integration points
nMax, & ! maximum number of ip neighbors eMax, & !< maximum number of elements
myNgrains, & ! number of grains in current IP nMax, & !< maximum number of ip neighbors
myNgrains, & !< number of grains in current IP
section, & section, &
j, & j, &
p, & p, &
@ -177,15 +198,14 @@ subroutine crystallite_init(Temperature)
mySize, & mySize, &
myPhase, & myPhase, &
myMat myMat
character(len=64) tag character(len=64) :: tag
character(len=1024) line character(len=1024) :: line
write(6,'(/,a)') ' <<<+- crystallite init -+>>>' write(6,'(/,a)') ' <<<+- crystallite init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
write(6,'(a16,a)') ' Current time : ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
gMax = homogenization_maxNgrains gMax = homogenization_maxNgrains
iMax = mesh_maxNips iMax = mesh_maxNips
eMax = mesh_NcpElems eMax = mesh_NcpElems
@ -276,9 +296,6 @@ subroutine crystallite_init(Temperature)
100 close(myFile) 100 close(myFile)
do i = 1_pInt,material_Ncrystallite ! sanity checks
enddo
do i = 1_pInt,material_Ncrystallite do i = 1_pInt,material_Ncrystallite
do j = 1_pInt,crystallite_Noutput(i) do j = 1_pInt,crystallite_Noutput(i)
select case(crystallite_output(j,i)) select case(crystallite_output(j,i))
@ -310,14 +327,13 @@ subroutine crystallite_init(Temperature)
crystallite_sizePostResults(microstructure_crystallite(j))) crystallite_sizePostResults(microstructure_crystallite(j)))
enddo enddo
!--------------------------------------------------------------------------------------------------
! write description file for crystallite output ! write description file for crystallite output
call IO_write_jobFile(myFile,'outputCrystallite') call IO_write_jobFile(myFile,'outputCrystallite')
do p = 1_pInt,material_Ncrystallite do p = 1_pInt,material_Ncrystallite
write(myFile,*) write(myFile,'(/,a,/)') '['//trim(crystallite_name(p))//']'
write(myFile,'(a)') '['//trim(crystallite_name(p))//']'
write(myFile,*)
do e = 1_pInt,crystallite_Noutput(p) do e = 1_pInt,crystallite_Noutput(p)
write(myFile,'(a,i4)') trim(crystallite_output(e,p))//char(9),crystallite_sizePostResult(e,p) write(myFile,'(a,i4)') trim(crystallite_output(e,p))//char(9),crystallite_sizePostResult(e,p)
enddo enddo
@ -325,8 +341,9 @@ subroutine crystallite_init(Temperature)
close(myFile) close(myFile)
!--------------------------------------------------------------------------------------------------
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements ! initialize
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements
myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1_pInt:myNgrains) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1_pInt:myNgrains)
crystallite_Fp0(1:3,1:3,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation crystallite_Fp0(1:3,1:3,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation
@ -336,19 +353,18 @@ do e = FEsolving_execElem(1),FEsolving_execElem(2)
crystallite_Fp(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) crystallite_Fp(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e)
crystallite_requested(g,i,e) = .true. crystallite_requested(g,i,e) = .true.
endforall endforall
enddo enddo
if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601) if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601_pInt)
crystallite_partionedTemperature0 = Temperature ! isothermal assumption crystallite_partionedTemperature0 = Temperature ! isothermal assumption
crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFp0 = crystallite_Fp0
crystallite_partionedF0 = crystallite_F0 crystallite_partionedF0 = crystallite_F0
crystallite_partionedF = crystallite_F0 crystallite_partionedF = crystallite_F0
!--------------------------------------------------------------------------------------------------
! Initialize crystallite_symmetryID(g,i,e) ! Initialize crystallite_symmetryID
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1_pInt,myNgrains do g = 1_pInt,myNgrains
@ -356,41 +372,46 @@ do e = FEsolving_execElem(1),FEsolving_execElem(2)
myMat = phase_plasticityInstance(myPhase) myMat = phase_plasticityInstance(myPhase)
select case (phase_plasticity(myPhase)) select case (phase_plasticity(myPhase))
case (constitutive_phenopowerlaw_label) case (constitutive_phenopowerlaw_label)
crystallite_symmetryID(g,i,e) = lattice_symmetryType(constitutive_phenopowerlaw_structureName(myMat)) crystallite_symmetryID(g,i,e) = &
lattice_symmetryType(constitutive_phenopowerlaw_structureName(myMat))
case (constitutive_titanmod_label) case (constitutive_titanmod_label)
crystallite_symmetryID(g,i,e) = lattice_symmetryType(constitutive_titanmod_structureName(myMat)) crystallite_symmetryID(g,i,e) = &
lattice_symmetryType(constitutive_titanmod_structureName(myMat))
case (constitutive_dislotwin_label) case (constitutive_dislotwin_label)
crystallite_symmetryID(g,i,e) = lattice_symmetryType(constitutive_dislotwin_structureName(myMat)) crystallite_symmetryID(g,i,e) = &
lattice_symmetryType(constitutive_dislotwin_structureName(myMat))
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
crystallite_symmetryID(g,i,e) = lattice_symmetryType(constitutive_nonlocal_structureName(myMat)) crystallite_symmetryID(g,i,e) = &
lattice_symmetryType(constitutive_nonlocal_structureName(myMat))
case default case default
crystallite_symmetryID(g,i,e) = 0_pInt ! does this happen for j2 material? crystallite_symmetryID(g,i,e) = 0_pInt !< ToDo: does this happen for j2 material?
end select end select
enddo enddo
enddo enddo
enddo enddo
call crystallite_orientations()
crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations
call crystallite_orientations() !$OMP PARALLEL DO PRIVATE(myNgrains)
crystallite_orientation0 = crystallite_orientation ! Store initial orientations for calculation of grain rotations
!$OMP PARALLEL DO PRIVATE(myNgrains)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1_pInt,myNgrains do g = 1_pInt,myNgrains
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states crystallite_Fe(1:3,1:3,g,i,e), crystallite_Fp(1:3,1:3,g,i,e),g,i,e) ! update dependent state variables to be consistent with basic states
enddo enddo
enddo enddo
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
call crystallite_stressAndItsTangent(.true.,.false.) ! request elastic answers call crystallite_stressAndItsTangent(.true.,.false.) ! request elastic answers
crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback
! *** Output *** !--------------------------------------------------------------------------------------------------
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then ! debug output
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Temperature: ', shape(crystallite_Temperature) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Temperature: ', shape(crystallite_Temperature)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_dotTemperature: ', shape(crystallite_dotTemperature) write(6,'(a35,1x,7(i8,1x))') 'crystallite_dotTemperature: ', shape(crystallite_dotTemperature)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe) write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe)
@ -437,7 +458,7 @@ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,*) write(6,*)
write(6,*) 'Number of nonlocal grains: ',count(.not. crystallite_localPlasticity) write(6,*) 'Number of nonlocal grains: ',count(.not. crystallite_localPlasticity)
flush(6) flush(6)
endif endif
call debug_info call debug_info
call debug_reset call debug_reset
@ -449,7 +470,8 @@ end subroutine crystallite_init
!> @brief calculate stress (P) and tangent (dPdF) for crystallites !> @brief calculate stress (P) and tangent (dPdF) for crystallites
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
use numerics, only: subStepMinCryst, & use numerics, only: &
subStepMinCryst, &
subStepSizeCryst, & subStepSizeCryst, &
stepIncreaseCryst, & stepIncreaseCryst, &
pert_Fg, & pert_Fg, &
@ -460,7 +482,8 @@ use numerics, only: subStepMinCryst, &
numerics_timeSyncing, & numerics_timeSyncing, &
relevantStrain, & relevantStrain, &
analyticJaco analyticJaco
use debug, only: debug_level, & use debug, only: &
debug_level, &
debug_crystallite, & debug_crystallite, &
debug_levelBasic, & debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
@ -469,8 +492,10 @@ use debug, only: debug_level, &
debug_i, & debug_i, &
debug_g, & debug_g, &
debug_CrystalliteLoopDistribution debug_CrystalliteLoopDistribution
use IO, only: IO_warning use IO, only: &
use math, only: math_inv33, & IO_warning
use math, only: &
math_inv33, &
math_identity2nd, & math_identity2nd, &
math_transpose33, & math_transpose33, &
math_mul33x33, & math_mul33x33, &
@ -479,18 +504,22 @@ use math, only: math_inv33, &
math_Mandel33to6, & math_Mandel33to6, &
math_I3, & math_I3, &
math_mul3333xx3333 math_mul3333xx3333
use FEsolving, only: FEsolving_execElem, & use FEsolving, only: &
FEsolving_execElem, &
FEsolving_execIP FEsolving_execIP
use mesh, only: mesh_element, & use mesh, only: &
mesh_element, &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips, & mesh_maxNips, &
mesh_ipNeighborhood, & mesh_ipNeighborhood, &
FE_NipNeighbors, & FE_NipNeighbors, &
FE_geomtype, & FE_geomtype, &
FE_celltype FE_cellType
use material, only: homogenization_Ngrains, & use material, only: &
homogenization_Ngrains, &
homogenization_maxNgrains homogenization_maxNgrains
use constitutive, only: constitutive_sizeState, & use constitutive, only: &
constitutive_sizeState, &
constitutive_sizeDotState, & constitutive_sizeDotState, &
constitutive_state, & constitutive_state, &
constitutive_state_backup, & constitutive_state_backup, &
@ -501,30 +530,29 @@ use constitutive, only: constitutive_sizeState, &
constitutive_dotState_backup, & constitutive_dotState_backup, &
constitutive_TandItsTangent constitutive_TandItsTangent
implicit none
implicit none logical, intent(in) :: updateJaco, rate_sensitivity ! flag indicating whether we want to update the Jacobian (stiffness) or not
logical, intent(in) :: updateJaco, rate_sensitivity ! flag indicating wehther we want to update the Jacobian (stiffness) or not real(pReal) myPert, & ! perturbation with correct sign
real(pReal) myPert, & ! perturbation with correct sign
formerSubStep, & formerSubStep, &
subFracIntermediate subFracIntermediate
real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient
Fe_guess, & ! guess for elastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient
Tstar ! 2nd Piola-Kirchhoff stress tensor Tstar ! 2nd Piola-Kirchhoff stress tensor
real(pReal), dimension(3,3,3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(3,3,3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
dPdF_perturbation1, & dPdF_perturbation1, &
dPdF_perturbation2 dPdF_perturbation2
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
F_backup, & F_backup, &
Fp_backup, & Fp_backup, &
InvFp_backup, & InvFp_backup, &
Fe_backup, & Fe_backup, &
Lp_backup, & Lp_backup, &
P_backup P_backup
real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
Tstar_v_backup Tstar_v_backup
real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
Temperature_backup Temperature_backup
integer(pInt) NiterationCrystallite, & ! number of iterations in crystallite loop integer(pInt) NiterationCrystallite, & ! number of iterations in crystallite loop
e, & ! element index e, & ! element index
i, & ! integration point index i, & ! integration point index
g, & ! grain index g, & ! grain index
@ -537,24 +565,23 @@ integer(pInt) NiterationCrystallite, &
p, & p, &
perturbation , & ! loop counter for forward,backward perturbation mode perturbation , & ! loop counter for forward,backward perturbation mode
myNgrains myNgrains
logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
convergenceFlag_backup convergenceFlag_backup
! local variables used for calculating analytic Jacobian ! local variables used for calculating analytic Jacobian
real(pReal), dimension(3,3):: Fpinv_rate, & real(pReal), dimension(3,3):: Fpinv_rate, &
FDot_inv, & FDot_inv, &
junk junk
real(pReal), dimension(3,3,3,3) :: dSdFe, & real(pReal), dimension(3,3,3,3) :: dSdFe, &
dFedF, & dFedF, &
dFedFdot, & dFedFdot, &
dSdF, & dSdF, &
dSdFdot, & dSdFdot, &
dFp_invdFdot, & dFp_invdFdot, &
junk2 junk2
real(pReal) :: counter real(pReal) :: counter
! --+>> INITIALIZE TO STARTING CONDITION <<+--
if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt& if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt&
.and. debug_e > 0 .and. debug_e <= mesh_NcpElems & .and. debug_e > 0 .and. debug_e <= mesh_NcpElems &
.and. debug_i > 0 .and. debug_i <= mesh_maxNips & .and. debug_i > 0 .and. debug_i <= mesh_maxNips &
.and. debug_g > 0 .and. debug_g <= homogenization_maxNgrains) then .and. debug_g > 0 .and. debug_g <= homogenization_maxNgrains) then
@ -569,13 +596,16 @@ if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt&
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp0', & write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp0', &
math_transpose33(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e)) math_transpose33(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e))
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
crystallite_subStep = 0.0_pReal !--------------------------------------------------------------------------------------------------
! initialize to starting condition
crystallite_subStep = 0.0_pReal
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1_pInt:myNgrains, crystallite_requested(g,i,e)) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), &
g = 1_pInt:myNgrains, crystallite_requested(g,i,e))
crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature
constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure
crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_partionedFp0(1:3,1:3,g,i,e) ! ...plastic def grad crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_partionedFp0(1:3,1:3,g,i,e) ! ...plastic def grad
@ -590,14 +620,14 @@ do e = FEsolving_execElem(1),FEsolving_execElem(2)
crystallite_todo(g,i,e) = .true. crystallite_todo(g,i,e) = .true.
crystallite_converged(g,i,e) = .false. ! pretend failed step of twice the required size crystallite_converged(g,i,e) = .false. ! pretend failed step of twice the required size
endforall endforall
enddo enddo elementLooping1
! --+>> CRYSTALLITE CUTBACK LOOP <<+-- ! --+>> CRYSTALLITE CUTBACK LOOP <<+--
NiterationCrystallite = 0_pInt NiterationCrystallite = 0_pInt
numerics_integrationMode = 1_pInt numerics_integrationMode = 1_pInt
do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2)))) ! cutback loop for crystallites cutbackLooping: do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2))))
if (any(.not. crystallite_localPlasticity) .and. numerics_timeSyncing) then if (any(.not. crystallite_localPlasticity) .and. numerics_timeSyncing) then
@ -873,7 +903,8 @@ do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2)))
crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad
!$OMP FLUSH(crystallite_subF0) !$OMP FLUSH(crystallite_subF0)
crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation crystallite_subFe0(1:3,1:3,g,i,e) = &
math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation
crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient
constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure
crystallite_subTstar0_v(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) ! ...2nd PK stress crystallite_subTstar0_v(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) ! ...2nd PK stress
@ -948,7 +979,8 @@ do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2)))
if (crystallite_todo(g,i,e) .and. (crystallite_clearToWindForward(i,e) .or. crystallite_clearToCutback(i,e))) then if (crystallite_todo(g,i,e) .and. (crystallite_clearToWindForward(i,e) .or. crystallite_clearToCutback(i,e))) then
crystallite_subF(1:3,1:3,g,i,e) = crystallite_subF0(1:3,1:3,g,i,e) & crystallite_subF(1:3,1:3,g,i,e) = crystallite_subF0(1:3,1:3,g,i,e) &
+ crystallite_subStep(g,i,e) & + crystallite_subStep(g,i,e) &
* (crystallite_partionedF(1:3,1:3,g,i,e) - crystallite_partionedF0(1:3,1:3,g,i,e)) * (crystallite_partionedF(1:3,1:3,g,i,e) &
- crystallite_partionedF0(1:3,1:3,g,i,e))
!$OMP FLUSH(crystallite_subF) !$OMP FLUSH(crystallite_subF)
crystallite_Fe(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) crystallite_Fe(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e))
crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e)
@ -987,12 +1019,10 @@ do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2)))
if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,'(/,a,e12.5)') '<< CRYST >> min(subStep) ',minval(crystallite_subStep)
write(6,'(a,e12.5)') '<< CRYST >> min(subStep) ',minval(crystallite_subStep)
write(6,'(a,e12.5)') '<< CRYST >> max(subStep) ',maxval(crystallite_subStep) write(6,'(a,e12.5)') '<< CRYST >> max(subStep) ',maxval(crystallite_subStep)
write(6,'(a,e12.5)') '<< CRYST >> min(subFrac) ',minval(crystallite_subFrac) write(6,'(a,e12.5)') '<< CRYST >> min(subFrac) ',minval(crystallite_subFrac)
write(6,'(a,e12.5)') '<< CRYST >> max(subFrac) ',maxval(crystallite_subFrac) write(6,'(a,e12.5,/)') '<< CRYST >> max(subFrac) ',maxval(crystallite_subFrac)
write(6,*)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -1018,7 +1048,7 @@ do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2)))
NiterationCrystallite = NiterationCrystallite + 1_pInt NiterationCrystallite = NiterationCrystallite + 1_pInt
enddo ! cutback loop enddo cutbackLooping
! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+-- ! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+--
@ -1044,11 +1074,12 @@ do e = FEsolving_execElem(1),FEsolving_execElem(2)
.or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g
write(6,*) write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', &
write(6,'(a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp', math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)) write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp', &
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp', math_transpose33(crystallite_Lp(1:3,1:3,g,i,e)) math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
write(6,*) write(6,'(a,/,3(12x,3(f14.9,1x)/),/)') '<< CRYST >> Lp', &
math_transpose33(crystallite_Lp(1:3,1:3,g,i,e))
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
enddo enddo
@ -1184,7 +1215,7 @@ if(updateJaco) then
! --- STIFFNESS ACCORDING TO PERTURBATION METHOD AND CONVERGENCE --- ! --- STIFFNESS ACCORDING TO PERTURBATION METHOD AND CONVERGENCE ---
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
select case(pert_method) select case(pert_method)
case(1_pInt) case(1_pInt)
@ -1204,7 +1235,7 @@ if(updateJaco) then
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, & forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
crystallite_requested(g,i,e) .and. .not. convergenceFlag_backup(g,i,e)) & ! for any pertubation mode: if central solution did not converge... crystallite_requested(g,i,e) .and. .not. convergenceFlag_backup(g,i,e)) & ! for any pertubation mode: if central solution did not converge...
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! ...use (elastic) fallback crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! ...use (elastic) fallback
enddo elementLooping enddo elementLooping2
! --- RESTORE --- ! --- RESTORE ---
@ -1592,8 +1623,8 @@ subroutine crystallite_integrateStateRKCK45()
constitutive_microstructure constitutive_microstructure
implicit none implicit none
!*** local variables ***! integer(pInt) :: &
integer(pInt) e, & ! element index in element loop e, & ! element index in element loop
i, & ! integration point index in ip loop i, & ! integration point index in ip loop
g, & ! grain index in grain loop g, & ! grain index in grain loop
n, & ! stage index in integration stage loop n, & ! stage index in integration stage loop
@ -1604,58 +1635,31 @@ subroutine crystallite_integrateStateRKCK45()
gIter ! bounds for grain iteration gIter ! bounds for grain iteration
real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
RKCK45dotTemperature ! evolution of Temperature of each grain for Runge Kutta Cash Karp integration RKCK45dotTemperature ! evolution of Temperature of each grain for Runge Kutta Cash Karp integration
real(pReal), dimension(5,5) :: a ! coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6) real(pReal), dimension(5,5), parameter :: A = reshape([&
real(pReal), dimension(6) :: b, db ! coefficients in Butcher tableau (used for final integration and error estimate) .2_pReal, .075_pReal, .3_pReal, -11.0_pReal/54.0_pReal, 1631.0_pReal/55296.0_pReal, &
real(pReal), dimension(5) :: c ! coefficients in Butcher tableau (fractions of original time step in stages 2 to 6) .0_pReal, .225_pReal, -.9_pReal, 2.5_pReal, 175.0_pReal/512.0_pReal, &
.0_pReal, .0_pReal, 1.2_pReal, -70.0_pReal/27.0_pReal, 575.0_pReal/13824.0_pReal, &
.0_pReal, .0_pReal, .0_pReal, 35.0_pReal/27.0_pReal, 44275.0_pReal/110592.0_pReal, &
.0_pReal, .0_pReal, .0_pReal, .0_pReal, 253.0_pReal/4096.0_pReal], &
[5,5], order=[2,1]) !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
real(pReal), dimension(6), parameter :: B = &
[37.0_pReal/378.0_pReal, .0_pReal, 250.0_pReal/621.0_pReal, &
125.0_pReal/594.0_pReal, .0_pReal, 512.0_pReal/1771.0_pReal], & !< coefficients in Butcher tableau (used for final integration and error estimate)
DB = B - &
[2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,&
13525.0_pReal/55296.0_pReal,277.0_pReal/14336.0_pReal, 0.25_pReal] !< coefficients in Butcher tableau (used for final integration and error estimate)
real(pReal), dimension(5), parameter :: C = &
[0.2_pReal, 0.3_pReal, 0.6_pReal, 1.0_pReal, 0.875_pReal] !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
stateResiduum, & ! residuum from evolution in micrstructure stateResiduum, & ! residuum from evolution in micrstructure
relStateResiduum ! relative residuum from evolution in microstructure relStateResiduum ! relative residuum from evolution in microstructure
real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
temperatureResiduum, & ! residuum from evolution in temperature temperatureResiduum, & ! residuum from evolution in temperature
relTemperatureResiduum ! relative residuum from evolution in temperature relTemperatureResiduum ! relative residuum from evolution in temperature
logical singleRun ! flag indicating computation for single (g,i,e) triple logical :: singleRun ! flag indicating computation for single (g,i,e) triple
! --- FILL BUTCHER TABLEAU ---
a = 0.0_pReal
b = 0.0_pReal
db = 0.0_pReal
c = 0.0_pReal
a(1,1) = 0.2_pReal
a(1,2) = 0.075_pReal
a(2,2) = 0.225_pReal
a(1,3) = 0.3_pReal
a(2,3) = -0.9_pReal
a(3,3) = 1.2_pReal
a(1,4) = -11.0_pReal / 54.0_pReal
a(2,4) = 2.5_pReal
a(3,4) = -70.0_pReal / 27.0_pReal
a(4,4) = 35.0_pReal / 27.0_pReal
a(1,5) = 1631.0_pReal / 55296.0_pReal
a(2,5) = 175.0_pReal / 512.0_pReal
a(3,5) = 575.0_pReal / 13824.0_pReal
a(4,5) = 44275.0_pReal / 110592.0_pReal
a(5,5) = 253.0_pReal / 4096.0_pReal
b(1) = 37.0_pReal / 378.0_pReal
b(3) = 250.0_pReal / 621.0_pReal
b(4) = 125.0_pReal / 594.0_pReal
b(6) = 512.0_pReal / 1771.0_pReal
db(1) = b(1) - 2825.0_pReal / 27648.0_pReal
db(3) = b(3) - 18575.0_pReal / 48384.0_pReal
db(4) = b(4) - 13525.0_pReal / 55296.0_pReal
db(5) = - 277.0_pReal / 14336.0_pReal
db(6) = b(6) - 0.25_pReal
c(1) = 0.2_pReal
c(2) = 0.3_pReal
c(3) = 0.6_pReal
c(4) = 1.0_pReal
c(5) = 0.875_pReal
! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP --- ! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP ---
eIter = FEsolving_execElem(1:2) eIter = FEsolving_execElem(1:2)
@ -2352,8 +2356,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged' write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged'
write(6,*)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
if (.not. singleRun) then ! if not requesting Integration of just a single IP if (.not. singleRun) then ! if not requesting Integration of just a single IP
@ -2637,7 +2640,7 @@ subroutine crystallite_integrateStateFPI()
!$OMP DO !$OMP 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 !< @ToDo: Put in loop above?
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e) crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
@ -2800,14 +2803,10 @@ subroutine crystallite_integrateStateFPI()
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> update state at el ip g ',e,i,g write(6,'(a,i8,1x,i2,1x,i3,/)') '<< CRYST >> update state at el ip g ',e,i,g
write(6,*) write(6,'(a,f6.1,/)') '<< CRYST >> statedamper ',statedamper
write(6,'(a,f6.1)') '<< CRYST >> statedamper ',statedamper write(6,'(a,/,(12x,12(e12.6,1x)),/)') '<< CRYST >> state residuum',stateResiduum(1:mySizeDotState)
write(6,*) write(6,'(a,/,(12x,12(e12.6,1x)),/)') '<< CRYST >> new state',tempState(1:mySizeDotState)
write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> state residuum',stateResiduum(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state',tempState(1:mySizeDotState)
write(6,*)
endif endif
#endif #endif
@ -2880,10 +2879,10 @@ subroutine crystallite_integrateStateFPI()
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL(write2out) !$OMP CRITICAL(write2out)
write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)),' grains converged after non-local check' write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after state integration no. ',& ' grains converged after non-local check'
NiterationState write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_todo(:,:,:)), &
write(6,*) ' grains todo after state integration no. ', NiterationState
!$OMP END CRITICAL(write2out) !$OMP END CRITICAL(write2out)
endif endif
@ -2895,8 +2894,9 @@ end subroutine crystallite_integrateStateFPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates a jump in the state according to the current state and the current stress !> @brief calculates a jump in the state according to the current state and the current stress
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function crystallite_stateJump(g,i,e) logical function crystallite_stateJump(g,i,e)
use debug, only: debug_level, & use debug, only: &
debug_level, &
debug_crystallite, & debug_crystallite, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective, &
@ -2919,11 +2919,7 @@ function crystallite_stateJump(g,i,e)
i, & ! integration point index i, & ! integration point index
g ! grain index g ! grain index
!*** output variables ***! integer(pInt) :: mySizeDotState
logical crystallite_stateJump
!*** local variables ***!
integer(pInt) mySizeDotState
crystallite_stateJump = .false. crystallite_stateJump = .false.
@ -2943,12 +2939,9 @@ function crystallite_stateJump(g,i,e)
.and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> update state at el ip g ',e,i,g write(6,'(a,i8,1x,i2,1x,i3, /)') '<< CRYST >> update state at el ip g ',e,i,g
write(6,*) write(6,'(a,/,(12x,12(e12.6,1x)),/)') '<< CRYST >> deltaState', constitutive_deltaState(g,i,e)%p(1:mySizeDotState)
write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> deltaState', constitutive_deltaState(g,i,e)%p(1:mySizeDotState) write(6,'(a,/,(12x,12(e12.6,1x)),/)') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
write(6,*)
endif endif
#endif #endif
@ -2957,12 +2950,11 @@ function crystallite_stateJump(g,i,e)
end function crystallite_stateJump end function crystallite_stateJump
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction !> intermediate acceleration of the Newton-Raphson correction
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function crystallite_integrateStress(& logical function crystallite_integrateStress(&
g,& ! grain number g,& ! grain number
i,& ! integration point number i,& ! integration point number
e,& ! element number e,& ! element number
@ -3013,9 +3005,6 @@ function crystallite_integrateStress(&
g ! grain index g ! grain index
real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep
!*** output variables ***!
logical crystallite_integrateStress ! flag indicating if integration suceeded
!*** local variables ***! !*** local variables ***!
real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep
Fp_current, & ! plastic deformation gradient at start of timestep Fp_current, & ! plastic deformation gradient at start of timestep
@ -3406,17 +3395,17 @@ subroutine crystallite_orientations
neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure
if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me
crystallite_disorientation(:,n,1,i,e) = & crystallite_disorientation(:,n,1,i,e) = &
math_qDisorientation(crystallite_orientation(1:4,1,i,e), & math_qDisorientation( crystallite_orientation(1:4,1,i,e), &
crystallite_orientation(1:4,1,neighboring_i,neighboring_e), & crystallite_orientation(1:4,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 endif
else ! for neighbor with local plasticity else ! for neighbor with local plasticity
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
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
enddo enddo
@ -3471,11 +3460,9 @@ function crystallite_postResults(&
g ! grain index g ! grain index
real(pReal), intent(in):: dt ! time increment real(pReal), intent(in):: dt ! time increment
!*** output variables ***!
real(pReal), dimension(1+crystallite_sizePostResults(microstructure_crystallite(mesh_element(4,e)))+ & real(pReal), dimension(1+crystallite_sizePostResults(microstructure_crystallite(mesh_element(4,e)))+ &
1+constitutive_sizePostResults(g,i,e)) :: crystallite_postResults 1+constitutive_sizePostResults(g,i,e)) :: crystallite_postResults
!*** local variables ***!
real(pReal), dimension(3,3) :: Ee real(pReal), dimension(3,3) :: Ee
real(pReal), dimension(4) :: rotation real(pReal), dimension(4) :: rotation
real(pReal) detF real(pReal) detF