2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
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
2013-04-29 16:47:30 +05:30
use prec , only : &
pReal , &
pInt
2012-08-31 01:56:28 +05:30
implicit none
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
private
2013-12-12 22:39:59 +05:30
character ( len = 64 ) , dimension ( : , : ) , allocatable , private :: &
2013-02-22 04:38:36 +05:30
crystallite_output !< name of each post result output
2013-12-12 22:39:59 +05:30
integer ( pInt ) , public , protected :: &
2013-05-17 23:22:46 +05:30
crystallite_maxSizePostResults !< description not available
2013-12-12 22:39:59 +05:30
integer ( pInt ) , dimension ( : ) , allocatable , public , protected :: &
2013-05-17 23:22:46 +05:30
crystallite_sizePostResults !< description not available
2013-12-12 22:39:59 +05:30
integer ( pInt ) , dimension ( : , : ) , allocatable , private :: &
2013-05-17 23:22:46 +05:30
crystallite_sizePostResult !< description not available
2014-08-26 20:14:32 +05:30
2013-12-12 22:39:59 +05:30
real ( pReal ) , dimension ( : , : , : ) , allocatable , public :: &
2013-10-16 18:34:59 +05:30
crystallite_dt !< requested time increment of each grain
2013-12-12 22:39:59 +05:30
real ( pReal ) , dimension ( : , : , : ) , allocatable , private :: &
2013-02-22 04:38:36 +05:30
crystallite_subdt , & !< substepped time increment of each grain
crystallite_subFrac , & !< already calculated fraction of increment
2013-10-16 18:34:59 +05:30
crystallite_subStep !< size of next integration step
2013-12-12 22:39:59 +05:30
real ( pReal ) , dimension ( : , : , : , : ) , allocatable , public :: &
2013-02-22 04:38:36 +05:30
crystallite_Tstar_v , & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
crystallite_Tstar0_v , & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_partionedTstar0_v !< 2nd Piola-Kirchhoff stress vector at start of homog inc
2013-12-12 22:39:59 +05:30
real ( pReal ) , dimension ( : , : , : , : ) , allocatable , private :: &
2013-02-22 04:38:36 +05:30
crystallite_subTstar0_v , & !< 2nd Piola-Kirchhoff stress vector at start of crystallite inc
crystallite_orientation , & !< orientation as quaternion
crystallite_orientation0 , & !< initial orientation as quaternion
2014-08-26 20:14:32 +05:30
crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees) in crystal reference frame
2018-09-19 20:34:12 +05:30
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
2013-12-12 22:39:59 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , public :: &
2013-02-22 04:38:36 +05:30
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
2015-03-06 18:39:00 +05:30
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
2013-02-22 04:38:36 +05:30
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
2018-09-19 20:34:12 +05:30
crystallite_partionedLp0 , & !< plastic velocity grad at start of homog inc
2015-03-06 18:39:00 +05:30
crystallite_Li , & !< current intermediate velocitiy grad (end of converged time step)
crystallite_Li0 , & !< intermediate velocitiy grad at start of FE inc
2018-09-19 20:34:12 +05:30
crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc
2013-12-12 22:39:59 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , private :: &
2013-02-22 04:38:36 +05:30
crystallite_invFp , & !< inverse of current plastic def grad (end of converged time step)
crystallite_subFp0 , & !< plastic def grad at start of crystallite inc
2015-03-06 18:39:00 +05:30
crystallite_invFi , & !< inverse of current intermediate def grad (end of converged time step)
crystallite_subFi0 , & !< intermediate def grad at start of crystallite inc
2018-02-16 20:06:18 +05:30
crystallite_subF , & !< def grad to be reached at end of crystallite inc
2013-02-22 04:38:36 +05:30
crystallite_subF0 , & !< def grad at start of crystallite inc
crystallite_subLp0 , & !< plastic velocity grad at start of crystallite inc
2019-01-14 12:14:36 +05:30
crystallite_subLi0 !< intermediate velocity grad at start of crystallite inc
2013-12-12 22:39:59 +05:30
real ( pReal ) , dimension ( : , : , : , : , : , : , : ) , allocatable , public :: &
2019-01-15 12:34:50 +05:30
crystallite_dPdF !< current individual dPdF per grain (end of converged time step)
2013-12-12 22:39:59 +05:30
logical , dimension ( : , : , : ) , allocatable , public :: &
2013-02-22 04:38:36 +05:30
crystallite_requested !< flag to request crystallite calculation
2013-12-12 22:39:59 +05:30
logical , dimension ( : , : , : ) , allocatable , public , protected :: &
2019-01-15 12:22:01 +05:30
crystallite_converged !< convergence flag
2013-12-12 22:39:59 +05:30
logical , dimension ( : , : , : ) , allocatable , private :: &
2019-01-15 12:22:01 +05:30
crystallite_localPlasticity , & !< indicates this grain to have purely local constitutive law
2013-02-22 04:38:36 +05:30
crystallite_todo !< flag to indicate need for further computation
2014-08-26 20:14:32 +05:30
enum , bind ( c )
2013-12-12 22:39:59 +05:30
enumerator :: undefined_ID , &
phase_ID , &
texture_ID , &
volume_ID , &
orientation_ID , &
grainrotation_ID , &
eulerangles_ID , &
defgrad_ID , &
fe_ID , &
fp_ID , &
2015-03-06 18:39:00 +05:30
fi_ID , &
2013-12-12 22:39:59 +05:30
lp_ID , &
2015-03-06 18:39:00 +05:30
li_ID , &
2013-12-12 22:39:59 +05:30
p_ID , &
s_ID , &
elasmatrix_ID , &
neighboringip_ID , &
neighboringelement_ID
end enum
2014-08-26 20:14:32 +05:30
integer ( kind ( undefined_ID ) ) , dimension ( : , : ) , allocatable , private :: &
2013-12-12 22:39:59 +05:30
crystallite_outputID !< ID of each post result output
2018-08-05 20:36:03 +05:30
procedure ( ) , pointer :: integrateState
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
public :: &
crystallite_init , &
crystallite_stressAndItsTangent , &
crystallite_orientations , &
2014-09-10 14:07:12 +05:30
crystallite_push33ToRef , &
2013-02-22 04:38:36 +05:30
crystallite_postResults
private :: &
2018-08-05 20:36:03 +05:30
integrateState , &
2018-09-20 09:57:53 +05:30
integrateStateFPI , &
integrateStateEuler , &
integrateStateAdaptiveEuler , &
integrateStateRK4 , &
integrateStateRKCK45 , &
integrateStress , &
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-02-02 17:06:09 +05:30
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
2017-10-05 20:05:34 +05:30
use , intrinsic :: iso_fortran_env , only : &
compiler_version , &
compiler_options
#endif
2018-09-20 09:39:02 +05:30
#ifdef DEBUG
2013-04-29 16:47:30 +05:30
use debug , only : &
debug_info , &
debug_reset , &
debug_level , &
debug_crystallite , &
debug_levelBasic
2018-09-20 09:39:02 +05:30
#endif
2013-04-16 22:37:27 +05:30
use numerics , only : &
2018-08-05 20:36:03 +05:30
numerics_integrator , &
2014-10-10 01:53:06 +05:30
worldrank , &
2014-08-26 20:14:32 +05:30
usePingPong
2013-05-17 23:22:46 +05:30
use math , only : &
2013-04-29 16:47:30 +05:30
math_I3 , &
math_EulerToR , &
math_inv33 , &
math_mul33xx33 , &
math_mul33x33
use FEsolving , only : &
FEsolving_execElem , &
FEsolving_execIP
use mesh , only : &
mesh_element , &
mesh_NcpElems , &
mesh_maxNips , &
mesh_maxNipNeighbors
use IO , only : &
IO_timeStamp , &
IO_stringValue , &
IO_write_jobFile , &
2018-06-02 22:58:08 +05:30
IO_error
2012-08-31 01:56:28 +05:30
use material
2018-06-27 00:03:41 +05:30
use config , only : &
2018-08-22 18:00:51 +05:30
config_deallocate , &
2018-06-27 00:24:54 +05:30
config_crystallite , &
2018-09-20 09:54:03 +05:30
crystallite_name
2013-04-29 16:47:30 +05:30
use constitutive , only : &
2015-07-24 20:17:18 +05:30
constitutive_initialFi , &
2018-09-07 02:13:29 +05:30
constitutive_microstructure ! derived (shortcut) quantities of given state
2014-08-26 20:14:32 +05:30
2012-08-31 01:56:28 +05:30
implicit none
2014-08-26 20:14:32 +05:30
2018-06-02 22:58:08 +05:30
integer ( pInt ) , parameter :: FILEUNIT = 434_pInt
2013-04-29 16:47:30 +05:30
integer ( pInt ) :: &
2016-01-17 23:26:24 +05:30
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e , & !< counter in element loop
2018-02-16 20:06:18 +05:30
o = 0_pInt , & !< counter in output loop
2018-09-07 02:13:29 +05:30
r , &
2016-01-17 23:26:24 +05:30
cMax , & !< maximum number of integration point components
2013-04-29 16:47:30 +05:30
iMax , & !< maximum number of integration points
eMax , & !< maximum number of elements
nMax , & !< maximum number of ip neighbors
2016-01-17 23:26:24 +05:30
myNcomponents , & !< number of components at current IP
2014-03-09 02:20:31 +05:30
mySize
2018-06-22 11:33:22 +05:30
character ( len = 65536 ) , dimension ( : ) , allocatable :: str
2014-08-26 20:14:32 +05:30
2016-06-30 02:57:22 +05:30
write ( 6 , '(/,a)' ) ' <<<+- crystallite init -+>>>'
write ( 6 , '(a15,a)' ) ' Current time: ' , IO_timeStamp ( )
2012-02-01 00:48:55 +05:30
#include "compilation_info.f90"
2014-08-26 20:14:32 +05:30
2016-01-17 20:20:33 +05:30
cMax = homogenization_maxNgrains
2012-08-31 01:56:28 +05:30
iMax = mesh_maxNips
eMax = mesh_NcpElems
nMax = mesh_maxNipNeighbors
2014-05-27 20:16:03 +05:30
2014-08-26 20:14:32 +05:30
2016-01-17 20:20:33 +05:30
allocate ( crystallite_Tstar0_v ( 6 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_partionedTstar0_v ( 6 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subTstar0_v ( 6 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_Tstar_v ( 6 , 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 ( 4 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_orientation0 ( 4 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_rotation ( 4 , cMax , iMax , eMax ) , source = 0.0_pReal )
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 . )
2012-08-31 01:56:28 +05:30
allocate ( crystallite_output ( maxval ( crystallite_Noutput ) , &
2018-06-27 00:24:54 +05:30
size ( config_crystallite ) ) ) ; crystallite_output = ''
2013-12-12 22:39:59 +05:30
allocate ( crystallite_outputID ( maxval ( crystallite_Noutput ) , &
2018-06-27 00:24:54 +05:30
size ( config_crystallite ) ) , source = undefined_ID )
allocate ( crystallite_sizePostResults ( size ( config_crystallite ) ) , source = 0_pInt )
2012-08-31 01:56:28 +05:30
allocate ( crystallite_sizePostResult ( maxval ( crystallite_Noutput ) , &
2018-06-27 00:24:54 +05:30
size ( config_crystallite ) ) , source = 0_pInt )
2016-01-17 23:26:24 +05:30
2018-08-05 20:36:03 +05:30
select case ( numerics_integrator ( 1 ) )
case ( 1_pInt )
2018-09-20 09:57:53 +05:30
integrateState = > integrateStateFPI
2018-08-05 20:36:03 +05:30
case ( 2_pInt )
2018-09-20 09:57:53 +05:30
integrateState = > integrateStateEuler
2018-08-05 20:36:03 +05:30
case ( 3_pInt )
2018-09-20 09:57:53 +05:30
integrateState = > integrateStateAdaptiveEuler
2018-08-05 20:36:03 +05:30
case ( 4_pInt )
2018-09-20 09:57:53 +05:30
integrateState = > integrateStateRK4
2018-08-05 20:36:03 +05:30
case ( 5_pInt )
2018-09-20 09:57:53 +05:30
integrateState = > integrateStateRKCK45
2018-08-05 20:36:03 +05:30
end select
2018-06-27 00:24:54 +05:30
do c = 1_pInt , size ( config_crystallite )
2018-06-22 11:33:22 +05:30
#if defined(__GFORTRAN__)
str = [ 'GfortranBug86277' ]
2018-06-27 00:24:54 +05:30
str = config_crystallite ( c ) % getStrings ( '(output)' , defaultVal = str )
2018-06-22 11:33:22 +05:30
if ( str ( 1 ) == 'GfortranBug86277' ) str = [ character ( len = 65536 ) :: ]
#else
2018-06-27 00:24:54 +05:30
str = config_crystallite ( c ) % getStrings ( '(output)' , defaultVal = [ character ( len = 65536 ) :: ] )
2018-06-22 11:33:22 +05:30
#endif
2018-06-02 22:58:08 +05:30
do o = 1_pInt , size ( str )
2018-06-03 00:30:47 +05:30
crystallite_output ( o , c ) = str ( o )
2018-06-02 22:58:08 +05:30
outputName : select case ( str ( o ) )
2016-01-17 23:26:24 +05:30
case ( 'phase' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = phase_ID
2016-01-17 23:26:24 +05:30
case ( 'texture' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = texture_ID
2016-01-17 23:26:24 +05:30
case ( 'volume' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = volume_ID
2016-01-17 23:26:24 +05:30
case ( 'orientation' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = orientation_ID
2016-01-17 23:26:24 +05:30
case ( 'grainrotation' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = grainrotation_ID
2016-01-17 23:26:24 +05:30
case ( 'eulerangles' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = eulerangles_ID
2016-01-17 23:26:24 +05:30
case ( 'defgrad' , 'f' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = defgrad_ID
2016-01-17 23:26:24 +05:30
case ( 'fe' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = fe_ID
2016-01-17 23:26:24 +05:30
case ( 'fp' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = fp_ID
2016-01-17 23:26:24 +05:30
case ( 'fi' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = fi_ID
2016-01-17 23:26:24 +05:30
case ( 'lp' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = lp_ID
2016-01-17 23:26:24 +05:30
case ( 'li' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = li_ID
2016-01-17 23:26:24 +05:30
case ( 'p' , 'firstpiola' , '1stpiola' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = p_ID
2016-01-17 23:26:24 +05:30
case ( 's' , 'tstar' , 'secondpiola' , '2ndpiola' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = s_ID
2016-01-17 23:26:24 +05:30
case ( 'elasmatrix' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = elasmatrix_ID
2016-01-17 23:26:24 +05:30
case ( 'neighboringip' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = neighboringip_ID
2016-01-17 23:26:24 +05:30
case ( 'neighboringelement' ) outputName
2018-06-02 22:58:08 +05:30
crystallite_outputID ( o , c ) = neighboringelement_ID
2016-01-17 23:26:24 +05:30
case default outputName
2018-09-20 10:28:31 +05:30
call IO_error ( 105_pInt , ext_msg = trim ( str ( o ) ) / / ' (Crystallite)' )
2016-01-17 23:26:24 +05:30
end select outputName
2018-06-20 02:28:46 +05:30
enddo
2012-08-31 01:56:28 +05:30
enddo
2014-08-26 20:14:32 +05:30
2018-06-27 00:24:54 +05:30
do r = 1_pInt , size ( config_crystallite )
2016-01-17 23:26:24 +05:30
do o = 1_pInt , crystallite_Noutput ( r )
select case ( crystallite_outputID ( o , r ) )
2018-09-20 10:05:30 +05:30
case ( phase_ID , texture_ID , volume_ID )
2012-08-31 01:56:28 +05:30
mySize = 1_pInt
2016-01-17 23:26:24 +05:30
case ( orientation_ID , grainrotation_ID )
2012-08-31 01:56:28 +05:30
mySize = 4_pInt
2013-12-12 22:39:59 +05:30
case ( eulerangles_ID )
2012-08-31 01:56:28 +05:30
mySize = 3_pInt
2018-09-20 10:05:30 +05:30
case ( defgrad_ID , fe_ID , fp_ID , fi_ID , lp_ID , li_ID , p_ID , s_ID )
2012-08-31 01:56:28 +05:30
mySize = 9_pInt
2013-12-12 22:39:59 +05:30
case ( elasmatrix_ID )
2014-08-26 20:14:32 +05:30
mySize = 36_pInt
2013-12-12 22:39:59 +05:30
case ( neighboringip_ID , neighboringelement_ID )
2013-05-08 17:32:30 +05:30
mySize = mesh_maxNipNeighbors
2012-08-31 01:56:28 +05:30
case default
2014-08-26 20:14:32 +05:30
mySize = 0_pInt
2012-08-31 01:56:28 +05:30
end select
2016-01-17 23:26:24 +05:30
crystallite_sizePostResult ( o , r ) = mySize
crystallite_sizePostResults ( r ) = crystallite_sizePostResults ( r ) + mySize
2012-08-31 01:56:28 +05:30
enddo
enddo
2014-08-26 20:14:32 +05:30
2016-01-17 23:26:24 +05:30
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
2015-03-25 21:32:30 +05:30
if ( worldrank == 0_pInt ) then
call IO_write_jobFile ( FILEUNIT , 'outputCrystallite' )
2014-08-26 20:14:32 +05:30
2018-06-27 00:24:54 +05:30
do r = 1_pInt , size ( config_crystallite )
2016-01-17 23:26:24 +05:30
if ( any ( microstructure_crystallite ( mesh_element ( 4 , : ) ) == r ) ) then
write ( FILEUNIT , '(/,a,/)' ) '[' / / trim ( crystallite_name ( r ) ) / / ']'
do o = 1_pInt , crystallite_Noutput ( r )
write ( FILEUNIT , '(a,i4)' ) trim ( crystallite_output ( o , r ) ) / / char ( 9 ) , crystallite_sizePostResult ( o , r )
2015-04-21 17:53:00 +05:30
enddo
endif
2012-08-31 01:56:28 +05:30
enddo
2014-08-26 20:14:32 +05:30
2015-03-25 21:32:30 +05:30
close ( FILEUNIT )
2015-10-14 00:22:01 +05:30
endif
2014-08-26 20:14:32 +05:30
2018-06-27 00:03:41 +05:30
call config_deallocate ( 'material.config/crystallite' )
2013-04-29 16:47:30 +05:30
!--------------------------------------------------------------------------------------------------
! initialize
2016-01-17 23:26:24 +05:30
!$OMP PARALLEL DO PRIVATE(myNcomponents)
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
myNcomponents = homogenization_Ngrains ( mesh_element ( 3 , e ) )
forall ( i = FEsolving_execIP ( 1 , e ) : FEsolving_execIP ( 2 , e ) , c = 1_pInt : myNcomponents )
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 ( math_mul33x33 ( crystallite_Fi0 ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e ) ) ) ! assuming that euler angles are given in internal strain free configuration
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 .
2014-05-27 20:16:03 +05:30
endforall
enddo
!$OMP END PARALLEL DO
2014-08-26 20:14:32 +05:30
2013-05-17 23:22:46 +05:30
if ( any ( . not . crystallite_localPlasticity ) . and . . not . usePingPong ) call IO_error ( 601_pInt ) ! exit if nonlocal but no ping-pong
2014-08-26 20:14:32 +05:30
2013-04-29 16:47:30 +05:30
crystallite_partionedFp0 = crystallite_Fp0
2015-03-06 18:39:00 +05:30
crystallite_partionedFi0 = crystallite_Fi0
2016-01-17 23:26:24 +05:30
crystallite_partionedF0 = crystallite_F0
2016-04-13 23:36:04 +05:30
crystallite_partionedF = crystallite_F0
2014-08-26 20:14:32 +05:30
2013-04-29 16:47:30 +05:30
call crystallite_orientations ( )
crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations
2009-12-22 17:58:02 +05:30
2018-09-19 23:15:57 +05:30
!$OMP PARALLEL DO
2013-04-29 16:47:30 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
2018-09-19 23:15:57 +05:30
do c = 1_pInt , homogenization_Ngrains ( mesh_element ( 3 , e ) )
2016-01-17 23:26:24 +05:30
call constitutive_microstructure ( crystallite_orientation , & ! pass orientation to constitutive module
crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) , &
2017-02-04 00:49:02 +05:30
c , i , e ) ! update dependent state variables to be consistent with basic states
2014-06-25 04:51:25 +05:30
enddo
2013-04-29 16:47:30 +05:30
enddo
enddo
!$OMP END PARALLEL DO
2014-06-17 00:49:38 +05:30
2014-08-08 02:38:34 +05:30
call crystallite_stressAndItsTangent ( . true . ) ! request elastic answers
2014-08-26 20:14:32 +05:30
2018-09-20 09:39:02 +05:30
#ifdef DEBUG
2013-04-29 16:47:30 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt ) then
2018-09-20 10:28:31 +05:30
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: ' , nMax
write ( 6 , '(a42,1x,i10)' ) ' # of nonlocal constituents: ' , count ( . not . crystallite_localPlasticity )
2013-04-29 16:47:30 +05:30
flush ( 6 )
endif
2010-11-03 22:52:48 +05:30
2012-08-31 01:56:28 +05:30
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
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P) and tangent (dPdF) for crystallites
!--------------------------------------------------------------------------------------------------
2014-08-08 02:38:34 +05:30
subroutine crystallite_stressAndItsTangent ( updateJaco )
2015-04-21 17:53:00 +05:30
use prec , only : &
2016-05-27 15:16:34 +05:30
tol_math_check , &
2016-10-29 13:09:08 +05:30
dNeq0
2013-04-29 16:47:30 +05:30
use numerics , only : &
subStepMinCryst , &
subStepSizeCryst , &
2019-01-14 17:26:46 +05:30
stepIncreaseCryst
2018-09-20 09:54:03 +05:30
#ifdef DEBUG
2013-04-29 16:47:30 +05:30
use debug , only : &
debug_level , &
debug_crystallite , &
debug_levelBasic , &
debug_levelExtensive , &
debug_levelSelective , &
debug_e , &
debug_i , &
2018-05-17 20:03:35 +05:30
debug_g
2018-09-20 09:54:03 +05:30
#endif
2013-04-29 16:47:30 +05:30
use IO , only : &
2014-08-08 02:46:37 +05:30
IO_warning , &
IO_error
2013-04-29 16:47:30 +05:30
use math , only : &
math_inv33 , &
math_identity2nd , &
math_mul33x33 , &
math_mul66x6 , &
math_Mandel6to33 , &
math_Mandel33to6 , &
2015-03-06 18:39:00 +05:30
math_Plain3333to99 , &
math_Plain99to3333 , &
2013-04-29 16:47:30 +05:30
math_I3 , &
2013-10-19 00:27:28 +05:30
math_mul3333xx3333 , &
2014-08-08 02:38:34 +05:30
math_mul33xx33 , &
2014-11-01 00:33:08 +05:30
math_invert , &
math_det33
2013-04-29 16:47:30 +05:30
use FEsolving , only : &
2014-08-26 20:14:32 +05:30
FEsolving_execElem , &
2013-04-29 16:47:30 +05:30
FEsolving_execIP
use mesh , only : &
mesh_element , &
mesh_maxNips , &
mesh_ipNeighborhood , &
FE_NipNeighbors , &
FE_geomtype , &
FE_cellType
use material , only : &
homogenization_Ngrains , &
2014-05-12 18:30:37 +05:30
plasticState , &
2015-05-28 22:32:23 +05:30
sourceState , &
phase_Nsources , &
2017-10-03 18:50:53 +05:30
phaseAt , phasememberAt
2013-04-29 16:47:30 +05:30
use constitutive , only : &
2018-08-28 18:24:36 +05:30
constitutive_SandItsTangents , &
2018-08-28 18:37:39 +05:30
constitutive_LpAndItsTangents , &
constitutive_LiAndItsTangents
2014-08-26 20:14:32 +05:30
2013-04-29 16:47:30 +05:30
implicit none
2013-05-17 23:22:46 +05:30
logical , intent ( in ) :: &
2014-08-08 02:38:34 +05:30
updateJaco !< whether to update the Jacobian (stiffness) or not
2013-05-17 23:22:46 +05:30
real ( pReal ) :: &
formerSubStep , &
subFracIntermediate
real ( pReal ) , dimension ( 3 , 3 ) :: &
invFp , & ! inverse of the plastic deformation gradient
Fe_guess , & ! guess for elastic deformation gradient
Tstar ! 2nd Piola-Kirchhoff stress tensor
integer ( pInt ) :: &
NiterationCrystallite , & ! number of iterations in crystallite loop
2016-01-17 23:26:24 +05:30
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e , & !< counter in element loop
2013-05-17 23:22:46 +05:30
n , startIP , endIP , &
neighboring_e , &
neighboring_i , &
o , &
p , &
2015-05-28 22:32:23 +05:30
mySource
2013-04-29 16:47:30 +05:30
! local variables used for calculating analytic Jacobian
2015-03-06 18:39:00 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: temp_33
2013-04-29 16:47:30 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dSdFe , &
dSdF , &
2015-03-06 18:39:00 +05:30
dSdFi , &
2014-11-01 00:33:08 +05:30
dLidS , &
2015-03-06 18:39:00 +05:30
dLidFi , &
2014-11-01 00:33:08 +05:30
dLpdS , &
2015-03-06 18:39:00 +05:30
dLpdFi , &
dFidS , &
2014-11-01 00:33:08 +05:30
dFpinvdF , &
rhs_3333 , &
lhs_3333 , &
temp_3333
2014-11-03 16:13:36 +05:30
real ( pReal ) , dimension ( 9 , 9 ) :: temp_99
2014-08-08 02:38:34 +05:30
logical :: error
2014-08-26 20:14:32 +05:30
2018-09-20 09:39:02 +05:30
#ifdef DEBUG
2014-05-19 22:13:32 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt &
. and . FEsolving_execElem ( 1 ) < = debug_e &
. and . debug_e < = FEsolving_execElem ( 2 ) ) then
2016-01-17 23:26:24 +05:30
write ( 6 , '(/,a,i8,1x,a,i8,a,1x,i2,1x,i3)' ) '<< CRYST >> boundary values at el ip ipc ' , &
2014-08-27 21:24:11 +05:30
debug_e , '(' , mesh_element ( 1 , debug_e ) , ')' , debug_i , debug_g
2015-10-15 00:06:19 +05:30
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST >> F ' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_partionedF ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
2014-05-27 20:16:03 +05:30
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST >> F0 ' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_partionedF0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
2014-05-27 20:16:03 +05:30
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST >> Fp0' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_partionedFp0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
2015-03-06 18:39:00 +05:30
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST >> Fi0' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_partionedFi0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
2014-05-27 20:16:03 +05:30
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST >> Lp0' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_partionedLp0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
2015-03-06 18:39:00 +05:30
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST >> Li0' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_partionedLi0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
2013-04-29 16:47:30 +05:30
endif
2018-09-20 09:39:02 +05:30
#endif
2013-04-29 16:47:30 +05:30
!--------------------------------------------------------------------------------------------------
! initialize to starting condition
crystallite_subStep = 0.0_pReal
2018-09-19 23:15:57 +05:30
!$OMP PARALLEL DO
2014-05-27 20:16:03 +05:30
elementLooping1 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2018-09-19 23:15:57 +05:30
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e ) ; do c = 1_pInt , homogenization_Ngrains ( mesh_element ( 3 , e ) )
2016-01-17 23:26:24 +05:30
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 ) )
2018-09-07 02:13:29 +05:30
2016-01-17 23:26:24 +05:30
do mySource = 1_pInt , phase_Nsources ( phaseAt ( c , i , e ) )
sourceState ( phaseAt ( c , i , e ) ) % p ( mySource ) % subState0 ( : , phasememberAt ( c , i , e ) ) = &
sourceState ( phaseAt ( c , i , e ) ) % p ( mySource ) % partionedState0 ( : , phasememberAt ( c , i , e ) )
2015-05-28 22:32:23 +05:30
enddo
2016-01-17 23:26:24 +05:30
crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedFp0 ( 1 : 3 , 1 : 3 , c , i , e ) ! ...plastic def grad
crystallite_subLp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedLp0 ( 1 : 3 , 1 : 3 , c , i , e ) ! ...plastic velocity grad
crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedFi0 ( 1 : 3 , 1 : 3 , c , i , e ) ! ...intermediate def grad
crystallite_subLi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedLi0 ( 1 : 3 , 1 : 3 , c , i , e ) ! ...intermediate velocity grad
crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedF0 ( 1 : 3 , 1 : 3 , c , i , e ) ! ...def grad
crystallite_subTstar0_v ( 1 : 6 , c , i , e ) = crystallite_partionedTstar0_v ( 1 : 6 , c , i , e ) !...2nd PK stress
crystallite_subFrac ( c , i , e ) = 0.0_pReal
crystallite_subStep ( c , i , e ) = 1.0_pReal / subStepSizeCryst
crystallite_todo ( c , i , e ) = . true .
crystallite_converged ( c , i , e ) = . false . ! pretend failed step of twice the required size
2015-05-28 22:32:23 +05:30
endif
enddo ; enddo
2014-05-27 20:16:03 +05:30
enddo elementLooping1
!$OMP END PARALLEL DO
singleRun : if ( FEsolving_execELem ( 1 ) == FEsolving_execElem ( 2 ) . and . &
2013-05-17 23:22:46 +05:30
FEsolving_execIP ( 1 , FEsolving_execELem ( 1 ) ) == FEsolving_execIP ( 2 , FEsolving_execELem ( 1 ) ) ) then
startIP = FEsolving_execIP ( 1 , FEsolving_execELem ( 1 ) )
endIP = startIP
2014-05-27 20:16:03 +05:30
else singleRun
2013-05-17 23:22:46 +05:30
startIP = 1_pInt
endIP = mesh_maxNips
2014-05-27 20:16:03 +05:30
endif singleRun
2009-12-22 17:58:02 +05:30
2013-04-29 16:47:30 +05:30
NiterationCrystallite = 0_pInt
2013-05-17 23:22:46 +05:30
cutbackLooping : do while ( any ( crystallite_todo ( : , startIP : endIP , FEsolving_execELem ( 1 ) : FEsolving_execElem ( 2 ) ) ) )
2018-09-20 09:54:03 +05:30
#ifdef DEBUG
2015-08-05 02:56:22 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,i6)' ) '<< CRYST >> crystallite iteration ' , NiterationCrystallite
2018-09-20 09:54:03 +05:30
#endif
2015-08-05 02:56:22 +05:30
2018-09-19 23:15:57 +05:30
!$OMP PARALLEL DO PRIVATE(formerSubStep)
2013-10-19 00:27:28 +05:30
elementLooping3 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2013-04-29 16:47:30 +05:30
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e ) ! iterate over IPs of this element to be processed
2018-09-19 23:15:57 +05:30
do c = 1 , homogenization_Ngrains ( mesh_element ( 3 , e ) )
2013-04-29 16:47:30 +05:30
! --- wind forward ---
2014-08-26 20:14:32 +05:30
2019-01-15 12:22:01 +05:30
if ( crystallite_converged ( c , i , e ) ) then
2016-01-17 23:26:24 +05:30
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 ) , &
stepIncreaseCryst * crystallite_subStep ( c , i , e ) )
2018-09-19 23:15:57 +05:30
2016-01-17 23:26:24 +05:30
if ( crystallite_subStep ( c , i , e ) > 0.0_pReal ) then
crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) ! ...def grad
crystallite_subLp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e ) ! ...plastic velocity gradient
crystallite_subLi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e ) ! ...intermediate velocity gradient
crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) ! ...plastic def grad
crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) ! ...intermediate def grad
2014-05-27 20:16:03 +05:30
!if abbrevation, make c and p private in omp
2016-01-17 23:26:24 +05:30
plasticState ( phaseAt ( c , i , e ) ) % subState0 ( : , phasememberAt ( c , i , e ) ) = &
plasticState ( phaseAt ( c , i , e ) ) % state ( : , phasememberAt ( c , i , e ) )
do mySource = 1_pInt , phase_Nsources ( phaseAt ( c , i , e ) )
sourceState ( phaseAt ( c , i , e ) ) % p ( mySource ) % subState0 ( : , phasememberAt ( c , i , e ) ) = &
sourceState ( phaseAt ( c , i , e ) ) % p ( mySource ) % state ( : , phasememberAt ( c , i , e ) )
2015-05-28 22:32:23 +05:30
enddo
2016-01-17 23:26:24 +05:30
crystallite_subTstar0_v ( 1 : 6 , c , i , e ) = crystallite_Tstar_v ( 1 : 6 , c , i , e ) ! ...2nd PK stress
2019-01-14 17:26:46 +05:30
crystallite_todo ( c , i , e ) = . true .
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2013-04-29 16:47:30 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt &
2016-01-17 23:26:24 +05:30
. and . ( ( e == debug_e . and . i == debug_i . and . c == debug_g ) &
2013-10-19 00:27:28 +05:30
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) &
write ( 6 , '(a,f12.8,a,f12.8,a,i8,1x,i2,1x,i3,/)' ) '<< CRYST >> winding forward from ' , &
2016-01-17 23:26:24 +05:30
crystallite_subFrac ( c , i , e ) - formerSubStep , ' to current crystallite_subfrac ' , &
crystallite_subFrac ( c , i , e ) , ' in crystallite_stressAndItsTangent at el ip ipc ' , e , i , c
2011-03-29 12:57:19 +05:30
#endif
2013-08-02 13:29:55 +05:30
else ! this crystallite just converged for the entire timestep
2016-01-17 23:26:24 +05:30
crystallite_todo ( c , i , e ) = . false . ! so done here
2013-04-29 16:47:30 +05:30
endif
2014-08-26 20:14:32 +05:30
2013-04-29 16:47:30 +05:30
! --- cutback ---
2014-08-26 20:14:32 +05:30
2019-01-15 12:22:01 +05:30
elseif ( . not . crystallite_converged ( c , i , e ) ) then
2019-01-14 17:26:46 +05:30
crystallite_subStep ( c , i , e ) = subStepSizeCryst * crystallite_subStep ( c , i , e ) ! cut step in half and restore...
2016-01-17 23:26:24 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) ! ...plastic def grad
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 ) ! ...intermediate def grad
crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) = math_inv33 ( crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) )
crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subLp0 ( 1 : 3 , 1 : 3 , c , i , e ) ! ...plastic velocity grad
crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subLi0 ( 1 : 3 , 1 : 3 , c , i , e ) ! ...intermediate velocity grad
plasticState ( phaseAt ( c , i , e ) ) % state ( : , phasememberAt ( c , i , e ) ) = &
plasticState ( phaseAt ( c , i , e ) ) % subState0 ( : , phasememberAt ( c , i , e ) )
do mySource = 1_pInt , phase_Nsources ( phaseAt ( c , i , e ) )
sourceState ( phaseAt ( c , i , e ) ) % p ( mySource ) % state ( : , phasememberAt ( c , i , e ) ) = &
sourceState ( phaseAt ( c , i , e ) ) % p ( mySource ) % subState0 ( : , phasememberAt ( c , i , e ) )
2015-05-28 22:32:23 +05:30
enddo
2016-01-17 23:26:24 +05:30
crystallite_Tstar_v ( 1 : 6 , c , i , e ) = crystallite_subTstar0_v ( 1 : 6 , c , i , e ) ! ...2nd PK stress
2014-05-27 20:16:03 +05:30
2013-04-29 16:47:30 +05:30
! cant restore dotState here, since not yet calculated in first cutback after initialization
2016-01-17 23:26:24 +05:30
crystallite_todo ( c , i , e ) = crystallite_subStep ( c , i , e ) > subStepMinCryst ! still on track or already done (beyond repair)
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2017-12-14 05:48:45 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( e == debug_e . and . i == debug_i . and . c == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
2016-01-17 23:26:24 +05:30
if ( crystallite_todo ( c , i , e ) ) then
2014-07-10 14:17:00 +05:30
write ( 6 , '(a,f12.8,a,i8,1x,i2,1x,i3,/)' ) ' < < CRYST > > cutback step in crystallite_stressAndItsTangent &
2013-04-29 16:47:30 +05:30
& with new crystallite_subStep : ' , &
2016-01-17 23:26:24 +05:30
crystallite_subStep ( c , i , e ) , ' at el ip ipc ' , e , i , c
2013-04-29 16:47:30 +05:30
else
2014-07-10 14:17:00 +05:30
write ( 6 , '(a,i8,1x,i2,1x,i3,/)' ) ' < < CRYST > > reached minimum step size &
2016-01-17 23:26:24 +05:30
& in crystallite_stressAndItsTangent at el ip ipc ' , e , i , c
2013-04-29 16:47:30 +05:30
endif
endif
2012-11-22 15:28:36 +05:30
#endif
2013-04-29 16:47:30 +05:30
endif
2014-08-26 20:14:32 +05:30
2013-04-29 16:47:30 +05:30
! --- prepare for integration ---
2014-08-26 20:14:32 +05:30
2019-01-15 12:22:01 +05:30
if ( crystallite_todo ( c , i , e ) ) then
2016-01-17 23:26:24 +05:30
crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) &
2018-09-19 23:15:57 +05:30
+ crystallite_subStep ( c , i , e ) * ( crystallite_partionedF ( 1 : 3 , 1 : 3 , c , i , e ) &
2016-01-17 23:26:24 +05:30
- crystallite_partionedF0 ( 1 : 3 , 1 : 3 , c , i , e ) )
crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) = math_mul33x33 ( math_mul33x33 ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) , &
crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) )
crystallite_subdt ( c , i , e ) = crystallite_subStep ( c , i , e ) * crystallite_dt ( c , i , e )
crystallite_converged ( c , i , e ) = . false . ! start out non-converged
2013-04-29 16:47:30 +05:30
endif
2014-08-26 20:14:32 +05:30
2013-04-29 16:47:30 +05:30
enddo ! grains
enddo ! IPs
2013-10-19 00:27:28 +05:30
enddo elementLooping3
2013-04-29 16:47:30 +05:30
!$OMP END PARALLEL DO
2014-08-26 20:14:32 +05:30
2018-09-20 09:39:02 +05:30
#ifdef DEBUG
2013-04-29 16:47:30 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) then
2017-11-07 04:39:04 +05:30
write ( 6 , '(/,a,f8.5)' ) '<< CRYST >> min(subStep) ' , minval ( crystallite_subStep )
write ( 6 , '(a,f8.5)' ) '<< CRYST >> max(subStep) ' , maxval ( crystallite_subStep )
write ( 6 , '(a,f8.5)' ) '<< CRYST >> min(subFrac) ' , minval ( crystallite_subFrac )
write ( 6 , '(a,f8.5,/)' ) '<< CRYST >> max(subFrac) ' , maxval ( crystallite_subFrac )
2014-05-27 20:16:03 +05:30
flush ( 6 )
2017-12-14 05:48:45 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) then
write ( 6 , '(/,a,f8.5,1x,a,1x,f8.5,1x,a)' ) '<< CRYST >> subFrac + subStep = ' , &
crystallite_subFrac ( debug_g , debug_i , debug_e ) , '+' , crystallite_subStep ( debug_g , debug_i , debug_e ) , '@selective'
flush ( 6 )
endif
2013-04-29 16:47:30 +05:30
endif
2018-09-20 09:39:02 +05:30
#endif
2014-08-26 20:14:32 +05:30
2013-04-29 16:47:30 +05:30
! --- integrate --- requires fully defined state array (basic + dependent state)
2018-08-05 20:36:03 +05:30
if ( any ( crystallite_todo ) ) call integrateState ( )
2013-04-29 16:47:30 +05:30
where ( . not . crystallite_converged . and . crystallite_subStep > subStepMinCryst ) & ! do not try non-converged & fully cutbacked any further
crystallite_todo = . true .
2014-08-26 20:14:32 +05:30
2013-04-29 16:47:30 +05:30
NiterationCrystallite = NiterationCrystallite + 1_pInt
2014-08-26 20:14:32 +05:30
2013-04-29 16:47:30 +05:30
enddo cutbackLooping
2010-11-03 22:52:48 +05:30
! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+--
2013-11-21 16:28:41 +05:30
2013-10-19 00:27:28 +05:30
elementLooping5 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2018-09-07 02:13:29 +05:30
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e ) ! iterate over IPs of this element to be processed
2018-09-19 23:15:57 +05:30
do c = 1 , homogenization_Ngrains ( mesh_element ( 3 , e ) )
2018-09-07 02:13:29 +05:30
if ( . not . crystallite_converged ( c , i , e ) ) then ! respond fully elastically (might be not required due to becoming terminally ill anyway)
2018-09-20 09:39:02 +05:30
#ifdef DEBUG
2014-05-27 20:16:03 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt ) &
2016-01-17 23:26:24 +05:30
write ( 6 , '(a,i8,1x,a,i8,a,1x,i2,1x,i3,/)' ) '<< CRYST >> no convergence: respond fully elastic at el (elFE) ip ipc ' , &
e , '(' , mesh_element ( 1 , e ) , ')' , i , c
2018-09-20 09:39:02 +05:30
#endif
2016-01-17 23:26:24 +05:30
invFp = math_inv33 ( crystallite_partionedFp0 ( 1 : 3 , 1 : 3 , c , i , e ) )
Fe_guess = math_mul33x33 ( math_mul33x33 ( crystallite_partionedF ( 1 : 3 , 1 : 3 , c , i , e ) , invFp ) , &
math_inv33 ( crystallite_partionedFi0 ( 1 : 3 , 1 : 3 , c , i , e ) ) )
2018-08-28 18:24:36 +05:30
call constitutive_SandItsTangents ( Tstar , dSdFe , dSdFi , Fe_guess , crystallite_partionedFi0 ( 1 : 3 , 1 : 3 , c , i , e ) , c , i , e )
2016-01-17 23:26:24 +05:30
crystallite_P ( 1 : 3 , 1 : 3 , c , i , e ) = math_mul33x33 ( math_mul33x33 ( crystallite_partionedF ( 1 : 3 , 1 : 3 , c , i , e ) , invFp ) , &
2015-10-14 00:22:01 +05:30
math_mul33x33 ( Tstar , transpose ( invFp ) ) )
2013-05-17 23:22:46 +05:30
endif
2018-09-20 09:39:02 +05:30
#ifdef DEBUG
2014-08-26 20:14:32 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
2016-01-17 23:26:24 +05:30
. and . ( ( e == debug_e . and . i == debug_i . and . c == debug_g ) &
2012-11-07 21:13:29 +05:30
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
2016-01-17 23:26:24 +05:30
write ( 6 , '(a,i8,1x,i2,1x,i3)' ) '<< CRYST >> central solution of cryst_StressAndTangent at el ip ipc ' , e , i , c
2013-05-17 23:22:46 +05:30
write ( 6 , '(/,a,/,3(12x,3(f12.4,1x)/))' ) '<< CRYST >> P / MPa' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_P ( 1 : 3 , 1 : 3 , c , i , e ) ) * 1.0e-6_pReal
2013-05-17 23:22:46 +05:30
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST >> Fp' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) )
2015-03-06 18:39:00 +05:30
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST >> Fi' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) )
2013-05-17 23:22:46 +05:30
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/),/)' ) '<< CRYST >> Lp' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e ) )
2015-03-06 18:39:00 +05:30
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/),/)' ) '<< CRYST >> Li' , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e ) )
2014-05-27 20:16:03 +05:30
flush ( 6 )
2013-05-17 23:22:46 +05:30
endif
2018-09-20 09:39:02 +05:30
#endif
2013-05-17 23:22:46 +05:30
enddo
enddo
2013-10-19 00:27:28 +05:30
enddo elementLooping5
2010-11-03 22:52:48 +05:30
2014-09-03 01:16:52 +05:30
2010-11-03 22:52:48 +05:30
! --+>> STIFFNESS CALCULATION <<+--
2009-06-16 14:33:30 +05:30
2013-10-19 00:27:28 +05:30
computeJacobian : if ( updateJaco ) then
2016-07-25 23:37:12 +05:30
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,&
2018-09-19 23:15:57 +05:30
!$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,error)
2016-07-25 23:37:12 +05:30
elementLooping6 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e ) ! iterate over IPs of this element to be processed
2018-09-19 23:15:57 +05:30
do c = 1_pInt , homogenization_Ngrains ( mesh_element ( 3 , e ) )
2018-08-28 18:24:36 +05:30
call constitutive_SandItsTangents ( temp_33 , dSdFe , dSdFi , crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) , &
2016-07-25 23:37:12 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) , c , i , e ) ! call constitutive law to calculate elastic stress tangent
2018-08-28 18:37:39 +05:30
call constitutive_LiAndItsTangents ( temp_33 , dLidS , dLidFi , crystallite_Tstar_v ( 1 : 6 , c , i , e ) , &
2016-07-25 23:37:12 +05:30
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
temp_33 = math_inv33 ( crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) )
lhs_3333 = 0.0_pReal ; rhs_3333 = 0.0_pReal
do o = 1_pInt , 3_pInt ; do p = 1_pInt , 3_pInt
lhs_3333 ( 1 : 3 , 1 : 3 , o , p ) = lhs_3333 ( 1 : 3 , 1 : 3 , o , p ) + &
crystallite_subdt ( c , i , e ) * math_mul33x33 ( temp_33 , 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 ) * math_mul33x33 ( temp_33 , dLidS ( 1 : 3 , 1 : 3 , o , p ) )
enddo ; enddo
call math_invert ( 9_pInt , math_Plain3333to99 ( lhs_3333 ) , temp_99 , error )
if ( error ) then
call IO_warning ( warning_ID = 600_pInt , el = e , ip = i , g = c , &
ext_msg = 'inversion error in analytic tangent calculation' )
dFidS = 0.0_pReal
else
dFidS = math_mul3333xx3333 ( math_Plain99to3333 ( temp_99 ) , rhs_3333 )
endif
dLidS = math_mul3333xx3333 ( dLidFi , dFidS ) + dLidS
endif
2016-01-17 23:26:24 +05:30
2018-08-28 18:37:39 +05:30
call constitutive_LpAndItsTangents ( temp_33 , dLpdS , dLpdFi , crystallite_Tstar_v ( 1 : 6 , c , i , e ) , &
2016-07-25 23:37:12 +05:30
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
2018-06-02 22:58:08 +05:30
temp_33 = transpose ( math_mul33x33 ( crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) , &
2016-07-25 23:37:12 +05:30
crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) ) )
rhs_3333 = 0.0_pReal
forall ( p = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
rhs_3333 ( p , o , 1 : 3 , 1 : 3 ) = math_mul33x33 ( dSdFe ( p , o , 1 : 3 , 1 : 3 ) , temp_33 )
temp_3333 = 0.0_pReal
temp_33 = math_mul33x33 ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
math_inv33 ( crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) ) )
forall ( p = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
temp_3333 ( 1 : 3 , 1 : 3 , p , o ) = math_mul33x33 ( math_mul33x33 ( temp_33 , dLpdS ( 1 : 3 , 1 : 3 , p , o ) ) , &
crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) )
temp_33 = math_mul33x33 ( math_mul33x33 ( 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 ) ) )
forall ( p = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
temp_3333 ( 1 : 3 , 1 : 3 , p , o ) = temp_3333 ( 1 : 3 , 1 : 3 , p , o ) + math_mul33x33 ( temp_33 , dLidS ( 1 : 3 , 1 : 3 , p , o ) )
lhs_3333 = crystallite_subdt ( c , i , e ) * math_mul3333xx3333 ( dSdFe , temp_3333 ) + &
math_mul3333xx3333 ( dSdFi , dFidS )
call math_invert ( 9_pInt , math_identity2nd ( 9_pInt ) + math_Plain3333to99 ( lhs_3333 ) , temp_99 , error )
if ( error ) then
call IO_warning ( warning_ID = 600_pInt , el = e , ip = i , g = c , &
ext_msg = 'inversion error in analytic tangent calculation' )
dSdF = rhs_3333
else
dSdF = math_mul3333xx3333 ( math_Plain99to3333 ( temp_99 ) , rhs_3333 )
endif
2015-05-28 22:32:23 +05:30
2016-07-25 23:37:12 +05:30
dFpinvdF = 0.0_pReal
temp_3333 = math_mul3333xx3333 ( dLpdS , dSdF )
forall ( p = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) = - crystallite_subdt ( c , i , e ) * &
math_mul33x33 ( math_inv33 ( crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) ) , &
math_mul33x33 ( temp_3333 ( 1 : 3 , 1 : 3 , p , o ) , &
crystallite_invFi ( 1 : 3 , 1 : 3 , c , i , e ) ) )
crystallite_dPdF ( 1 : 3 , 1 : 3 , 1 : 3 , 1 : 3 , c , i , e ) = 0.0_pReal
temp_33 = math_mul33x33 ( crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) , &
math_mul33x33 ( math_Mandel6to33 ( crystallite_Tstar_v ( 1 : 6 , c , i , e ) ) , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) ) )
2016-07-25 23:37:12 +05:30
forall ( p = 1_pInt : 3_pInt ) &
2018-06-02 22:58:08 +05:30
crystallite_dPdF ( p , 1 : 3 , p , 1 : 3 , c , i , e ) = transpose ( temp_33 )
2016-07-25 23:37:12 +05:30
temp_33 = math_mul33x33 ( math_Mandel6to33 ( crystallite_Tstar_v ( 1 : 6 , c , i , e ) ) , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) )
2016-07-25 23:37:12 +05:30
forall ( p = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
crystallite_dPdF ( 1 : 3 , 1 : 3 , p , o , c , i , e ) = crystallite_dPdF ( 1 : 3 , 1 : 3 , p , o , c , i , e ) + &
math_mul33x33 ( math_mul33x33 ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) , temp_33 )
temp_33 = math_mul33x33 ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) )
forall ( p = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
crystallite_dPdF ( 1 : 3 , 1 : 3 , p , o , c , i , e ) = crystallite_dPdF ( 1 : 3 , 1 : 3 , p , o , c , i , e ) + &
math_mul33x33 ( math_mul33x33 ( temp_33 , dSdF ( 1 : 3 , 1 : 3 , p , o ) ) , &
2018-06-02 22:58:08 +05:30
transpose ( crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) )
2016-07-25 23:37:12 +05:30
temp_33 = math_mul33x33 ( math_mul33x33 ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_invFp ( 1 : 3 , 1 : 3 , c , i , e ) ) , &
math_Mandel6to33 ( crystallite_Tstar_v ( 1 : 6 , c , i , e ) ) )
forall ( p = 1_pInt : 3_pInt , o = 1_pInt : 3_pInt ) &
crystallite_dPdF ( 1 : 3 , 1 : 3 , p , o , c , i , e ) = crystallite_dPdF ( 1 : 3 , 1 : 3 , p , o , c , i , e ) + &
2018-06-02 22:58:08 +05:30
math_mul33x33 ( temp_33 , transpose ( dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) )
2014-09-03 01:16:52 +05:30
2016-07-25 23:37:12 +05:30
enddo ; enddo
enddo elementLooping6
!$OMP END PARALLEL DO
2013-10-19 00:27:28 +05:30
endif computeJacobian
2012-03-09 01:55:28 +05:30
end subroutine crystallite_stressAndItsTangent
2009-05-07 21:57:36 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +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
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
logical function stateJump ( ipc , ip , el )
2017-05-04 04:02:44 +05:30
use , intrinsic :: &
IEEE_arithmetic
2019-01-15 08:57:57 +05:30
use prec , only : &
dNeq0
2018-02-16 20:06:18 +05:30
#ifdef DEBUG
2018-09-20 09:54:03 +05:30
use debug , only : &
2018-02-16 20:06:18 +05:30
debug_e , &
debug_i , &
debug_g , &
2013-10-19 00:27:28 +05:30
debug_level , &
debug_crystallite , &
debug_levelExtensive , &
2018-05-17 20:03:35 +05:30
debug_levelSelective
2018-09-20 09:54:03 +05:30
#endif
2013-10-19 00:27:28 +05:30
use material , only : &
2014-05-12 18:30:37 +05:30
plasticState , &
2015-05-28 22:32:23 +05:30
sourceState , &
phase_Nsources , &
2016-01-15 05:49:44 +05:30
phaseAt , phasememberAt
2013-10-19 00:27:28 +05:30
use constitutive , only : &
2019-01-15 08:57:57 +05:30
constitutive_collectDeltaState
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
implicit none
2019-01-15 08:57:57 +05:30
integer ( pInt ) , intent ( in ) :: &
el , & ! element index
ip , & ! integration point index
ipc ! grain index
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
integer ( pInt ) :: &
c , &
p , &
mySource , &
myOffsetPlasticDeltaState , &
myOffsetSourceDeltaState , &
mySizePlasticDeltaState , &
mySizeSourceDeltaState
2014-05-27 20:16:03 +05:30
2019-01-15 08:57:57 +05:30
c = phasememberAt ( ipc , ip , el )
p = phaseAt ( ipc , ip , el )
2014-05-27 20:16:03 +05:30
2019-01-15 08:57:57 +05:30
call constitutive_collectDeltaState ( crystallite_Tstar_v ( 1 : 6 , ipc , ip , el ) , &
crystallite_Fe ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
ipc , ip , el )
2014-05-27 20:16:03 +05:30
2019-01-15 08:57:57 +05:30
myOffsetPlasticDeltaState = plasticState ( p ) % offsetDeltaState
mySizePlasticDeltaState = plasticState ( p ) % sizeDeltaState
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
if ( any ( IEEE_is_NaN ( plasticState ( p ) % deltaState ( 1 : mySizePlasticDeltaState , c ) ) ) ) then ! NaN occured in deltaState
stateJump = . false .
return
endif
2014-05-27 20:16:03 +05:30
2019-01-15 08:57:57 +05:30
plasticState ( p ) % state ( myOffsetPlasticDeltaState + 1_pInt : &
myOffsetPlasticDeltaState + mySizePlasticDeltaState , c ) = &
plasticState ( p ) % state ( myOffsetPlasticDeltaState + 1_pInt : &
myOffsetPlasticDeltaState + mySizePlasticDeltaState , c ) + &
plasticState ( p ) % deltaState ( 1 : mySizePlasticDeltaState , c )
2014-07-02 17:57:39 +05:30
2019-01-15 08:57:57 +05:30
do mySource = 1_pInt , phase_Nsources ( p )
myOffsetSourceDeltaState = sourceState ( p ) % p ( mySource ) % offsetDeltaState
mySizeSourceDeltaState = sourceState ( p ) % p ( mySource ) % sizeDeltaState
if ( any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % deltaState ( 1 : mySizeSourceDeltaState , c ) ) ) ) then ! NaN occured in deltaState
stateJump = . false .
return
endif
sourceState ( p ) % p ( mySource ) % state ( myOffsetSourceDeltaState + 1_pInt : &
myOffsetSourceDeltaState + mySizeSourceDeltaState , c ) = &
sourceState ( p ) % p ( mySource ) % state ( myOffsetSourceDeltaState + 1_pInt : &
myOffsetSourceDeltaState + mySizeSourceDeltaState , c ) + &
sourceState ( p ) % p ( mySource ) % deltaState ( 1 : mySizeSourceDeltaState , c )
enddo
2014-09-03 01:16:52 +05:30
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
if ( any ( dNeq0 ( plasticState ( p ) % deltaState ( 1 : mySizePlasticDeltaState , c ) ) ) &
. and . iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,i8,1x,i2,1x,i3, /)' ) '<< CRYST >> update state at el ip ipc ' , el , ip , ipc
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> deltaState' , plasticState ( p ) % deltaState ( 1 : mySizePlasticDeltaState , c )
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> new state' , &
plasticState ( p ) % state ( myOffsetPlasticDeltaState + 1_pInt : &
myOffsetPlasticDeltaState + mySizePlasticDeltaState , c )
endif
2014-09-03 01:16:52 +05:30
#endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
stateJump = . true .
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
end function stateJump
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 )
use math , only : &
math_mul33x33 , &
math_inv33 , &
math_EulerToR
use material , only : &
material_EulerAngles
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 ) :: crystallite_push33ToRef
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: tensor33
real ( pReal ) , dimension ( 3 , 3 ) :: T
integer ( pInt ) , intent ( in ) :: &
el , & ! element index
ip , & ! integration point index
ipc ! grain index
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
T = math_mul33x33 ( math_EulerToR ( material_EulerAngles ( 1 : 3 , ipc , ip , el ) ) , &
transpose ( math_inv33 ( crystallite_subF ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) ) )
crystallite_push33ToRef = math_mul33x33 ( transpose ( T ) , math_mul33x33 ( 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-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
!--------------------------------------------------------------------------------------------------
logical function integrateStress ( &
ipc , & ! grain number
ip , & ! integration point number
el , & ! element number
timeFraction &
)
use , intrinsic :: &
IEEE_arithmetic
use prec , only : pLongInt , &
tol_math_check , &
dEq0
use numerics , only : nStress , &
aTol_crystalliteStress , &
rTol_crystalliteStress , &
iJacoLpresiduum , &
subStepSizeLp , &
subStepSizeLi
#ifdef DEBUG
use debug , only : debug_level , &
debug_e , &
debug_i , &
debug_g , &
debug_crystallite , &
debug_levelBasic , &
debug_levelExtensive , &
debug_levelSelective
#endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
use constitutive , only : constitutive_LpAndItsTangents , &
constitutive_LiAndItsTangents , &
constitutive_SandItsTangents
use math , only : math_mul33x33 , &
math_mul33xx33 , &
math_mul3333xx3333 , &
math_mul66x6 , &
math_mul99x99 , &
math_inv33 , &
math_invert , &
math_det33 , &
math_I3 , &
math_identity2nd , &
math_Mandel66to3333 , &
math_Mandel6to33 , &
math_Mandel33to6 , &
math_Plain3333to99 , &
math_Plain33to9 , &
math_Plain9to33 , &
math_Plain99to3333
#ifdef DEBUG
use mesh , only : mesh_element
#endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
implicit none
integer ( pInt ) , intent ( in ) :: el , & ! element index
ip , & ! integration point index
ipc ! grain index
real ( pReal ) , optional , intent ( in ) :: timeFraction ! fraction of timestep
2015-10-14 00:22:01 +05:30
2019-01-15 08:57:57 +05:30
!*** local variables ***!
real ( pReal ) , dimension ( 3 , 3 ) :: Fg_new , & ! deformation gradient at end of timestep
Fp_current , & ! plastic deformation gradient at start of timestep
Fi_current , & ! intermediate deformation gradient at start 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
Tstar , & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
A , &
B , &
Fe , & ! elastic deformation gradient
temp_33
real ( pReal ) , dimension ( 6 ) :: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation
real ( pReal ) , dimension ( 9 ) :: work ! needed for matrix inversion by LAPACK
integer ( pInt ) , dimension ( 9 ) :: ipiv ! 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 ( pInt ) 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
2014-05-27 20:16:03 +05:30
2019-01-15 08:57:57 +05:30
!* be pessimistic
integrateStress = . false .
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) &
write ( 6 , '(a,i8,1x,i2,1x,i3)' ) '<< CRYST >> integrateStress at el ip ipc ' , el , ip , ipc
#endif
2014-07-02 17:57:39 +05:30
2019-01-15 08:57:57 +05:30
!* only integrate over fraction of timestep?
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +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
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
Fp_current = crystallite_subFp0 ( 1 : 3 , 1 : 3 , ipc , ip , el ) ! "Fp_current" is only used as temp var here...
Lpguess = crystallite_Lp ( 1 : 3 , 1 : 3 , ipc , ip , el ) ! ... and take it as first guess
Fi_current = crystallite_subFi0 ( 1 : 3 , 1 : 3 , ipc , ip , el ) ! intermediate configuration, assume decomposition as F = Fe Fi Fp
Liguess = crystallite_Li ( 1 : 3 , 1 : 3 , ipc , ip , el ) ! ... and take it as first guess
Liguess_old = Liguess
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!* inversion of Fp_current...
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
invFp_current = math_inv33 ( Fp_current )
failedInversionFp : if ( all ( dEq0 ( invFp_current ) ) ) then
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt ) then
write ( 6 , '(a,i8,1x,a,i8,a,1x,i2,1x,i3)' ) '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip ipc ' , &
el , '(' , mesh_element ( 1 , el ) , ')' , ip , ipc
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) > 0_pInt ) &
write ( 6 , '(/,a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST >> Fp_current' , transpose ( Fp_current ( 1 : 3 , 1 : 3 ) )
2013-02-22 04:38:36 +05:30
endif
2019-01-15 08:57:57 +05:30
#endif
return
endif failedInversionFp
A = math_mul33x33 ( Fg_new , invFp_current ) ! intermediate tensor needed later to calculate dFe_dLp
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!* inversion of Fi_current...
invFi_current = math_inv33 ( Fi_current )
failedInversionFi : if ( all ( dEq0 ( invFi_current ) ) ) then
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt ) then
write ( 6 , '(a,i8,1x,a,i8,a,1x,i2,1x,i3)' ) '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip ipc ' , &
el , '(' , mesh_element ( 1 , el ) , ')' , ip , ipc
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) > 0_pInt ) &
write ( 6 , '(/,a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST >> Fp_current' , transpose ( Fi_current ( 1 : 3 , 1 : 3 ) )
endif
#endif
return
endif failedInversionFi
2010-10-01 17:48:49 +05:30
2019-01-15 08:57:57 +05:30
!* start LpLoop with normal step length
2010-10-01 17:48:49 +05:30
2019-01-15 08:57:57 +05:30
NiterationStressLi = 0_pInt
jacoCounterLi = 0_pInt
steplengthLi = 1.0_pReal
residuumLi_old = 0.0_pReal
LiLoop : do
NiterationStressLi = NiterationStressLi + 1_pInt
IloopsExeced : if ( NiterationStressLi > nStress ) then
2018-02-16 20:06:18 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt ) &
write ( 6 , '(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)' ) '<< CRYST >> integrateStress reached inelastic loop limit' , nStress , &
' at el (elFE) ip ipc ' , el , '(' , mesh_element ( 1 , el ) , ')' , ip , ipc
2018-09-20 09:54:03 +05:30
#endif
2019-01-15 08:57:57 +05:30
return
endif IloopsExeced
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
invFi_new = math_mul33x33 ( invFi_current , math_I3 - dt * Liguess )
Fi_new = math_inv33 ( invFi_new )
detInvFi = math_det33 ( invFi_new )
2013-11-21 16:28:41 +05:30
2019-01-15 08:57:57 +05:30
NiterationStressLp = 0_pInt
jacoCounterLp = 0_pInt
steplengthLp = 1.0_pReal
residuumLp_old = 0.0_pReal
Lpguess_old = Lpguess
2013-11-21 16:28:41 +05:30
2019-01-15 08:57:57 +05:30
LpLoop : do ! inner stress integration loop for consistency with Fi
NiterationStressLp = NiterationStressLp + 1_pInt
loopsExeced : if ( NiterationStressLp > nStress ) then
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt ) &
write ( 6 , '(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)' ) '<< CRYST >> integrateStress reached loop limit' , nStress , &
' at el (elFE) ip ipc ' , el , '(' , mesh_element ( 1 , el ) , ')' , ip , ipc
#endif
return
endif loopsExeced
2013-11-21 16:28:41 +05:30
2019-01-15 08:57:57 +05:30
!* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
B = math_I3 - dt * Lpguess
Fe = math_mul33x33 ( math_mul33x33 ( A , B ) , invFi_new ) ! current elastic deformation tensor
call constitutive_SandItsTangents ( Tstar , 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
Tstar_v = math_Mandel33to6 ( Tstar )
!* calculate plastic velocity gradient and its tangent from constitutive law
2014-08-26 20:14:32 +05:30
2018-09-20 09:54:03 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,i3,/)' ) '<< CRYST >> stress iteration ' , NiterationStressLp
write ( 6 , '(a,/,3(12x,3(e20.10,1x)/))' ) '<< CRYST >> Lpguess' , transpose ( Lpguess )
write ( 6 , '(a,/,3(12x,3(e20.10,1x)/))' ) '<< CRYST >> Fi' , transpose ( Fi_new )
write ( 6 , '(a,/,3(12x,3(e20.10,1x)/))' ) '<< CRYST >> Fe' , transpose ( Fe )
write ( 6 , '(a,/,6(e20.10,1x))' ) '<< CRYST >> Tstar' , Tstar_v
endif
2018-09-20 09:54:03 +05:30
#endif
2019-01-15 08:57:57 +05:30
call constitutive_LpAndItsTangents ( Lp_constitutive , dLp_dS , dLp_dFi , &
Tstar_v , Fi_new , ipc , ip , el )
2014-05-27 20:16:03 +05:30
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,/,3(12x,3(e20.10,1x)/))' ) '<< CRYST >> Lp_constitutive' , transpose ( Lp_constitutive )
endif
#endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!* update current residuum and check for convergence of loop
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
aTolLp = max ( rTol_crystalliteStress * max ( norm2 ( Lpguess ) , norm2 ( Lp_constitutive ) ) , & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress ) ! minimum lower cutoff
residuumLp = Lpguess - Lp_constitutive
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
if ( any ( IEEE_is_NaN ( residuumLp ) ) ) then ! NaN in residuum...
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt ) &
write ( 6 , '(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)' ) '<< CRYST >> integrateStress encountered NaN at el (elFE) ip ipc ' , &
el , '(' , mesh_element ( 1 , el ) , ')' , ip , ipc , &
' ; iteration ' , NiterationStressLp , &
' >> returning..!'
#endif
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_pInt &
. 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...
steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction
Lpguess = Lpguess_old + steplengthLp * deltaLp
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,1x,f7.4)' ) '<< CRYST >> linear search for Lpguess with step' , steplengthLp
2014-05-27 20:16:03 +05:30
endif
2019-01-15 08:57:57 +05:30
#endif
cycle LpLoop
2013-02-22 04:38:36 +05:30
endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!* calculate Jacobian for correction term
2014-05-27 20:16:03 +05:30
2019-01-15 08:57:57 +05:30
if ( mod ( jacoCounterLp , iJacoLpresiduum ) == 0_pInt ) then
dFe_dLp = 0.0_pReal
forall ( o = 1_pInt : 3_pInt , p = 1_pInt : 3_pInt ) &
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)
dFe_dLp = - dt * dFe_dLp
dRLp_dLp = math_identity2nd ( 9_pInt ) &
- math_Plain3333to99 ( math_mul3333xx3333 ( math_mul3333xx3333 ( dLp_dS , dS_dFe ) , dFe_dLp ) )
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,/,9(12x,9(e12.4,1x)/))' ) '<< CRYST >> dLp_dS' , math_Plain3333to99 ( dLp_dS )
write ( 6 , '(a,1x,e20.10)' ) '<< CRYST >> dLp_dS norm' , norm2 ( math_Plain3333to99 ( dLp_dS ) )
write ( 6 , '(a,/,9(12x,9(e12.4,1x)/))' ) '<< CRYST >> dRLp_dLp' , dRLp_dLp - math_identity2nd ( 9_pInt )
write ( 6 , '(a,1x,e20.10)' ) '<< CRYST >> dRLp_dLp norm' , norm2 ( dRLp_dLp - math_identity2nd ( 9_pInt ) )
endif
#endif
dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine
work = math_plain33to9 ( residuumLp )
call dgesv ( 9 , 1 , dRLp_dLp2 , 9 , ipiv , work , 9 , ierr ) ! solve dRLp/dLp * delta Lp = -res for delta Lp
if ( ierr / = 0_pInt ) then
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt ) then
write ( 6 , '(a,i8,1x,a,i8,a,1x,i2,1x,i3)' ) '<< CRYST >> integrateStress failed on dR/dLp inversion at el (elFE) ip ipc ' , &
el , '(' , mesh_element ( 1 , el ) , ')' , ip , ipc
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , * )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST >> dR_dLp' , transpose ( dRLp_dLp )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST >> dFe_dLp' , transpose ( math_Plain3333to99 ( dFe_dLp ) )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST >> dS_dFe_constitutive' , transpose ( math_Plain3333to99 ( dS_dFe ) )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST >> dLp_dS_constitutive' , transpose ( math_Plain3333to99 ( dLp_dS ) )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST >> A' , transpose ( A )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST >> B' , transpose ( B )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST >> Lp_constitutive' , transpose ( Lp_constitutive )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST >> Lpguess' , transpose ( Lpguess )
endif
endif
#endif
return
endif
deltaLp = - math_plain9to33 ( work )
endif
jacoCounterLp = jacoCounterLp + 1_pInt ! increase counter for jaco update
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
Lpguess = Lpguess + steplengthLp * deltaLp
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +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
2014-07-02 17:57:39 +05:30
2019-01-15 08:57:57 +05:30
call constitutive_LiAndItsTangents ( Li_constitutive , dLi_dS , dLi_dFi , &
Tstar_v , Fi_new , ipc , ip , el )
2014-09-03 01:41:57 +05:30
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST >> Li_constitutive' , transpose ( Li_constitutive )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST >> Liguess' , transpose ( Liguess )
endif
#endif
!* update current residuum and check for convergence of loop
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
aTolLi = max ( rTol_crystalliteStress * max ( norm2 ( Liguess ) , norm2 ( Li_constitutive ) ) , & ! absolute tolerance from largest acceptable relative error
aTol_crystalliteStress ) ! minimum lower cutoff
residuumLi = Liguess - Li_constitutive
if ( any ( IEEE_is_NaN ( residuumLi ) ) ) then ! NaN in residuum...
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_pInt &
. 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...
steplengthLi = subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction
Liguess = Liguess_old + steplengthLi * deltaLi
cycle LiLoop
endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!* calculate Jacobian for correction term
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
if ( mod ( jacoCounterLi , iJacoLpresiduum ) == 0_pInt ) then
temp_33 = math_mul33x33 ( math_mul33x33 ( A , B ) , invFi_current )
dFe_dLi = 0.0_pReal
dFi_dLi = 0.0_pReal
forall ( o = 1_pInt : 3_pInt , p = 1_pInt : 3_pInt )
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
end forall
forall ( o = 1_pInt : 3_pInt , p = 1_pInt : 3_pInt ) &
dFi_dLi ( 1 : 3 , 1 : 3 , o , p ) = math_mul33x33 ( math_mul33x33 ( Fi_new , dFi_dLi ( 1 : 3 , 1 : 3 , o , p ) ) , Fi_new )
dRLi_dLi = math_identity2nd ( 9_pInt ) &
- math_Plain3333to99 ( math_mul3333xx3333 ( dLi_dS , math_mul3333xx3333 ( dS_dFe , dFe_dLi ) + &
math_mul3333xx3333 ( dS_dFi , dFi_dLi ) ) ) &
- math_Plain3333to99 ( math_mul3333xx3333 ( dLi_dFi , dFi_dLi ) )
work = math_plain33to9 ( residuumLi )
call dgesv ( 9 , 1 , dRLi_dLi , 9 , ipiv , work , 9 , ierr ) ! solve dRLi/dLp * delta Li = -res for delta Li
if ( ierr / = 0_pInt ) then
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt ) then
write ( 6 , '(a,i8,1x,a,i8,a,1x,i2,1x,i3)' ) '<< CRYST >> integrateStress failed on dR/dLi inversion at el (elFE) ip ipc ' , &
el , '(' , mesh_element ( 1 , el ) , ')' , ip , ipc
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , * )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST >> dR_dLi' , transpose ( dRLi_dLi )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST >> dFe_dLi' , transpose ( math_Plain3333to99 ( dFe_dLi ) )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST >> dS_dFi_constitutive' , transpose ( math_Plain3333to99 ( dS_dFi ) )
write ( 6 , '(a,/,9(12x,9(e15.3,1x)/))' ) '<< CRYST >> dLi_dS_constitutive' , transpose ( math_Plain3333to99 ( dLi_dS ) )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST >> Li_constitutive' , transpose ( Li_constitutive )
write ( 6 , '(a,/,3(12x,3(e20.7,1x)/))' ) '<< CRYST >> Liguess' , transpose ( Liguess )
2013-02-22 04:38:36 +05:30
endif
endif
2019-01-15 08:57:57 +05:30
#endif
return
endif
deltaLi = - math_plain9to33 ( work )
endif
jacoCounterLi = jacoCounterLi + 1_pInt ! increase counter for jaco update
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
Liguess = Liguess + steplengthLi * deltaLi
enddo LiLoop
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!* calculate new plastic and elastic deformation gradient
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
invFp_new = math_mul33x33 ( invFp_current , B )
invFp_new = invFp_new / math_det33 ( invFp_new ) ** ( 1.0_pReal / 3.0_pReal ) ! regularize by det
Fp_new = math_inv33 ( invFp_new )
failedInversionInvFp : if ( all ( dEq0 ( Fp_new ) ) ) then
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0_pInt ) then
write ( 6 , '(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)' ) '<< CRYST >> integrateStress failed on invFp_new inversion at el (elFE) ip ipc ' , &
el , '(' , mesh_element ( 1 , el ) , ')' , ip , ipc , ' ; iteration ' , NiterationStressLp
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) &
write ( 6 , '(/,a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST >> invFp_new' , transpose ( invFp_new )
endif
#endif
return
endif failedInversionInvFp
Fe_new = math_mul33x33 ( math_mul33x33 ( Fg_new , invFp_new ) , invFi_new ) ! calc resulting Fe
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!* calculate 1st Piola-Kirchhoff stress
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
crystallite_P ( 1 : 3 , 1 : 3 , ipc , ip , el ) = math_mul33x33 ( math_mul33x33 ( Fg_new , invFp_new ) , &
math_mul33x33 ( math_Mandel6to33 ( Tstar_v ) , &
transpose ( invFp_new ) ) )
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!* store local values in global variables
crystallite_Lp ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Lpguess
crystallite_Li ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Liguess
crystallite_Tstar_v ( 1 : 6 , ipc , ip , el ) = Tstar_v
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-15 08:57:57 +05:30
!* set return flag to true
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
integrateStress = . true .
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST >> 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 >> Cauchy / MPa' , &
math_mul33x33 ( 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 >> Fe Lp Fe^-1' , &
transpose ( math_mul33x33 ( Fe_new , math_mul33x33 ( crystallite_Lp ( 1 : 3 , 1 : 3 , ipc , ip , el ) , math_inv33 ( Fe_new ) ) ) ) ! transpose to get correct print out order
write ( 6 , '(a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST >> Fp' , transpose ( crystallite_Fp ( 1 : 3 , 1 : 3 , ipc , ip , el ) )
write ( 6 , '(a,/,3(12x,3(f12.7,1x)/))' ) '<< CRYST >> Fi' , transpose ( crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) )
endif
2011-03-29 12:57:19 +05:30
#endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
end function integrateStress
2014-08-26 20:14:32 +05:30
2013-05-17 23:22:46 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
!> @brief calculates orientations
!--------------------------------------------------------------------------------------------------
subroutine crystallite_orientations
use math , only : &
math_rotationalPart33 , &
math_RtoQ
use FEsolving , only : &
FEsolving_execElem , &
FEsolving_execIP
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
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
implicit none
integer ( pInt ) &
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e , & !< counter in element loop
myPhase ! phase
2015-10-14 00:22:01 +05:30
2019-01-15 08:57:57 +05:30
! --- CALCULATE ORIENTATION AND LATTICE ROTATION ---
2014-07-02 17:57:39 +05:30
2019-01-15 08:57:57 +05:30
!$OMP PARALLEL DO
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
do c = 1_pInt , homogenization_Ngrains ( mesh_element ( 3 , e ) )
! somehow this subroutine is not threadsafe, so need critical statement here; not clear, what exactly the problem is
!$OMP CRITICAL (polarDecomp)
crystallite_orientation ( 1 : 4 , c , i , e ) = math_RtoQ ( transpose ( math_rotationalPart33 ( crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) ) ) )
!$OMP END CRITICAL (polarDecomp)
crystallite_rotation ( 1 : 4 , c , i , e ) = lattice_qDisorientation ( crystallite_orientation0 ( 1 : 4 , c , i , e ) , & ! active rotation from initial
crystallite_orientation ( 1 : 4 , c , i , e ) ) ! to current orientation (with no symmetry)
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL ---
! --- we use crystallite_orientation from above, so need a separate loop
nonlocalPresent : if ( any ( plasticState % nonLocal ) ) then
!$OMP PARALLEL DO PRIVATE(myPhase)
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 , e ) , FEsolving_execIP ( 2 , e )
myPhase = material_phase ( 1 , i , e ) ! get my phase (non-local models make no sense with more than one grain per material point)
if ( plasticState ( myPhase ) % nonLocal ) then ! if nonlocal model
! --- calculate compatibility and transmissivity between me and my neighbor ---
call plastic_nonlocal_updateCompatibility ( crystallite_orientation , i , e )
2014-07-02 17:57:39 +05:30
endif
2019-01-15 08:57:57 +05:30
enddo ; enddo
!$OMP END PARALLEL DO
endif nonlocalPresent
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
end subroutine crystallite_orientations
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief return results of particular grain
!--------------------------------------------------------------------------------------------------
function crystallite_postResults ( ipc , ip , el )
use math , only : &
math_qToEuler , &
math_qToEulerAxisAngle , &
math_mul33x33 , &
math_det33 , &
math_I3 , &
inDeg , &
math_Mandel6to33
use mesh , only : &
mesh_element , &
mesh_ipVolume , &
mesh_maxNipNeighbors , &
mesh_ipNeighborhood , &
FE_NipNeighbors , &
FE_geomtype , &
FE_celltype
use material , only : &
plasticState , &
sourceState , &
microstructure_crystallite , &
crystallite_Noutput , &
material_phase , &
material_texture , &
homogenization_Ngrains
use constitutive , only : &
constitutive_homogenizedC , &
constitutive_postResults
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
implicit none
integer ( pInt ) , intent ( in ) :: &
el , & !< element index
ip , & !< integration point index
ipc !< grain index
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
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
integer ( pInt ) :: &
o , &
c , &
crystID , &
mySize , &
n
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
crystID = microstructure_crystallite ( mesh_element ( 4 , el ) )
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
crystallite_postResults = 0.0_pReal
c = 0_pInt
crystallite_postResults ( c + 1 ) = real ( crystallite_sizePostResults ( crystID ) , pReal ) ! size of results from cryst
c = c + 1_pInt
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
do o = 1_pInt , crystallite_Noutput ( crystID )
mySize = 0_pInt
select case ( crystallite_outputID ( o , crystID ) )
case ( phase_ID )
mySize = 1_pInt
crystallite_postResults ( c + 1 ) = real ( material_phase ( ipc , ip , el ) , pReal ) ! phaseID of grain
case ( texture_ID )
mySize = 1_pInt
crystallite_postResults ( c + 1 ) = real ( material_texture ( ipc , ip , el ) , pReal ) ! textureID of grain
case ( volume_ID )
mySize = 1_pInt
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 )
mySize = 4_pInt
crystallite_postResults ( c + 1 : c + mySize ) = crystallite_orientation ( 1 : 4 , ipc , ip , el ) ! grain orientation as quaternion
case ( eulerangles_ID )
mySize = 3_pInt
crystallite_postResults ( c + 1 : c + mySize ) = inDeg &
* math_qToEuler ( crystallite_orientation ( 1 : 4 , ipc , ip , el ) ) ! grain orientation as Euler angles in degree
case ( grainrotation_ID )
mySize = 4_pInt
crystallite_postResults ( c + 1 : c + mySize ) = &
math_qToEulerAxisAngle ( crystallite_rotation ( 1 : 4 , ipc , ip , el ) ) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates
crystallite_postResults ( c + 4 ) = inDeg * crystallite_postResults ( c + 4 ) ! angle in degree
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! 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"
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
case ( defgrad_ID )
mySize = 9_pInt
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_partionedF ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( fe_ID )
mySize = 9_pInt
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_Fe ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( fp_ID )
mySize = 9_pInt
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_Fp ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( fi_ID )
mySize = 9_pInt
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( lp_ID )
mySize = 9_pInt
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_Lp ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( li_ID )
mySize = 9_pInt
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_Li ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( p_ID )
mySize = 9_pInt
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( transpose ( crystallite_P ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) , [ mySize ] )
case ( s_ID )
mySize = 9_pInt
crystallite_postResults ( c + 1 : c + mySize ) = &
reshape ( math_Mandel6to33 ( crystallite_Tstar_v ( 1 : 6 , ipc , ip , el ) ) , [ mySize ] )
case ( elasmatrix_ID )
mySize = 36_pInt
crystallite_postResults ( c + 1 : c + mySize ) = reshape ( constitutive_homogenizedC ( ipc , ip , el ) , [ mySize ] )
case ( neighboringelement_ID )
mySize = mesh_maxNipNeighbors
crystallite_postResults ( c + 1 : c + mySize ) = 0.0_pReal
forall ( n = 1_pInt : FE_NipNeighbors ( FE_celltype ( FE_geomtype ( mesh_element ( 2 , el ) ) ) ) ) &
crystallite_postResults ( c + n ) = real ( mesh_ipNeighborhood ( 1 , n , ip , el ) , pReal )
case ( neighboringip_ID )
mySize = mesh_maxNipNeighbors
crystallite_postResults ( c + 1 : c + mySize ) = 0.0_pReal
forall ( n = 1_pInt : FE_NipNeighbors ( FE_celltype ( FE_geomtype ( mesh_element ( 2 , el ) ) ) ) ) &
crystallite_postResults ( c + n ) = real ( mesh_ipNeighborhood ( 2 , n , ip , el ) , pReal )
end select
c = c + mySize
enddo
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
crystallite_postResults ( c + 1 ) = real ( plasticState ( material_phase ( ipc , ip , el ) ) % sizePostResults , pReal ) ! size of constitutive results
c = c + 1_pInt
if ( size ( crystallite_postResults ) - c > 0_pInt ) &
crystallite_postResults ( c + 1 : size ( crystallite_postResults ) ) = &
constitutive_postResults ( crystallite_Tstar_v ( 1 : 6 , ipc , ip , el ) , crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
crystallite_Fe , ipc , ip , el )
2015-10-14 00:22:01 +05:30
2019-01-15 08:57:57 +05:30
end function crystallite_postResults
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-01-15 08:57:57 +05:30
subroutine integrateStateFPI ( )
2017-05-04 04:02:44 +05:30
use , intrinsic :: &
IEEE_arithmetic
2018-02-16 20:06:18 +05:30
#ifdef DEBUG
2018-09-20 09:54:03 +05:30
use debug , only : &
2018-02-16 20:06:18 +05:30
debug_e , &
debug_i , &
debug_g , &
2019-01-15 08:57:57 +05:30
debug_level , &
2013-11-21 16:28:41 +05:30
debug_crystallite , &
debug_levelBasic , &
debug_levelExtensive , &
2018-05-17 20:03:35 +05:30
debug_levelSelective
2018-09-20 09:54:03 +05:30
#endif
2013-11-21 16:28:41 +05:30
use numerics , only : &
2019-01-15 08:57:57 +05:30
nState , &
2018-08-05 20:36:03 +05:30
rTol_crystalliteState
2014-08-26 20:14:32 +05:30
use FEsolving , only : &
FEsolving_execElem , &
2013-11-21 16:28:41 +05:30
FEsolving_execIP
use mesh , only : &
mesh_element , &
2019-01-15 08:57:57 +05:30
mesh_NcpElems
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 : &
2014-05-27 20:16:03 +05:30
constitutive_collectDotState , &
constitutive_microstructure , &
2015-05-28 22:32:23 +05:30
constitutive_plasticity_maxSizeDotState , &
constitutive_source_maxSizeDotState
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
implicit none
2019-01-15 08:57:57 +05:30
2013-11-21 16:28:41 +05:30
integer ( pInt ) :: &
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 , &
2015-05-28 22:32:23 +05:30
mySource , &
2014-09-03 01:16:52 +05:30
mySizePlasticDotState , & ! size of dot states
2015-05-28 22:32:23 +05:30
mySizeSourceDotState
2013-11-21 16:28:41 +05:30
integer ( pInt ) , dimension ( 2 ) :: &
eIter ! bounds for element iteration
integer ( pInt ) , dimension ( 2 , mesh_NcpElems ) :: &
iIter , & ! bounds for ip iteration
gIter ! bounds for grain iteration
2019-01-15 08:57:57 +05:30
real ( pReal ) :: &
dot_prod12 , &
dot_prod22 , &
plasticStateDamper , & ! damper for integration of state
sourceStateDamper
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState ) :: &
plasticStateResiduum , &
tempPlasticState
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , maxval ( phase_Nsources ) ) :: &
2015-05-28 22:32:23 +05:30
sourceStateResiduum , & ! residuum from evolution in micrstructure
2019-01-15 08:57:57 +05:30
tempSourceState
2013-11-21 16:28:41 +05:30
logical :: &
2015-05-28 22:32:23 +05:30
converged , &
2015-09-10 03:22:00 +05:30
NaN , &
2019-01-15 08:57:57 +05:30
singleRun , & ! flag indicating computation for single (g,i,e) triple
doneWithIntegration
2014-08-26 20:14:32 +05:30
2013-04-26 18:53:36 +05:30
eIter = FEsolving_execElem ( 1 : 2 )
do e = eIter ( 1 ) , eIter ( 2 )
iIter ( 1 : 2 , e ) = FEsolving_execIP ( 1 : 2 , e )
gIter ( 1 : 2 , e ) = [ 1_pInt , homogenization_Ngrains ( mesh_element ( 3 , e ) ) ]
enddo
2014-08-26 20:14:32 +05:30
2013-04-30 17:44:07 +05:30
singleRun = ( eIter ( 1 ) == eIter ( 2 ) . and . iIter ( 1 , eIter ( 1 ) ) == iIter ( 2 , eIter ( 2 ) ) )
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,i8,a)' ) '<< CRYST >> ' , count ( crystallite_todo ( : , : , : ) ) , ' grains todo at start of state integration'
#endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
! initialize dotState
if ( . not . singleRun ) then
forall ( p = 1_pInt : size ( plasticState ) )
plasticState ( p ) % previousDotState = 0.0_pReal
plasticState ( p ) % previousDotState2 = 0.0_pReal
end forall
do p = 1_pInt , size ( sourceState ) ; do mySource = 1_pInt , phase_Nsources ( p )
sourceState ( p ) % p ( mySource ) % previousDotState = 0.0_pReal
sourceState ( p ) % p ( mySource ) % previousDotState2 = 0.0_pReal
enddo ; enddo
else
e = eIter ( 1 )
i = iIter ( 1 , e )
do g = gIter ( 1 , e ) , gIter ( 2 , e )
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
plasticState ( p ) % previousDotState ( : , c ) = 0.0_pReal
plasticState ( p ) % previousDotState2 ( : , c ) = 0.0_pReal
do mySource = 1_pInt , phase_Nsources ( p )
sourceState ( p ) % p ( mySource ) % previousDotState ( : , c ) = 0.0_pReal
sourceState ( p ) % p ( mySource ) % previousDotState2 ( : , c ) = 0.0_pReal
enddo
enddo
endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --+>> PREGUESS FOR STATE <<+--
! --- DOT STATES ---
2014-05-27 20:16:03 +05:30
2014-08-26 20:14:32 +05:30
!$OMP PARALLEL
2019-01-15 08:57:57 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) &
call constitutive_collectDotState ( crystallite_Tstar_v ( 1 : 6 , g , i , e ) , &
crystallite_Fe , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp , &
crystallite_subdt ( g , i , e ) , crystallite_subFrac , g , i , e )
enddo ; enddo ; enddo
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP ENDDO
!$OMP DO PRIVATE(p,c,NaN)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % dotState ( : , c ) ) )
do mySource = 1_pInt , phase_Nsources ( p )
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % dotState ( : , c ) ) )
enddo
if ( NaN ) then ! NaN occured in any dotState
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , * ) '<< CRYST >> dotstate ' , plasticState ( p ) % dotState ( : , c )
#endif
if ( . not . crystallite_localPlasticity ( g , i , e ) ) then ! if broken is a non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo . and . crystallite_localPlasticity ! ...all non-locals done (and broken)
!$OMP END CRITICAL (checkTodo)
else ! broken one was local...
crystallite_todo ( g , i , e ) = . false . ! ... done (and broken)
2014-05-27 20:16:03 +05:30
endif
2013-02-22 04:38:36 +05:30
endif
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
!$OMP ENDDO
! --- UPDATE STATE ---
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,i8,a)' ) '<< CRYST >> ' , count ( crystallite_todo ( : , : , : ) ) , ' grains todo after preguess of state'
#endif
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
mySizePlasticDotState = plasticState ( p ) % sizeDotState
plasticState ( p ) % state ( 1 : mySizePlasticDotState , c ) = &
plasticState ( p ) % subState0 ( 1 : mySizePlasticDotState , c ) &
+ plasticState ( p ) % dotState ( 1 : mySizePlasticDotState , c ) &
* crystallite_subdt ( g , i , e )
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , c ) = &
sourceState ( p ) % p ( mySource ) % subState0 ( 1 : mySizeSourceDotState , c ) &
+ sourceState ( p ) % p ( mySource ) % dotState ( 1 : mySizeSourceDotState , c ) &
* crystallite_subdt ( g , i , e )
enddo
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP END PARALLEL
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --+>> STATE LOOP <<+--
2015-10-14 00:22:01 +05:30
2019-01-15 08:57:57 +05:30
NiterationState = 0_pInt
doneWithIntegration = . false .
crystalliteLooping : do while ( . not . doneWithIntegration . and . NiterationState < nState )
NiterationState = NiterationState + 1_pInt
!$OMP PARALLEL
! --- UPDATE DEPENDENT STATES ---
!$OMP DO PRIVATE(p,c)
2013-02-22 04:38:36 +05:30
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) &
call constitutive_microstructure ( crystallite_orientation , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp ( 1 : 3 , 1 : 3 , g , i , e ) , &
g , i , e ) ! update dependent state variables to be consistent with basic states
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
plasticState ( p ) % previousDotState2 ( : , c ) = plasticState ( p ) % previousDotState ( : , c )
plasticState ( p ) % previousDotState ( : , c ) = plasticState ( p ) % dotState ( : , c )
do mySource = 1_pInt , phase_Nsources ( p )
sourceState ( p ) % p ( mySource ) % previousDotState2 ( : , c ) = sourceState ( p ) % p ( mySource ) % previousDotState ( : , c )
sourceState ( p ) % p ( mySource ) % previousDotState ( : , c ) = sourceState ( p ) % p ( mySource ) % dotState ( : , c )
enddo
2013-02-22 04:38:36 +05:30
enddo ; enddo ; enddo
2014-08-26 20:14:32 +05:30
!$OMP ENDDO
2019-01-15 08:57:57 +05:30
! --- STRESS INTEGRATION ---
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,i8,a)' ) '<< CRYST >> ' , count ( crystallite_todo ( : , : , : ) ) , ' grains todo before stress integration'
#endif
2013-02-22 04:38:36 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
crystallite_todo ( g , i , e ) = integrateStress ( g , i , e )
2013-02-22 04:38:36 +05:30
!$OMP FLUSH(crystallite_todo)
2019-01-15 08:57:57 +05:30
if ( . not . crystallite_todo ( g , i , e ) . and . . not . crystallite_localPlasticity ( g , i , e ) ) then ! broken non-local...
2013-02-22 04:38:36 +05:30
!$OMP CRITICAL (checkTodo)
2019-01-15 08:57:57 +05:30
crystallite_todo = crystallite_todo . and . crystallite_localPlasticity ! ... then all non-locals skipped
2013-02-22 04:38:36 +05:30
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,i8,a)' ) '<< CRYST >> ' , count ( crystallite_todo ( : , : , : ) ) , ' grains todo after stress integration'
#endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- DOT STATE ---
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
!$OMP DO
2019-01-15 08:57:57 +05:30
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) &
2015-05-28 22:32:23 +05:30
call constitutive_collectDotState ( crystallite_Tstar_v ( 1 : 6 , g , i , e ) , &
2014-09-27 02:19:25 +05:30
crystallite_Fe , &
2018-08-25 19:29:34 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
2014-10-10 17:58:57 +05:30
crystallite_Fp , &
2013-11-21 16:28:41 +05:30
crystallite_subdt ( g , i , e ) , crystallite_subFrac , g , i , e )
2013-02-22 04:38:36 +05:30
enddo ; enddo ; enddo
!$OMP ENDDO
2019-01-15 08:57:57 +05:30
!$OMP DO PRIVATE(p,c)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
2013-02-22 04:38:36 +05:30
!$OMP FLUSH(crystallite_todo)
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
2016-01-15 05:49:44 +05:30
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
2017-05-04 04:02:44 +05:30
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % dotState ( : , c ) ) )
2015-05-28 22:32:23 +05:30
do mySource = 1_pInt , phase_Nsources ( p )
2017-05-04 04:02:44 +05:30
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % dotState ( : , c ) ) )
2015-10-14 00:22:01 +05:30
enddo
2019-01-15 08:57:57 +05:30
if ( NaN ) then ! NaN occured in any dotState
crystallite_todo ( g , i , e ) = . false . ! ... skip me next time
if ( . not . crystallite_localPlasticity ( g , i , e ) ) then ! if me is non-local...
2014-05-27 20:16:03 +05:30
!$OMP CRITICAL (checkTodo)
2019-01-15 08:57:57 +05:30
crystallite_todo = crystallite_todo . and . crystallite_localPlasticity ! ...all non-locals skipped
2014-05-27 20:16:03 +05:30
!$OMP END CRITICAL (checkTodo)
endif
endif
2019-01-15 08:57:57 +05:30
2013-02-22 04:38:36 +05:30
endif
2019-01-15 08:57:57 +05:30
2013-02-22 04:38:36 +05:30
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- UPDATE STATE ---
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO PRIVATE(dot_prod12,dot_prod22, &
!$OMP& mySizePlasticDotState,mySizeSourceDotState, &
!$OMP& plasticStateResiduum,sourceStateResiduum, &
!$OMP& plasticStatedamper,sourceStateDamper, &
!$OMP& tempPlasticState,tempSourceState,converged,p,c)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
2014-05-27 20:16:03 +05:30
2016-01-15 05:49:44 +05:30
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
2019-01-15 08:57:57 +05:30
dot_prod12 = dot_product ( plasticState ( p ) % dotState ( : , c ) &
- plasticState ( p ) % previousDotState ( : , c ) , &
plasticState ( p ) % previousDotState ( : , c ) &
- plasticState ( p ) % previousDotState2 ( : , c ) )
dot_prod22 = dot_product ( plasticState ( p ) % previousDotState ( : , c ) &
- plasticState ( p ) % previousDotState2 ( : , c ) , &
plasticState ( p ) % previousDotState ( : , c ) &
- plasticState ( p ) % previousDotState2 ( : , c ) )
if ( dot_prod22 > 0.0_pReal &
. and . ( dot_prod12 < 0.0_pReal &
. or . dot_product ( plasticState ( p ) % dotState ( : , c ) , &
plasticState ( p ) % previousDotState ( : , c ) ) < 0.0_pReal ) ) then
plasticStateDamper = 0.75_pReal + 0.25_pReal * tanh ( 2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22 )
else
plasticStateDamper = 1.0_pReal
endif
! --- get residui ---
2015-05-28 22:32:23 +05:30
mySizePlasticDotState = plasticState ( p ) % sizeDotState
2019-01-15 08:57:57 +05:30
plasticStateResiduum ( 1 : mySizePlasticDotState ) = &
plasticState ( p ) % state ( 1 : mySizePlasticDotState , c ) &
- plasticState ( p ) % subState0 ( 1 : mySizePlasticDotState , c ) &
- ( plasticState ( p ) % dotState ( 1 : mySizePlasticDotState , c ) * plasticStateDamper &
+ plasticState ( p ) % previousDotState ( 1 : mySizePlasticDotState , c ) &
* ( 1.0_pReal - plasticStateDamper ) ) * crystallite_subdt ( g , i , e )
! --- correct state with residuum ---
tempPlasticState ( 1 : mySizePlasticDotState ) = &
plasticState ( p ) % state ( 1 : mySizePlasticDotState , c ) &
- plasticStateResiduum ( 1 : mySizePlasticDotState ) ! need to copy to local variable, since we cant flush a pointer in openmp
! --- store corrected dotState --- (cannot do this before state update, because not sure how to flush pointers in openmp)
plasticState ( p ) % dotState ( : , c ) = plasticState ( p ) % dotState ( : , c ) * plasticStateDamper &
+ plasticState ( p ) % previousDotState ( : , c ) &
* ( 1.0_pReal - plasticStateDamper )
2015-10-14 00:22:01 +05:30
2015-05-28 22:32:23 +05:30
do mySource = 1_pInt , phase_Nsources ( p )
2019-01-15 08:57:57 +05:30
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
dot_prod12 = dot_product ( sourceState ( p ) % p ( mySource ) % dotState ( : , c ) &
- sourceState ( p ) % p ( mySource ) % previousDotState ( : , c ) , &
sourceState ( p ) % p ( mySource ) % previousDotState ( : , c ) &
- sourceState ( p ) % p ( mySource ) % previousDotState2 ( : , c ) )
dot_prod22 = dot_product ( sourceState ( p ) % p ( mySource ) % previousDotState ( : , c ) &
- sourceState ( p ) % p ( mySource ) % previousDotState2 ( : , c ) , &
sourceState ( p ) % p ( mySource ) % previousDotState ( : , c ) &
- sourceState ( p ) % p ( mySource ) % previousDotState2 ( : , c ) )
2015-10-14 00:22:01 +05:30
2019-01-15 08:57:57 +05:30
if ( dot_prod22 > 0.0_pReal &
. and . ( dot_prod12 < 0.0_pReal &
. or . dot_product ( sourceState ( p ) % p ( mySource ) % dotState ( : , c ) , &
sourceState ( p ) % p ( mySource ) % previousDotState ( : , c ) ) < 0.0_pReal ) ) then
sourceStateDamper = 0.75_pReal + 0.25_pReal * tanh ( 2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22 )
else
sourceStateDamper = 1.0_pReal
endif
! --- get residui ---
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceStateResiduum ( 1 : mySizeSourceDotState , mySource ) = &
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , c ) &
- sourceState ( p ) % p ( mySource ) % subState0 ( 1 : mySizeSourceDotState , c ) &
- ( sourceState ( p ) % p ( mySource ) % dotState ( 1 : mySizeSourceDotState , c ) * sourceStateDamper &
+ sourceState ( p ) % p ( mySource ) % previousDotState ( 1 : mySizeSourceDotState , c ) &
* ( 1.0_pReal - sourceStateDamper ) ) * crystallite_subdt ( g , i , e )
! --- correct state with residuum ---
tempSourceState ( 1 : mySizeSourceDotState , mySource ) = &
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , c ) &
- sourceStateResiduum ( 1 : mySizeSourceDotState , mySource ) ! need to copy to local variable, since we cant flush a pointer in openmp
! --- store corrected dotState --- (cannot do this before state update, because not sure how to flush pointers in openmp)
sourceState ( p ) % p ( mySource ) % dotState ( : , c ) = &
sourceState ( p ) % p ( mySource ) % dotState ( : , c ) * sourceStateDamper &
+ sourceState ( p ) % p ( mySource ) % previousDotState ( : , c ) &
* ( 1.0_pReal - sourceStateDamper )
enddo
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( e == debug_e . and . i == debug_i . and . g == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,i8,1x,i2,1x,i3,/)' ) '<< CRYST >> update state at el ip g ' , e , i , g
write ( 6 , '(a,f6.1,/)' ) '<< CRYST >> plasticstatedamper ' , plasticStatedamper
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> plastic state residuum' , &
abs ( plasticStateResiduum ( 1 : mySizePlasticDotState ) )
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> abstol dotstate' , plasticState ( p ) % aTolState ( 1 : mySizePlasticDotState )
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> reltol dotstate' , rTol_crystalliteState * &
abs ( tempPlasticState ( 1 : mySizePlasticDotState ) )
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> new state' , tempPlasticState ( 1 : mySizePlasticDotState )
2014-05-27 20:16:03 +05:30
endif
2012-03-09 01:55:28 +05:30
#endif
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
! --- converged ? ---
2019-01-15 08:57:57 +05:30
converged = all ( abs ( plasticStateResiduum ( 1 : mySizePlasticDotState ) ) < &
plasticState ( p ) % aTolState ( 1 : mySizePlasticDotState ) &
. or . abs ( plasticStateResiduum ( 1 : mySizePlasticDotState ) ) < &
rTol_crystalliteState * abs ( tempPlasticState ( 1 : mySizePlasticDotState ) ) )
2015-10-14 00:22:01 +05:30
do mySource = 1_pInt , phase_Nsources ( p )
2015-05-28 22:32:23 +05:30
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
converged = converged . and . &
2019-01-15 08:57:57 +05:30
all ( abs ( sourceStateResiduum ( 1 : mySizeSourceDotState , mySource ) ) < &
sourceState ( p ) % p ( mySource ) % aTolState ( 1 : mySizeSourceDotState ) &
. or . abs ( sourceStateResiduum ( 1 : mySizeSourceDotState , mySource ) ) < &
rTol_crystalliteState * abs ( tempSourceState ( 1 : mySizeSourceDotState , mySource ) ) )
2015-05-28 22:32:23 +05:30
enddo
2019-01-15 08:57:57 +05:30
if ( converged ) crystallite_converged ( g , i , e ) = . true . ! ... converged per definition
plasticState ( p ) % state ( 1 : mySizePlasticDotState , c ) = &
tempPlasticState ( 1 : mySizePlasticDotState )
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , c ) = &
tempSourceState ( 1 : mySizeSourceDotState , mySource )
enddo
endif
enddo ; enddo ; enddo
!$OMP ENDDO
! --- STATE JUMP ---
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$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
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,i8,a,i2)' ) '<< CRYST >> ' , count ( crystallite_converged ( : , : , : ) ) , &
' grains converged after state integration #' , NiterationState
#endif
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- NON-LOCAL CONVERGENCE CHECK ---
if ( . not . singleRun ) then ! if not requesting Integration of just a single IP
if ( any ( . not . crystallite_converged . and . . not . crystallite_localPlasticity ) ) & ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged . and . crystallite_localPlasticity ! ...restart all non-local as not converged
endif
2014-08-26 20:14:32 +05:30
2018-09-20 09:54:03 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) then
write ( 6 , '(a,i8,a)' ) '<< CRYST >> ' , count ( crystallite_converged ( : , : , : ) ) , &
' grains converged after non-local check'
write ( 6 , '(a,i8,a,i2,/)' ) '<< CRYST >> ' , count ( crystallite_todo ( : , : , : ) ) , &
' grains todo after state integration #' , NiterationState
endif
2018-09-20 09:54:03 +05:30
#endif
2014-05-27 20:16:03 +05:30
2019-01-15 08:57:57 +05:30
! --- CHECK IF DONE WITH INTEGRATION ---
2014-05-27 20:16:03 +05:30
2019-01-15 08:57:57 +05:30
doneWithIntegration = . true .
elemLoop : do e = eIter ( 1 ) , eIter ( 2 )
do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
doneWithIntegration = . false .
exit elemLoop
endif
enddo ; enddo
enddo elemLoop
enddo crystalliteLooping
end subroutine integrateStateFPI
2010-10-01 17:48:49 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-10-16 18:34:59 +05:30
!> @brief integrate stress, and state with 1st order explicit Euler method
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2018-09-20 09:57:53 +05:30
subroutine integrateStateEuler ( )
2017-05-04 04:02:44 +05:30
use , intrinsic :: &
IEEE_arithmetic
2018-02-16 20:06:18 +05:30
#ifdef DEBUG
2018-09-20 09:54:03 +05:30
use debug , only : &
2018-02-16 20:06:18 +05:30
debug_e , &
debug_i , &
debug_g , &
2013-11-21 16:28:41 +05:30
debug_level , &
debug_crystallite , &
debug_levelBasic , &
debug_levelExtensive , &
2018-05-17 20:03:35 +05:30
debug_levelSelective
2018-09-20 09:54:03 +05:30
#endif
2014-08-26 20:14:32 +05:30
use FEsolving , only : &
FEsolving_execElem , &
2013-11-21 16:28:41 +05:30
FEsolving_execIP
use mesh , only : &
mesh_element , &
mesh_NcpElems
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 , &
2013-11-21 16:28:41 +05:30
homogenization_Ngrains
use constitutive , only : &
2014-05-27 20:16:03 +05:30
constitutive_collectDotState , &
2014-06-23 00:28:29 +05:30
constitutive_microstructure
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
implicit none
2013-11-21 16:28:41 +05:30
integer ( pInt ) :: &
e , & ! element index in element loop
i , & ! integration point index in ip loop
g , & ! grain index in grain loop
2015-05-28 22:32:23 +05:30
p , & ! phase loop
2014-06-25 04:51:25 +05:30
c , &
2015-05-28 22:32:23 +05:30
mySource , &
2014-06-25 04:51:25 +05:30
mySizePlasticDotState , &
2015-05-28 22:32:23 +05:30
mySizeSourceDotState
2013-11-21 16:28:41 +05:30
integer ( pInt ) , dimension ( 2 ) :: &
eIter ! bounds for element iteration
integer ( pInt ) , dimension ( 2 , mesh_NcpElems ) :: &
iIter , & ! bounds for ip iteration
gIter ! bounds for grain iteration
logical :: &
2015-09-10 03:22:00 +05:30
NaN , &
2013-11-21 16:28:41 +05:30
singleRun ! flag indicating computation for single (g,i,e) triple
2013-04-26 18:53:36 +05:30
eIter = FEsolving_execElem ( 1 : 2 )
do e = eIter ( 1 ) , eIter ( 2 )
iIter ( 1 : 2 , e ) = FEsolving_execIP ( 1 : 2 , e )
gIter ( 1 : 2 , e ) = [ 1_pInt , homogenization_Ngrains ( mesh_element ( 3 , e ) ) ]
enddo
2014-08-26 20:14:32 +05:30
2013-04-30 17:44:07 +05:30
singleRun = ( eIter ( 1 ) == eIter ( 2 ) . and . iIter ( 1 , eIter ( 1 ) ) == iIter ( 2 , eIter ( 2 ) ) )
2014-08-26 20:14:32 +05:30
2014-05-27 20:16:03 +05:30
!$OMP PARALLEL
2014-08-26 20:14:32 +05:30
2013-10-16 18:34:59 +05:30
! --- DOT STATE ---
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
2013-10-19 00:27:28 +05:30
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) &
2015-05-28 22:32:23 +05:30
call constitutive_collectDotState ( crystallite_Tstar_v ( 1 : 6 , g , i , e ) , &
2014-09-27 02:19:25 +05:30
crystallite_Fe , &
2018-08-25 19:29:34 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
2014-10-10 17:58:57 +05:30
crystallite_Fp , &
2013-11-21 16:28:41 +05:30
crystallite_subdt ( g , i , e ) , crystallite_subFrac , g , i , e )
2013-02-22 04:38:36 +05:30
enddo ; enddo ; enddo
!$OMP ENDDO
2015-09-10 03:22:00 +05:30
!$OMP DO PRIVATE(p,c,NaN)
2013-02-22 04:38:36 +05:30
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
2016-01-15 05:49:44 +05:30
c = phasememberAt ( g , i , e )
p = phaseAt ( g , i , e )
2017-05-04 04:02:44 +05:30
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % dotState ( : , c ) ) )
2015-05-28 22:32:23 +05:30
do mySource = 1_pInt , phase_Nsources ( p )
2017-05-04 04:02:44 +05:30
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % dotState ( : , c ) ) )
2015-10-14 00:22:01 +05:30
enddo
2015-09-10 03:22:00 +05:30
if ( NaN ) then ! NaN occured in any dotState
2019-01-14 17:26:46 +05:30
if ( . not . crystallite_localPlasticity ( g , i , e ) ) then ! if broken non-local...
2014-05-27 20:16:03 +05:30
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo . and . crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
else ! if broken local...
crystallite_todo ( g , i , e ) = . false . ! ... skip this one next time
endif
endif
2013-02-22 04:38:36 +05:30
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2013-10-16 18:34:59 +05:30
! --- UPDATE STATE ---
2015-10-14 00:22:01 +05:30
2015-05-28 22:32:23 +05:30
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c)
2013-02-22 04:38:36 +05:30
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
2016-01-15 05:49:44 +05:30
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
2014-06-25 04:51:25 +05:30
mySizePlasticDotState = plasticState ( p ) % sizeDotState
2015-05-28 22:32:23 +05:30
plasticState ( p ) % state ( 1 : mySizePlasticDotState , c ) = &
plasticState ( p ) % state ( 1 : mySizePlasticDotState , c ) &
+ plasticState ( p ) % dotState ( 1 : mySizePlasticDotState , c ) &
2015-10-14 00:22:01 +05:30
* crystallite_subdt ( g , i , e )
2015-05-28 22:32:23 +05:30
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , c ) = &
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , c ) &
+ sourceState ( p ) % p ( mySource ) % dotState ( 1 : mySizeSourceDotState , c ) &
* crystallite_subdt ( g , i , e )
enddo
2014-08-26 20:14:32 +05:30
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2014-07-02 17:57:39 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( e == debug_e . and . i == debug_i . and . g == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
2016-01-15 05:49:44 +05:30
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
2014-07-02 17:57:39 +05:30
write ( 6 , '(a,i8,1x,i2,1x,i3,/)' ) '<< CRYST >> update state at el ip g ' , e , i , g
2015-05-28 22:32:23 +05:30
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> dotState' , plasticState ( p ) % dotState ( 1 : mySizePlasticDotState , c )
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> new state' , plasticState ( p ) % state ( 1 : mySizePlasticDotState , c )
2014-07-02 17:57:39 +05:30
endif
#endif
2013-02-22 04:38:36 +05:30
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
! --- STATE JUMP ---
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
!$OMP DO
2015-09-10 03:22:00 +05:30
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
2013-02-22 04:38:36 +05:30
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
2018-09-20 09:57:53 +05:30
crystallite_todo ( g , i , e ) = stateJump ( g , i , e )
2013-02-22 04:38:36 +05:30
!$OMP FLUSH(crystallite_todo)
2019-01-14 17:26:46 +05:30
if ( . not . crystallite_todo ( g , i , e ) . and . . not . crystallite_localPlasticity ( g , i , e ) ) then ! if broken non-local...
2013-02-22 04:38:36 +05:30
!$OMP CRITICAL (checkTodo)
2015-09-10 03:22:00 +05:30
crystallite_todo = crystallite_todo . and . crystallite_localPlasticity ! ...all non-locals skipped
2013-02-22 04:38:36 +05:30
!$OMP END CRITICAL (checkTodo)
endif
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
! --- UPDATE DEPENDENT STATES ---
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
2013-10-19 00:27:28 +05:30
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) &
2015-10-14 00:22:01 +05:30
call constitutive_microstructure ( crystallite_orientation , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
2014-11-29 01:36:24 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , g , i , e ) , &
2017-02-04 00:49:02 +05:30
g , i , e ) ! update dependent state variables to be consistent with basic states
2014-06-25 04:51:25 +05:30
enddo ; enddo ; enddo
2013-02-22 04:38:36 +05:30
!$OMP ENDDO
2014-05-27 20:16:03 +05:30
!$OMP END PARALLEL
2014-08-26 20:14:32 +05:30
2014-05-27 20:16:03 +05:30
!$OMP PARALLEL
2013-02-22 04:38:36 +05:30
! --- STRESS INTEGRATION ---
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
2018-09-20 09:57:53 +05:30
crystallite_todo ( g , i , e ) = integrateStress ( g , i , e )
2013-02-22 04:38:36 +05:30
!$OMP FLUSH(crystallite_todo)
2019-01-14 17:26:46 +05:30
if ( . not . crystallite_todo ( g , i , e ) . and . . not . crystallite_localPlasticity ( g , i , e ) ) then ! if broken non-local...
2013-02-22 04:38:36 +05:30
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo . and . crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
2011-11-04 18:14:50 +05:30
enddo ; enddo ; enddo
2013-02-22 04:38:36 +05:30
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
! --- SET CONVERGENCE FLAG ---
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
2018-07-04 01:42:25 +05:30
crystallite_converged ( g , i , e ) = crystallite_todo ( g , i , e ) . or . crystallite_converged ( g , i , e ) ! if still "to do" then converged per definitionem
2013-02-22 04:38:36 +05:30
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
!$OMP END PARALLEL
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
! --- CHECK NON-LOCAL CONVERGENCE ---
2014-08-26 20:14:32 +05:30
if ( . not . singleRun ) then ! if not requesting Integration of just a single IP
2019-01-14 17:26:46 +05:30
if ( any ( . not . crystallite_converged . and . . not . crystallite_localPlasticity ) ) & ! any non-local not yet converged (or broken)...
2013-02-22 04:38:36 +05:30
crystallite_converged = crystallite_converged . and . crystallite_localPlasticity ! ...restart all non-local as not converged
endif
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-01-15 08:57:57 +05:30
subroutine integrateStateAdaptiveEuler ( )
2017-05-04 04:02:44 +05:30
use , intrinsic :: &
IEEE_arithmetic
2018-02-16 20:06:18 +05:30
#ifdef DEBUG
2018-09-20 09:54:03 +05:30
use debug , only : &
2015-04-21 20:03:38 +05:30
debug_e , &
debug_i , &
debug_g , &
2019-01-15 08:57:57 +05:30
debug_level , &
2013-11-21 16:28:41 +05:30
debug_crystallite , &
debug_levelBasic , &
debug_levelExtensive , &
2018-05-17 20:03:35 +05:30
debug_levelSelective
2018-09-20 09:54:03 +05:30
#endif
2013-11-21 16:28:41 +05:30
use numerics , only : &
rTol_crystalliteState
use FEsolving , only : &
2014-08-26 20:14:32 +05:30
FEsolving_execElem , &
2013-11-21 16:28:41 +05:30
FEsolving_execIP
use mesh , only : &
mesh_element , &
2019-01-15 08:57:57 +05:30
mesh_NcpElems , &
mesh_maxNips
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 : &
2014-05-27 20:16:03 +05:30
constitutive_collectDotState , &
constitutive_microstructure , &
2015-05-28 22:32:23 +05:30
constitutive_plasticity_maxSizeDotState , &
constitutive_source_maxSizeDotState
2014-07-02 17:57:39 +05:30
2014-08-26 20:14:32 +05:30
implicit none
2013-05-17 23:22:46 +05:30
integer ( pInt ) :: &
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
s , & ! state index
2014-06-25 04:51:25 +05:30
p , &
c , &
2015-05-28 22:32:23 +05:30
mySource , &
2014-09-03 01:16:52 +05:30
mySizePlasticDotState , & ! size of dot states
2015-05-28 22:32:23 +05:30
mySizeSourceDotState
2019-01-15 08:57:57 +05:30
integer ( pInt ) , dimension ( 2 ) :: &
eIter ! bounds for element iteration
integer ( pInt ) , dimension ( 2 , mesh_NcpElems ) :: &
iIter , & ! bounds for ip iteration
gIter ! bounds for grain iteration
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState , &
homogenization_maxNgrains , mesh_maxNips , mesh_NcpElems ) :: &
plasticStateResiduum , & ! residuum from evolution in micrstructure
relPlasticStateResiduum ! relative residuum from evolution in microstructure
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , &
maxval ( phase_Nsources ) , &
homogenization_maxNgrains , mesh_maxNips , mesh_NcpElems ) :: &
sourceStateResiduum , & ! residuum from evolution in micrstructure
relSourceStateResiduum ! relative residuum from evolution in microstructure
2014-07-09 19:04:33 +05:30
2019-01-15 08:57:57 +05:30
logical :: &
converged , &
NaN , &
singleRun ! flag indicating computation for single (g,i,e) triple
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP ---
eIter = FEsolving_execElem ( 1 : 2 )
do e = eIter ( 1 ) , eIter ( 2 )
iIter ( 1 : 2 , e ) = FEsolving_execIP ( 1 : 2 , e )
gIter ( 1 : 2 , e ) = [ 1_pInt , homogenization_Ngrains ( mesh_element ( 3 , e ) ) ]
enddo
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
singleRun = ( eIter ( 1 ) == eIter ( 2 ) . and . iIter ( 1 , eIter ( 1 ) ) == iIter ( 2 , eIter ( 2 ) ) )
2014-08-26 20:14:32 +05:30
2014-07-09 19:04:33 +05:30
2019-01-15 08:57:57 +05:30
plasticStateResiduum = 0.0_pReal
relPlasticStateResiduum = 0.0_pReal
sourceStateResiduum = 0.0_pReal
relSourceStateResiduum = 0.0_pReal
2014-05-27 20:16:03 +05:30
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP PARALLEL
! --- DOT STATE (EULER INTEGRATION) ---
2014-08-26 20:14:32 +05:30
2013-02-22 04:38:36 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) &
2015-05-28 22:32:23 +05:30
call constitutive_collectDotState ( crystallite_Tstar_v ( 1 : 6 , g , i , e ) , &
2014-09-27 02:19:25 +05:30
crystallite_Fe , &
2018-08-25 19:29:34 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
2014-10-10 17:58:57 +05:30
crystallite_Fp , &
2013-11-21 16:28:41 +05:30
crystallite_subdt ( g , i , e ) , crystallite_subFrac , g , i , e )
2019-01-15 08:57:57 +05:30
enddo ; enddo ; enddo
2013-02-22 04:38:36 +05:30
!$OMP ENDDO
2019-01-15 08:57:57 +05:30
!$OMP DO PRIVATE(p,c,NaN)
2013-02-22 04:38:36 +05:30
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2016-01-15 05:49:44 +05:30
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
2017-05-04 04:02:44 +05:30
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % dotState ( : , c ) ) )
2015-05-28 22:32:23 +05:30
do mySource = 1_pInt , phase_Nsources ( p )
2017-05-04 04:02:44 +05:30
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % dotState ( : , c ) ) )
2015-10-14 00:22:01 +05:30
enddo
2015-09-10 03:22:00 +05:30
if ( NaN ) then ! NaN occured in any dotState
2019-01-15 08:57:57 +05:30
if ( . not . crystallite_localPlasticity ( g , i , e ) ) then ! if broken non-local...
2014-05-27 20:16:03 +05:30
!$OMP CRITICAL (checkTodo)
2014-09-03 01:16:52 +05:30
crystallite_todo = crystallite_todo . and . crystallite_localPlasticity ! ...all non-locals skipped
2014-05-27 20:16:03 +05:30
!$OMP END CRITICAL (checkTodo)
2019-01-15 08:57:57 +05:30
else ! if broken local...
crystallite_todo ( g , i , e ) = . false . ! ... skip this one next time
2014-08-26 20:14:32 +05:30
endif
endif
2013-02-22 04:38:36 +05:30
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-05-27 20:16:03 +05:30
2014-07-09 19:04:33 +05:30
2019-01-15 08:57:57 +05:30
! --- STATE UPDATE (EULER INTEGRATION) ---
2013-10-16 18:34:59 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
2016-01-15 05:49:44 +05:30
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
2019-01-15 08:57:57 +05:30
mySizePlasticDotState = plasticState ( p ) % sizeDotState
plasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) = &
- 0.5_pReal &
* plasticState ( p ) % dotstate ( 1 : mySizePlasticDotState , c ) &
* crystallite_subdt ( g , i , e ) ! contribution to absolute residuum in state
plasticState ( p ) % state ( 1 : mySizePlasticDotState , c ) = &
plasticState ( p ) % state ( 1 : mySizePlasticDotState , c ) &
+ plasticState ( p ) % dotstate ( 1 : mySizePlasticDotState , c ) &
* crystallite_subdt ( g , i , e )
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceStateResiduum ( 1 : mySizeSourceDotState , mySource , g , i , e ) = &
- 0.5_pReal &
* sourceState ( p ) % p ( mySource ) % dotstate ( 1 : mySizeSourceDotState , c ) &
* crystallite_subdt ( g , i , e ) ! contribution to absolute residuum in state
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , c ) = &
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , c ) &
+ sourceState ( p ) % p ( mySource ) % dotstate ( 1 : mySizeSourceDotState , c ) &
* crystallite_subdt ( g , i , e )
enddo
endif
enddo ; enddo ; enddo
!$OMP ENDDO
! --- STATE JUMP ---
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
crystallite_todo ( g , i , e ) = stateJump ( g , i , e )
!$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)
2014-07-07 19:51:58 +05:30
endif
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-09-03 01:16:52 +05:30
2019-01-15 08:57:57 +05:30
! --- UPDATE DEPENDENT STATES (EULER INTEGRATION) ---
2015-10-14 00:22:01 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) &
call constitutive_microstructure ( crystallite_orientation , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp ( 1 : 3 , 1 : 3 , g , i , e ) , &
g , i , e ) ! update dependent state variables to be consistent with basic states
enddo ; enddo ; enddo
!$OMP ENDDO
!$OMP END PARALLEL
2014-12-08 16:15:12 +05:30
2015-10-14 00:22:01 +05:30
2019-01-15 08:57:57 +05:30
! --- STRESS INTEGRATION (EULER INTEGRATION) ---
2015-05-28 22:32:23 +05:30
2019-01-15 08:57:57 +05:30
!$OMP PARALLEL DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
crystallite_todo ( g , i , e ) = integrateStress ( g , i , e )
!$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
!$OMP PARALLEL
! --- DOT STATE (HEUN METHOD) ---
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) &
call constitutive_collectDotState ( crystallite_Tstar_v ( 1 : 6 , g , i , e ) , &
crystallite_Fe , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp , &
crystallite_subdt ( g , i , e ) , crystallite_subFrac , g , i , e )
enddo ; enddo ; enddo
!$OMP ENDDO
!$OMP DO PRIVATE(p,c,NaN)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % dotState ( : , c ) ) )
do mySource = 1_pInt , phase_Nsources ( p )
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % dotState ( : , c ) ) )
enddo
if ( NaN ) then ! NaN occured in any dotState
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)
else ! if broken local...
crystallite_todo ( g , i , e ) = . false . ! ... skip this one next time
endif
endif
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2015-05-28 22:32:23 +05:30
2019-01-15 08:57:57 +05:30
! --- ERROR ESTIMATE FOR STATE (HEUN METHOD) ---
!$OMP SINGLE
relPlasticStateResiduum = 0.0_pReal
relSourceStateResiduum = 0.0_pReal
!$OMP END SINGLE
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,converged,p,c,s)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
! --- contribution of heun step to absolute residui ---
mySizePlasticDotState = plasticState ( p ) % sizeDotState
plasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) = &
plasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) &
+ 0.5_pReal * plasticState ( p ) % dotState ( : , c ) &
* crystallite_subdt ( g , i , e ) ! contribution to absolute residuum in state
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceStateResiduum ( 1 : mySizeSourceDotState , mySource , g , i , e ) = &
sourceStateResiduum ( 1 : mySizeSourceDotState , mySource , g , i , e ) &
+ 0.5_pReal * sourceState ( p ) % p ( mySource ) % dotState ( : , c ) &
* crystallite_subdt ( g , i , e ) ! contribution to absolute residuum in state
enddo
! --- relative residui ---
forall ( s = 1_pInt : mySizePlasticDotState , abs ( plasticState ( p ) % dotState ( s , c ) ) > 0.0_pReal ) &
relPlasticStateResiduum ( s , g , i , e ) = &
plasticStateResiduum ( s , g , i , e ) / plasticState ( p ) % dotState ( s , c )
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
forall ( s = 1_pInt : mySizeSourceDotState , abs ( sourceState ( p ) % p ( mySource ) % dotState ( s , c ) ) > 0.0_pReal ) &
relSourceStateResiduum ( s , mySource , g , i , e ) = &
sourceStateResiduum ( s , mySource , g , i , e ) / sourceState ( p ) % p ( mySource ) % dotState ( s , c )
2015-10-14 00:22:01 +05:30
enddo
2015-05-28 22:32:23 +05:30
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2014-05-15 18:31:54 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
2019-01-15 08:57:57 +05:30
. and . ( ( e == debug_e . and . i == debug_i . and . g == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,i8,1x,i2,1x,i3,/)' ) '<< CRYST >> updateState at el ip g ' , e , i , g
write ( 6 , '(a,/,(12x,12(f12.1,1x)),/)' ) '<< CRYST >> absolute residuum tolerance' , &
plasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) / plasticState ( p ) % aTolState ( 1 : mySizePlasticDotState )
write ( 6 , '(a,/,(12x,12(f12.1,1x)),/)' ) '<< CRYST >> relative residuum tolerance' , &
relPlasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) / rTol_crystalliteState
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> dotState' , plasticState ( p ) % dotState ( 1 : mySizePlasticDotState , c ) &
- 2.0_pReal * plasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) / crystallite_subdt ( g , i , e ) ! calculate former dotstate from higher order solution and state residuum
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> new state' , plasticState ( p ) % state ( 1 : mySizePlasticDotState , c )
2013-02-22 04:38:36 +05:30
endif
2012-05-17 20:55:21 +05:30
#endif
2014-05-27 20:16:03 +05:30
! --- converged ? ---
2019-01-15 08:57:57 +05:30
converged = all ( abs ( relPlasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) ) < &
rTol_crystalliteState . or . &
abs ( plasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) ) < &
plasticState ( p ) % aTolState ( 1 : mySizePlasticDotState ) )
2015-10-14 00:22:01 +05:30
do mySource = 1_pInt , phase_Nsources ( p )
2015-05-28 22:32:23 +05:30
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
converged = converged . and . &
2019-01-15 08:57:57 +05:30
all ( abs ( relSourceStateResiduum ( 1 : mySizeSourceDotState , mySource , g , i , e ) ) < &
rTol_crystalliteState . or . &
abs ( sourceStateResiduum ( 1 : mySizeSourceDotState , mySource , g , i , e ) ) < &
sourceState ( p ) % p ( mySource ) % aTolState ( 1 : mySizeSourceDotState ) )
2015-10-14 00:22:01 +05:30
enddo
2019-01-15 08:57:57 +05:30
if ( converged ) crystallite_converged ( g , i , e ) = . true . ! ... converged per definitionem
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-15 08:57:57 +05:30
! --- NONLOCAL CONVERGENCE CHECK ---
2018-09-20 09:54:03 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,i8,a,i2,/)' ) '<< CRYST >> ' , count ( crystallite_converged ( : , : , : ) ) , ' grains converged'
2018-09-20 09:54:03 +05:30
#endif
2019-01-15 08:57:57 +05:30
if ( ( . not . singleRun ) . and . any ( . not . crystallite_converged . and . . not . crystallite_localPlasticity ) ) & ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged . and . crystallite_localPlasticity ! ...restart all non-local as not converged
2018-09-20 09:54:03 +05:30
2014-01-22 00:15:41 +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
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
subroutine integrateStateRK4 ( )
2017-05-04 04:02:44 +05:30
use , intrinsic :: &
IEEE_arithmetic
2018-02-16 20:06:18 +05:30
#ifdef DEBUG
2013-04-29 16:47:30 +05:30
use debug , only : &
2018-02-16 20:06:18 +05:30
debug_e , &
debug_i , &
debug_g , &
2015-04-21 20:03:38 +05:30
debug_level , &
2013-04-29 16:47:30 +05:30
debug_crystallite , &
2019-01-15 08:57:57 +05:30
debug_levelBasic , &
2013-04-29 16:47:30 +05:30
debug_levelExtensive , &
2018-02-16 20:06:18 +05:30
debug_levelSelective
#endif
2019-01-15 08:57:57 +05:30
use FEsolving , only : &
FEsolving_execElem , &
FEsolving_execIP
use mesh , only : &
mesh_element , &
mesh_NcpElems
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
use config , only : &
material_Nphase
use constitutive , only : &
constitutive_collectDotState , &
constitutive_microstructure
2014-09-10 14:07:12 +05:30
implicit none
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-01-15 08:57:57 +05:30
integer ( pInt ) :: e , & ! element index in element loop
i , & ! integration point index in ip loop
g , & ! grain index in grain loop
p , & ! phase loop
c , &
n , &
mySource , &
mySizePlasticDotState , &
mySizeSourceDotState
integer ( pInt ) , dimension ( 2 ) :: eIter ! bounds for element iteration
integer ( pInt ) , dimension ( 2 , mesh_NcpElems ) :: iIter , & ! bounds for ip iteration
gIter ! bounds for grain iteration
logical :: NaN , &
singleRun ! flag indicating computation for single (g,i,e) triple
2014-09-10 14:07:12 +05:30
2019-01-15 08:57:57 +05:30
eIter = FEsolving_execElem ( 1 : 2 )
do e = eIter ( 1 ) , eIter ( 2 )
iIter ( 1 : 2 , e ) = FEsolving_execIP ( 1 : 2 , e )
gIter ( 1 : 2 , e ) = [ 1_pInt , homogenization_Ngrains ( mesh_element ( 3 , e ) ) ]
enddo
2014-11-06 17:19:37 +05:30
2019-01-15 08:57:57 +05:30
singleRun = ( eIter ( 1 ) == eIter ( 2 ) . and . iIter ( 1 , eIter ( 1 ) ) == iIter ( 2 , eIter ( 2 ) ) )
2014-11-06 17:19:37 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
! initialize dotState
if ( . not . singleRun ) then
do p = 1_pInt , material_Nphase
plasticState ( p ) % RK4dotState = 0.0_pReal
do mySource = 1_pInt , phase_Nsources ( p )
sourceState ( p ) % p ( mySource ) % RK4dotState = 0.0_pReal
enddo
enddo
else
e = eIter ( 1 )
i = iIter ( 1 , e )
do g = gIter ( 1 , e ) , gIter ( 2 , e )
plasticState ( phaseAt ( g , i , e ) ) % RK4dotState ( : , phasememberAt ( g , i , e ) ) = 0.0_pReal
do mySource = 1_pInt , phase_Nsources ( phaseAt ( g , i , e ) )
sourceState ( phaseAt ( g , i , e ) ) % p ( mySource ) % RK4dotState ( : , phasememberAt ( g , i , e ) ) = 0.0_pReal
enddo
enddo
endif
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
! first Runge-Kutta step
!$OMP PARALLEL
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) &
call constitutive_collectDotState ( crystallite_Tstar_v ( 1 : 6 , g , i , e ) , &
crystallite_Fe , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp , &
crystallite_subdt ( g , i , e ) , crystallite_subFrac , g , i , e )
enddo ; enddo ; enddo
!$OMP ENDDO
2018-05-09 20:24:06 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO PRIVATE(p,c,NaN)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
c = phasememberAt ( g , i , e )
p = phaseAt ( g , i , e )
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % dotState ( : , c ) ) )
do mySource = 1_pInt , phase_Nsources ( p )
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % dotState ( : , c ) ) )
enddo
if ( NaN ) then ! NaN occured in any dotState
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)
else ! if broken local...
crystallite_todo ( g , i , e ) = . false . ! ... skip this one next time
endif
endif
endif
enddo ; enddo ; enddo
!$OMP ENDDO
!$OMP END PARALLEL
!--------------------------------------------------------------------------------------------------
! --- SECOND TO FOURTH RUNGE KUTTA STEP PLUS FINAL INTEGRATION ---
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
do n = 1_pInt , 4_pInt
! --- state update ---
2013-04-29 16:47:30 +05:30
2019-01-15 08:57:57 +05:30
!$OMP PARALLEL
!$OMP DO PRIVATE(p,c)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
plasticState ( p ) % RK4dotState ( : , c ) = plasticState ( p ) % RK4dotState ( : , c ) &
+ weight ( n ) * plasticState ( p ) % dotState ( : , c )
do mySource = 1_pInt , phase_Nsources ( p )
sourceState ( p ) % p ( mySource ) % RK4dotState ( : , c ) = sourceState ( p ) % p ( mySource ) % RK4dotState ( : , c ) &
+ weight ( n ) * sourceState ( p ) % p ( mySource ) % dotState ( : , c )
enddo
endif
enddo ; enddo ; enddo
!$OMP ENDDO
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
mySizePlasticDotState = plasticState ( p ) % sizeDotState
plasticState ( p ) % state ( 1 : mySizePlasticDotState , c ) = &
plasticState ( p ) % subState0 ( 1 : mySizePlasticDotState , c ) &
+ plasticState ( p ) % dotState ( 1 : mySizePlasticDotState , c ) &
* crystallite_subdt ( g , i , e ) * timeStepFraction ( n )
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , c ) = &
sourceState ( p ) % p ( mySource ) % subState0 ( 1 : mySizeSourceDotState , c ) &
+ sourceState ( p ) % p ( mySource ) % dotState ( 1 : mySizeSourceDotState , c ) &
* crystallite_subdt ( g , i , e ) * timeStepFraction ( n )
enddo
2013-06-11 22:05:04 +05:30
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
if ( n == 4 &
. and . iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( e == debug_e . and . i == debug_i . and . g == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then ! final integration step
write ( 6 , '(a,i8,1x,i2,1x,i3,/)' ) '<< CRYST >> updateState at el ip g ' , e , i , g
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> dotState' , plasticState ( p ) % dotState ( 1 : mySizePlasticDotState , c )
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> new state' , plasticState ( p ) % state ( 1 : mySizePlasticDotState , c )
endif
2011-03-29 12:57:19 +05:30
#endif
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- state jump ---
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
crystallite_todo ( g , i , e ) = stateJump ( g , i , e )
!$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 ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- update dependent states ---
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) &
call constitutive_microstructure ( crystallite_orientation , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp ( 1 : 3 , 1 : 3 , g , i , e ) , &
g , i , e ) ! update dependent state variables to be consistent with basic states
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- stress integration ---
2014-11-13 18:23:20 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
crystallite_todo ( g , i , e ) = integrateStress ( g , i , e , timeStepFraction ( n ) ) ! fraction of original times step
!$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 ENDDO
2014-11-13 18:23:20 +05:30
2019-01-15 08:57:57 +05:30
! --- dot state and RK dot state---
first3steps : if ( n < 4 ) then
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) &
call constitutive_collectDotState ( crystallite_Tstar_v ( 1 : 6 , g , i , e ) , &
crystallite_Fe , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp , &
timeStepFraction ( n ) * crystallite_subdt ( g , i , e ) , & ! fraction of original timestep
crystallite_subFrac , g , i , e )
enddo ; enddo ; enddo
!$OMP ENDDO
!$OMP DO PRIVATE(p,c,NaN)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
c = phasememberAt ( g , i , e )
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % dotState ( : , c ) ) )
do mySource = 1_pInt , phase_Nsources ( p )
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % dotState ( : , c ) ) )
enddo
if ( NaN ) then ! NaN occured in any dotState
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)
else ! if broken local...
crystallite_todo ( g , i , e ) = . false . ! ... skip this one next time
endif
endif
endif
enddo ; enddo ; enddo
!$OMP ENDDO
endif first3steps
!$OMP END PARALLEL
enddo
! --- SET CONVERGENCE FLAG ---
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
crystallite_converged ( g , i , e ) = crystallite_todo ( g , i , e ) . or . crystallite_converged ( g , i , e ) ! if still "to do" then converged per definitionem
enddo ; enddo ; enddo
! --- CHECK NONLOCAL CONVERGENCE ---
if ( . not . singleRun ) then ! if not requesting Integration of just a single IP
if ( any ( . not . crystallite_converged . and . . not . crystallite_localPlasticity ) ) then ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged . and . crystallite_localPlasticity ! ...restart all non-local as not converged
2014-11-01 00:33:08 +05:30
endif
2019-01-15 08:57:57 +05:30
endif
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")
!--------------------------------------------------------------------------------------------------
subroutine integrateStateRKCK45 ( )
use , intrinsic :: &
IEEE_arithmetic
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
use debug , only : &
debug_e , &
debug_i , &
debug_g , &
debug_level , &
debug_crystallite , &
debug_levelBasic , &
debug_levelExtensive , &
debug_levelSelective
2011-03-29 12:57:19 +05:30
#endif
2019-01-15 08:57:57 +05:30
use numerics , only : &
rTol_crystalliteState
use FEsolving , only : &
FEsolving_execElem , &
FEsolving_execIP
use mesh , only : &
mesh_element , &
mesh_NcpElems , &
mesh_maxNips
use material , only : &
homogenization_Ngrains , &
plasticState , &
sourceState , &
phase_Nsources , &
phaseAt , phasememberAt , &
homogenization_maxNgrains
use constitutive , only : &
constitutive_collectDotState , &
constitutive_plasticity_maxSizeDotState , &
constitutive_source_maxSizeDotState , &
constitutive_microstructure
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
implicit none
real ( pReal ) , dimension ( 5 , 5 ) , parameter :: &
A = reshape ( [ &
. 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 ] , &
[ 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 - &
[ 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)
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)
integer ( pInt ) :: &
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
s , & ! state index
n , &
p , &
cc , &
mySource , &
mySizePlasticDotState , & ! size of dot States
mySizeSourceDotState
integer ( pInt ) , dimension ( 2 ) :: &
eIter ! bounds for element iteration
integer ( pInt ) , dimension ( 2 , mesh_NcpElems ) :: &
iIter , & ! bounds for ip iteration
gIter ! bounds for grain iteration
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState , &
homogenization_maxNgrains , mesh_maxNips , mesh_NcpElems ) :: &
plasticStateResiduum , & ! residuum from evolution in microstructure
relPlasticStateResiduum ! relative residuum from evolution in microstructure
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , &
maxval ( phase_Nsources ) , &
homogenization_maxNgrains , mesh_maxNips , mesh_NcpElems ) :: &
sourceStateResiduum , & ! residuum from evolution in microstructure
relSourceStateResiduum ! relative residuum from evolution in microstructure
logical :: &
NaN , &
singleRun ! flag indicating computation for single (g,i,e) triple
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
eIter = FEsolving_execElem ( 1 : 2 )
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,1x,i1)' ) '<< CRYST >> Runge--Kutta step' , 1
2014-11-01 00:33:08 +05:30
#endif
2015-10-14 00:22:01 +05:30
2019-01-15 08:57:57 +05:30
! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP ---
do e = eIter ( 1 ) , eIter ( 2 )
iIter ( 1 : 2 , e ) = FEsolving_execIP ( 1 : 2 , e )
gIter ( 1 : 2 , e ) = [ 1_pInt , homogenization_Ngrains ( mesh_element ( 3 , e ) ) ]
enddo
2014-11-01 00:33:08 +05:30
2019-01-15 08:57:57 +05:30
singleRun = ( eIter ( 1 ) == eIter ( 2 ) . and . iIter ( 1 , eIter ( 1 ) ) == iIter ( 2 , eIter ( 2 ) ) )
2014-11-01 00:33:08 +05:30
2014-11-12 22:10:50 +05:30
2014-11-01 00:33:08 +05:30
2019-01-15 08:57:57 +05:30
! --- FIRST RUNGE KUTTA STEP ---
!$OMP PARALLEL
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) &
call constitutive_collectDotState ( crystallite_Tstar_v ( 1 : 6 , g , i , e ) , &
crystallite_Fe , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp , &
crystallite_subdt ( g , i , e ) , crystallite_subFrac , g , i , e )
enddo ; enddo ; enddo
!$OMP ENDDO
!$OMP DO PRIVATE(p,cc,NaN)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
cc = phasememberAt ( g , i , e )
p = phaseAt ( g , i , e )
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % dotState ( : , cc ) ) )
do mySource = 1_pInt , phase_Nsources ( p )
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % dotState ( : , cc ) ) )
enddo
if ( NaN ) then ! NaN occured in any dotState
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)
else ! if broken local...
crystallite_todo ( g , i , e ) = . false . ! ... skip this one next time
endif
endif
2014-11-12 22:10:50 +05:30
endif
2019-01-15 08:57:57 +05:30
enddo ; enddo ; enddo
!$OMP ENDDO
!$OMP END PARALLEL
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-01-15 08:57:57 +05:30
do stage = 1_pInt , 5_pInt
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
! --- state update ---
!$OMP PARALLEL
!$OMP DO PRIVATE(p,cc)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
cc = phasememberAt ( g , i , e )
plasticState ( p ) % RKCK45dotState ( stage , : , cc ) = plasticState ( p ) % dotState ( : , cc ) ! store Runge-Kutta dotState
do mySource = 1_pInt , phase_Nsources ( p )
sourceState ( p ) % p ( mySource ) % RKCK45dotState ( stage , : , cc ) = sourceState ( p ) % p ( mySource ) % dotState ( : , cc )
enddo
2017-12-14 05:48:45 +05:30
endif
2019-01-15 08:57:57 +05:30
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO PRIVATE(p,cc,n)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
cc = phasememberAt ( g , i , e )
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
plasticState ( p ) % dotState ( : , cc ) = A ( 1 , stage ) * plasticState ( p ) % RKCK45dotState ( 1 , : , cc )
do mySource = 1_pInt , phase_Nsources ( p )
sourceState ( p ) % p ( mySource ) % dotState ( : , cc ) = A ( 1 , stage ) * sourceState ( p ) % p ( mySource ) % RKCK45dotState ( 1 , : , cc )
enddo
do n = 2_pInt , stage
plasticState ( p ) % dotState ( : , cc ) = &
plasticState ( p ) % dotState ( : , cc ) + A ( n , stage ) * plasticState ( p ) % RKCK45dotState ( n , : , cc )
do mySource = 1_pInt , phase_Nsources ( p )
sourceState ( p ) % p ( mySource ) % dotState ( : , cc ) = &
sourceState ( p ) % p ( mySource ) % dotState ( : , cc ) + A ( n , stage ) * sourceState ( p ) % p ( mySource ) % RKCK45dotState ( n , : , cc )
enddo
enddo
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-11-01 00:33:08 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
cc = phasememberAt ( g , i , e )
mySizePlasticDotState = plasticState ( p ) % sizeDotState
plasticState ( p ) % state ( 1 : mySizePlasticDotState , cc ) = &
plasticState ( p ) % subState0 ( 1 : mySizePlasticDotState , cc ) &
+ plasticState ( p ) % dotState ( 1 : mySizePlasticDotState , cc ) &
* crystallite_subdt ( g , i , e )
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , cc ) = &
sourceState ( p ) % p ( mySource ) % subState0 ( 1 : mySizeSourceDotState , cc ) &
+ sourceState ( p ) % p ( mySource ) % dotState ( 1 : mySizeSourceDotState , cc ) &
* crystallite_subdt ( g , i , e )
enddo
2017-12-14 05:48:45 +05:30
endif
2019-01-15 08:57:57 +05:30
enddo ; enddo ; enddo
!$OMP ENDDO
! --- state jump ---
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
crystallite_todo ( g , i , e ) = stateJump ( g , i , e )
!$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)
2014-11-12 22:10:50 +05:30
endif
2014-11-01 00:33:08 +05:30
endif
2019-01-15 08:57:57 +05:30
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- update dependent states ---
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) &
call constitutive_microstructure ( crystallite_orientation , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp ( 1 : 3 , 1 : 3 , g , i , e ) , &
g , i , e ) ! update dependent state variables to be consistent with basic states
enddo ; enddo ; enddo
!$OMP ENDDO
! --- stress integration ---
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
crystallite_todo ( g , i , e ) = integrateStress ( g , i , e , C ( stage ) ) ! fraction of original time step
!$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 ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- dot state and RK dot state---
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,1x,i1)' ) '<< CRYST >> Runge--Kutta step' , stage + 1_pInt
2015-08-05 02:56:22 +05:30
#endif
2019-01-15 08:57:57 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) &
call constitutive_collectDotState ( crystallite_Tstar_v ( 1 : 6 , g , i , e ) , &
crystallite_Fe , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp , &
C ( stage ) * crystallite_subdt ( g , i , e ) , & ! fraction of original timestep
crystallite_subFrac , g , i , e )
enddo ; enddo ; enddo
!$OMP ENDDO
!$OMP DO PRIVATE(p,cc,NaN)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
p = phaseAt ( g , i , e )
cc = phasememberAt ( g , i , e )
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % dotState ( : , cc ) ) )
do mySource = 1_pInt , phase_Nsources ( p )
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % dotState ( : , cc ) ) )
enddo
if ( NaN ) then ! NaN occured in any dotState
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)
else ! if broken local...
crystallite_todo ( g , i , e ) = . false . ! ... skip this one next time
endif
endif
endif
enddo ; enddo ; enddo
!$OMP ENDDO
!$OMP END PARALLEL
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 ---
relPlasticStateResiduum = 0.0_pReal
relSourceStateResiduum = 0.0_pReal
!$OMP PARALLEL
!$OMP DO PRIVATE(p,cc)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
cc = phasememberAt ( g , i , e )
plasticState ( p ) % RKCK45dotState ( 6 , : , cc ) = plasticState ( p ) % dotState ( : , cc ) ! store Runge-Kutta dotState
do mySource = 1_pInt , phase_Nsources ( p )
sourceState ( p ) % p ( mySource ) % RKCK45dotState ( 6 , : , cc ) = sourceState ( p ) % p ( mySource ) % dotState ( : , cc ) ! store Runge-Kutta dotState
enddo
2016-03-21 03:50:58 +05:30
endif
2019-01-15 08:57:57 +05:30
enddo ; enddo ; enddo
!$OMP ENDDO
2015-08-05 02:56:22 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
cc = phasememberAt ( g , i , e )
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- absolute residuum in state ---
mySizePlasticDotState = plasticState ( p ) % sizeDotState
plasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) = &
matmul ( transpose ( plasticState ( p ) % RKCK45dotState ( 1 : 6 , 1 : mySizePlasticDotState , cc ) ) , DB ) &
* crystallite_subdt ( g , i , e )
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceStateResiduum ( 1 : mySizeSourceDotState , mySource , g , i , e ) = &
matmul ( transpose ( sourceState ( p ) % p ( mySource ) % RKCK45dotState ( 1 : 6 , 1 : mySizeSourceDotState , cc ) ) , DB ) &
* crystallite_subdt ( g , i , e )
enddo
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- dot state ---
plasticState ( p ) % dotState ( : , cc ) = &
matmul ( transpose ( plasticState ( p ) % RKCK45dotState ( 1 : 6 , 1 : mySizePlasticDotState , cc ) ) , B )
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceState ( p ) % p ( mySource ) % dotState ( : , cc ) = &
matmul ( transpose ( sourceState ( p ) % p ( mySource ) % RKCK45dotState ( 1 : 6 , 1 : mySizeSourceDotState , cc ) ) , B )
enddo
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- state and update ---
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
p = phaseAt ( g , i , e )
cc = phasememberAt ( g , i , e )
mySizePlasticDotState = plasticState ( p ) % sizeDotState
plasticState ( p ) % state ( 1 : mySizePlasticDotState , cc ) = &
plasticState ( p ) % subState0 ( 1 : mySizePlasticDotState , cc ) &
+ plasticState ( p ) % dotState ( 1 : mySizePlasticDotState , cc ) &
* crystallite_subdt ( g , i , e )
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
sourceState ( p ) % p ( mySource ) % state ( 1 : mySizeSourceDotState , cc ) = &
sourceState ( p ) % p ( mySource ) % subState0 ( 1 : mySizeSourceDotState , cc ) &
+ sourceState ( p ) % p ( mySource ) % dotState ( 1 : mySizeSourceDotState , cc ) &
* crystallite_subdt ( g , i , e )
enddo
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- relative residui and state convergence ---
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc,s)
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) then
p = phaseAt ( g , i , e )
cc = phasememberAt ( g , i , e )
mySizePlasticDotState = plasticState ( p ) % sizeDotState
forall ( s = 1_pInt : mySizePlasticDotState , abs ( plasticState ( p ) % state ( s , cc ) ) > 0.0_pReal ) &
relPlasticStateResiduum ( s , g , i , e ) = &
plasticStateResiduum ( s , g , i , e ) / plasticState ( p ) % state ( s , cc )
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
forall ( s = 1_pInt : mySizeSourceDotState , abs ( sourceState ( p ) % p ( mySource ) % state ( s , cc ) ) > 0.0_pReal ) &
relSourceStateResiduum ( s , mySource , g , i , e ) = &
sourceStateResiduum ( s , mySource , g , i , e ) / sourceState ( p ) % p ( mySource ) % state ( s , cc )
enddo
crystallite_todo ( g , i , e ) = all ( abs ( relPlasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) ) < &
rTol_crystalliteState . or . &
abs ( plasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) ) < &
plasticState ( p ) % aTolState ( 1 : mySizePlasticDotState ) )
do mySource = 1_pInt , phase_Nsources ( p )
mySizeSourceDotState = sourceState ( p ) % p ( mySource ) % sizeDotState
crystallite_todo ( g , i , e ) = crystallite_todo ( g , i , e ) . and . &
all ( abs ( relSourceStateResiduum ( 1 : mySizeSourceDotState , mySource , g , i , e ) ) < &
rTol_crystalliteState . or . &
abs ( sourceStateResiduum ( 1 : mySizeSourceDotState , mySource , g , i , e ) ) < &
sourceState ( p ) % p ( mySource ) % aTolState ( 1 : mySizeSourceDotState ) )
enddo
2014-08-26 20:14:32 +05:30
2017-09-30 04:02:07 +05:30
#ifdef DEBUG
2019-01-15 08:57:57 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( e == debug_e . and . i == debug_i . and . g == debug_g ) &
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,i8,1x,i3,1x,i3,/)' ) '<< CRYST >> updateState at el ip ipc ' , e , i , g
write ( 6 , '(a,/,(12x,12(f12.1,1x)),/)' ) '<< CRYST >> absolute residuum tolerance' , &
plasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) / plasticState ( p ) % aTolState ( 1 : mySizePlasticDotState )
write ( 6 , '(a,/,(12x,12(f12.1,1x)),/)' ) '<< CRYST >> relative residuum tolerance' , &
relPlasticStateResiduum ( 1 : mySizePlasticDotState , g , i , e ) / rTol_crystalliteState
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> dotState' , &
plasticState ( p ) % dotState ( 1 : mySizePlasticDotState , cc )
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> new state' , &
plasticState ( p ) % state ( 1 : mySizePlasticDotState , cc )
endif
2011-03-29 12:57:19 +05:30
#endif
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
! --- STATE JUMP ---
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
crystallite_todo ( g , i , e ) = stateJump ( g , i , e )
!$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)
2013-02-22 04:38:36 +05:30
endif
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2014-08-26 20:14:32 +05:30
2009-12-18 21:16:33 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
! --- UPDATE DEPENDENT STATES IF RESIDUUM BELOW TOLERANCE ---
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
if ( crystallite_todo ( g , i , e ) ) &
call constitutive_microstructure ( crystallite_orientation , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp ( 1 : 3 , 1 : 3 , g , i , e ) , &
g , i , e ) ! update dependent state variables to be consistent with basic states
enddo ; enddo ; enddo
!$OMP ENDDO
2015-10-14 00:22:01 +05:30
2010-02-25 23:09:11 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
! --- FINAL STRESS INTEGRATION STEP IF RESIDUUM BELOW TOLERANCE ---
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) ) then
crystallite_todo ( g , i , e ) = integrateStress ( g , i , e )
!$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 ENDDO
2010-02-25 23:09:11 +05:30
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
! --- SET CONVERGENCE FLAG ---
!$OMP DO
do e = eIter ( 1 ) , eIter ( 2 ) ; do i = iIter ( 1 , e ) , iIter ( 2 , e ) ; do g = gIter ( 1 , e ) , gIter ( 2 , e ) ! iterate over elements, ips and grains
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
!$OMP ENDDO
2011-05-11 22:08:45 +05:30
2019-01-15 08:57:57 +05:30
!$OMP END PARALLEL
2014-08-26 20:14:32 +05:30
2011-05-11 22:08:45 +05:30
2019-01-15 08:57:57 +05:30
! --- nonlocal convergence check ---
#ifdef DEBUG
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,i8,a,i2,/)' ) '<< CRYST >> ' , count ( crystallite_converged ( : , : , : ) ) , ' grains converged' ! if not requesting Integration of just a single IP
#endif
if ( ( . not . singleRun ) . and . any ( . not . crystallite_converged . and . . not . crystallite_localPlasticity ) ) & ! any non-local not yet converged (or broken)...
crystallite_converged = crystallite_converged . and . crystallite_localPlasticity ! ...restart all non-local as not converged
2014-06-25 04:51:25 +05:30
2019-01-15 08:57:57 +05:30
end subroutine integrateStateRKCK45
2009-05-07 21:57:36 +05:30
2013-02-22 04:38:36 +05:30
end module crystallite