2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-18 16:46:26 +05:30
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
2016-02-29 18:56:06 +05:30
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
2013-02-22 04:38:36 +05:30
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
2016-02-29 18:56:06 +05:30
!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Chen Zhang, Michigan State University
2013-02-22 04:38:36 +05:30
!> @brief crystallite state integration functions and reporting of results
!--------------------------------------------------------------------------------------------------
2010-10-01 17:48:49 +05:30
2012-08-31 01:56:28 +05:30
module crystallite
2019-04-03 16:24:07 +05:30
use prec , only : &
2019-04-06 10:01:02 +05:30
pReal , &
pStringLen
2019-04-03 16:24:07 +05:30
use rotations , only : &
rotation
use FEsolving , only : &
FEsolving_execElem , &
FEsolving_execIP
use material , only : &
homogenization_Ngrains
2019-04-03 21:54:15 +05:30
use future
2019-04-03 16:24:07 +05:30
2019-04-11 11:05:58 +05:30
implicit none
2019-04-03 16:24:07 +05:30
private
character ( len = 64 ) , dimension ( : , : ) , allocatable , private :: &
crystallite_output !< name of each post result output
integer , public , protected :: &
crystallite_maxSizePostResults !< description not available
integer , dimension ( : ) , allocatable , public , protected :: &
crystallite_sizePostResults !< description not available
integer , dimension ( : , : ) , allocatable , private :: &
crystallite_sizePostResult !< description not available
real ( pReal ) , dimension ( : , : , : ) , allocatable , public :: &
crystallite_dt !< requested time increment of each grain
real ( pReal ) , dimension ( : , : , : ) , allocatable , private :: &
crystallite_subdt , & !< substepped time increment of each grain
crystallite_subFrac , & !< already calculated fraction of increment
crystallite_subStep !< size of next integration step
type ( rotation ) , dimension ( : , : , : ) , allocatable , private :: &
crystallite_orientation , & !< orientation
crystallite_orientation0 !< initial orientation
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , public , protected :: &
crystallite_Fe , & !< current "elastic" def grad (end of converged time step)
crystallite_P !< 1st Piola-Kirchhoff stress per grain
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , public :: &
crystallite_S , & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
crystallite_S0 , & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_partionedS0 , & !< 2nd Piola-Kirchhoff stress vector at start of homog inc
crystallite_Fp , & !< current plastic def grad (end of converged time step)
crystallite_Fp0 , & !< plastic def grad at start of FE inc
crystallite_partionedFp0 , & !< plastic def grad at start of homog inc
crystallite_Fi , & !< current intermediate def grad (end of converged time step)
crystallite_Fi0 , & !< intermediate def grad at start of FE inc
crystallite_partionedFi0 , & !< intermediate def grad at start of homog inc
crystallite_F0 , & !< def grad at start of FE inc
crystallite_partionedF , & !< def grad to be reached at end of homog inc
crystallite_partionedF0 , & !< def grad at start of homog inc
crystallite_Lp , & !< current plastic velocitiy grad (end of converged time step)
crystallite_Lp0 , & !< plastic velocitiy grad at start of FE inc
crystallite_partionedLp0 , & !< plastic velocity grad at start of homog inc
crystallite_Li , & !< current intermediate velocitiy grad (end of converged time step)
crystallite_Li0 , & !< intermediate velocitiy grad at start of FE inc
crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , private :: &
crystallite_subS0 , & !< 2nd Piola-Kirchhoff stress vector at start of crystallite inc
crystallite_invFp , & !< inverse of current plastic def grad (end of converged time step)
crystallite_subFp0 , & !< plastic def grad at start of crystallite inc
crystallite_invFi , & !< inverse of current intermediate def grad (end of converged time step)
crystallite_subFi0 , & !< intermediate def grad at start of crystallite inc
crystallite_subF , & !< def grad to be reached at end of crystallite inc
crystallite_subF0 , & !< def grad at start of crystallite inc
crystallite_subLp0 , & !< plastic velocity grad at start of crystallite inc
crystallite_subLi0 !< intermediate velocity grad at start of crystallite inc
real ( pReal ) , dimension ( : , : , : , : , : , : , : ) , allocatable , public :: &
crystallite_dPdF !< current individual dPdF per grain (end of converged time step)
logical , dimension ( : , : , : ) , allocatable , public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
logical , dimension ( : , : , : ) , allocatable , private :: &
crystallite_converged , & !< convergence flag
crystallite_todo , & !< flag to indicate need for further computation
crystallite_localPlasticity !< indicates this grain to have purely local constitutive law
enum , bind ( c )
enumerator :: undefined_ID , &
phase_ID , &
texture_ID , &
volume_ID , &
orientation_ID , &
grainrotation_ID , &
defgrad_ID , &
fe_ID , &
fp_ID , &
fi_ID , &
lp_ID , &
li_ID , &
p_ID , &
s_ID , &
elasmatrix_ID , &
neighboringip_ID , &
neighboringelement_ID
end enum
integer ( kind ( undefined_ID ) ) , dimension ( : , : ) , allocatable , private :: &
crystallite_outputID !< ID of each post result output
2019-04-06 15:36:34 +05:30
type , private :: tOutput !< new requested output (per phase)
character ( len = 65536 ) , allocatable , dimension ( : ) :: &
label
end type tOutput
type ( tOutput ) , allocatable , dimension ( : ) , private :: output_constituent
2019-04-11 14:57:03 +05:30
type , private :: tNumerics
integer :: &
iJacoLpresiduum , & !< frequency of Jacobian update of residuum in Lp
nState , & !< state loop limit
nStress !< stress loop limit
2019-04-11 10:54:04 +05:30
real ( pReal ) :: &
subStepMinCryst , & !< minimum (relative) size of sub-step allowed during cutback
subStepSizeCryst , & !< size of first substep when cutback
subStepSizeLp , & !< size of first substep when cutback in Lp calculation
subStepSizeLi , & !< size of first substep when cutback in Li calculation
stepIncreaseCryst , & !< increase of next substep size when previous substep converged
rTol_crystalliteState , & !< relative tolerance in state loop
rTol_crystalliteStress , & !< relative tolerance in stress loop
aTol_crystalliteStress !< absolute tolerance in stress loop
end type tNumerics
type ( tNumerics ) :: num ! numerics parameters. Better name?
2019-04-03 16:24:07 +05:30
procedure ( ) , pointer :: integrateState
public :: &
crystallite_init , &
crystallite_stress , &
crystallite_stressTangent , &
crystallite_orientations , &
crystallite_push33ToRef , &
2019-04-06 10:01:02 +05:30
crystallite_postResults , &
crystallite_results
2019-04-03 16:24:07 +05:30
private :: &
integrateStress , &
integrateState , &
integrateStateFPI , &
integrateStateEuler , &
integrateStateAdaptiveEuler , &
integrateStateRK4 , &
integrateStateRKCK45 , &
stateJump
2014-08-26 20:14:32 +05:30
2012-08-31 01:56:28 +05:30
contains
2009-05-07 21:57:36 +05:30
2011-03-29 12:57:19 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief allocates and initialize per grain variables
!--------------------------------------------------------------------------------------------------
2014-10-10 17:58:57 +05:30
subroutine crystallite_init
2018-09-20 09:39:02 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
use debug , only : &
debug_info , &
debug_reset , &
debug_level , &
debug_crystallite , &
debug_levelBasic
2018-09-20 09:39:02 +05:30
#endif
2019-04-03 16:24:07 +05:30
use numerics , only : &
numerics_integrator , &
worldrank , &
usePingPong
use math , only : &
math_I3 , &
math_EulerToR , &
math_inv33
use mesh , only : &
theMesh , &
mesh_element
use IO , only : &
IO_stringValue , &
IO_write_jobFile , &
IO_error
use material
use config , only : &
config_deallocate , &
config_crystallite , &
2019-04-11 10:54:04 +05:30
config_numerics , &
2019-04-06 10:01:02 +05:30
config_phase , &
2019-04-03 16:24:07 +05:30
crystallite_name
use constitutive , only : &
constitutive_initialFi , &
constitutive_microstructure ! derived (shortcut) quantities of given state
integer , parameter :: FILEUNIT = 434
logical , dimension ( : , : ) , allocatable :: devNull
integer :: &
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e , & !< counter in element loop
o = 0 , & !< counter in output loop
r , &
cMax , & !< maximum number of integration point components
iMax , & !< maximum number of integration points
eMax , & !< maximum number of elements
myNcomponents , & !< number of components at current IP
mySize
character ( len = 65536 ) , dimension ( : ) , allocatable :: str
write ( 6 , '(/,a)' ) ' <<<+- crystallite init -+>>>'
cMax = homogenization_maxNgrains
iMax = theMesh % elem % nIPs
eMax = theMesh % nElems
allocate ( crystallite_S0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_partionedS0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_S ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subS0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_P ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_F0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_partionedF0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_partionedF ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subF0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subF ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_Fp0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_partionedFp0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subFp0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_Fp ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_invFp ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_Fi0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_partionedFi0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subFi0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_Fi ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_invFi ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_Fe ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_Lp0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_partionedLp0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subLp0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_Lp ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_Li0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_partionedLi0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subLi0 ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_Li ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_dPdF ( 3 , 3 , 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_dt ( cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subdt ( cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subFrac ( cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subStep ( cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_orientation ( cMax , iMax , eMax ) )
allocate ( crystallite_orientation0 ( cMax , iMax , eMax ) )
allocate ( crystallite_localPlasticity ( cMax , iMax , eMax ) , source = . true . )
allocate ( crystallite_requested ( cMax , iMax , eMax ) , source = . false . )
allocate ( crystallite_todo ( cMax , iMax , eMax ) , source = . false . )
allocate ( crystallite_converged ( cMax , iMax , eMax ) , source = . true . )
allocate ( crystallite_output ( maxval ( crystallite_Noutput ) , &
size ( config_crystallite ) ) ) ; crystallite_output = ''
allocate ( crystallite_outputID ( maxval ( crystallite_Noutput ) , &
size ( config_crystallite ) ) , source = undefined_ID )
allocate ( crystallite_sizePostResults ( size ( config_crystallite ) ) , source = 0 )
allocate ( crystallite_sizePostResult ( maxval ( crystallite_Noutput ) , &
size ( config_crystallite ) ) , source = 0 )
2019-04-11 10:54:04 +05:30
2019-04-11 14:57:03 +05:30
num % subStepMinCryst = config_numerics % getFloat ( 'substepmincryst' , defaultVal = 1.0e-3_pReal )
num % subStepSizeCryst = config_numerics % getFloat ( 'substepsizecryst' , defaultVal = 0.25_pReal )
num % stepIncreaseCryst = config_numerics % getFloat ( 'stepincreasecryst' , defaultVal = 1.5_pReal )
2019-04-11 10:54:04 +05:30
2019-04-11 14:57:03 +05:30
num % subStepSizeLp = config_numerics % getFloat ( 'substepsizelp' , defaultVal = 0.5_pReal )
num % subStepSizeLi = config_numerics % getFloat ( 'substepsizeli' , defaultVal = 0.5_pReal )
num % rTol_crystalliteState = config_numerics % getFloat ( 'rtol_crystallitestate' , defaultVal = 1.0e-6_pReal )
num % rTol_crystalliteStress = config_numerics % getFloat ( 'rtol_crystallitestress' , defaultVal = 1.0e-6_pReal )
num % aTol_crystalliteStress = config_numerics % getFloat ( 'atol_crystallitestress' , defaultVal = 1.0e-8_pReal )
num % iJacoLpresiduum = config_numerics % getInt ( 'ijacolpresiduum' , defaultVal = 1 )
num % nState = config_numerics % getInt ( 'nstate' , defaultVal = 20 )
num % nStress = config_numerics % getInt ( 'nstress' , defaultVal = 40 )
if ( num % subStepMinCryst < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepMinCryst' )
if ( num % subStepSizeCryst < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepSizeCryst' )
if ( num % stepIncreaseCryst < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'stepIncreaseCryst' )
if ( num % subStepSizeLp < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepSizeLp' )
if ( num % subStepSizeLi < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepSizeLi' )
if ( num % rTol_crystalliteState < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'rTol_crystalliteState' )
if ( num % rTol_crystalliteStress < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'rTol_crystalliteStress' )
if ( num % aTol_crystalliteStress < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'aTol_crystalliteStress' )
if ( num % iJacoLpresiduum < 1 ) call IO_error ( 301 , ext_msg = 'iJacoLpresiduum' )
if ( num % nState < 1 ) call IO_error ( 301 , ext_msg = 'nState' )
if ( num % nStress < 1 ) call IO_error ( 301 , ext_msg = 'nStress' )
2019-04-03 16:24:07 +05:30
select case ( numerics_integrator )
case ( 1 )
integrateState = > integrateStateFPI
case ( 2 )
integrateState = > integrateStateEuler
case ( 3 )
integrateState = > integrateStateAdaptiveEuler
case ( 4 )
integrateState = > integrateStateRK4
case ( 5 )
integrateState = > integrateStateRKCK45
end select
do c = 1 , size ( config_crystallite )
2018-06-22 11:33:22 +05:30
#if defined(__GFORTRAN__)
2019-04-03 16:24:07 +05:30
str = [ 'GfortranBug86277' ]
str = config_crystallite ( c ) % getStrings ( '(output)' , defaultVal = str )
if ( str ( 1 ) == 'GfortranBug86277' ) str = [ character ( len = 65536 ) :: ]
2018-06-22 11:33:22 +05:30
#else
2019-04-03 16:24:07 +05:30
str = config_crystallite ( c ) % getStrings ( '(output)' , defaultVal = [ character ( len = 65536 ) :: ] )
2018-06-22 11:33:22 +05:30
#endif
2019-04-03 16:24:07 +05:30
do o = 1 , size ( str )
crystallite_output ( o , c ) = str ( o )
outputName : select case ( str ( o ) )
case ( 'phase' ) outputName
crystallite_outputID ( o , c ) = phase_ID
case ( 'texture' ) outputName
crystallite_outputID ( o , c ) = texture_ID
case ( 'volume' ) outputName
crystallite_outputID ( o , c ) = volume_ID
case ( 'orientation' ) outputName
crystallite_outputID ( o , c ) = orientation_ID
case ( 'grainrotation' ) outputName
crystallite_outputID ( o , c ) = grainrotation_ID
case ( 'defgrad' , 'f' ) outputName ! ToDo: no alias (f only)
crystallite_outputID ( o , c ) = defgrad_ID
case ( 'fe' ) outputName
crystallite_outputID ( o , c ) = fe_ID
case ( 'fp' ) outputName
crystallite_outputID ( o , c ) = fp_ID
case ( 'fi' ) outputName
crystallite_outputID ( o , c ) = fi_ID
case ( 'lp' ) outputName
crystallite_outputID ( o , c ) = lp_ID
case ( 'li' ) outputName
crystallite_outputID ( o , c ) = li_ID
case ( 'p' , 'firstpiola' , '1stpiola' ) outputName ! ToDo: no alias (p only)
crystallite_outputID ( o , c ) = p_ID
case ( 's' , 'tstar' , 'secondpiola' , '2ndpiola' ) outputName ! ToDo: no alias (s only)
crystallite_outputID ( o , c ) = s_ID
case ( 'elasmatrix' ) outputName
crystallite_outputID ( o , c ) = elasmatrix_ID
case ( 'neighboringip' ) outputName ! ToDo: this is not a result, it is static. Should be written out by mesh
crystallite_outputID ( o , c ) = neighboringip_ID
case ( 'neighboringelement' ) outputName ! ToDo: this is not a result, it is static. Should be written out by mesh
crystallite_outputID ( o , c ) = neighboringelement_ID
case default outputName
call IO_error ( 105 , ext_msg = trim ( str ( o ) ) / / ' (Crystallite)' )
end select outputName
enddo
enddo
2019-04-06 10:01:02 +05:30
2019-04-06 15:36:34 +05:30
allocate ( output_constituent ( size ( config_phase ) ) )
2019-04-06 10:01:02 +05:30
do c = 1 , size ( config_phase )
#if defined(__GFORTRAN__)
2019-04-06 15:36:34 +05:30
allocate ( output_constituent ( c ) % label ( 1 ) )
output_constituent ( c ) % label ( 1 ) = 'GfortranBug86277'
output_constituent ( c ) % label = config_phase ( c ) % getStrings ( '(output)' , defaultVal = output_constituent ( c ) % label )
if ( output_constituent ( c ) % label ( 1 ) == 'GfortranBug86277' ) output_constituent ( c ) % label = [ character ( len = pStringLen ) :: ]
2019-04-06 10:01:02 +05:30
#else
2019-04-06 15:36:34 +05:30
output_constituent ( c ) % label = config_phase ( c ) % getStrings ( '(output)' , defaultVal = [ character ( len = pStringLen ) :: ] )
2019-04-06 10:01:02 +05:30
#endif
enddo
2019-04-03 16:24:07 +05:30
do r = 1 , size ( config_crystallite )
do o = 1 , crystallite_Noutput ( r )
select case ( crystallite_outputID ( o , r ) )
case ( phase_ID , texture_ID , volume_ID )
mySize = 1
case ( orientation_ID , grainrotation_ID )
mySize = 4
case ( defgrad_ID , fe_ID , fp_ID , fi_ID , lp_ID , li_ID , p_ID , s_ID )
mySize = 9
case ( elasmatrix_ID )
mySize = 36
case ( neighboringip_ID , neighboringelement_ID )
mySize = theMesh % elem % nIPneighbors
case default
mySize = 0
end select
crystallite_sizePostResult ( o , r ) = mySize
crystallite_sizePostResults ( r ) = crystallite_sizePostResults ( r ) + mySize
enddo
enddo
crystallite_maxSizePostResults = &
maxval ( crystallite_sizePostResults ( microstructure_crystallite ) , microstructure_active )
2018-06-27 00:03:41 +05:30
2009-06-16 14:33:30 +05:30
2013-04-29 16:47:30 +05:30
!--------------------------------------------------------------------------------------------------
2010-02-25 23:09:11 +05:30
! write description file for crystallite output
2019-04-03 16:24:07 +05:30
if ( worldrank == 0 ) then
call IO_write_jobFile ( FILEUNIT , 'outputCrystallite' )
do r = 1 , size ( config_crystallite )
if ( any ( microstructure_crystallite ( mesh_element ( 4 , : ) ) == r ) ) then
write ( FILEUNIT , '(/,a,/)' ) '[' / / trim ( crystallite_name ( r ) ) / / ']'
do o = 1 , crystallite_Noutput ( r )
write ( FILEUNIT , '(a,i4)' ) trim ( crystallite_output ( o , r ) ) / / char ( 9 ) , crystallite_sizePostResult ( o , r )
enddo
endif
enddo
close ( FILEUNIT )
endif
2019-04-06 10:01:02 +05:30
call config_deallocate ( 'material.config/phase' )
2019-04-03 16:24:07 +05:30
call config_deallocate ( 'material.config/crystallite' )
2018-06-27 00:03:41 +05:30
2013-04-29 16:47:30 +05:30
!--------------------------------------------------------------------------------------------------
! initialize
2019-01-18 16:46:26 +05:30
!$OMP PARALLEL DO PRIVATE(myNcomponents,i,c)
2019-04-03 16:24:07 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
myNcomponents = homogenization_Ngrains ( mesh_element ( 3 , e ) )
2019-04-03 16:47:21 +05:30
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e ) ; do c = 1 , myNcomponents
2019-04-03 16:24:07 +05:30
crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e ) = math_EulerToR ( material_EulerAngles ( 1 : 3 , c , i , e ) ) ! plastic def gradient reflects init orientation
crystallite_Fi0 ( 1 : 3 , 1 : 3 , c , i , e ) = constitutive_initialFi ( c , i , e )
crystallite_F0 ( 1 : 3 , 1 : 3 , c , i , e ) = math_I3
crystallite_localPlasticity ( c , i , e ) = phase_localPlasticity ( material_phase ( c , i , e ) )
crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) = math_inv33 ( matmul ( crystallite_Fi0 ( 1 : 3 , 1 : 3 , c , i , e ) , &
2019-04-03 21:54:15 +05:30
crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e ) ) ) ! assuming that euler angles are given in internal strain free configuration
2019-04-03 16:24:07 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_requested ( c , i , e ) = . true .
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-04-03 16:24:07 +05:30
enddo
!$OMP END PARALLEL DO
if ( any ( . not . crystallite_localPlasticity ) . and . . not . usePingPong ) call IO_error ( 601 ) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal?
crystallite_partionedFp0 = crystallite_Fp0
crystallite_partionedFi0 = crystallite_Fi0
crystallite_partionedF0 = crystallite_F0
crystallite_partionedF = crystallite_F0
call crystallite_orientations ( )
crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations
!$OMP PARALLEL DO
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
do c = 1 , homogenization_Ngrains ( mesh_element ( 3 , e ) )
call constitutive_microstructure ( crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) , &
c , i , e ) ! update dependent state variables to be consistent with basic states
enddo
2019-01-19 03:39:46 +05:30
enddo
2019-04-03 16:24:07 +05:30
enddo
!$OMP END PARALLEL DO
devNull = crystallite_stress ( )
call crystallite_stressTangent
2019-01-18 19:12:44 +05:30
2018-09-20 09:39:02 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) then
write ( 6 , '(a42,1x,i10)' ) ' # of elements: ' , eMax
write ( 6 , '(a42,1x,i10)' ) 'max # of integration points/element: ' , iMax
write ( 6 , '(a42,1x,i10)' ) 'max # of constituents/integration point: ' , cMax
write ( 6 , '(a42,1x,i10)' ) 'max # of neigbours/integration point: ' , theMesh % elem % nIPneighbors
write ( 6 , '(a42,1x,i10)' ) ' # of nonlocal constituents: ' , count ( . not . crystallite_localPlasticity )
flush ( 6 )
endif
call debug_info
call debug_reset
2018-09-20 09:39:02 +05:30
#endif
2010-11-03 22:52:48 +05:30
2012-03-09 01:55:28 +05:30
end subroutine crystallite_init
2009-05-07 21:57:36 +05:30
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
2019-02-17 16:45:46 +05:30
function crystallite_stress ( dummyArgumentToPreventInternalCompilerErrorWithGCC )
2019-04-03 16:24:07 +05:30
use prec , only : &
tol_math_check , &
dNeq0
2019-01-18 19:12:44 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
use debug , only : &
debug_level , &
debug_crystallite , &
debug_levelBasic , &
debug_levelExtensive , &
debug_levelSelective , &
debug_e , &
debug_i , &
debug_g
2019-01-18 19:12:44 +05:30
#endif
2019-04-03 16:24:07 +05:30
use IO , only : &
IO_warning , &
IO_error
use math , only : &
math_inv33
use mesh , only : &
theMesh , &
mesh_element
use material , only : &
homogenization_Ngrains , &
plasticState , &
sourceState , &
phase_Nsources , &
phaseAt , phasememberAt
logical , dimension ( theMesh % elem % nIPs , theMesh % Nelems ) :: crystallite_stress
real ( pReal ) , intent ( in ) , optional :: &
dummyArgumentToPreventInternalCompilerErrorWithGCC
real ( pReal ) :: &
formerSubStep
integer :: &
NiterationCrystallite , & ! number of iterations in crystallite loop
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e , & !< counter in element loop
startIP , endIP , &
s
2019-01-18 19:12:44 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 &
. and . FEsolving_execElem ( 1 ) < = debug_e &
. and . debug_e < = FEsolving_execElem ( 2 ) ) then
write ( 6 , '(/,a,i8,1x,i2,1x,i3)' ) '<< CRYST stress >> boundary and initial values at el ip ipc ' , &
debug_e , debug_i , debug_g
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> F ' , &
transpose ( crystallite_partionedF ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> F0 ' , &
transpose ( crystallite_partionedF0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> Fp0' , &
transpose ( crystallite_partionedFp0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> Fi0' , &
transpose ( crystallite_partionedFi0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> Lp0' , &
transpose ( crystallite_partionedLp0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> Li0' , &
transpose ( crystallite_partionedLi0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
endif
2019-01-18 19:12:44 +05:30
#endif
!--------------------------------------------------------------------------------------------------
! initialize to starting condition
2019-04-03 16:24:07 +05:30
crystallite_subStep = 0.0_pReal
!$OMP PARALLEL DO
elementLooping1 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e ) ; do c = 1 , homogenization_Ngrains ( mesh_element ( 3 , e ) )
homogenizationRequestsCalculation : if ( crystallite_requested ( c , i , e ) ) then
plasticState ( phaseAt ( c , i , e ) ) % subState0 ( : , phasememberAt ( c , i , e ) ) = &
plasticState ( phaseAt ( c , i , e ) ) % partionedState0 ( : , phasememberAt ( c , i , e ) )
do s = 1 , phase_Nsources ( phaseAt ( c , i , e ) )
sourceState ( phaseAt ( c , i , e ) ) % p ( s ) % subState0 ( : , phasememberAt ( c , i , e ) ) = &
sourceState ( phaseAt ( c , i , e ) ) % p ( s ) % partionedState0 ( : , phasememberAt ( c , i , e ) )
enddo
crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedFp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subLp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedLp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedFi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subLi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedLi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedF0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subS0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedS0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subFrac ( c , i , e ) = 0.0_pReal
2019-04-11 10:54:04 +05:30
crystallite_subStep ( c , i , e ) = 1.0_pReal / num % subStepSizeCryst
2019-04-03 16:24:07 +05:30
crystallite_todo ( c , i , e ) = . true .
crystallite_converged ( c , i , e ) = . false . ! pretend failed step of 1/subStepSizeCryst
endif homogenizationRequestsCalculation
enddo ; enddo
enddo elementLooping1
!$OMP END PARALLEL DO
singleRun : if ( FEsolving_execELem ( 1 ) == FEsolving_execElem ( 2 ) . and . &
FEsolving_execIP ( 1 , FEsolving_execELem ( 1 ) ) == FEsolving_execIP ( 2 , FEsolving_execELem ( 1 ) ) ) then
startIP = FEsolving_execIP ( 1 , FEsolving_execELem ( 1 ) )
endIP = startIP
else singleRun
startIP = 1
endIP = theMesh % elem % nIPs
endif singleRun
NiterationCrystallite = 0
cutbackLooping : do while ( any ( crystallite_todo ( : , startIP : endIP , FEsolving_execELem ( 1 ) : FEsolving_execElem ( 2 ) ) ) )
NiterationCrystallite = NiterationCrystallite + 1
2019-02-27 00:17:46 +05:30
2019-01-18 19:12:44 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 ) &
write ( 6 , '(a,i6)' ) '<< CRYST stress >> crystallite iteration ' , NiterationCrystallite
2019-01-18 19:12:44 +05:30
#endif
2019-04-03 16:24:07 +05:30
!$OMP PARALLEL DO PRIVATE(formerSubStep)
elementLooping3 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
do c = 1 , homogenization_Ngrains ( mesh_element ( 3 , e ) )
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
! wind forward
2019-04-03 16:24:07 +05:30
if ( crystallite_converged ( c , i , e ) ) then
formerSubStep = crystallite_subStep ( c , i , e )
crystallite_subFrac ( c , i , e ) = crystallite_subFrac ( c , i , e ) + crystallite_subStep ( c , i , e )
crystallite_subStep ( c , i , e ) = min ( 1.0_pReal - crystallite_subFrac ( c , i , e ) , &
2019-04-11 10:54:04 +05:30
num % stepIncreaseCryst * crystallite_subStep ( c , i , e ) )
2019-04-03 16:24:07 +05:30
crystallite_todo ( c , i , e ) = crystallite_subStep ( c , i , e ) > 0.0_pReal ! still time left to integrate on?
if ( crystallite_todo ( c , i , e ) ) then
crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subLp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subLi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subS0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_S ( 1 : 3 , 1 : 3 , c , i , e )
!if abbrevation, make c and p private in omp
plasticState ( phaseAt ( c , i , e ) ) % subState0 ( : , phasememberAt ( c , i , e ) ) &
= plasticState ( phaseAt ( c , i , e ) ) % state ( : , phasememberAt ( c , i , e ) )
do s = 1 , phase_Nsources ( phaseAt ( c , i , e ) )
sourceState ( phaseAt ( c , i , e ) ) % p ( s ) % subState0 ( : , phasememberAt ( c , i , e ) ) &
= sourceState ( phaseAt ( c , i , e ) ) % p ( s ) % state ( : , phasememberAt ( c , i , e ) )
enddo
2019-01-18 19:12:44 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 &
. and . ( ( e == debug_e . and . i == debug_i . and . c == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) &
write ( 6 , '(a,f12.8,a,f12.8,a,i8,1x,i2,1x,i3,/)' ) '<< CRYST stress >> winding forward from ' , &
crystallite_subFrac ( c , i , e ) - formerSubStep , ' to current crystallite_subfrac ' , &
crystallite_subFrac ( c , i , e ) , ' in crystallite_stress at el ip ipc ' , e , i , c
2019-01-18 19:12:44 +05:30
#endif
2019-04-03 16:24:07 +05:30
endif
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore)
2019-04-03 16:24:07 +05:30
else
2019-04-11 10:54:04 +05:30
crystallite_subStep ( c , i , e ) = num % subStepSizeCryst * crystallite_subStep ( c , i , e )
2019-04-03 16:24:07 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) = math_inv33 ( crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) )
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) = math_inv33 ( crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) )
crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_S0 ( 1 : 3 , 1 : 3 , c , i , e )
if ( crystallite_subStep ( c , i , e ) < 1.0_pReal ) then ! actual (not initial) cutback
crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subLp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subLi0 ( 1 : 3 , 1 : 3 , c , i , e )
endif
plasticState ( phaseAt ( c , i , e ) ) % state ( : , phasememberAt ( c , i , e ) ) &
= plasticState ( phaseAt ( c , i , e ) ) % subState0 ( : , phasememberAt ( c , i , e ) )
do s = 1 , phase_Nsources ( phaseAt ( c , i , e ) )
sourceState ( phaseAt ( c , i , e ) ) % p ( s ) % state ( : , phasememberAt ( c , i , e ) ) &
= sourceState ( phaseAt ( c , i , e ) ) % p ( s ) % subState0 ( : , phasememberAt ( c , i , e ) )
enddo
! cant restore dotState here, since not yet calculated in first cutback after initialization
2019-04-11 10:54:04 +05:30
crystallite_todo ( c , i , e ) = crystallite_subStep ( c , i , e ) > num % subStepMinCryst ! still on track or already done (beyond repair)
2019-01-18 19:12:44 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( e == debug_e . and . i == debug_i . and . c == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
if ( crystallite_todo ( c , i , e ) ) then
write ( 6 , '(a,f12.8,a,i8,1x,i2,1x,i3,/)' ) '<< CRYST stress >> cutback with new crystallite_subStep: ' , &
crystallite_subStep ( c , i , e ) , ' at el ip ipc ' , e , i , c
else
write ( 6 , '(a,i8,1x,i2,1x,i3,/)' ) '<< CRYST stress >> reached minimum step size at el ip ipc ' , e , i , c
endif
endif
2019-01-18 19:12:44 +05:30
#endif
2019-04-03 16:24:07 +05:30
endif
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
! prepare for integration
2019-04-03 16:24:07 +05:30
if ( crystallite_todo ( c , i , e ) ) then
crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) &
+ crystallite_subStep ( c , i , e ) * ( crystallite_partionedF ( 1 : 3 , 1 : 3 , c , i , e ) &
- crystallite_partionedF0 ( 1 : 3 , 1 : 3 , c , i , e ) )
crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) = matmul ( matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
2019-04-03 21:54:15 +05:30
crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) , &
crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) )
2019-04-03 16:24:07 +05:30
crystallite_subdt ( c , i , e ) = crystallite_subStep ( c , i , e ) * crystallite_dt ( c , i , e )
crystallite_converged ( c , i , e ) = . false .
endif
enddo
enddo
enddo elementLooping3
!$OMP END PARALLEL DO
2019-01-18 19:12:44 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 ) then
write ( 6 , '(/,a,f8.5,a,f8.5,/)' ) '<< CRYST stress >> ' , minval ( crystallite_subStep ) , &
' ≤ subStep ≤ ' , maxval ( crystallite_subStep )
write ( 6 , '(/,a,f8.5,a,f8.5,/)' ) '<< CRYST stress >> ' , minval ( crystallite_subFrac ) , &
' ≤ subFrac ≤ ' , maxval ( crystallite_subFrac )
flush ( 6 )
if ( iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) then
write ( 6 , '(/,a,f8.5,1x,a,1x,f8.5,1x,a)' ) '<< CRYST stress >> subFrac + subStep = ' , &
crystallite_subFrac ( debug_g , debug_i , debug_e ) , '+' , crystallite_subStep ( debug_g , debug_i , debug_e ) , '@selective'
flush ( 6 )
endif
endif
2019-01-18 19:12:44 +05:30
#endif
!--------------------------------------------------------------------------------------------------
! integrate --- requires fully defined state array (basic + dependent state)
2019-04-03 16:24:07 +05:30
if ( any ( crystallite_todo ) ) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation
2019-04-11 10:54:04 +05:30
where ( . not . crystallite_converged . and . crystallite_subStep > num % subStepMinCryst ) & ! do not try non-converged but fully cutbacked any further
2019-04-03 16:24:07 +05:30
crystallite_todo = . true . ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation
2019-01-18 19:12:44 +05:30
2019-04-03 16:24:07 +05:30
enddo cutbackLooping
2019-01-18 19:12:44 +05:30
! return whether converged or not
2019-04-03 16:24:07 +05:30
crystallite_stress = . false .
elementLooping5 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
crystallite_stress ( i , e ) = all ( crystallite_converged ( : , i , e ) )
enddo
enddo elementLooping5
2019-01-18 19:12:44 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
elementLooping6 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
do c = 1 , homogenization_Ngrains ( mesh_element ( 3 , e ) )
if ( . not . crystallite_converged ( c , i , e ) ) then
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) &
write ( 6 , '(a,i8,1x,i2,1x,i3,/)' ) '<< CRYST stress >> no convergence at el ip ipc ' , &
e , i , c
endif
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( e == debug_e . and . i == debug_i . and . c == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
write ( 6 , '(a,i8,1x,i2,1x,i3)' ) '<< CRYST stress >> solution at el ip ipc ' , e , i , c
write ( 6 , '(/,a,/,3(12x,3(f12.4,1x)/))' ) '<< CRYST stress >> P / MPa' , &
transpose ( crystallite_P ( 1 : 3 , 1 : 3 , c , i , e ) ) * 1.0e-6_pReal
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> Fp' , &
transpose ( crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> Fi' , &
transpose ( crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/),/)' ) '<< CRYST stress >> Lp' , &
transpose ( crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/),/)' ) '<< CRYST stress >> Li' , &
transpose ( crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e ) )
flush ( 6 )
endif
enddo
enddo
enddo elementLooping6
2019-01-18 19:12:44 +05:30
#endif
end function crystallite_stress
2019-01-18 16:46:26 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate tangent (dPdF)
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine crystallite_stressTangent
2019-04-03 16:24:07 +05:30
use prec , only : &
tol_math_check , &
dNeq0
use IO , only : &
IO_warning , &
IO_error
use math , only : &
math_inv33 , &
math_identity2nd , &
math_3333to99 , &
math_99to3333 , &
math_I3 , &
math_mul3333xx3333 , &
math_mul33xx33 , &
math_invert2 , &
math_det33
use mesh , only : &
mesh_element
use material , only : &
homogenization_Ngrains
use constitutive , only : &
constitutive_SandItsTangents , &
constitutive_LpAndItsTangents , &
constitutive_LiAndItsTangents
integer :: &
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e , & !< counter in element loop
o , &
p
real ( pReal ) , dimension ( 3 , 3 ) :: temp_33_1 , devNull , invSubFi0 , temp_33_2 , temp_33_3 , temp_33_4
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dSdFe , &
dSdF , &
dSdFi , &
dLidS , &
dLidFi , &
dLpdS , &
dLpdFi , &
dFidS , &
dFpinvdF , &
rhs_3333 , &
lhs_3333 , &
temp_3333
real ( pReal ) , dimension ( 9 , 9 ) :: temp_99
logical :: error
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,invSubFi0,o,p, &
!$OMP rhs_3333,lhs_3333,temp_99,temp_33_1,temp_33_2,temp_33_3,temp_33_4,temp_3333,error)
elementLooping : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
do c = 1 , homogenization_Ngrains ( mesh_element ( 3 , e ) )
call constitutive_SandItsTangents ( devNull , dSdFe , dSdFi , &
crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) , c , i , e ) ! call constitutive law to calculate elastic stress tangent
call constitutive_LiAndItsTangents ( devNull , dLidS , dLidFi , &
crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) , &
c , i , e ) ! call constitutive law to calculate Li tangent in lattice configuration
if ( sum ( abs ( dLidS ) ) < tol_math_check ) then
dFidS = 0.0_pReal
else
invSubFi0 = math_inv33 ( crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) )
lhs_3333 = 0.0_pReal ; rhs_3333 = 0.0_pReal
do o = 1 , 3 ; do p = 1 , 3
lhs_3333 ( 1 : 3 , 1 : 3 , o , p ) = lhs_3333 ( 1 : 3 , 1 : 3 , o , p ) &
+ crystallite_subdt ( c , i , e ) * matmul ( invSubFi0 , dLidFi ( 1 : 3 , 1 : 3 , o , p ) )
lhs_3333 ( 1 : 3 , o , 1 : 3 , p ) = lhs_3333 ( 1 : 3 , o , 1 : 3 , p ) &
+ crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) * crystallite_invFi ( p , o , c , i , e )
rhs_3333 ( 1 : 3 , 1 : 3 , o , p ) = rhs_3333 ( 1 : 3 , 1 : 3 , o , p ) &
- crystallite_subdt ( c , i , e ) * matmul ( invSubFi0 , dLidS ( 1 : 3 , 1 : 3 , o , p ) )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-04-03 16:24:07 +05:30
call math_invert2 ( temp_99 , error , math_3333to99 ( lhs_3333 ) )
if ( error ) then
call IO_warning ( warning_ID = 600 , el = e , ip = i , g = c , &
ext_msg = 'inversion error in analytic tangent calculation' )
dFidS = 0.0_pReal
else
dFidS = math_mul3333xx3333 ( math_99to3333 ( temp_99 ) , rhs_3333 )
endif
dLidS = math_mul3333xx3333 ( dLidFi , dFidS ) + dLidS
endif
call constitutive_LpAndItsTangents ( devNull , dLpdS , dLpdFi , &
crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) , c , i , e ) ! call constitutive law to calculate Lp tangent in lattice configuration
dLpdS = math_mul3333xx3333 ( dLpdFi , dFidS ) + dLpdS
2019-01-18 16:46:26 +05:30
!--------------------------------------------------------------------------------------------------
! calculate dSdF
2019-04-03 16:24:07 +05:30
temp_33_1 = transpose ( matmul ( crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) ) )
temp_33_2 = matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
math_inv33 ( crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) ) )
temp_33_3 = matmul ( matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) , &
math_inv33 ( crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) ) )
2019-04-03 16:47:21 +05:30
do o = 1 , 3 ; do p = 1 , 3
2019-04-03 16:24:07 +05:30
rhs_3333 ( p , o , 1 : 3 , 1 : 3 ) = matmul ( dSdFe ( p , o , 1 : 3 , 1 : 3 ) , temp_33_1 )
temp_3333 ( 1 : 3 , 1 : 3 , p , o ) = matmul ( matmul ( temp_33_2 , dLpdS ( 1 : 3 , 1 : 3 , p , o ) ) , &
crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) ) &
+ matmul ( temp_33_3 , dLidS ( 1 : 3 , 1 : 3 , p , o ) )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-04-03 16:24:07 +05:30
lhs_3333 = crystallite_subdt ( c , i , e ) * math_mul3333xx3333 ( dSdFe , temp_3333 ) &
+ math_mul3333xx3333 ( dSdFi , dFidS )
call math_invert2 ( temp_99 , error , math_identity2nd ( 9 ) + math_3333to99 ( lhs_3333 ) )
if ( error ) then
call IO_warning ( warning_ID = 600 , el = e , ip = i , g = c , &
ext_msg = 'inversion error in analytic tangent calculation' )
dSdF = rhs_3333
else
dSdF = math_mul3333xx3333 ( math_99to3333 ( temp_99 ) , rhs_3333 )
endif
2019-01-18 16:46:26 +05:30
!--------------------------------------------------------------------------------------------------
! calculate dFpinvdF
2019-04-03 16:24:07 +05:30
temp_3333 = math_mul3333xx3333 ( dLpdS , dSdF )
2019-04-03 16:47:21 +05:30
do o = 1 , 3 ; do p = 1 , 3
2019-04-03 16:24:07 +05:30
dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) &
= - crystallite_subdt ( c , i , e ) &
* matmul ( math_inv33 ( crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) ) , &
matmul ( temp_3333 ( 1 : 3 , 1 : 3 , p , o ) , crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) ) )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-01-18 16:46:26 +05:30
!--------------------------------------------------------------------------------------------------
! assemble dPdF
2019-04-03 16:24:07 +05:30
temp_33_1 = matmul ( crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) , &
matmul ( crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) , &
transpose ( crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) ) )
temp_33_2 = matmul ( crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) , &
transpose ( crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) )
temp_33_3 = matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) )
temp_33_4 = matmul ( matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) , &
crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) )
crystallite_dPdF ( 1 : 3 , 1 : 3 , 1 : 3 , 1 : 3 , c , i , e ) = 0.0_pReal
2019-04-03 16:47:21 +05:30
do p = 1 , 3
2019-04-03 16:24:07 +05:30
crystallite_dPdF ( p , 1 : 3 , p , 1 : 3 , c , i , e ) = transpose ( temp_33_1 )
enddo
2019-04-03 16:47:21 +05:30
do o = 1 , 3 ; do p = 1 , 3
2019-04-03 16:24:07 +05:30
crystallite_dPdF ( 1 : 3 , 1 : 3 , p , o , c , i , e ) = crystallite_dPdF ( 1 : 3 , 1 : 3 , p , o , c , i , e ) + &
matmul ( matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) , temp_33_2 ) + &
matmul ( matmul ( temp_33_3 , dSdF ( 1 : 3 , 1 : 3 , p , o ) ) , transpose ( crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) ) + &
matmul ( temp_33_4 , transpose ( dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-04-03 16:24:07 +05:30
enddo ; enddo
enddo elementLooping
!$OMP END PARALLEL DO
2019-01-18 16:46:26 +05:30
end subroutine crystallite_stressTangent
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-19 14:05:45 +05:30
!> @brief calculates orientations
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-19 14:05:45 +05:30
subroutine crystallite_orientations
2019-04-03 16:24:07 +05:30
use math , only : &
math_rotationalPart33 , &
math_RtoQ
use material , only : &
plasticState , &
material_phase , &
homogenization_Ngrains
use mesh , only : &
mesh_element
use lattice , only : &
lattice_qDisorientation
use plastic_nonlocal , only : &
plastic_nonlocal_updateCompatibility
2019-01-19 14:05:45 +05:30
2019-04-03 16:24:07 +05:30
integer &
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e !< counter in element loop
!$OMP PARALLEL DO
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
do c = 1 , homogenization_Ngrains ( mesh_element ( 3 , e ) )
call crystallite_orientation ( c , i , e ) % fromRotationMatrix ( transpose ( math_rotationalPart33 ( crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) ) ) )
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
nonlocalPresent : if ( any ( plasticState % nonLocal ) ) then
!$OMP PARALLEL DO
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
if ( plasticState ( material_phase ( 1 , i , e ) ) % nonLocal ) & ! if nonlocal model
call plastic_nonlocal_updateCompatibility ( crystallite_orientation , i , e )
enddo ; enddo
!$OMP END PARALLEL DO
endif nonlocalPresent
2014-08-26 20:14:32 +05:30
2019-01-19 14:05:45 +05:30
end subroutine crystallite_orientations
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Map 2nd order tensor to reference config
!--------------------------------------------------------------------------------------------------
function crystallite_push33ToRef ( ipc , ip , el , tensor33 )
2019-04-03 16:24:07 +05:30
use math , only : &
math_inv33 , &
math_EulerToR
use material , only : &
material_EulerAngles ! ToDo: Why stored? We also have crystallite_orientation0
real ( pReal ) , dimension ( 3 , 3 ) :: crystallite_push33ToRef
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: tensor33
real ( pReal ) , dimension ( 3 , 3 ) :: T
integer , intent ( in ) :: &
el , &
ip , &
ipc
T = matmul ( math_EulerToR ( material_EulerAngles ( 1 : 3 , ipc , ip , el ) ) , &
transpose ( math_inv33 ( crystallite_subF ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) ) )
crystallite_push33ToRef = matmul ( transpose ( T ) , matmul ( tensor33 , T ) )
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
end function crystallite_push33ToRef
2014-08-26 20:14:32 +05:30
2019-01-19 14:05:45 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief return results of particular grain
!--------------------------------------------------------------------------------------------------
function crystallite_postResults ( ipc , ip , el )
use math , only : &
math_det33 , &
math_I3 , &
2019-03-09 05:03:11 +05:30
inDeg
2019-01-19 14:05:45 +05:30
use mesh , only : &
2019-02-02 15:05:10 +05:30
theMesh , &
2019-01-19 14:05:45 +05:30
mesh_element , &
mesh_ipVolume , &
2019-02-02 15:05:10 +05:30
mesh_ipNeighborhood
2019-01-19 14:05:45 +05:30
use material , only : &
plasticState , &
sourceState , &
microstructure_crystallite , &
crystallite_Noutput , &
material_phase , &
material_texture , &
homogenization_Ngrains
use constitutive , only : &
constitutive_homogenizedC , &
constitutive_postResults
2019-02-02 00:58:07 +05:30
use rotations , only : &
rotation
2019-01-19 14:05:45 +05:30
2019-04-03 16:02:30 +05:30
integer , intent ( in ) :: &
2019-01-19 14:05:45 +05:30
el , & !< element index
ip , & !< integration point index
ipc !< grain index
real ( pReal ) , dimension ( 1 + crystallite_sizePostResults ( microstructure_crystallite ( mesh_element ( 4 , el ) ) ) + &
1 + plasticState ( material_phase ( ipc , ip , el ) ) % sizePostResults + &
sum ( sourceState ( material_phase ( ipc , ip , el ) ) % p ( : ) % sizePostResults ) ) :: &
crystallite_postResults
real ( pReal ) :: &
detF
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-19 14:05:45 +05:30
o , &
c , &
crystID , &
mySize , &
n
2019-02-02 00:58:07 +05:30
type ( rotation ) :: rot
2019-01-19 14:05:45 +05:30
crystID = microstructure_crystallite ( mesh_element ( 4 , el ) )
crystallite_postResults = 0.0_pReal
2019-01-28 15:58:46 +05:30
crystallite_postResults ( 1 ) = real ( crystallite_sizePostResults ( crystID ) , pReal ) ! header-like information (length)
2019-04-03 16:02:30 +05:30
c = 1
2019-01-19 14:05:45 +05:30
2019-04-03 16:02:30 +05:30
do o = 1 , crystallite_Noutput ( crystID )
mySize = 0
2019-01-19 14:05:45 +05:30
select case ( crystallite_outputID ( o , crystID ) )
case ( phase_ID )
2019-04-03 16:02:30 +05:30
mySize = 1
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 ) = real ( material_phase ( ipc , ip , el ) , pReal ) ! phaseID of grain
case ( texture_ID )
2019-04-03 16:02:30 +05:30
mySize = 1
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 ) = real ( material_texture ( ipc , ip , el ) , pReal ) ! textureID of grain
case ( volume_ID )
2019-04-03 16:02:30 +05:30
mySize = 1
2019-01-19 14:05:45 +05:30
detF = math_det33 ( crystallite_partionedF ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) ! V_current = det(F) * V_reference
crystallite_postResults ( c + 1 ) = detF * mesh_ipVolume ( ip , el ) &
/ real ( homogenization_Ngrains ( mesh_element ( 3 , el ) ) , pReal ) ! grain volume (not fraction but absolute)
case ( orientation_ID )
2019-04-03 16:02:30 +05:30
mySize = 4
2019-02-02 00:58:07 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = crystallite_orientation ( ipc , ip , el ) % asQuaternion ( )
2019-01-19 14:05:45 +05:30
case ( grainrotation_ID )
2019-02-02 00:58:07 +05:30
rot = crystallite_orientation0 ( ipc , ip , el ) % misorientation ( crystallite_orientation ( ipc , ip , el ) )
2019-04-03 16:02:30 +05:30
mySize = 4
2019-02-02 00:58:07 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = rot % asAxisAnglePair ( )
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 4 ) = inDeg * crystallite_postResults ( c + 4 ) ! angle in degree
! remark: tensor output is of the form 11,12,13, 21,22,23, 31,32,33
! thus row index i is slow, while column index j is fast. reminder: "row is slow"
case ( defgrad_ID )
2019-04-03 16:02:30 +05:30
mySize = 9
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_partionedF ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( fe_ID )
2019-04-03 16:02:30 +05:30
mySize = 9
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_Fe ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( fp_ID )
2019-04-03 16:02:30 +05:30
mySize = 9
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_Fp ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( fi_ID )
2019-04-03 16:02:30 +05:30
mySize = 9
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( lp_ID )
2019-04-03 16:02:30 +05:30
mySize = 9
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_Lp ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( li_ID )
2019-04-03 16:02:30 +05:30
mySize = 9
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_Li ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( p_ID )
2019-04-03 16:02:30 +05:30
mySize = 9
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_P ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( s_ID )
2019-04-03 16:02:30 +05:30
mySize = 9
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = &
2019-03-09 05:03:11 +05:30
reshape ( crystallite_S ( 1 : 3 , 1 : 3 , ipc , ip , el ) , [ mySize ] )
2019-01-19 14:05:45 +05:30
case ( elasmatrix_ID )
2019-04-03 16:02:30 +05:30
mySize = 36
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = reshape ( constitutive_homogenizedC ( ipc , ip , el ) , [ mySize ] )
case ( neighboringelement_ID )
2019-02-02 15:05:10 +05:30
mySize = theMesh % elem % nIPneighbors
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = 0.0_pReal
2019-04-03 16:02:30 +05:30
forall ( n = 1 : mySize ) &
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + n ) = real ( mesh_ipNeighborhood ( 1 , n , ip , el ) , pReal )
case ( neighboringip_ID )
2019-02-02 15:05:10 +05:30
mySize = theMesh % elem % nIPneighbors
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : c + mySize ) = 0.0_pReal
2019-04-03 16:02:30 +05:30
forall ( n = 1 : mySize ) &
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + n ) = real ( mesh_ipNeighborhood ( 2 , n , ip , el ) , pReal )
end select
c = c + mySize
enddo
crystallite_postResults ( c + 1 ) = real ( plasticState ( material_phase ( ipc , ip , el ) ) % sizePostResults , pReal ) ! size of constitutive results
2019-04-03 16:02:30 +05:30
c = c + 1
if ( size ( crystallite_postResults ) - c > 0 ) &
2019-01-19 14:05:45 +05:30
crystallite_postResults ( c + 1 : size ( crystallite_postResults ) ) = &
2019-03-09 05:03:11 +05:30
constitutive_postResults ( crystallite_S ( 1 : 3 , 1 : 3 , ipc , ip , el ) , crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
2019-02-26 12:24:03 +05:30
ipc , ip , el )
2019-01-19 14:05:45 +05:30
end function crystallite_postResults
2019-04-06 10:01:02 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-13 13:17:56 +05:30
!> @brief writes crystallite results to HDF5 output file
2019-04-06 10:01:02 +05:30
!--------------------------------------------------------------------------------------------------
subroutine crystallite_results
#if defined(PETSc) || defined(DAMASK_HDF5)
2019-04-07 18:17:21 +05:30
use lattice
2019-04-06 10:01:02 +05:30
use results
use HDF5_utilities
2019-04-07 17:58:08 +05:30
use rotations
2019-04-06 10:01:02 +05:30
use config , only : &
config_name_phase = > phase_name ! anticipate logical name
2019-04-06 15:36:34 +05:30
integer :: p , o
2019-04-07 17:58:08 +05:30
real ( pReal ) , allocatable , dimension ( : , : , : ) :: selected_tensors
type ( rotation ) , allocatable , dimension ( : ) :: selected_rotations
2019-04-07 18:17:21 +05:30
character ( len = 256 ) :: group , lattice_label
2019-04-06 10:01:02 +05:30
do p = 1 , size ( config_name_phase )
2019-04-18 15:25:50 +05:30
group = trim ( 'current/constituent' ) / / '/' / / trim ( config_name_phase ( p ) ) / / '/generic'
call HDF5_closeGroup ( results_addGroup ( group ) )
2019-04-13 13:17:56 +05:30
2019-04-06 15:36:34 +05:30
do o = 1 , size ( output_constituent ( p ) % label )
select case ( output_constituent ( p ) % label ( o ) )
case ( 'f' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_partionedF , p )
call results_writeDataset ( group , selected_tensors , 'F' , &
2019-04-06 15:36:34 +05:30
'deformation gradient' , '1' )
case ( 'fe' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Fe , p )
call results_writeDataset ( group , selected_tensors , 'Fe' , &
2019-04-06 15:36:34 +05:30
'elastic deformation gradient' , '1' )
case ( 'fp' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Fp , p )
call results_writeDataset ( group , selected_tensors , 'Fp' , &
2019-04-06 15:36:34 +05:30
'plastic deformation gradient' , '1' )
case ( 'fi' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Fi , p )
call results_writeDataset ( group , selected_tensors , 'Fi' , &
2019-04-06 15:36:34 +05:30
'inelastic deformation gradient' , '1' )
case ( 'lp' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Lp , p )
call results_writeDataset ( group , selected_tensors , 'Lp' , &
2019-04-06 15:36:34 +05:30
'plastic velocity gradient' , '1/s' )
case ( 'li' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Li , p )
call results_writeDataset ( group , selected_tensors , 'Li' , &
2019-04-06 15:36:34 +05:30
'inelastic velocity gradient' , '1/s' )
case ( 'p' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_P , p )
call results_writeDataset ( group , selected_tensors , 'P' , &
2019-04-06 15:36:34 +05:30
'1st Piola-Kirchoff stress' , 'Pa' )
case ( 's' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_S , p )
call results_writeDataset ( group , selected_tensors , 'S' , &
'2nd Piola-Kirchoff stress' , 'Pa' )
case ( 'orientation' )
2019-04-07 18:17:21 +05:30
select case ( lattice_structure ( p ) )
case ( LATTICE_iso_ID )
lattice_label = 'iso'
case ( LATTICE_fcc_ID )
lattice_label = 'fcc'
case ( LATTICE_bcc_ID )
lattice_label = 'bcc'
case ( LATTICE_bct_ID )
lattice_label = 'bct'
case ( LATTICE_hex_ID )
lattice_label = 'hex'
case ( LATTICE_ort_ID )
lattice_label = 'ort'
end select
2019-04-07 17:58:08 +05:30
selected_rotations = select_rotations ( crystallite_orientation , p )
call results_writeDataset ( group , selected_rotations , 'orientation' , &
2019-04-07 18:17:21 +05:30
'crystal orientation as quaternion' , lattice_label )
2019-04-06 15:36:34 +05:30
end select
enddo
2019-04-06 10:01:02 +05:30
enddo
contains
2019-04-07 17:58:08 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief select tensors for output
!--------------------------------------------------------------------------------------------------
2019-04-13 19:04:51 +05:30
function select_tensors ( dataset , instance )
2019-04-07 17:58:08 +05:30
use material , only : &
homogenization_maxNgrains , &
material_phaseAt
2019-04-13 19:04:51 +05:30
integer , intent ( in ) :: instance
real ( pReal ) , dimension ( : , : , : , : , : ) , intent ( in ) :: dataset
real ( pReal ) , allocatable , dimension ( : , : , : ) :: select_tensors
integer :: e , i , c , j
allocate ( select_tensors ( 3 , 3 , count ( material_phaseAt == instance ) * homogenization_maxNgrains ) )
j = 0
do e = 1 , size ( material_phaseAt , 2 )
do i = 1 , homogenization_maxNgrains !ToDo: this needs to be changed for varying Ngrains
do c = 1 , size ( material_phaseAt , 1 )
if ( material_phaseAt ( c , e ) == instance ) then
2019-04-07 17:58:08 +05:30
j = j + 1
2019-04-13 19:04:51 +05:30
select_tensors ( 1 : 3 , 1 : 3 , j ) = dataset ( 1 : 3 , 1 : 3 , c , i , e )
endif
enddo
enddo
2019-04-07 17:58:08 +05:30
enddo
2019-04-13 19:04:51 +05:30
end function select_tensors
2019-04-06 10:01:02 +05:30
2019-04-07 17:58:08 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief select rotations for output
!--------------------------------------------------------------------------------------------------
2019-04-13 19:04:51 +05:30
function select_rotations ( dataset , instance )
2019-04-06 15:36:34 +05:30
use material , only : &
homogenization_maxNgrains , &
material_phaseAt
2019-04-06 10:01:02 +05:30
2019-04-13 19:04:51 +05:30
integer , intent ( in ) :: instance
type ( rotation ) , dimension ( : , : , : ) , intent ( in ) :: dataset
type ( rotation ) , allocatable , dimension ( : ) :: select_rotations
integer :: e , i , c , j
allocate ( select_rotations ( count ( material_phaseAt == instance ) * homogenization_maxNgrains ) )
j = 0
do e = 1 , size ( material_phaseAt , 2 )
do i = 1 , homogenization_maxNgrains !ToDo: this needs to be changed for varying Ngrains
do c = 1 , size ( material_phaseAt , 1 )
if ( material_phaseAt ( c , e ) == instance ) then
2019-04-06 15:36:34 +05:30
j = j + 1
2019-04-13 19:04:51 +05:30
select_rotations ( j ) = dataset ( c , i , e )
endif
enddo
enddo
enddo
2019-04-06 10:01:02 +05:30
2019-04-07 17:58:08 +05:30
end function select_rotations
2019-04-06 10:01:02 +05:30
#endif
end subroutine crystallite_results
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction
!--------------------------------------------------------------------------------------------------
2019-04-03 16:24:07 +05:30
logical function integrateStress ( ipc , ip , el , timeFraction )
use , intrinsic :: &
IEEE_arithmetic
use prec , only : tol_math_check , &
dEq0
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
use debug , only : debug_level , &
debug_e , &
debug_i , &
debug_g , &
debug_crystallite , &
debug_levelBasic , &
debug_levelExtensive , &
debug_levelSelective
2019-01-15 08:57:57 +05:30
#endif
2014-08-26 20:14:32 +05:30
2019-04-03 16:24:07 +05:30
use constitutive , only : constitutive_LpAndItsTangents , &
constitutive_LiAndItsTangents , &
constitutive_SandItsTangents
use math , only : math_mul33xx33 , &
math_mul3333xx3333 , &
math_inv33 , &
math_det33 , &
math_I3 , &
math_identity2nd , &
math_3333to99 , &
math_33to9 , &
math_9to33
integer , intent ( in ) :: el , & ! element index
ip , & ! integration point index
ipc ! grain index
real ( pReal ) , optional , intent ( in ) :: timeFraction ! fraction of timestep
real ( pReal ) , dimension ( 3 , 3 ) :: Fg_new , & ! deformation gradient at end of timestep
Fp_new , & ! plastic deformation gradient at end of timestep
Fe_new , & ! elastic deformation gradient at end of timestep
invFp_new , & ! inverse of Fp_new
Fi_new , & ! gradient of intermediate deformation stages
invFi_new , &
invFp_current , & ! inverse of Fp_current
invFi_current , & ! inverse of Fp_current
Lpguess , & ! current guess for plastic velocity gradient
Lpguess_old , & ! known last good guess for plastic velocity gradient
Lp_constitutive , & ! plastic velocity gradient resulting from constitutive law
residuumLp , & ! current residuum of plastic velocity gradient
residuumLp_old , & ! last residuum of plastic velocity gradient
deltaLp , & ! direction of next guess
Liguess , & ! current guess for intermediate velocity gradient
Liguess_old , & ! known last good guess for intermediate velocity gradient
Li_constitutive , & ! intermediate velocity gradient resulting from constitutive law
residuumLi , & ! current residuum of intermediate velocity gradient
residuumLi_old , & ! last residuum of intermediate velocity gradient
deltaLi , & ! direction of next guess
S , & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
A , &
B , &
Fe , & ! elastic deformation gradient
temp_33
real ( pReal ) , dimension ( 9 ) :: work ! needed for matrix inversion by LAPACK
integer , dimension ( 9 ) :: devNull ! needed for matrix inversion by LAPACK
real ( pReal ) , dimension ( 9 , 9 ) :: dRLp_dLp , & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme)
dRLp_dLp2 , & ! working copy of dRdLp
dRLi_dLi ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme)
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dS_dFe , & ! partial derivative of 2nd Piola-Kirchhoff stress
dS_dFi , &
dFe_dLp , & ! partial derivative of elastic deformation gradient
dFe_dLi , &
dFi_dLi , &
dLp_dFi , &
dLi_dFi , &
dLp_dS , &
dLi_dS
real ( pReal ) detInvFi , & ! determinant of InvFi
steplengthLp , &
steplengthLi , &
dt , & ! time increment
aTolLp , &
aTolLi
integer NiterationStressLp , & ! number of stress integrations
NiterationStressLi , & ! number of inner stress integrations
ierr , & ! error indicator for LAPACK
o , &
p , &
jacoCounterLp , &
jacoCounterLi ! counters to check for Jacobian update
external :: &
dgesv
!* be pessimistic
integrateStress = . false .
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) &
write ( 6 , '(a,i8,1x,i2,1x,i3)' ) '<< CRYST integrateStress >> at el ip ipc ' , el , ip , ipc
2019-01-15 08:57:57 +05:30
#endif
2014-07-02 17:57:39 +05:30
2019-04-03 16:24:07 +05:30
if ( present ( timeFraction ) ) then
dt = crystallite_subdt ( ipc , ip , el ) * timeFraction
Fg_new = crystallite_subF0 ( 1 : 3 , 1 : 3 , ipc , ip , el ) &
+ ( crystallite_subF ( 1 : 3 , 1 : 3 , ipc , ip , el ) - crystallite_subF0 ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) * timeFraction
else
dt = crystallite_subdt ( ipc , ip , el )
Fg_new = crystallite_subF ( 1 : 3 , 1 : 3 , ipc , ip , el )
endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!* feed local variables
2019-04-03 16:24:07 +05:30
Lpguess = crystallite_Lp ( 1 : 3 , 1 : 3 , ipc , ip , el ) ! ... and take it as first guess
Liguess = crystallite_Li ( 1 : 3 , 1 : 3 , ipc , ip , el ) ! ... and take it as first guess
Liguess_old = Liguess
invFp_current = math_inv33 ( crystallite_subFp0 ( 1 : 3 , 1 : 3 , ipc , ip , el ) )
failedInversionFp : if ( all ( dEq0 ( invFp_current ) ) ) then
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) &
write ( 6 , '(a,i8,1x,i2,1x,i3)' ) '<< CRYST integrateStress >> failed on inversion of current Fp at el ip ipc ' , &
el , ip , ipc
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) > 0 ) &
write ( 6 , '(/,a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST >> current Fp ' , transpose ( crystallite_subFp0 ( 1 : 3 , 1 : 3 , ipc , ip , el ) )
2019-01-15 08:57:57 +05:30
#endif
2019-04-03 16:24:07 +05:30
return
endif failedInversionFp
A = matmul ( Fg_new , invFp_current ) ! intermediate tensor needed later to calculate dFe_dLp
2014-08-26 20:14:32 +05:30
2019-04-03 16:24:07 +05:30
invFi_current = math_inv33 ( crystallite_subFi0 ( 1 : 3 , 1 : 3 , ipc , ip , el ) )
failedInversionFi : if ( all ( dEq0 ( invFi_current ) ) ) then
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) &
write ( 6 , '(a,i8,1x,i2,1x,i3)' ) '<< CRYST integrateStress >> failed on inversion of current Fi at el ip ipc ' , &
el , ip , ipc
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) > 0 ) &
write ( 6 , '(/,a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST integrateStress >> current Fi ' , &
transpose ( crystallite_subFi0 ( 1 : 3 , 1 : 3 , ipc , ip , el ) )
2019-01-15 08:57:57 +05:30
#endif
2019-04-03 16:24:07 +05:30
return
endif failedInversionFi
!* start Li loop with normal step length
NiterationStressLi = 0
jacoCounterLi = 0
steplengthLi = 1.0_pReal
residuumLi_old = 0.0_pReal
LiLoop : do
NiterationStressLi = NiterationStressLi + 1
2019-04-11 14:57:03 +05:30
LiLoopLimit : if ( NiterationStressLi > num % nStress ) then
2018-02-16 20:06:18 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) &
2019-04-11 14:57:03 +05:30
write ( 6 , '(a,i3,a,i8,1x,i2,1x,i3,/)' ) '<< CRYST integrateStress >> reached Li loop limit' , num % nStress , &
2019-01-19 03:39:46 +05:30
' at el ip ipc ' , el , ip , ipc
2018-09-20 09:54:03 +05:30
#endif
2019-04-03 16:24:07 +05:30
return
endif LiLoopLimit
invFi_new = matmul ( invFi_current , math_I3 - dt * Liguess )
Fi_new = math_inv33 ( invFi_new )
detInvFi = math_det33 ( invFi_new )
!* start Lp loop with normal step length
NiterationStressLp = 0
jacoCounterLp = 0
steplengthLp = 1.0_pReal
residuumLp_old = 0.0_pReal
Lpguess_old = Lpguess
LpLoop : do
NiterationStressLp = NiterationStressLp + 1
2019-04-11 14:57:03 +05:30
LpLoopLimit : if ( NiterationStressLp > num % nStress ) then
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) &
2019-04-11 14:57:03 +05:30
write ( 6 , '(a,i3,a,i8,1x,i2,1x,i3,/)' ) '<< CRYST integrateStress >> reached Lp loop limit' , num % nStress , &
2019-04-03 16:24:07 +05:30
' at el ip ipc ' , el , ip , ipc
2019-01-15 08:57:57 +05:30
#endif
2019-04-03 16:24:07 +05:30
return
endif LpLoopLimit
!* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
B = math_I3 - dt * Lpguess
Fe = matmul ( matmul ( A , B ) , invFi_new )
call constitutive_SandItsTangents ( S , dS_dFe , dS_dFi , &
Fe , Fi_new , ipc , ip , el ) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
!* calculate plastic velocity gradient and its tangent from constitutive law
call constitutive_LpAndItsTangents ( Lp_constitutive , dLp_dS , dLp_dFi , &
S , Fi_new , ipc , ip , el )
2018-09-20 09:54:03 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
write ( 6 , '(a,i3,/)' ) '<< CRYST integrateStress >> iteration ' , NiterationStressLp
write ( 6 , '(a,/,3(12x,3(e20.10,1x)/))' ) '<< CRYST integrateStress >> Lpguess' , transpose ( Lpguess )
write ( 6 , '(a,/,3(12x,3(e20.10,1x)/))' ) '<< CRYST integrateStress >> Lp_constitutive' , transpose ( Lp_constitutive )
write ( 6 , '(a,/,3(12x,3(e20.10,1x)/))' ) '<< CRYST integrateStress >> Fi' , transpose ( Fi_new )
write ( 6 , '(a,/,3(12x,3(e20.10,1x)/))' ) '<< CRYST integrateStress >> Fe' , transpose ( Fe )
write ( 6 , '(a,/,3(12x,3(e20.10,1x)/))' ) '<< CRYST integrateStress >> S' , transpose ( S )
endif
2019-01-15 08:57:57 +05:30
#endif
2014-08-26 20:14:32 +05:30
2019-04-03 16:24:07 +05:30
!* update current residuum and check for convergence of loop
2019-04-11 10:54:04 +05:30
aTolLp = max ( num % rTol_crystalliteStress * max ( norm2 ( Lpguess ) , norm2 ( Lp_constitutive ) ) , & ! absolute tolerance from largest acceptable relative error
num % aTol_crystalliteStress ) ! minimum lower cutoff
2019-04-03 16:24:07 +05:30
residuumLp = Lpguess - Lp_constitutive
if ( any ( IEEE_is_NaN ( residuumLp ) ) ) then
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) &
write ( 6 , '(a,i8,1x,i2,1x,i3,a,i3,a)' ) '<< CRYST integrateStress >> encountered NaN for Lp-residuum at el ip ipc ' , &
el , ip , ipc , &
' ; iteration ' , NiterationStressLp , &
' >> returning..!'
2019-01-15 08:57:57 +05:30
#endif
2019-04-03 16:24:07 +05:30
return ! ...me = .false. to inform integrator about problem
elseif ( norm2 ( residuumLp ) < aTolLp ) then ! converged if below absolute tolerance
exit LpLoop ! ...leave iteration loop
elseif ( NiterationStressLp == 1 &
. or . norm2 ( residuumLp ) < norm2 ( residuumLp_old ) ) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLp_old = residuumLp ! ...remember old values and...
Lpguess_old = Lpguess
steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
2019-04-11 10:54:04 +05:30
steplengthLp = num % subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction
2019-04-03 16:24:07 +05:30
Lpguess = Lpguess_old + steplengthLp * deltaLp
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
write ( 6 , '(a,1x,f7.4)' ) '<< CRYST integrateStress >> linear search for Lpguess with step' , steplengthLp
endif
2019-01-15 08:57:57 +05:30
#endif
2019-04-03 16:24:07 +05:30
cycle LpLoop
endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!* calculate Jacobian for correction term
2019-04-11 14:57:03 +05:30
if ( mod ( jacoCounterLp , num % iJacoLpresiduum ) == 0 ) then
2019-04-03 16:47:21 +05:30
do o = 1 , 3 ; do p = 1 , 3
dFe_dLp ( o , 1 : 3 , p , 1 : 3 ) = A ( o , p ) * transpose ( invFi_new ) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
enddo ; enddo
2019-04-03 16:24:07 +05:30
dFe_dLp = - dt * dFe_dLp
dRLp_dLp = math_identity2nd ( 9 ) &
- math_3333to99 ( math_mul3333xx3333 ( math_mul3333xx3333 ( dLp_dS , dS_dFe ) , dFe_dLp ) )
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
write ( 6 , '(a,/,9(12x,9(e12.4,1x)/))' ) '<< CRYST integrateStress >> dLp_dS' , math_3333to99 ( dLp_dS )
write ( 6 , '(a,1x,e20.10)' ) '<< CRYST integrateStress >> dLp_dS norm' , norm2 ( math_3333to99 ( dLp_dS ) )
write ( 6 , '(a,/,9(12x,9(e12.4,1x)/))' ) '<< CRYST integrateStress >> dRLp_dLp' , dRLp_dLp - math_identity2nd ( 9 )
write ( 6 , '(a,1x,e20.10)' ) '<< CRYST integrateStress >> dRLp_dLp norm' , norm2 ( dRLp_dLp - math_identity2nd ( 9 ) )
endif
2019-01-15 08:57:57 +05:30
#endif
2019-04-03 16:24:07 +05:30
dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine
work = math_33to9 ( residuumLp )
call dgesv ( 9 , 1 , dRLp_dLp2 , 9 , devNull , work , 9 , ierr ) ! solve dRLp/dLp * delta Lp = -res for delta Lp
if ( ierr / = 0 ) then
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) then
write ( 6 , '(a,i8,1x,i2,1x,i3)' ) '<< CRYST integrateStress >> failed on dR/dLp inversion at el ip ipc ' , &
el , ip , ipc
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
write ( 6 , * )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST integrateStress >> dR_dLp' , transpose ( dRLp_dLp )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST integrateStress >> dFe_dLp' , transpose ( math_3333to99 ( dFe_dLp ) )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST integrateStress >> dS_dFe (cnst)' , transpose ( math_3333to99 ( dS_dFe ) )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST integrateStress >> dLp_dS (cnst)' , transpose ( math_3333to99 ( dLp_dS ) )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST integrateStress >> A' , transpose ( A )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST integrateStress >> B' , transpose ( B )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST integrateStress >> Lp_constitutive' , transpose ( Lp_constitutive )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST integrateStress >> Lpguess' , transpose ( Lpguess )
endif
endif
2019-01-15 08:57:57 +05:30
#endif
2019-04-03 16:24:07 +05:30
return
endif
deltaLp = - math_9to33 ( work )
endif
jacoCounterLp = jacoCounterLp + 1
2014-08-26 20:14:32 +05:30
2019-04-03 16:24:07 +05:30
Lpguess = Lpguess + steplengthLp * deltaLp
2014-08-26 20:14:32 +05:30
2019-04-03 16:24:07 +05:30
enddo LpLoop
2014-09-03 01:41:57 +05:30
2019-01-15 08:57:57 +05:30
!* calculate intermediate velocity gradient and its tangent from constitutive law
2019-04-03 16:24:07 +05:30
call constitutive_LiAndItsTangents ( Li_constitutive , dLi_dS , dLi_dFi , &
S , Fi_new , ipc , ip , el )
2014-09-03 01:41:57 +05:30
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST integrateStress >> Li_constitutive' , transpose ( Li_constitutive )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST integrateStress >> Liguess' , transpose ( Liguess )
endif
2019-01-15 08:57:57 +05:30
#endif
2014-08-26 20:14:32 +05:30
2019-04-03 16:24:07 +05:30
!* update current residuum and check for convergence of loop
2019-04-11 10:54:04 +05:30
aTolLi = max ( num % rTol_crystalliteStress * max ( norm2 ( Liguess ) , norm2 ( Li_constitutive ) ) , & ! absolute tolerance from largest acceptable relative error
num % aTol_crystalliteStress ) ! minimum lower cutoff
2019-04-03 16:24:07 +05:30
residuumLi = Liguess - Li_constitutive
if ( any ( IEEE_is_NaN ( residuumLi ) ) ) then ! NaN in residuum...
2019-01-19 03:39:46 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) &
write ( 6 , '(a,i8,1x,i2,1x,i3,a,i3,a)' ) '<< CRYST integrateStress >> encountered NaN for Li-residuum at el ip ipc ' , &
el , ip , ipc , &
' ; iteration ' , NiterationStressLi , &
' >> returning..!'
2019-01-19 03:39:46 +05:30
#endif
2019-04-03 16:24:07 +05:30
return ! ...me = .false. to inform integrator about problem
elseif ( norm2 ( residuumLi ) < aTolLi ) then ! converged if below absolute tolerance
exit LiLoop ! ...leave iteration loop
elseif ( NiterationStressLi == 1 &
. or . norm2 ( residuumLi ) < norm2 ( residuumLi_old ) ) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLi_old = residuumLi ! ...remember old values and...
Liguess_old = Liguess
steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
2019-04-11 10:54:04 +05:30
steplengthLi = num % subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction
2019-04-03 16:24:07 +05:30
Liguess = Liguess_old + steplengthLi * deltaLi
cycle LiLoop
endif
!* calculate Jacobian for correction term
2019-04-11 14:57:03 +05:30
if ( mod ( jacoCounterLi , num % iJacoLpresiduum ) == 0 ) then
2019-04-03 16:24:07 +05:30
temp_33 = matmul ( matmul ( A , B ) , invFi_current )
2019-04-03 16:47:21 +05:30
do o = 1 , 3 ; do p = 1 , 3
2019-04-03 16:24:07 +05:30
dFe_dLi ( 1 : 3 , o , 1 : 3 , p ) = - dt * math_I3 ( o , p ) * temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
dFi_dLi ( 1 : 3 , o , 1 : 3 , p ) = - dt * math_I3 ( o , p ) * invFi_current
2019-04-03 16:47:21 +05:30
enddo ; enddo
do o = 1 , 3 ; do p = 1 , 3
2019-04-03 16:24:07 +05:30
dFi_dLi ( 1 : 3 , 1 : 3 , o , p ) = matmul ( matmul ( Fi_new , dFi_dLi ( 1 : 3 , 1 : 3 , o , p ) ) , Fi_new )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-04-03 16:24:07 +05:30
dRLi_dLi = math_identity2nd ( 9 ) &
- math_3333to99 ( math_mul3333xx3333 ( dLi_dS , math_mul3333xx3333 ( dS_dFe , dFe_dLi ) + &
math_mul3333xx3333 ( dS_dFi , dFi_dLi ) ) ) &
- math_3333to99 ( math_mul3333xx3333 ( dLi_dFi , dFi_dLi ) )
work = math_33to9 ( residuumLi )
call dgesv ( 9 , 1 , dRLi_dLi , 9 , devNull , work , 9 , ierr ) ! solve dRLi/dLp * delta Li = -res for delta Li
if ( ierr / = 0 ) then
2019-01-19 14:05:45 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) then
write ( 6 , '(a,i8,1x,i2,1x,i3)' ) '<< CRYST integrateStress >> failed on dR/dLi inversion at el ip ipc ' , &
el , ip , ipc
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
write ( 6 , * )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST integrateStress >> dR_dLi' , transpose ( dRLi_dLi )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST integrateStress >> dFe_dLi' , transpose ( math_3333to99 ( dFe_dLi ) )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST integrateStress >> dS_dFi (cnst)' , transpose ( math_3333to99 ( dS_dFi ) )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST integrateStress >> dLi_dS (cnst)' , transpose ( math_3333to99 ( dLi_dS ) )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST integrateStress >> Li_constitutive' , transpose ( Li_constitutive )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST integrateStress >> Liguess' , transpose ( Liguess )
endif
endif
2019-01-19 14:05:45 +05:30
#endif
2019-04-03 16:24:07 +05:30
return
endif
deltaLi = - math_9to33 ( work )
endif
jacoCounterLi = jacoCounterLi + 1
Liguess = Liguess + steplengthLi * deltaLi
enddo LiLoop
!* calculate new plastic and elastic deformation gradient
invFp_new = matmul ( invFp_current , B )
invFp_new = invFp_new / math_det33 ( invFp_new ) ** ( 1.0_pReal / 3.0_pReal ) ! regularize
Fp_new = math_inv33 ( invFp_new )
failedInversionInvFp : if ( all ( dEq0 ( Fp_new ) ) ) then
2019-01-19 14:05:45 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) then
write ( 6 , '(a,i8,1x,i2,1x,i3)' ) '<< CRYST integrateStress >> failed on invFp_new inversion at el ip ipc ' , &
el , ip , ipc
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) &
write ( 6 , '(/,a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST integrateStress >> invFp_new' , transpose ( invFp_new )
endif
2019-01-19 14:05:45 +05:30
#endif
2019-04-03 16:24:07 +05:30
return
endif failedInversionInvFp
Fe_new = matmul ( matmul ( Fg_new , invFp_new ) , invFi_new )
2014-08-26 20:14:32 +05:30
2019-01-19 14:05:45 +05:30
!--------------------------------------------------------------------------------------------------
! stress integration was successful
2019-04-03 16:24:07 +05:30
integrateStress = . true .
crystallite_P ( 1 : 3 , 1 : 3 , ipc , ip , el ) = matmul ( matmul ( Fg_new , invFp_new ) , matmul ( S , transpose ( invFp_new ) ) )
crystallite_S ( 1 : 3 , 1 : 3 , ipc , ip , el ) = S
crystallite_Lp ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Lpguess
crystallite_Li ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Liguess
crystallite_Fp ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Fp_new
crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Fi_new
crystallite_Fe ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Fe_new
crystallite_invFp ( 1 : 3 , 1 : 3 , ipc , ip , el ) = invFp_new
crystallite_invFi ( 1 : 3 , 1 : 3 , ipc , ip , el ) = invFi_new
2014-08-26 20:14:32 +05:30
2019-01-19 14:05:45 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
write ( 6 , '(a,/)' ) '<< CRYST integrateStress >> successful integration'
write ( 6 , '(a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST integrateStress >> P / MPa' , &
transpose ( crystallite_P ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) * 1.0e-6_pReal
write ( 6 , '(a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST integrateStress >> Cauchy / MPa' , &
matmul ( crystallite_P ( 1 : 3 , 1 : 3 , ipc , ip , el ) , transpose ( Fg_new ) ) * 1.0e-6_pReal / math_det33 ( Fg_new )
write ( 6 , '(a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST integrateStress >> Fe Lp Fe^-1' , &
transpose ( matmul ( Fe_new , matmul ( crystallite_Lp ( 1 : 3 , 1 : 3 , ipc , ip , el ) , math_inv33 ( Fe_new ) ) ) )
write ( 6 , '(a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST integrateStress >> Fp' , transpose ( crystallite_Fp ( 1 : 3 , 1 : 3 , ipc , ip , el ) )
write ( 6 , '(a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST integrateStress >> Fi' , transpose ( crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) )
endif
2019-01-19 14:05:45 +05:30
#endif
2015-10-14 00:22:01 +05:30
2019-01-19 14:05:45 +05:30
end function integrateStress
2010-10-01 17:48:49 +05:30
2012-05-17 20:55:21 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-11 11:05:58 +05:30
subroutine integrateStateFPI
2019-02-27 00:17:46 +05:30
#ifdef DEBUG
use debug , only : debug_level , &
debug_e , &
debug_i , &
debug_g , &
debug_crystallite , &
debug_levelBasic , &
debug_levelExtensive , &
debug_levelSelective
#endif
2013-11-21 16:28:41 +05:30
use mesh , only : &
2019-01-29 09:41:29 +05:30
mesh_element
2013-11-21 16:28:41 +05:30
use material , only : &
2014-05-27 20:16:03 +05:30
plasticState , &
2015-05-28 22:32:23 +05:30
sourceState , &
2016-01-15 05:49:44 +05:30
phaseAt , phasememberAt , &
2015-05-28 22:32:23 +05:30
phase_Nsources , &
2019-01-15 08:57:57 +05:30
homogenization_Ngrains
2013-11-21 16:28:41 +05:30
use constitutive , only : &
2015-05-28 22:32:23 +05:30
constitutive_plasticity_maxSizeDotState , &
constitutive_source_maxSizeDotState
2014-08-26 20:14:32 +05:30
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-15 08:57:57 +05:30
NiterationState , & !< number of iterations in state loop
e , & !< element index in element loop
i , & !< integration point index in ip loop
g , & !< grain index in grain loop
2014-06-25 04:51:25 +05:30
p , &
c , &
2019-01-24 03:32:21 +05:30
s , &
2019-01-29 10:09:01 +05:30
sizeDotState
2019-01-15 08:57:57 +05:30
real ( pReal ) :: &
2019-01-30 02:59:36 +05:30
zeta
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState ) :: &
2019-01-30 02:49:38 +05:30
residuum_plastic ! residuum for plastic state
2019-01-30 02:59:36 +05:30
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState ) :: &
2019-01-30 02:49:38 +05:30
residuum_source ! residuum for source state
2013-11-21 16:28:41 +05:30
logical :: &
2019-01-15 08:57:57 +05:30
doneWithIntegration
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --+>> PREGUESS FOR STATE <<+--
2019-01-28 16:19:24 +05:30
call update_dotState ( 1.0_pReal )
call update_state ( 1.0_pReal )
2015-10-14 00:22:01 +05:30
2019-04-03 16:02:30 +05:30
NiterationState = 0
2019-01-15 08:57:57 +05:30
doneWithIntegration = . false .
2019-04-11 14:57:03 +05:30
crystalliteLooping : do while ( . not . doneWithIntegration . and . NiterationState < num % nState )
2019-04-03 16:02:30 +05:30
NiterationState = NiterationState + 1
2019-01-15 08:57:57 +05:30
2019-02-27 02:01:47 +05:30
#ifdef DEBUG
2019-04-03 16:02:30 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 ) &
2019-02-27 00:17:46 +05:30
write ( 6 , '(a,i6)' ) '<< CRYST stateFPI >> state iteration ' , NiterationState
2019-02-27 02:01:47 +05:30
#endif
2019-02-27 00:17:46 +05:30
2019-01-24 03:32:21 +05:30
! store previousDotState and previousDotState2
2019-01-29 12:59:19 +05:30
2019-01-24 03:32:21 +05:30
!$OMP PARALLEL DO PRIVATE(p,c)
2019-01-29 12:59:19 +05:30
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 ) )
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
p = phaseAt ( g , i , e ) ; c = phasememberAt ( g , i , e )
plasticState ( p ) % previousDotState2 ( : , c ) = merge ( plasticState ( p ) % previousDotState ( : , c ) , &
0.0_pReal , &
2019-04-03 16:02:30 +05:30
NiterationState > 1 )
2019-01-29 12:59:19 +05:30
plasticState ( p ) % previousDotState ( : , c ) = plasticState ( p ) % dotState ( : , c )
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-29 12:59:19 +05:30
sourceState ( p ) % p ( s ) % previousDotState2 ( : , c ) = merge ( sourceState ( p ) % p ( s ) % previousDotState ( : , c ) , &
0.0_pReal , &
2019-04-03 16:02:30 +05:30
NiterationState > 1 )
2019-01-29 12:59:19 +05:30
sourceState ( p ) % p ( s ) % previousDotState ( : , c ) = sourceState ( p ) % p ( s ) % dotState ( : , c )
enddo
endif
2019-01-15 08:57:57 +05:30
enddo
2019-01-24 03:32:21 +05:30
enddo
enddo
!$OMP END PARALLEL DO
call update_dependentState
2019-01-24 22:29:38 +05:30
call update_stress ( 1.0_pReal )
2019-01-23 10:44:19 +05:30
call update_dotState ( 1.0_pReal )
2019-01-29 12:59:19 +05:30
!$OMP PARALLEL
2019-01-30 02:59:36 +05:30
!$OMP DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c)
2019-02-27 00:17:46 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2019-01-28 16:19:24 +05:30
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
do g = 1 , homogenization_Ngrains ( mesh_element ( 3 , e ) )
2019-01-29 12:59:19 +05:30
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
p = phaseAt ( g , i , e ) ; c = phasememberAt ( g , i , e )
2019-01-30 04:01:26 +05:30
sizeDotState = plasticState ( p ) % sizeDotState
2014-05-27 20:16:03 +05:30
2019-01-30 02:59:36 +05:30
zeta = damper ( plasticState ( p ) % dotState ( : , c ) , &
plasticState ( p ) % previousDotState ( : , c ) , &
plasticState ( p ) % previousDotState2 ( : , c ) )
2019-01-30 04:01:26 +05:30
2019-01-30 02:49:38 +05:30
residuum_plastic ( 1 : SizeDotState ) = plasticState ( p ) % state ( 1 : sizeDotState , c ) &
2019-01-29 12:59:19 +05:30
- plasticState ( p ) % subState0 ( 1 : sizeDotState , c ) &
2019-01-30 02:59:36 +05:30
- ( plasticState ( p ) % dotState ( : , c ) * zeta &
+ plasticState ( p ) % previousDotState ( : , c ) * ( 1.0_pReal - zeta ) &
2019-01-29 12:59:19 +05:30
) * crystallite_subdt ( g , i , e )
2019-01-15 08:57:57 +05:30
2019-01-29 12:59:19 +05:30
plasticState ( p ) % state ( 1 : sizeDotState , c ) = plasticState ( p ) % state ( 1 : sizeDotState , c ) &
2019-01-30 02:49:38 +05:30
- residuum_plastic ( 1 : sizeDotState )
2019-01-30 02:59:36 +05:30
plasticState ( p ) % dotState ( : , c ) = plasticState ( p ) % dotState ( : , c ) * zeta &
+ plasticState ( p ) % previousDotState ( : , c ) * ( 1.0_pReal - zeta )
2019-01-31 13:42:44 +05:30
crystallite_converged ( g , i , e ) = converged ( residuum_plastic ( 1 : sizeDotState ) , &
plasticState ( p ) % state ( 1 : sizeDotState , c ) , &
plasticState ( p ) % aTolState ( 1 : sizeDotState ) )
2019-01-29 09:50:16 +05:30
2014-05-27 20:16:03 +05:30
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 04:01:26 +05:30
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
2019-01-30 02:59:36 +05:30
zeta = damper ( sourceState ( p ) % p ( s ) % dotState ( : , c ) , &
sourceState ( p ) % p ( s ) % previousDotState ( : , c ) , &
sourceState ( p ) % p ( s ) % previousDotState2 ( : , c ) )
2019-01-30 04:01:26 +05:30
2019-01-30 02:59:36 +05:30
residuum_source ( 1 : sizeDotState ) = sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) &
- sourceState ( p ) % p ( s ) % subState0 ( 1 : sizeDotState , c ) &
- ( sourceState ( p ) % p ( s ) % dotState ( : , c ) * zeta &
+ sourceState ( p ) % p ( s ) % previousDotState ( : , c ) * ( 1.0_pReal - zeta ) &
) * crystallite_subdt ( g , i , e )
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) = sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) &
- residuum_source ( 1 : sizeDotState )
sourceState ( p ) % p ( s ) % dotState ( : , c ) = sourceState ( p ) % p ( s ) % dotState ( : , c ) * zeta &
+ sourceState ( p ) % p ( s ) % previousDotState ( : , c ) * ( 1.0_pReal - zeta )
2019-02-27 00:26:49 +05:30
crystallite_converged ( g , i , e ) = &
crystallite_converged ( g , i , e ) . and . converged ( residuum_source ( 1 : sizeDotState ) , &
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) , &
sourceState ( p ) % p ( s ) % aTolState ( 1 : sizeDotState ) )
2019-01-30 02:59:36 +05:30
enddo
2019-01-15 08:57:57 +05:30
endif
2019-01-30 02:59:36 +05:30
enddo ; enddo ; enddo
2019-01-15 08:57:57 +05:30
!$OMP ENDDO
!$OMP DO
2019-01-28 16:19:24 +05:30
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 ) )
2019-01-15 08:57:57 +05:30
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) . and . crystallite_converged ( g , i , e ) ) then ! converged and still alive...
crystallite_todo ( g , i , e ) = stateJump ( g , i , e )
!$OMP FLUSH(crystallite_todo)
if ( . not . crystallite_todo ( g , i , e ) ) then ! if state jump fails, then convergence is broken
crystallite_converged ( g , i , e ) = . false .
if ( . not . crystallite_localPlasticity ( g , i , e ) ) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo . and . crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
2013-02-22 04:38:36 +05:30
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2019-01-15 08:57:57 +05:30
!$OMP END PARALLEL
2014-08-26 20:14:32 +05:30
2019-01-29 09:41:29 +05:30
if ( any ( plasticState ( : ) % nonlocal ) ) call nonlocalConvergenceCheck
2014-08-26 20:14:32 +05:30
2014-05-27 20:16:03 +05:30
2019-01-15 08:57:57 +05:30
! --- CHECK IF DONE WITH INTEGRATION ---
doneWithIntegration = . true .
2019-01-28 16:19:24 +05:30
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 ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
doneWithIntegration = . false .
2019-01-28 16:19:24 +05:30
exit
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo
2019-01-28 16:19:24 +05:30
enddo
2019-01-15 08:57:57 +05:30
enddo crystalliteLooping
2019-01-28 16:19:24 +05:30
2019-01-28 15:58:46 +05:30
contains
2019-01-29 04:57:58 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate the damping for correction of state and dot state
!--------------------------------------------------------------------------------------------------
2019-01-28 15:58:46 +05:30
real ( pReal ) pure function damper ( current , previous , previous2 )
real ( pReal ) , dimension ( : ) , intent ( in ) :: &
current , previous , previous2
real ( pReal ) :: dot_prod12 , dot_prod22
2019-01-29 15:22:00 +05:30
dot_prod12 = dot_product ( current - previous , previous - previous2 )
dot_prod22 = dot_product ( previous - previous2 , previous - previous2 )
2019-01-30 11:17:36 +05:30
if ( ( dot_product ( current , previous ) < 0.0_pReal . or . dot_prod12 < 0.0_pReal ) . and . dot_prod22 > 0.0_pReal ) then
2019-01-28 15:58:46 +05:30
damper = 0.75_pReal + 0.25_pReal * tanh ( 2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22 )
else
damper = 1.0_pReal
endif
end function damper
2019-01-15 08:57:57 +05:30
end subroutine integrateStateFPI
2010-10-01 17:48:49 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-29 09:41:29 +05:30
!> @brief integrate state with 1st order explicit Euler method
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine integrateStateEuler
2013-11-21 16:28:41 +05:30
use material , only : &
2019-01-28 16:19:24 +05:30
plasticState
2019-01-29 09:41:29 +05:30
2019-01-23 10:44:19 +05:30
call update_dotState ( 1.0_pReal )
2019-01-29 09:41:29 +05:30
call update_state ( 1.0_pReal )
2019-01-24 16:03:04 +05:30
call update_deltaState
call update_dependentState
2019-01-24 22:29:38 +05:30
call update_stress ( 1.0_pReal )
2019-01-25 11:50:05 +05:30
call setConvergenceFlag
2019-01-29 09:41:29 +05:30
if ( any ( plasticState ( : ) % nonlocal ) ) call nonlocalConvergenceCheck
2014-08-26 20:14:32 +05:30
2018-09-20 09:57:53 +05:30
end subroutine integrateStateEuler
2010-10-01 17:48:49 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine integrateStateAdaptiveEuler
2013-11-21 16:28:41 +05:30
use mesh , only : &
2019-02-02 16:45:05 +05:30
theMesh , &
mesh_element
2013-11-21 16:28:41 +05:30
use material , only : &
2019-01-15 08:57:57 +05:30
homogenization_Ngrains , &
2014-05-27 20:16:03 +05:30
plasticState , &
2015-05-28 22:32:23 +05:30
sourceState , &
2016-01-15 05:49:44 +05:30
phaseAt , phasememberAt , &
2015-05-28 22:32:23 +05:30
phase_Nsources , &
2019-01-15 08:57:57 +05:30
homogenization_maxNgrains
2013-11-21 16:28:41 +05:30
use constitutive , only : &
2015-05-28 22:32:23 +05:30
constitutive_plasticity_maxSizeDotState , &
constitutive_source_maxSizeDotState
2014-07-02 17:57:39 +05:30
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-15 08:57:57 +05:30
e , & ! element index in element loop
i , & ! integration point index in ip loop
g , & ! grain index in grain loop
2014-06-25 04:51:25 +05:30
p , &
c , &
2019-01-30 03:34:50 +05:30
s , &
2019-01-30 03:16:21 +05:30
sizeDotState
2019-01-30 04:01:26 +05:30
2019-01-30 14:58:47 +05:30
! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState , &
2019-02-02 16:45:05 +05:30
homogenization_maxNgrains , theMesh % elem % nIPs , theMesh % Nelems ) :: &
2019-01-30 14:58:47 +05:30
residuum_plastic
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , &
maxval ( phase_Nsources ) , &
2019-02-02 16:45:05 +05:30
homogenization_maxNgrains , theMesh % elem % nIPs , theMesh % Nelems ) :: &
2019-01-30 04:15:41 +05:30
residuum_source
2014-05-27 20:16:03 +05:30
2019-01-23 10:44:19 +05:30
!--------------------------------------------------------------------------------------------------
! contribution to state and relative residui and from Euler integration
call update_dotState ( 1.0_pReal )
2014-08-26 20:14:32 +05:30
2019-01-30 03:16:21 +05:30
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
2019-01-30 04:01:26 +05:30
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 ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-01-30 04:01:26 +05:30
p = phaseAt ( g , i , e ) ; c = phasememberAt ( g , i , e )
2019-01-30 03:16:21 +05:30
sizeDotState = plasticState ( p ) % sizeDotState
2019-01-30 04:01:26 +05:30
2019-01-30 04:15:41 +05:30
residuum_plastic ( 1 : sizeDotState , g , i , e ) = plasticState ( p ) % dotstate ( 1 : sizeDotState , c ) &
* ( - 0.5_pReal * crystallite_subdt ( g , i , e ) )
2019-02-27 02:01:47 +05:30
plasticState ( p ) % state ( 1 : sizeDotState , c ) = &
plasticState ( p ) % state ( 1 : sizeDotState , c ) + plasticState ( p ) % dotstate ( 1 : sizeDotState , c ) * crystallite_subdt ( g , i , e ) !ToDo: state, partitioned state?
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 03:34:50 +05:30
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
2019-01-30 04:01:26 +05:30
2019-01-30 04:15:41 +05:30
residuum_source ( 1 : sizeDotState , s , g , i , e ) = sourceState ( p ) % p ( s ) % dotstate ( 1 : sizeDotState , c ) &
* ( - 0.5_pReal * crystallite_subdt ( g , i , e ) )
2019-02-27 02:01:47 +05:30
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) = &
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) + sourceState ( p ) % p ( s ) % dotstate ( 1 : sizeDotState , c ) * crystallite_subdt ( g , i , e ) !ToDo: state, partitioned state?
2019-01-15 08:57:57 +05:30
enddo
endif
enddo ; enddo ; enddo
2019-01-29 05:24:02 +05:30
!$OMP END PARALLEL DO
2019-01-15 08:57:57 +05:30
2019-01-30 14:58:47 +05:30
call update_deltaState
call update_dependentState
call update_stress ( 1.0_pReal )
call update_dotState ( 1.0_pReal )
2019-01-15 08:57:57 +05:30
2019-01-30 04:15:41 +05:30
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
2019-01-30 04:01:26 +05:30
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 ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-01-29 05:24:02 +05:30
p = phaseAt ( g , i , e ) ; c = phasememberAt ( g , i , e )
2019-01-30 04:01:26 +05:30
sizeDotState = plasticState ( p ) % sizeDotState
2019-01-30 04:15:41 +05:30
residuum_plastic ( 1 : sizeDotState , g , i , e ) = residuum_plastic ( 1 : sizeDotState , g , i , e ) &
+ 0.5_pReal * plasticState ( p ) % dotState ( : , c ) * crystallite_subdt ( g , i , e )
2019-01-31 13:42:44 +05:30
crystallite_converged ( g , i , e ) = converged ( residuum_plastic ( 1 : sizeDotState , g , i , e ) , &
2019-02-27 02:01:47 +05:30
plasticState ( p ) % state ( 1 : sizeDotState , c ) , &
plasticState ( p ) % aTolState ( 1 : sizeDotState ) )
2015-05-28 22:32:23 +05:30
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 03:34:50 +05:30
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
2019-01-30 04:01:26 +05:30
2019-02-27 02:01:47 +05:30
residuum_source ( 1 : sizeDotState , s , g , i , e ) = &
residuum_source ( 1 : sizeDotState , s , g , i , e ) + 0.5_pReal * sourceState ( p ) % p ( s ) % dotState ( : , c ) * crystallite_subdt ( g , i , e )
2019-01-30 17:06:02 +05:30
2019-02-27 02:01:47 +05:30
crystallite_converged ( g , i , e ) = &
crystallite_converged ( g , i , e ) . and . converged ( residuum_source ( 1 : sizeDotState , s , g , i , e ) , &
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) , &
sourceState ( p ) % p ( s ) % aTolState ( 1 : sizeDotState ) )
2019-01-30 17:06:02 +05:30
enddo
2013-02-22 04:38:36 +05:30
endif
2019-01-30 04:01:26 +05:30
enddo ; enddo ; enddo
2019-01-29 05:24:02 +05:30
!$OMP END PARALLEL DO
2014-01-22 00:15:41 +05:30
2019-01-29 09:41:29 +05:30
if ( any ( plasticState ( : ) % nonlocal ) ) call nonlocalConvergenceCheck
2019-01-30 17:06:02 +05:30
2019-01-15 08:57:57 +05:30
end subroutine integrateStateAdaptiveEuler
2010-10-01 17:48:49 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
!> @brief integrate stress, state with 4th order explicit Runge Kutta method
2019-01-30 04:47:04 +05:30
! ToDo: This is totally BROKEN: RK4dotState is never used!!!
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine integrateStateRK4
2019-01-15 08:57:57 +05:30
use mesh , only : &
2019-01-30 04:41:10 +05:30
mesh_element
2013-11-21 16:28:41 +05:30
use material , only : &
2019-01-15 08:57:57 +05:30
homogenization_Ngrains , &
2014-05-27 20:16:03 +05:30
plasticState , &
2015-06-01 21:32:27 +05:30
sourceState , &
2019-01-15 08:57:57 +05:30
phase_Nsources , &
phaseAt , phasememberAt
2014-09-10 14:07:12 +05:30
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( 4 ) , parameter :: &
TIMESTEPFRACTION = [ 0.5_pReal , 0.5_pReal , 1.0_pReal , 1.0_pReal ] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
real ( pReal ) , dimension ( 4 ) , parameter :: &
WEIGHT = [ 1.0_pReal , 2.0_pReal , 2.0_pReal , 1.0_pReal / 6.0_pReal ] ! weight of slope used for Runge Kutta integration (final weight divided by 6)
2014-09-10 14:07:12 +05:30
2019-04-03 16:02:30 +05:30
integer :: e , & ! element index in element loop
2019-01-15 08:57:57 +05:30
i , & ! integration point index in ip loop
g , & ! grain index in grain loop
p , & ! phase loop
c , &
n , &
2019-01-30 04:31:40 +05:30
s
2019-01-15 08:57:57 +05:30
2019-01-23 10:44:19 +05:30
call update_dotState ( 1.0_pReal )
2018-05-09 20:24:06 +05:30
2014-08-26 20:14:32 +05:30
2019-04-03 16:02:30 +05:30
do n = 1 , 4
2013-04-29 16:47:30 +05:30
2019-01-30 04:41:10 +05:30
!$OMP PARALLEL DO PRIVATE(p,c)
2019-01-29 09:41:29 +05:30
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 ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-01-30 04:31:40 +05:30
p = phaseAt ( g , i , e ) ; c = phasememberAt ( g , i , e )
2019-01-30 04:41:10 +05:30
plasticState ( p ) % RK4dotState ( : , c ) = WEIGHT ( n ) * plasticState ( p ) % dotState ( : , c ) &
2019-04-03 16:02:30 +05:30
+ merge ( plasticState ( p ) % RK4dotState ( : , c ) , 0.0_pReal , n > 1 )
do s = 1 , phase_Nsources ( p )
2019-01-30 04:41:10 +05:30
sourceState ( p ) % p ( s ) % RK4dotState ( : , c ) = WEIGHT ( n ) * sourceState ( p ) % p ( s ) % dotState ( : , c ) &
2019-04-03 16:02:30 +05:30
+ merge ( sourceState ( p ) % p ( s ) % RK4dotState ( : , c ) , 0.0_pReal , n > 1 )
2019-01-15 08:57:57 +05:30
enddo
endif
enddo ; enddo ; enddo
2019-01-30 04:41:10 +05:30
!$OMP END PARALLEL DO
2019-01-15 08:57:57 +05:30
2019-01-23 16:21:43 +05:30
call update_state ( TIMESTEPFRACTION ( n ) )
2019-01-24 16:03:04 +05:30
call update_deltaState
2019-01-23 22:34:19 +05:30
call update_dependentState
2019-01-24 22:29:38 +05:30
call update_stress ( TIMESTEPFRACTION ( n ) )
2019-01-15 08:57:57 +05:30
! --- dot state and RK dot state---
first3steps : if ( n < 4 ) then
2019-01-30 04:41:10 +05:30
call update_dotState ( TIMESTEPFRACTION ( n ) )
2019-01-15 08:57:57 +05:30
endif first3steps
enddo
2019-01-29 09:41:29 +05:30
call setConvergenceFlag
if ( any ( plasticState ( : ) % nonlocal ) ) call nonlocalConvergenceCheck
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
end subroutine integrateStateRK4
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with
!> adaptive step size (use 5th order solution to advance = "local extrapolation")
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine integrateStateRKCK45
2019-01-15 08:57:57 +05:30
use mesh , only : &
mesh_element , &
2019-02-02 16:45:05 +05:30
theMesh
2019-01-15 08:57:57 +05:30
use material , only : &
homogenization_Ngrains , &
plasticState , &
sourceState , &
phase_Nsources , &
phaseAt , phasememberAt , &
homogenization_maxNgrains
use constitutive , only : &
constitutive_plasticity_maxSizeDotState , &
2019-01-30 15:16:53 +05:30
constitutive_source_maxSizeDotState
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( 5 , 5 ) , parameter :: &
A = reshape ( [ &
2019-02-27 02:01:47 +05:30
. 2_pReal , . 075_pReal , . 3_pReal , - 1 1.0_pReal / 5 4.0_pReal , 163 1.0_pReal / 5529 6.0_pReal , &
. 0_pReal , . 225_pReal , - . 9_pReal , 2.5_pReal , 17 5.0_pReal / 51 2.0_pReal , &
. 0_pReal , . 0_pReal , 1.2_pReal , - 7 0.0_pReal / 2 7.0_pReal , 57 5.0_pReal / 1382 4.0_pReal , &
. 0_pReal , . 0_pReal , . 0_pReal , 3 5.0_pReal / 2 7.0_pReal , 4427 5.0_pReal / 11059 2.0_pReal , &
. 0_pReal , . 0_pReal , . 0_pReal , . 0_pReal , 25 3.0_pReal / 409 6.0_pReal ] , &
2019-01-15 08:57:57 +05:30
[ 5 , 5 ] , order = [ 2 , 1 ] ) !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( 6 ) , parameter :: &
B = &
[ 3 7.0_pReal / 37 8.0_pReal , . 0_pReal , 25 0.0_pReal / 62 1.0_pReal , &
12 5.0_pReal / 59 4.0_pReal , . 0_pReal , 51 2.0_pReal / 177 1.0_pReal ] , & !< coefficients in Butcher tableau (used for final integration and error estimate)
DB = B - &
2019-02-27 02:01:47 +05:30
[ 282 5.0_pReal / 2764 8.0_pReal , . 0_pReal , 1857 5.0_pReal / 4838 4.0_pReal , &
1352 5.0_pReal / 5529 6.0_pReal , 27 7.0_pReal / 1433 6.0_pReal , 0.25_pReal ] !< coefficients in Butcher tableau (used for final integration and error estimate)
2019-01-15 08:57:57 +05:30
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)
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-15 08:57:57 +05:30
e , & ! element index in element loop
i , & ! integration point index in ip loop
g , & ! grain index in grain loop
stage , & ! stage index in integration stage loop
n , &
p , &
cc , &
2019-01-30 15:07:18 +05:30
s , &
2019-01-30 15:18:59 +05:30
sizeDotState
2019-01-29 09:41:29 +05:30
2019-01-30 15:34:49 +05:30
! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of RKCK45
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState , &
2019-02-02 16:45:05 +05:30
homogenization_maxNgrains , theMesh % elem % nIPs , theMesh % Nelems ) :: &
2019-01-30 15:34:49 +05:30
residuum_plastic ! relative residuum from evolution in microstructure
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , &
maxval ( phase_Nsources ) , &
2019-02-02 16:45:05 +05:30
homogenization_maxNgrains , theMesh % elem % nIPs , theMesh % Nelems ) :: &
2019-01-30 15:34:49 +05:30
residuum_source ! relative residuum from evolution in microstructure
2014-11-01 00:33:08 +05:30
2014-11-12 22:10:50 +05:30
2019-01-23 10:44:19 +05:30
call update_dotState ( 1.0_pReal )
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
! --- SECOND TO SIXTH RUNGE KUTTA STEP ---
2014-11-01 00:33:08 +05:30
2019-04-03 16:02:30 +05:30
do stage = 1 , 5
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
! --- state update ---
2019-01-30 13:26:16 +05:30
!$OMP PARALLEL DO PRIVATE(p,cc)
2019-01-29 09:41:29 +05:30
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 ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-01-30 13:41:12 +05:30
p = phaseAt ( g , i , e ) ; cc = phasememberAt ( g , i , e )
2014-11-12 22:10:50 +05:30
2019-01-30 13:26:16 +05:30
plasticState ( p ) % RKCK45dotState ( stage , : , cc ) = plasticState ( p ) % dotState ( : , cc )
2019-01-15 08:57:57 +05:30
plasticState ( p ) % dotState ( : , cc ) = A ( 1 , stage ) * plasticState ( p ) % RKCK45dotState ( 1 , : , cc )
2019-01-30 13:26:16 +05:30
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 15:07:18 +05:30
sourceState ( p ) % p ( s ) % RKCK45dotState ( stage , : , cc ) = sourceState ( p ) % p ( s ) % dotState ( : , cc )
sourceState ( p ) % p ( s ) % dotState ( : , cc ) = A ( 1 , stage ) * sourceState ( p ) % p ( s ) % RKCK45dotState ( 1 , : , cc )
2019-01-15 08:57:57 +05:30
enddo
2019-01-30 13:26:16 +05:30
2019-04-03 16:02:30 +05:30
do n = 2 , stage
2019-01-30 13:26:16 +05:30
plasticState ( p ) % dotState ( : , cc ) = plasticState ( p ) % dotState ( : , cc ) &
+ A ( n , stage ) * plasticState ( p ) % RKCK45dotState ( n , : , cc )
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 15:07:18 +05:30
sourceState ( p ) % p ( s ) % dotState ( : , cc ) = sourceState ( p ) % p ( s ) % dotState ( : , cc ) &
2019-01-30 15:16:53 +05:30
+ A ( n , stage ) * sourceState ( p ) % p ( s ) % RKCK45dotState ( n , : , cc )
2019-01-15 08:57:57 +05:30
enddo
enddo
2019-01-30 13:26:16 +05:30
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
2019-01-30 13:26:16 +05:30
!$OMP END PARALLEL DO
2014-11-01 00:33:08 +05:30
2019-01-23 16:21:43 +05:30
call update_state ( 1.0_pReal ) !MD: 1.0 correct?
2019-01-24 16:03:04 +05:30
call update_deltaState
call update_dependentState
2019-01-24 22:29:38 +05:30
call update_stress ( C ( stage ) )
call update_dotState ( C ( stage ) )
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
enddo
2014-08-26 20:14:32 +05:30
2015-03-06 18:39:00 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE ---
2019-01-30 15:18:59 +05:30
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc)
2019-01-29 09:41:29 +05:30
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 ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-01-30 13:51:33 +05:30
p = phaseAt ( g , i , e ) ; cc = phasememberAt ( g , i , e )
2019-01-30 15:18:59 +05:30
sizeDotState = plasticState ( p ) % sizeDotState
2019-01-30 13:51:33 +05:30
plasticState ( p ) % RKCK45dotState ( 6 , : , cc ) = plasticState ( p ) % dotState ( : , cc )
2019-01-30 15:21:24 +05:30
residuum_plastic ( 1 : sizeDotState , g , i , e ) = &
2019-01-31 09:47:46 +05:30
matmul ( transpose ( plasticState ( p ) % RKCK45dotState ( 1 : 6 , 1 : sizeDotState , cc ) ) , DB ) & ! why transpose? Better to transpose constant DB
2019-01-15 08:57:57 +05:30
* crystallite_subdt ( g , i , e )
2019-01-30 13:51:33 +05:30
plasticState ( p ) % dotState ( : , cc ) = &
2019-01-31 09:47:46 +05:30
matmul ( transpose ( plasticState ( p ) % RKCK45dotState ( 1 : 6 , 1 : sizeDotState , cc ) ) , B ) ! why transpose? Better to transpose constant B
2019-01-30 13:51:33 +05:30
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 15:18:59 +05:30
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
2019-01-30 13:51:33 +05:30
2019-01-30 19:22:12 +05:30
sourceState ( p ) % p ( s ) % RKCK45dotState ( 6 , : , cc ) = sourceState ( p ) % p ( s ) % dotState ( : , cc )
2019-01-30 15:21:24 +05:30
residuum_source ( 1 : sizeDotState , s , g , i , e ) = &
2019-01-30 15:18:59 +05:30
matmul ( transpose ( sourceState ( p ) % p ( s ) % RKCK45dotState ( 1 : 6 , 1 : sizeDotState , cc ) ) , DB ) &
2019-01-15 08:57:57 +05:30
* crystallite_subdt ( g , i , e )
2014-08-26 20:14:32 +05:30
2019-01-30 15:07:18 +05:30
sourceState ( p ) % p ( s ) % dotState ( : , cc ) = &
2019-01-30 15:18:59 +05:30
matmul ( transpose ( sourceState ( p ) % p ( s ) % RKCK45dotState ( 1 : 6 , 1 : sizeDotState , cc ) ) , B )
2019-01-15 08:57:57 +05:30
enddo
2019-01-30 13:51:33 +05:30
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
2019-01-30 13:51:33 +05:30
!$OMP END PARALLEL DO
2014-08-26 20:14:32 +05:30
2019-01-23 16:21:43 +05:30
call update_state ( 1.0_pReal )
2019-01-15 08:57:57 +05:30
! --- relative residui and state convergence ---
2014-08-26 20:14:32 +05:30
2019-01-30 15:34:49 +05:30
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc)
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 ) )
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e ) ; cc = phasememberAt ( g , i , e )
2019-01-30 15:07:18 +05:30
2019-01-30 15:34:49 +05:30
sizeDotState = plasticState ( p ) % sizeDotState
2019-01-30 17:06:02 +05:30
2019-01-31 13:42:44 +05:30
crystallite_todo ( g , i , e ) = converged ( residuum_plastic ( 1 : sizeDotState , g , i , e ) , &
2019-01-30 21:34:58 +05:30
plasticState ( p ) % state ( 1 : sizeDotState , cc ) , &
2019-01-31 13:42:44 +05:30
plasticState ( p ) % aTolState ( 1 : sizeDotState ) )
2014-08-26 20:14:32 +05:30
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 15:34:49 +05:30
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
2019-02-27 02:01:47 +05:30
crystallite_todo ( g , i , e ) = &
crystallite_todo ( g , i , e ) . and . converged ( residuum_source ( 1 : sizeDotState , s , g , i , e ) , &
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , cc ) , &
sourceState ( p ) % p ( s ) % aTolState ( 1 : sizeDotState ) )
2019-01-30 19:22:12 +05:30
enddo
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
2019-01-30 15:34:49 +05:30
!$OMP END PARALLEL DO
2014-08-26 20:14:32 +05:30
2019-01-24 16:03:04 +05:30
call update_deltaState
2019-01-23 22:34:19 +05:30
call update_dependentState
2019-01-24 22:29:38 +05:30
call update_stress ( 1.0_pReal )
2019-01-25 11:50:05 +05:30
call setConvergenceFlag
2019-01-29 09:41:29 +05:30
if ( any ( plasticState ( : ) % nonlocal ) ) call nonlocalConvergenceCheck
2019-01-30 17:06:02 +05:30
2019-01-29 09:41:29 +05:30
end subroutine integrateStateRKCK45
2011-05-11 22:08:45 +05:30
2014-08-26 20:14:32 +05:30
2019-01-29 09:41:29 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief sets convergence flag for nonlocal calculations
!> @detail one non-converged nonlocal sets all other nonlocals to non-converged to trigger cut back
!--------------------------------------------------------------------------------------------------
2019-04-11 11:05:58 +05:30
subroutine nonlocalConvergenceCheck
2019-01-29 09:41:29 +05:30
if ( any ( . not . crystallite_converged . and . . not . crystallite_localPlasticity ) ) & ! any non-local not yet converged (or broken)...
where ( . not . crystallite_localPlasticity ) crystallite_converged = . false .
2014-06-25 04:51:25 +05:30
2019-01-29 09:41:29 +05:30
end subroutine nonlocalConvergenceCheck
2009-05-07 21:57:36 +05:30
2019-01-19 14:05:45 +05:30
2019-01-25 11:50:05 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Sets convergence flag based on "todo": every point that survived the integration (todo is
! still .true. is considered as converged
!> @details: For explicitEuler, RK4 and RKCK45, adaptive Euler and FPI have their on criteria
!--------------------------------------------------------------------------------------------------
2019-04-11 11:05:58 +05:30
subroutine setConvergenceFlag
2019-02-02 15:05:10 +05:30
use mesh , only : &
mesh_element
2019-04-11 11:05:58 +05:30
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-25 11:50:05 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
g !< grain index in grain loop
2019-04-03 16:47:21 +05:30
!OMP DO PARALLEL PRIVATE
2019-01-25 11:50:05 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2019-04-03 16:47:21 +05:30
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
do g = 1 , homogenization_Ngrains ( mesh_element ( 3 , e ) )
crystallite_converged ( g , i , e ) = crystallite_todo ( g , i , e ) . or . crystallite_converged ( g , i , e ) ! if still "to do" then converged per definition
enddo ; enddo ; enddo
2019-01-25 11:50:05 +05:30
!OMP END DO PARALLEL
end subroutine setConvergenceFlag
2019-01-31 13:44:02 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief determines whether a point is converged
!--------------------------------------------------------------------------------------------------
logical pure function converged ( residuum , state , aTol )
real ( pReal ) , intent ( in ) , dimension ( : ) :: &
residuum , state , aTol
2019-04-11 10:54:04 +05:30
real ( pReal ) :: &
rTol
rTol = num % rTol_crystalliteState
2019-01-31 13:44:02 +05:30
converged = all ( abs ( residuum ) < = max ( aTol , rTol * abs ( state ) ) )
end function converged
2019-01-24 22:29:38 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
!--------------------------------------------------------------------------------------------------
subroutine update_stress ( timeFraction )
2019-02-02 15:05:10 +05:30
use mesh , only : &
mesh_element
2019-04-11 11:05:58 +05:30
2019-01-24 22:29:38 +05:30
real ( pReal ) , intent ( in ) :: &
timeFraction
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-24 22:29:38 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
g
!$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 ) )
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
crystallite_todo ( g , i , e ) = integrateStress ( g , i , e , timeFraction )
!$OMP FLUSH(crystallite_todo)
if ( . not . crystallite_todo ( g , i , e ) . and . . not . crystallite_localPlasticity ( g , i , e ) ) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo . and . crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
end subroutine update_stress
2019-01-23 22:34:19 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief tbd
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine update_dependentState
2019-02-02 15:05:10 +05:30
use mesh , only : &
mesh_element
2019-01-23 22:34:19 +05:30
use constitutive , only : &
constitutive_dependentState = > constitutive_microstructure
2019-04-03 16:02:30 +05:30
integer :: e , & ! element index in element loop
2019-01-23 22:34:19 +05:30
i , & ! integration point index in ip loop
g ! grain index in grain loop
2019-01-24 03:48:14 +05:30
!$OMP PARALLEL DO
2019-01-23 22:34:19 +05:30
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 ) )
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) &
2019-02-02 00:58:07 +05:30
call constitutive_dependentState ( crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
2019-01-23 22:34:19 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , g , i , e ) , &
g , i , e )
enddo ; enddo ; enddo
2019-01-24 03:48:14 +05:30
!$OMP END PARALLEL DO
2019-01-23 22:34:19 +05:30
end subroutine update_dependentState
2019-01-24 22:29:38 +05:30
2019-01-23 16:21:43 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
!--------------------------------------------------------------------------------------------------
subroutine update_state ( timeFraction )
use material , only : &
plasticState , &
sourceState , &
phase_Nsources , &
phaseAt , phasememberAt
2019-02-02 15:05:10 +05:30
use mesh , only : &
mesh_element
2019-01-23 16:21:43 +05:30
real ( pReal ) , intent ( in ) :: &
timeFraction
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-23 16:21:43 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
g , & !< grain index in grain loop
p , &
c , &
s , &
mySize
!$OMP PARALLEL DO PRIVATE(mySize,p,c)
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 ) )
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
p = phaseAt ( g , i , e ) ; c = phasememberAt ( g , i , e )
mySize = plasticState ( p ) % sizeDotState
plasticState ( p ) % state ( 1 : mySize , c ) = plasticState ( p ) % subState0 ( 1 : mySize , c ) &
+ plasticState ( p ) % dotState ( 1 : mySize , c ) &
* crystallite_subdt ( g , i , e ) * timeFraction
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-23 16:21:43 +05:30
mySize = sourceState ( p ) % p ( s ) % sizeDotState
sourceState ( p ) % p ( s ) % state ( 1 : mySize , c ) = sourceState ( p ) % p ( s ) % subState0 ( 1 : mySize , c ) &
+ sourceState ( p ) % p ( s ) % dotState ( 1 : mySize , c ) &
* crystallite_subdt ( g , i , e ) * timeFraction
enddo
endif
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
end subroutine update_state
2019-01-23 10:44:19 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief triggers calculation of all new rates
!> if NaN occurs, crystallite_todo is set to FALSE. Any NaN in a nonlocal propagates to all others
!--------------------------------------------------------------------------------------------------
subroutine update_dotState ( timeFraction )
use , intrinsic :: &
IEEE_arithmetic
use material , only : &
plasticState , &
sourceState , &
phaseAt , phasememberAt , &
phase_Nsources
2019-02-02 15:05:10 +05:30
use mesh , only : &
mesh_element
2019-01-23 10:44:19 +05:30
use constitutive , only : &
constitutive_collectDotState
real ( pReal ) , intent ( in ) :: &
timeFraction
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-23 10:44:19 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
g , & !< grain index in grain loop
p , &
c , &
s
logical :: &
2019-01-24 18:45:26 +05:30
NaN , &
nonlocalStop
nonlocalStop = . false .
2019-01-23 10:44:19 +05:30
2019-01-24 18:45:26 +05:30
!$OMP PARALLEL DO PRIVATE (p,c,NaN)
2019-01-23 10:44:19 +05:30
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 ) )
2019-01-24 18:45:26 +05:30
!$OMP FLUSH(nonlocalStop)
2019-01-30 20:36:14 +05:30
if ( ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) . and . . not . nonlocalStop ) then
2019-03-09 05:03:11 +05:30
call constitutive_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
2019-01-23 10:44:19 +05:30
crystallite_Fe , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp , &
2019-02-26 12:24:03 +05:30
crystallite_subdt ( g , i , e ) * timeFraction , g , i , e )
2019-01-23 10:44:19 +05:30
p = phaseAt ( g , i , e ) ; c = phasememberAt ( g , i , e )
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % dotState ( : , c ) ) )
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-23 10:44:19 +05:30
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( s ) % dotState ( : , c ) ) )
enddo
if ( NaN ) then
crystallite_todo ( g , i , e ) = . false . ! this one done (and broken)
2019-01-24 18:45:26 +05:30
if ( . not . crystallite_localPlasticity ( g , i , e ) ) nonlocalStop = . True .
2019-01-23 10:44:19 +05:30
endif
endif
enddo ; enddo ; enddo
2019-01-24 18:45:26 +05:30
!$OMP END PARALLEL DO
if ( nonlocalStop ) crystallite_todo = crystallite_todo . and . crystallite_localPlasticity
2019-01-23 10:44:19 +05:30
end subroutine update_DotState
2019-01-24 11:42:20 +05:30
subroutine update_deltaState
2019-01-24 12:04:30 +05:30
use , intrinsic :: &
IEEE_arithmetic
use prec , only : &
dNeq0
2019-02-02 15:05:10 +05:30
use mesh , only : &
mesh_element
2019-01-24 12:04:30 +05:30
use material , only : &
plasticState , &
sourceState , &
phase_Nsources , &
phaseAt , phasememberAt
use constitutive , only : &
constitutive_collectDeltaState
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-24 11:42:20 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
g , & !< grain index in grain loop
p , &
2019-01-30 04:41:10 +05:30
mySize , &
2019-01-24 12:04:30 +05:30
myOffset , &
2019-01-24 11:42:20 +05:30
c , &
s
2019-01-24 18:45:26 +05:30
logical :: &
NaN , &
nonlocalStop
nonlocalStop = . false .
2019-01-24 11:42:20 +05:30
2019-01-30 04:41:10 +05:30
!$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize,NaN)
2019-01-24 11:42:20 +05:30
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 ) )
2019-01-24 18:45:26 +05:30
!$OMP FLUSH(nonlocalStop)
2019-01-30 20:36:14 +05:30
if ( ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) . and . . not . nonlocalStop ) then
2019-03-09 05:03:11 +05:30
call constitutive_collectDeltaState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
g , i , e )
2019-01-24 12:04:30 +05:30
p = phaseAt ( g , i , e ) ; c = phasememberAt ( g , i , e )
myOffset = plasticState ( p ) % offsetDeltaState
mySize = plasticState ( p ) % sizeDeltaState
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % deltaState ( 1 : mySize , c ) ) )
if ( . not . NaN ) then
2019-04-03 16:02:30 +05:30
plasticState ( p ) % state ( myOffset + 1 : myOffset + mySize , c ) = &
plasticState ( p ) % state ( myOffset + 1 : myOffset + mySize , c ) + plasticState ( p ) % deltaState ( 1 : mySize , c )
do s = 1 , phase_Nsources ( p )
2019-01-30 04:41:10 +05:30
myOffset = sourceState ( p ) % p ( s ) % offsetDeltaState
mySize = sourceState ( p ) % p ( s ) % sizeDeltaState
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( s ) % deltaState ( 1 : mySize , c ) ) )
2019-01-24 16:03:04 +05:30
if ( . not . NaN ) then
2019-04-03 16:02:30 +05:30
sourceState ( p ) % p ( s ) % state ( myOffset + 1 : myOffset + mySize , c ) = &
sourceState ( p ) % p ( s ) % state ( myOffset + 1 : myOffset + mySize , c ) + sourceState ( p ) % p ( s ) % deltaState ( 1 : mySize , c )
2019-01-24 16:03:04 +05:30
endif
enddo
2019-01-24 12:04:30 +05:30
endif
2019-01-24 16:03:04 +05:30
crystallite_todo ( g , i , e ) = . not . NaN
2019-01-24 11:42:20 +05:30
if ( . not . crystallite_todo ( g , i , e ) ) then ! if state jump fails, then convergence is broken
crystallite_converged ( g , i , e ) = . false .
2019-01-24 18:45:26 +05:30
if ( . not . crystallite_localPlasticity ( g , i , e ) ) nonlocalStop = . true .
2019-01-24 11:42:20 +05:30
endif
endif
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
2019-01-24 18:45:26 +05:30
if ( nonlocalStop ) crystallite_todo = crystallite_todo . and . crystallite_localPlasticity
2019-01-24 11:42:20 +05:30
end subroutine update_deltaState
2019-01-19 14:05:45 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates a jump in the state according to the current state and the current stress
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
!--------------------------------------------------------------------------------------------------
logical function stateJump ( ipc , ip , el )
use , intrinsic :: &
IEEE_arithmetic
use prec , only : &
dNeq0
#ifdef DEBUG
use debug , only : &
debug_e , &
debug_i , &
debug_g , &
debug_level , &
debug_crystallite , &
debug_levelExtensive , &
debug_levelSelective
#endif
use material , only : &
plasticState , &
sourceState , &
phase_Nsources , &
phaseAt , phasememberAt
use constitutive , only : &
constitutive_collectDeltaState
2019-04-03 16:02:30 +05:30
integer , intent ( in ) :: &
2019-01-19 14:05:45 +05:30
el , & ! element index
ip , & ! integration point index
ipc ! grain index
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-19 14:05:45 +05:30
c , &
p , &
mySource , &
2019-01-24 11:42:20 +05:30
myOffset , &
mySize
2019-01-19 14:05:45 +05:30
c = phasememberAt ( ipc , ip , el )
p = phaseAt ( ipc , ip , el )
2019-03-09 05:03:11 +05:30
call constitutive_collectDeltaState ( crystallite_S ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
crystallite_Fe ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
ipc , ip , el )
2019-01-19 14:05:45 +05:30
2019-01-24 11:42:20 +05:30
myOffset = plasticState ( p ) % offsetDeltaState
mySize = plasticState ( p ) % sizeDeltaState
2019-01-19 14:05:45 +05:30
2019-01-24 11:42:20 +05:30
if ( any ( IEEE_is_NaN ( plasticState ( p ) % deltaState ( 1 : mySize , c ) ) ) ) then ! NaN occured in deltaState
2019-01-19 14:05:45 +05:30
stateJump = . false .
return
endif
2019-04-03 16:02:30 +05:30
plasticState ( p ) % state ( myOffset + 1 : myOffset + mySize , c ) = &
plasticState ( p ) % state ( myOffset + 1 : myOffset + mySize , c ) + plasticState ( p ) % deltaState ( 1 : mySize , c )
2019-01-19 14:05:45 +05:30
2019-04-03 16:02:30 +05:30
do mySource = 1 , phase_Nsources ( p )
2019-01-24 11:42:20 +05:30
myOffset = sourceState ( p ) % p ( mySource ) % offsetDeltaState
2019-01-25 11:50:05 +05:30
mySize = sourceState ( p ) % p ( mySource ) % sizeDeltaState
2019-01-24 11:42:20 +05:30
if ( any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % deltaState ( 1 : mySize , c ) ) ) ) then ! NaN occured in deltaState
2019-01-19 14:05:45 +05:30
stateJump = . false .
return
endif
2019-04-03 16:02:30 +05:30
sourceState ( p ) % p ( mySource ) % state ( myOffset + 1 : myOffset + mySize , c ) = &
sourceState ( p ) % p ( mySource ) % state ( myOffset + 1 : myOffset + mySize , c ) + sourceState ( p ) % p ( mySource ) % deltaState ( 1 : mySize , c )
2019-01-19 14:05:45 +05:30
enddo
#ifdef DEBUG
2019-01-24 11:42:20 +05:30
if ( any ( dNeq0 ( plasticState ( p ) % deltaState ( 1 : mySize , c ) ) ) &
2019-04-03 16:02:30 +05:30
. and . iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
2019-01-19 14:05:45 +05:30
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
2019-04-03 16:02:30 +05:30
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
2019-01-19 14:05:45 +05:30
write ( 6 , '(a,i8,1x,i2,1x,i3, /)' ) '<< CRYST >> update state at el ip ipc ' , el , ip , ipc
2019-01-24 11:42:20 +05:30
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> deltaState' , plasticState ( p ) % deltaState ( 1 : mySize , c )
2019-01-19 14:05:45 +05:30
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> new state' , &
2019-04-03 16:02:30 +05:30
plasticState ( p ) % state ( myOffset + 1 : &
2019-01-24 11:42:20 +05:30
myOffset + mySize , c )
2019-01-19 14:05:45 +05:30
endif
#endif
stateJump = . true .
end function stateJump
2013-02-22 04:38:36 +05:30
end module crystallite
