2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-18 16:46:26 +05:30
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
2016-02-29 18:56:06 +05:30
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
2013-02-22 04:38:36 +05:30
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
2016-02-29 18:56:06 +05:30
!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Chen Zhang, Michigan State University
2013-02-22 04:38:36 +05:30
!> @brief crystallite state integration functions and reporting of results
!--------------------------------------------------------------------------------------------------
2010-10-01 17:48:49 +05:30
2012-08-31 01:56:28 +05:30
module crystallite
2019-05-17 02:26:48 +05:30
use prec
2020-09-13 13:49:38 +05:30
use parallelization
2019-05-17 02:26:48 +05:30
use IO
2020-02-25 14:10:38 +05:30
use HDF5_utilities
use DAMASK_interface
2019-05-17 02:26:48 +05:30
use config
use rotations
use math
use FEsolving
use material
use constitutive
2019-06-07 02:19:17 +05:30
use discretization
2019-05-17 02:26:48 +05:30
use lattice
use results
2020-02-13 21:40:27 +05:30
implicit none
2019-04-03 16:24:07 +05:30
private
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
real ( pReal ) , dimension ( : , : , : ) , allocatable , public :: &
crystallite_dt !< requested time increment of each grain
2019-05-17 02:26:48 +05:30
real ( pReal ) , dimension ( : , : , : ) , allocatable :: &
2019-04-03 16:24:07 +05:30
crystallite_subdt , & !< substepped time increment of each grain
crystallite_subFrac , & !< already calculated fraction of increment
crystallite_subStep !< size of next integration step
2019-05-17 02:26:48 +05:30
type ( rotation ) , dimension ( : , : , : ) , allocatable :: &
2020-02-13 21:40:27 +05:30
crystallite_orientation !< current orientation
2020-09-28 21:31:43 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: &
2020-02-25 14:20:21 +05:30
crystallite_F0 , & !< def grad at start of FE inc
2020-09-29 12:33:58 +05:30
crystallite_subF , & !< def grad to be reached at end of crystallite inc
crystallite_subF0 , & !< def grad at start of crystallite inc
!
crystallite_Fe , & !< current "elastic" def grad (end of converged time step)
!
2019-04-03 16:24:07 +05:30
crystallite_Fp , & !< current plastic def grad (end of converged time step)
2020-09-29 12:33:58 +05:30
crystallite_Fp0 , & !< plastic def grad at start of FE inc
2020-10-08 01:45:13 +05:30
crystallite_partitionedFp0 , & !< plastic def grad at start of homog inc
2020-09-29 12:33:58 +05:30
crystallite_subFp0 , & !< plastic def grad at start of crystallite inc
!
2019-04-03 16:24:07 +05:30
crystallite_Fi , & !< current intermediate def grad (end of converged time step)
2020-09-29 12:33:58 +05:30
crystallite_Fi0 , & !< intermediate def grad at start of FE inc
2020-10-08 01:45:13 +05:30
crystallite_partitionedFi0 , & !< intermediate def grad at start of homog inc
2020-09-29 12:33:58 +05:30
crystallite_subFi0 , & !< intermediate def grad at start of crystallite inc
!
crystallite_Lp0 , & !< plastic velocitiy grad at start of FE inc
2020-10-08 01:45:13 +05:30
crystallite_partitionedLp0 , & !< plastic velocity grad at start of homog inc
2020-09-29 12:33:58 +05:30
!
2019-04-03 16:24:07 +05:30
crystallite_Li , & !< current intermediate velocitiy grad (end of converged time step)
2020-09-29 12:33:58 +05:30
crystallite_Li0 , & !< intermediate velocitiy grad at start of FE inc
2020-10-08 01:45:13 +05:30
crystallite_partitionedLi0 , & !< intermediate velocity grad at start of homog inc
2020-09-29 12:33:58 +05:30
!
crystallite_S0 , & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
2020-10-08 01:45:13 +05:30
crystallite_partitionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc
2020-09-28 21:31:43 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , public , protected :: &
crystallite_P , & !< 1st Piola-Kirchhoff stress per grain
crystallite_Lp , & !< current plastic velocitiy grad (end of converged time step)
crystallite_S , & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
2020-10-08 01:45:13 +05:30
crystallite_partitionedF0 !< def grad at start of homog inc
2020-09-28 21:31:43 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , public :: &
2020-10-08 01:45:13 +05:30
crystallite_partitionedF !< def grad to be reached at end of homog inc
2020-09-28 21:31:43 +05:30
2019-04-03 16:24:07 +05:30
logical , dimension ( : , : , : ) , allocatable , public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
2019-05-17 02:26:48 +05:30
logical , dimension ( : , : , : ) , allocatable :: &
2020-04-01 22:28:48 +05:30
crystallite_converged !< convergence flag
2020-02-13 21:40:27 +05:30
2019-05-17 02:26:48 +05:30
type :: tOutput !< new requested output (per phase)
2019-12-21 17:07:02 +05:30
character ( len = pStringLen ) , allocatable , dimension ( : ) :: &
2019-04-06 15:36:34 +05:30
label
end type tOutput
2019-05-17 02:26:48 +05:30
type ( tOutput ) , allocatable , dimension ( : ) :: output_constituent
2020-02-13 21:40:27 +05:30
2019-05-17 02:26:48 +05:30
type :: tNumerics
2019-04-11 14:57:03 +05:30
integer :: &
iJacoLpresiduum , & !< frequency of Jacobian update of residuum in Lp
nState , & !< state loop limit
2020-06-16 21:23:14 +05:30
nStress !< stress loop limit
2019-04-11 10:54:04 +05:30
real ( pReal ) :: &
subStepMinCryst , & !< minimum (relative) size of sub-step allowed during cutback
subStepSizeCryst , & !< size of first substep when cutback
subStepSizeLp , & !< size of first substep when cutback in Lp calculation
subStepSizeLi , & !< size of first substep when cutback in Li calculation
stepIncreaseCryst , & !< increase of next substep size when previous substep converged
2020-03-14 23:41:26 +05:30
rtol_crystalliteState , & !< relative tolerance in state loop
rtol_crystalliteStress , & !< relative tolerance in stress loop
atol_crystalliteStress !< absolute tolerance in stress loop
2019-04-11 10:54:04 +05:30
end type tNumerics
2020-02-13 21:40:27 +05:30
2019-04-11 10:54:04 +05:30
type ( tNumerics ) :: num ! numerics parameters. Better name?
2020-02-13 21:40:27 +05:30
2020-07-01 23:24:14 +05:30
type :: tDebugOptions
logical :: &
basic , &
extensive , &
selective
integer :: &
element , &
ip , &
grain
end type tDebugOptions
2020-07-02 04:55:24 +05:30
type ( tDebugOptions ) :: debugCrystallite
2020-07-01 23:24:14 +05:30
2020-04-01 12:24:20 +05:30
procedure ( integrateStateFPI ) , pointer :: integrateState
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
public :: &
crystallite_init , &
crystallite_stress , &
crystallite_stressTangent , &
crystallite_orientations , &
crystallite_push33ToRef , &
2020-02-25 14:10:38 +05:30
crystallite_results , &
2020-02-25 22:23:15 +05:30
crystallite_restartWrite , &
crystallite_restartRead , &
2020-09-28 21:26:48 +05:30
crystallite_forward , &
crystallite_initializeRestorationPoints , &
crystallite_windForward , &
crystallite_restore
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
2020-02-13 21:40:27 +05:30
2020-10-28 01:57:26 +05:30
logical , dimension ( discretization_nIPs , discretization_Nelems ) :: devNull
2019-04-03 16:24:07 +05:30
integer :: &
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e , & !< counter in element loop
cMax , & !< maximum number of integration point components
iMax , & !< maximum number of integration points
2020-09-28 21:31:43 +05:30
eMax !< maximum number of elements
2020-08-13 00:44:00 +05:30
2020-07-03 20:15:11 +05:30
class ( tNode ) , pointer :: &
2020-06-18 18:51:52 +05:30
num_crystallite , &
2020-08-15 19:32:10 +05:30
debug_crystallite , & ! pointer to debug options for crystallite
phases , &
phase , &
2020-11-03 03:16:46 +05:30
mech
2020-06-18 18:51:52 +05:30
2020-09-17 22:58:41 +05:30
print '(/,a)' , ' <<<+- crystallite init -+>>>'
2020-02-13 21:40:27 +05:30
2020-09-13 14:09:17 +05:30
debug_crystallite = > config_debug % get ( 'crystallite' , defaultVal = emptyList )
2020-07-02 04:55:24 +05:30
debugCrystallite % basic = debug_crystallite % contains ( 'basic' )
debugCrystallite % extensive = debug_crystallite % contains ( 'extensive' )
debugCrystallite % selective = debug_crystallite % contains ( 'selective' )
2020-09-13 14:09:17 +05:30
debugCrystallite % element = config_debug % get_asInt ( 'element' , defaultVal = 1 )
debugCrystallite % ip = config_debug % get_asInt ( 'integrationpoint' , defaultVal = 1 )
debugCrystallite % grain = config_debug % get_asInt ( 'grain' , defaultVal = 1 )
2020-07-01 23:24:14 +05:30
2020-10-28 01:57:26 +05:30
cMax = homogenization_maxNconstituents
iMax = discretization_nIPs
eMax = discretization_Nelems
2020-02-13 21:40:27 +05:30
2020-10-08 01:45:13 +05:30
allocate ( crystallite_partitionedF ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
2020-02-29 18:23:04 +05:30
allocate ( crystallite_S0 , &
crystallite_F0 , crystallite_Fi0 , crystallite_Fp0 , &
crystallite_Li0 , crystallite_Lp0 , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedS0 , &
crystallite_partitionedF0 , crystallite_partitionedFp0 , crystallite_partitionedFi0 , &
crystallite_partitionedLp0 , crystallite_partitionedLi0 , &
2020-02-29 18:23:04 +05:30
crystallite_S , crystallite_P , &
crystallite_Fe , crystallite_Fi , crystallite_Fp , &
crystallite_Li , crystallite_Lp , &
crystallite_subF , crystallite_subF0 , &
crystallite_subFp0 , crystallite_subFi0 , &
2020-10-08 01:45:13 +05:30
source = crystallite_partitionedF )
2020-02-29 18:23:04 +05:30
allocate ( crystallite_dt ( cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subdt , crystallite_subFrac , crystallite_subStep , &
source = crystallite_dt )
2019-04-03 16:24:07 +05:30
allocate ( crystallite_orientation ( cMax , iMax , eMax ) )
2020-02-29 18:23:04 +05:30
2019-04-03 16:24:07 +05:30
allocate ( crystallite_requested ( cMax , iMax , eMax ) , source = . false . )
allocate ( crystallite_converged ( cMax , iMax , eMax ) , source = . true . )
2020-02-13 21:40:27 +05:30
2020-09-13 14:09:17 +05:30
num_crystallite = > config_numerics % get ( 'crystallite' , defaultVal = emptyDict )
2020-06-29 18:39:13 +05:30
2020-06-16 21:23:14 +05:30
num % subStepMinCryst = num_crystallite % get_asFloat ( 'subStepMin' , defaultVal = 1.0e-3_pReal )
num % subStepSizeCryst = num_crystallite % get_asFloat ( 'subStepSize' , defaultVal = 0.25_pReal )
num % stepIncreaseCryst = num_crystallite % get_asFloat ( 'stepIncrease' , defaultVal = 1.5_pReal )
num % subStepSizeLp = num_crystallite % get_asFloat ( 'subStepSizeLp' , defaultVal = 0.5_pReal )
num % subStepSizeLi = num_crystallite % get_asFloat ( 'subStepSizeLi' , defaultVal = 0.5_pReal )
num % rtol_crystalliteState = num_crystallite % get_asFloat ( 'rtol_State' , defaultVal = 1.0e-6_pReal )
num % rtol_crystalliteStress = num_crystallite % get_asFloat ( 'rtol_Stress' , defaultVal = 1.0e-6_pReal )
num % atol_crystalliteStress = num_crystallite % get_asFloat ( 'atol_Stress' , defaultVal = 1.0e-8_pReal )
num % iJacoLpresiduum = num_crystallite % get_asInt ( 'iJacoLpresiduum' , defaultVal = 1 )
num % nState = num_crystallite % get_asInt ( 'nState' , defaultVal = 20 )
num % nStress = num_crystallite % get_asInt ( 'nStress' , defaultVal = 40 )
2020-02-13 21:40:27 +05:30
2019-04-11 14:57:03 +05:30
if ( num % subStepMinCryst < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepMinCryst' )
if ( num % subStepSizeCryst < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepSizeCryst' )
if ( num % stepIncreaseCryst < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'stepIncreaseCryst' )
if ( num % subStepSizeLp < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepSizeLp' )
if ( num % subStepSizeLi < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepSizeLi' )
2020-03-14 23:41:26 +05:30
if ( num % rtol_crystalliteState < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'rtol_crystalliteState' )
if ( num % rtol_crystalliteStress < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'rtol_crystalliteStress' )
if ( num % atol_crystalliteStress < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'atol_crystalliteStress' )
2020-02-13 21:40:27 +05:30
2019-04-11 14:57:03 +05:30
if ( num % iJacoLpresiduum < 1 ) call IO_error ( 301 , ext_msg = 'iJacoLpresiduum' )
2020-02-13 21:40:27 +05:30
2019-04-11 14:57:03 +05:30
if ( num % nState < 1 ) call IO_error ( 301 , ext_msg = 'nState' )
if ( num % nStress < 1 ) call IO_error ( 301 , ext_msg = 'nStress' )
2020-02-13 21:40:27 +05:30
2020-09-13 16:01:01 +05:30
select case ( num_crystallite % get_asString ( 'integrator' , defaultVal = 'FPI' ) )
2020-06-16 21:23:14 +05:30
case ( 'FPI' )
2019-04-03 16:24:07 +05:30
integrateState = > integrateStateFPI
2020-06-16 21:23:14 +05:30
case ( 'Euler' )
2019-04-03 16:24:07 +05:30
integrateState = > integrateStateEuler
2020-06-16 21:23:14 +05:30
case ( 'AdaptiveEuler' )
2019-04-03 16:24:07 +05:30
integrateState = > integrateStateAdaptiveEuler
2020-06-16 21:23:14 +05:30
case ( 'RK4' )
2019-04-03 16:24:07 +05:30
integrateState = > integrateStateRK4
2020-06-16 21:23:14 +05:30
case ( 'RKCK45' )
2019-04-03 16:24:07 +05:30
integrateState = > integrateStateRKCK45
2020-06-16 21:23:14 +05:30
case default
call IO_error ( 301 , ext_msg = 'integrator' )
2019-04-03 16:24:07 +05:30
end select
2020-02-13 21:40:27 +05:30
2020-09-13 14:09:17 +05:30
phases = > config_material % get ( 'phase' )
2020-08-15 19:32:10 +05:30
allocate ( output_constituent ( phases % length ) )
do c = 1 , phases % length
phase = > phases % get ( c )
2020-11-18 01:54:40 +05:30
mech = > phase % get ( 'mechanics' , defaultVal = emptyDict )
2019-04-06 10:01:02 +05:30
#if defined(__GFORTRAN__)
2020-11-03 03:16:46 +05:30
output_constituent ( c ) % label = output_asStrings ( mech )
2019-04-06 10:01:02 +05:30
#else
2020-11-03 03:16:46 +05:30
output_constituent ( c ) % label = mech % get_asStrings ( 'output' , defaultVal = emptyStringArray )
2019-04-06 10:01:02 +05:30
#endif
enddo
2019-04-03 16:24:07 +05:30
2018-06-27 00:03:41 +05:30
2013-04-29 16:47:30 +05:30
!--------------------------------------------------------------------------------------------------
! initialize
2020-09-28 21:31:43 +05:30
!$OMP PARALLEL DO PRIVATE(i,c)
2019-04-03 16:24:07 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-10-28 01:57:26 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 ) ; do c = 1 , homogenization_Nconstituents ( material_homogenizationAt ( e ) )
2020-09-14 18:41:48 +05:30
crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e ) = material_orientation0 ( c , i , e ) % asMatrix ( ) ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
2020-02-13 21:40:27 +05:30
crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e ) &
/ math_det33 ( crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e ) ) ** ( 1.0_pReal / 3.0_pReal )
2019-04-03 16:24:07 +05:30
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_Fe ( 1 : 3 , 1 : 3 , c , i , e ) = math_inv33 ( matmul ( crystallite_Fi0 ( 1 : 3 , 1 : 3 , c , i , e ) , &
2019-04-03 21:54:15 +05:30
crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e ) ) ) ! assuming that euler angles are given in internal strain free configuration
2019-04-03 16:24:07 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_requested ( c , i , e ) = . true .
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-04-03 16:24:07 +05:30
enddo
!$OMP END PARALLEL DO
2020-02-13 21:40:27 +05:30
2020-10-08 01:45:13 +05:30
crystallite_partitionedFp0 = crystallite_Fp0
crystallite_partitionedFi0 = crystallite_Fi0
crystallite_partitionedF0 = crystallite_F0
crystallite_partitionedF = crystallite_F0
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
call crystallite_orientations ( )
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
!$OMP PARALLEL DO
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2020-10-28 01:57:26 +05:30
do c = 1 , homogenization_Nconstituents ( material_homogenizationAt ( e ) )
2020-10-08 01:45:13 +05:30
call constitutive_dependentState ( crystallite_partitionedF0 ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_partitionedFp0 ( 1 : 3 , 1 : 3 , c , i , e ) , &
2019-04-03 16:24:07 +05:30
c , i , e ) ! update dependent state variables to be consistent with basic states
enddo
2019-01-19 03:39:46 +05:30
enddo
2019-04-03 16:24:07 +05:30
enddo
!$OMP END PARALLEL DO
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
devNull = crystallite_stress ( )
2019-01-18 19:12:44 +05:30
2018-09-20 09:39:02 +05:30
#ifdef DEBUG
2020-07-02 04:55:24 +05:30
if ( debugCrystallite % basic ) then
2020-09-17 22:58:41 +05:30
print '(a42,1x,i10)' , ' # of elements: ' , eMax
print '(a42,1x,i10)' , ' # of integration points/element: ' , iMax
print '(a42,1x,i10)' , 'max # of constituents/integration point: ' , cMax
2020-09-22 16:39:12 +05:30
flush ( IO_STDOUT )
2019-04-03 16:24:07 +05:30
endif
2018-09-20 09:39:02 +05:30
#endif
2010-11-03 22:52:48 +05:30
2012-03-09 01:55:28 +05:30
end subroutine crystallite_init
2009-05-07 21:57:36 +05:30
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
2020-06-17 03:07:24 +05:30
function crystallite_stress ( )
2020-02-13 21:40:27 +05:30
2020-10-28 01:57:26 +05:30
logical , dimension ( discretization_nIPs , discretization_Nelems ) :: crystallite_stress
2019-04-03 16:24:07 +05:30
real ( pReal ) :: &
formerSubStep
integer :: &
NiterationCrystallite , & ! number of iterations in crystallite loop
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e , & !< counter in element loop
2020-07-01 23:24:14 +05:30
s
2020-10-28 01:57:26 +05:30
logical , dimension ( homogenization_maxNconstituents , discretization_nIPs , discretization_Nelems ) :: todo !ToDo: need to set some values to false for different Ngrains
2020-09-30 03:02:17 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: &
subLp0 , & !< plastic velocity grad at start of crystallite inc
subLi0 !< intermediate velocity grad at start of crystallite inc
todo = . false .
2020-10-08 01:45:13 +05:30
subLp0 = crystallite_partitionedLp0
subLi0 = crystallite_partitionedLi0
2020-09-30 03:02:17 +05:30
2020-07-01 23:24:14 +05:30
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
! initialize to starting condition
2019-04-03 16:24:07 +05:30
crystallite_subStep = 0.0_pReal
!$OMP PARALLEL DO
elementLooping1 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-10-28 01:57:26 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 ) ; do c = 1 , homogenization_Nconstituents ( material_homogenizationAt ( e ) )
2019-04-03 16:24:07 +05:30
homogenizationRequestsCalculation : if ( crystallite_requested ( c , i , e ) ) then
2019-06-14 12:19:05 +05:30
plasticState ( material_phaseAt ( c , e ) ) % subState0 ( : , material_phaseMemberAt ( c , i , e ) ) = &
2020-10-08 01:45:13 +05:30
plasticState ( material_phaseAt ( c , e ) ) % partitionedState0 ( : , material_phaseMemberAt ( c , i , e ) )
2019-04-03 16:24:07 +05:30
2019-06-14 12:19:05 +05:30
do s = 1 , phase_Nsources ( material_phaseAt ( c , e ) )
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % subState0 ( : , material_phaseMemberAt ( c , i , e ) ) = &
2020-10-08 01:45:13 +05:30
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % partitionedState0 ( : , material_phaseMemberAt ( c , i , e ) )
2019-04-03 16:24:07 +05:30
enddo
2020-10-08 01:45:13 +05:30
crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partitionedFp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partitionedFi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partitionedF0 ( 1 : 3 , 1 : 3 , c , i , e )
2019-04-03 16:24:07 +05:30
crystallite_subFrac ( c , i , e ) = 0.0_pReal
2019-04-11 10:54:04 +05:30
crystallite_subStep ( c , i , e ) = 1.0_pReal / num % subStepSizeCryst
2020-04-01 12:24:20 +05:30
todo ( c , i , e ) = . true .
2019-04-03 16:24:07 +05:30
crystallite_converged ( c , i , e ) = . false . ! pretend failed step of 1/subStepSizeCryst
endif homogenizationRequestsCalculation
enddo ; enddo
enddo elementLooping1
!$OMP END PARALLEL DO
NiterationCrystallite = 0
2020-10-23 19:03:01 +05:30
cutbackLooping : do while ( any ( todo ( : , FEsolving_execIP ( 1 ) : FEsolving_execIP ( 2 ) , FEsolving_execELem ( 1 ) : FEsolving_execElem ( 2 ) ) ) )
2019-04-03 16:24:07 +05:30
NiterationCrystallite = NiterationCrystallite + 1
2019-02-27 00:17:46 +05:30
2019-01-18 19:12:44 +05:30
#ifdef DEBUG
2020-07-02 04:55:24 +05:30
if ( debugCrystallite % extensive ) &
2020-09-17 22:58:41 +05:30
print '(a,i6)' , '<< CRYST stress >> crystallite iteration ' , NiterationCrystallite
2019-01-18 19:12:44 +05:30
#endif
2019-04-03 16:24:07 +05:30
!$OMP PARALLEL DO PRIVATE(formerSubStep)
elementLooping3 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2020-10-28 01:57:26 +05:30
do c = 1 , homogenization_Nconstituents ( material_homogenizationAt ( e ) )
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
! wind forward
2019-04-03 16:24:07 +05:30
if ( crystallite_converged ( c , i , e ) ) then
formerSubStep = crystallite_subStep ( c , i , e )
crystallite_subFrac ( c , i , e ) = crystallite_subFrac ( c , i , e ) + crystallite_subStep ( c , i , e )
crystallite_subStep ( c , i , e ) = min ( 1.0_pReal - crystallite_subFrac ( c , i , e ) , &
2019-04-11 10:54:04 +05:30
num % stepIncreaseCryst * crystallite_subStep ( c , i , e ) )
2019-04-03 16:24:07 +05:30
2020-04-01 12:24:20 +05:30
todo ( c , i , e ) = crystallite_subStep ( c , i , e ) > 0.0_pReal ! still time left to integrate on?
if ( todo ( c , i , e ) ) then
2019-04-03 16:24:07 +05:30
crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e )
2020-09-30 03:02:17 +05:30
subLp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e )
subLi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e )
2019-04-03 16:24:07 +05:30
crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e )
2019-06-14 12:19:05 +05:30
plasticState ( material_phaseAt ( c , e ) ) % subState0 ( : , material_phaseMemberAt ( c , i , e ) ) &
= plasticState ( material_phaseAt ( c , e ) ) % state ( : , material_phaseMemberAt ( c , i , e ) )
do s = 1 , phase_Nsources ( material_phaseAt ( c , e ) )
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % subState0 ( : , material_phaseMemberAt ( c , i , e ) ) &
= sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % state ( : , material_phaseMemberAt ( c , i , e ) )
2019-04-03 16:24:07 +05:30
enddo
endif
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore)
2019-04-03 16:24:07 +05:30
else
2019-04-11 10:54:04 +05:30
crystallite_subStep ( c , i , e ) = num % subStepSizeCryst * crystallite_subStep ( c , i , e )
2019-04-03 16:24:07 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_S0 ( 1 : 3 , 1 : 3 , c , i , e )
if ( crystallite_subStep ( c , i , e ) < 1.0_pReal ) then ! actual (not initial) cutback
2020-09-30 03:02:17 +05:30
crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e ) = subLp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e ) = subLi0 ( 1 : 3 , 1 : 3 , c , i , e )
2019-04-03 16:24:07 +05:30
endif
2019-06-14 12:19:05 +05:30
plasticState ( material_phaseAt ( c , e ) ) % state ( : , material_phaseMemberAt ( c , i , e ) ) &
= plasticState ( material_phaseAt ( c , e ) ) % subState0 ( : , material_phaseMemberAt ( c , i , e ) )
do s = 1 , phase_Nsources ( material_phaseAt ( c , e ) )
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % state ( : , material_phaseMemberAt ( c , i , e ) ) &
= sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % subState0 ( : , material_phaseMemberAt ( c , i , e ) )
2019-04-03 16:24:07 +05:30
enddo
! cant restore dotState here, since not yet calculated in first cutback after initialization
2020-04-01 12:24:20 +05:30
todo ( c , i , e ) = crystallite_subStep ( c , i , e ) > num % subStepMinCryst ! still on track or already done (beyond repair)
2019-04-03 16:24:07 +05:30
endif
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
! prepare for integration
2020-04-01 12:24:20 +05:30
if ( todo ( c , i , e ) ) then
2019-04-03 16:24:07 +05:30
crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) &
2020-10-08 01:45:13 +05:30
+ crystallite_subStep ( c , i , e ) * ( crystallite_partitionedF ( 1 : 3 , 1 : 3 , c , i , e ) &
- crystallite_partitionedF0 ( 1 : 3 , 1 : 3 , c , i , e ) )
2020-10-23 19:11:51 +05:30
crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) = matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
math_inv33 ( matmul ( crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) , &
2020-10-23 21:06:41 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) ) ) )
2019-04-03 16:24:07 +05:30
crystallite_subdt ( c , i , e ) = crystallite_subStep ( c , i , e ) * crystallite_dt ( c , i , e )
crystallite_converged ( c , i , e ) = . false .
2020-09-29 16:43:53 +05:30
call integrateState ( c , i , e )
2019-04-03 16:24:07 +05:30
endif
enddo
enddo
enddo elementLooping3
!$OMP END PARALLEL DO
2020-09-29 17:06:15 +05:30
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
! integrate --- requires fully defined state array (basic + dependent state)
2019-04-11 10:54:04 +05:30
where ( . not . crystallite_converged . and . crystallite_subStep > num % subStepMinCryst ) & ! do not try non-converged but fully cutbacked any further
2020-04-01 12:24:20 +05:30
todo = . true . ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation
2019-01-18 19:12:44 +05:30
2019-04-03 16:24:07 +05:30
enddo cutbackLooping
2019-01-18 19:12:44 +05:30
! return whether converged or not
2019-04-03 16:24:07 +05:30
crystallite_stress = . false .
elementLooping5 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2020-02-13 21:40:27 +05:30
crystallite_stress ( i , e ) = all ( crystallite_converged ( : , i , e ) )
2019-04-03 16:24:07 +05:30
enddo
enddo elementLooping5
2019-01-18 19:12:44 +05:30
end function crystallite_stress
2020-09-28 21:26:48 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-29 12:33:58 +05:30
!> @brief Backup data for homog cutback.
2020-09-28 21:26:48 +05:30
!--------------------------------------------------------------------------------------------------
subroutine crystallite_initializeRestorationPoints ( i , e )
integer , intent ( in ) :: &
i , & !< integration point number
e !< element number
integer :: &
2020-09-28 21:31:43 +05:30
c , & !< constituent number
2020-09-28 21:26:48 +05:30
s
2020-10-28 01:57:26 +05:30
do c = 1 , homogenization_Nconstituents ( material_homogenizationAt ( e ) )
2020-10-08 01:45:13 +05:30
crystallite_partitionedFp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_partitionedLp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Lp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_partitionedFi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_partitionedLi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Li0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_partitionedF0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_F0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_partitionedS0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_S0 ( 1 : 3 , 1 : 3 , c , i , e )
plasticState ( material_phaseAt ( c , e ) ) % partitionedState0 ( : , material_phasememberAt ( c , i , e ) ) = &
2020-09-28 21:26:48 +05:30
plasticState ( material_phaseAt ( c , e ) ) % state0 ( : , material_phasememberAt ( c , i , e ) )
do s = 1 , phase_Nsources ( material_phaseAt ( c , e ) )
2020-10-08 01:45:13 +05:30
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % partitionedState0 ( : , material_phasememberAt ( c , i , e ) ) = &
2020-09-28 21:26:48 +05:30
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % state0 ( : , material_phasememberAt ( c , i , e ) )
enddo
enddo
end subroutine crystallite_initializeRestorationPoints
!--------------------------------------------------------------------------------------------------
2020-09-29 12:33:58 +05:30
!> @brief Wind homog inc forward.
2020-09-28 21:26:48 +05:30
!--------------------------------------------------------------------------------------------------
subroutine crystallite_windForward ( i , e )
integer , intent ( in ) :: &
i , & !< integration point number
e !< element number
integer :: &
2020-09-28 21:31:43 +05:30
c , & !< constituent number
2020-09-28 21:26:48 +05:30
s
2020-10-28 01:57:26 +05:30
do c = 1 , homogenization_Nconstituents ( material_homogenizationAt ( e ) )
2020-10-08 01:45:13 +05:30
crystallite_partitionedF0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partitionedF ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_partitionedFp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_partitionedLp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_partitionedFi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_partitionedLi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_partitionedS0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_S ( 1 : 3 , 1 : 3 , c , i , e )
plasticState ( material_phaseAt ( c , e ) ) % partitionedState0 ( : , material_phasememberAt ( c , i , e ) ) = &
2020-09-28 21:26:48 +05:30
plasticState ( material_phaseAt ( c , e ) ) % state ( : , material_phasememberAt ( c , i , e ) )
do s = 1 , phase_Nsources ( material_phaseAt ( c , e ) )
2020-10-08 01:45:13 +05:30
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % partitionedState0 ( : , material_phasememberAt ( c , i , e ) ) = &
2020-09-28 21:26:48 +05:30
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % state ( : , material_phasememberAt ( c , i , e ) )
enddo
enddo
end subroutine crystallite_windForward
!--------------------------------------------------------------------------------------------------
2020-09-29 12:33:58 +05:30
!> @brief Restore data after homog cutback.
2020-09-28 21:26:48 +05:30
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restore ( i , e , includeL )
integer , intent ( in ) :: &
i , & !< integration point number
e !< element number
logical , intent ( in ) :: &
includeL !< protect agains fake cutback
integer :: &
2020-09-28 21:31:43 +05:30
c , & !< constituent number
2020-09-28 21:26:48 +05:30
s
2020-10-28 01:57:26 +05:30
do c = 1 , homogenization_Nconstituents ( material_homogenizationAt ( e ) )
2020-09-28 21:26:48 +05:30
if ( includeL ) then
2020-10-08 01:45:13 +05:30
crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partitionedLp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partitionedLi0 ( 1 : 3 , 1 : 3 , c , i , e )
2020-09-28 21:26:48 +05:30
endif ! maybe protecting everything from overwriting makes more sense
2020-10-08 01:45:13 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partitionedFp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partitionedFi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partitionedS0 ( 1 : 3 , 1 : 3 , c , i , e )
2020-09-28 21:26:48 +05:30
plasticState ( material_phaseAt ( c , e ) ) % state ( : , material_phasememberAt ( c , i , e ) ) = &
2020-10-08 01:45:13 +05:30
plasticState ( material_phaseAt ( c , e ) ) % partitionedState0 ( : , material_phasememberAt ( c , i , e ) )
2020-09-28 21:26:48 +05:30
do s = 1 , phase_Nsources ( material_phaseAt ( c , e ) )
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % state ( : , material_phasememberAt ( c , i , e ) ) = &
2020-10-08 01:45:13 +05:30
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % partitionedState0 ( : , material_phasememberAt ( c , i , e ) )
2020-09-28 21:26:48 +05:30
enddo
enddo
end subroutine crystallite_restore
2019-01-18 16:46:26 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-29 12:33:58 +05:30
!> @brief Calculate tangent (dPdF).
2019-01-18 16:46:26 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-30 14:23:05 +05:30
function crystallite_stressTangent ( c , i , e ) result ( dPdF )
2019-04-03 16:24:07 +05:30
2020-09-30 14:23:05 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dPdF
integer , intent ( in ) :: &
2020-09-28 21:31:43 +05:30
c , & !< counter in constituent loop
2019-04-03 16:24:07 +05:30
i , & !< counter in integration point loop
2020-09-30 14:23:05 +05:30
e !< counter in element loop
integer :: &
2019-04-03 16:24:07 +05:30
o , &
p
2020-02-21 13:11:08 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: devNull , &
invSubFp0 , invSubFi0 , invFp , invFi , &
temp_33_1 , temp_33_2 , temp_33_3 , temp_33_4
2019-04-03 16:24:07 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dSdFe , &
dSdF , &
dSdFi , &
2020-02-14 10:52:38 +05:30
dLidS , & ! tangent in lattice configuration
2019-04-03 16:24:07 +05:30
dLidFi , &
dLpdS , &
dLpdFi , &
dFidS , &
dFpinvdF , &
rhs_3333 , &
lhs_3333 , &
temp_3333
real ( pReal ) , dimension ( 9 , 9 ) :: temp_99
logical :: error
call constitutive_SandItsTangents ( devNull , dSdFe , dSdFi , &
crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) , &
2020-02-14 10:52:38 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) , c , i , e )
2019-04-03 16:24:07 +05:30
call constitutive_LiAndItsTangents ( devNull , dLidS , dLidFi , &
crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) , &
2020-02-14 10:52:38 +05:30
c , i , e )
2019-04-03 16:24:07 +05:30
2020-02-21 13:11:08 +05:30
invFp = math_inv33 ( crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) )
invFi = math_inv33 ( crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) )
invSubFp0 = math_inv33 ( crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) )
invSubFi0 = math_inv33 ( crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) )
2019-04-03 16:24:07 +05:30
if ( sum ( abs ( dLidS ) ) < tol_math_check ) then
dFidS = 0.0_pReal
else
lhs_3333 = 0.0_pReal ; rhs_3333 = 0.0_pReal
do o = 1 , 3 ; do p = 1 , 3
lhs_3333 ( 1 : 3 , 1 : 3 , o , p ) = lhs_3333 ( 1 : 3 , 1 : 3 , o , p ) &
+ crystallite_subdt ( c , i , e ) * matmul ( invSubFi0 , dLidFi ( 1 : 3 , 1 : 3 , o , p ) )
lhs_3333 ( 1 : 3 , o , 1 : 3 , p ) = lhs_3333 ( 1 : 3 , o , 1 : 3 , p ) &
2020-02-21 13:11:08 +05:30
+ invFi * invFi ( p , o )
2019-04-03 16:24:07 +05:30
rhs_3333 ( 1 : 3 , 1 : 3 , o , p ) = rhs_3333 ( 1 : 3 , 1 : 3 , o , p ) &
- crystallite_subdt ( c , i , e ) * matmul ( invSubFi0 , dLidS ( 1 : 3 , 1 : 3 , o , p ) )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-09-21 06:50:33 +05:30
call math_invert ( temp_99 , error , math_3333to99 ( lhs_3333 ) )
2019-04-03 16:24:07 +05:30
if ( error ) then
call IO_warning ( warning_ID = 600 , el = e , ip = i , g = c , &
ext_msg = 'inversion error in analytic tangent calculation' )
dFidS = 0.0_pReal
else
dFidS = math_mul3333xx3333 ( math_99to3333 ( temp_99 ) , rhs_3333 )
endif
dLidS = math_mul3333xx3333 ( dLidFi , dFidS ) + dLidS
endif
call constitutive_LpAndItsTangents ( devNull , dLpdS , dLpdFi , &
crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) , c , i , e ) ! call constitutive law to calculate Lp tangent in lattice configuration
dLpdS = math_mul3333xx3333 ( dLpdFi , dFidS ) + dLpdS
2019-01-18 16:46:26 +05:30
!--------------------------------------------------------------------------------------------------
! calculate dSdF
2020-02-21 13:11:08 +05:30
temp_33_1 = transpose ( matmul ( invFp , invFi ) )
2020-02-21 13:15:11 +05:30
temp_33_2 = matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , invSubFp0 )
temp_33_3 = matmul ( matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , invFp ) , invSubFi0 )
2019-04-03 16:24:07 +05:30
2019-04-03 16:47:21 +05:30
do o = 1 , 3 ; do p = 1 , 3
2019-04-03 16:24:07 +05:30
rhs_3333 ( p , o , 1 : 3 , 1 : 3 ) = matmul ( dSdFe ( p , o , 1 : 3 , 1 : 3 ) , temp_33_1 )
2020-02-21 13:15:11 +05:30
temp_3333 ( 1 : 3 , 1 : 3 , p , o ) = matmul ( matmul ( temp_33_2 , dLpdS ( 1 : 3 , 1 : 3 , p , o ) ) , invFi ) &
2019-04-03 16:24:07 +05:30
+ matmul ( temp_33_3 , dLidS ( 1 : 3 , 1 : 3 , p , o ) )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-04-03 16:24:07 +05:30
lhs_3333 = crystallite_subdt ( c , i , e ) * math_mul3333xx3333 ( dSdFe , temp_3333 ) &
+ math_mul3333xx3333 ( dSdFi , dFidS )
2020-09-13 14:48:09 +05:30
call math_invert ( temp_99 , error , math_eye ( 9 ) + math_3333to99 ( lhs_3333 ) )
2019-04-03 16:24:07 +05:30
if ( error ) then
call IO_warning ( warning_ID = 600 , el = e , ip = i , g = c , &
ext_msg = 'inversion error in analytic tangent calculation' )
dSdF = rhs_3333
else
dSdF = math_mul3333xx3333 ( math_99to3333 ( temp_99 ) , rhs_3333 )
endif
2019-01-18 16:46:26 +05:30
!--------------------------------------------------------------------------------------------------
! calculate dFpinvdF
2019-04-03 16:24:07 +05:30
temp_3333 = math_mul3333xx3333 ( dLpdS , dSdF )
2019-04-03 16:47:21 +05:30
do o = 1 , 3 ; do p = 1 , 3
2020-02-14 10:52:38 +05:30
dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) = - crystallite_subdt ( c , i , e ) &
2020-02-21 13:15:11 +05:30
* matmul ( invSubFp0 , matmul ( temp_3333 ( 1 : 3 , 1 : 3 , p , o ) , invFi ) )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-01-18 16:46:26 +05:30
!--------------------------------------------------------------------------------------------------
! assemble dPdF
2020-02-21 13:11:08 +05:30
temp_33_1 = matmul ( crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) , transpose ( invFp ) )
temp_33_2 = matmul ( invFp , temp_33_1 )
temp_33_3 = matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , invFp )
2020-02-11 22:11:30 +05:30
temp_33_4 = matmul ( temp_33_3 , crystallite_S ( 1 : 3 , 1 : 3 , c , i , e ) )
2019-04-03 16:24:07 +05:30
2020-09-30 14:23:05 +05:30
dPdF = 0.0_pReal
2019-04-03 16:47:21 +05:30
do p = 1 , 3
2020-09-30 14:23:05 +05:30
dPdF ( p , 1 : 3 , p , 1 : 3 ) = transpose ( temp_33_2 )
2019-04-03 16:24:07 +05:30
enddo
2019-04-03 16:47:21 +05:30
do o = 1 , 3 ; do p = 1 , 3
2020-09-30 14:23:05 +05:30
dPdF ( 1 : 3 , 1 : 3 , p , o ) = dPdF ( 1 : 3 , 1 : 3 , p , o ) &
+ matmul ( matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) , temp_33_1 ) &
+ matmul ( matmul ( temp_33_3 , dSdF ( 1 : 3 , 1 : 3 , p , o ) ) , &
transpose ( invFp ) ) &
+ matmul ( temp_33_4 , transpose ( dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-04-03 16:24:07 +05:30
2020-09-30 14:23:05 +05:30
end function crystallite_stressTangent
2019-01-18 16:46:26 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-19 14:05:45 +05:30
!> @brief calculates orientations
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-19 14:05:45 +05:30
subroutine crystallite_orientations
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
integer &
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e !< counter in element loop
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
!$OMP PARALLEL DO
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2020-10-28 01:57:26 +05:30
do c = 1 , homogenization_Nconstituents ( material_homogenizationAt ( e ) )
2020-03-15 18:51:11 +05:30
call crystallite_orientation ( c , i , e ) % fromMatrix ( transpose ( math_rotationalPart ( crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) ) ) )
2019-04-03 16:24:07 +05:30
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
2020-02-13 21:40:27 +05:30
2020-04-01 21:23:07 +05:30
nonlocalPresent : if ( any ( plasticState % nonlocal ) ) then
2019-04-03 16:24:07 +05:30
!$OMP PARALLEL DO
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-04-01 21:23:07 +05:30
if ( plasticState ( material_phaseAt ( 1 , e ) ) % nonlocal ) then
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2020-03-16 14:19:59 +05:30
call plastic_nonlocal_updateCompatibility ( crystallite_orientation , &
2020-07-17 03:08:37 +05:30
phase_plasticityInstance ( material_phaseAt ( 1 , e ) ) , i , e )
2020-04-01 21:23:07 +05:30
enddo
endif
enddo
2019-04-03 16:24:07 +05:30
!$OMP END PARALLEL DO
endif nonlocalPresent
2014-08-26 20:14:32 +05:30
2019-01-19 14:05:45 +05:30
end subroutine crystallite_orientations
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Map 2nd order tensor to reference config
!--------------------------------------------------------------------------------------------------
function crystallite_push33ToRef ( ipc , ip , el , tensor33 )
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: crystallite_push33ToRef
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: tensor33
real ( pReal ) , dimension ( 3 , 3 ) :: T
integer , intent ( in ) :: &
el , &
ip , &
ipc
2020-02-13 21:40:27 +05:30
2019-09-21 05:48:09 +05:30
T = matmul ( material_orientation0 ( ipc , ip , el ) % asMatrix ( ) , & ! ToDo: initial orientation correct?
2019-09-20 08:12:28 +05:30
transpose ( math_inv33 ( crystallite_subF ( 1 : 3 , 1 : 3 , ipc , ip , el ) ) ) )
2019-04-03 16:24:07 +05:30
crystallite_push33ToRef = matmul ( transpose ( T ) , matmul ( tensor33 , T ) )
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
end function crystallite_push33ToRef
2014-08-26 20:14:32 +05:30
2019-04-06 10:01:02 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-13 13:17:56 +05:30
!> @brief writes crystallite results to HDF5 output file
2019-04-06 10:01:02 +05:30
!--------------------------------------------------------------------------------------------------
subroutine crystallite_results
2019-12-19 00:35:51 +05:30
2019-04-06 15:36:34 +05:30
integer :: p , o
2019-04-07 17:58:08 +05:30
real ( pReal ) , allocatable , dimension ( : , : , : ) :: selected_tensors
type ( rotation ) , allocatable , dimension ( : ) :: selected_rotations
2020-11-06 17:55:00 +05:30
character ( len = : ) , allocatable :: group , structureLabel
2020-02-13 21:40:27 +05:30
2020-08-15 19:32:10 +05:30
do p = 1 , size ( material_name_phase )
2020-11-18 19:22:16 +05:30
group = trim ( 'current/phase' ) / / '/' / / trim ( material_name_phase ( p ) ) / / '/mechanics'
2020-02-13 21:40:27 +05:30
call results_closeGroup ( results_addGroup ( group ) )
2019-04-13 13:17:56 +05:30
2019-04-06 15:36:34 +05:30
do o = 1 , size ( output_constituent ( p ) % label )
select case ( output_constituent ( p ) % label ( o ) )
2020-08-15 19:32:10 +05:30
case ( 'F' )
2020-10-08 01:45:13 +05:30
selected_tensors = select_tensors ( crystallite_partitionedF , p )
2020-08-15 19:32:10 +05:30
call results_writeDataset ( group , selected_tensors , output_constituent ( p ) % label ( o ) , &
2019-04-06 15:36:34 +05:30
'deformation gradient' , '1' )
2020-10-29 18:20:12 +05:30
case ( 'F_e' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Fe , p )
2020-08-15 19:32:10 +05:30
call results_writeDataset ( group , selected_tensors , output_constituent ( p ) % label ( o ) , &
2019-04-06 15:36:34 +05:30
'elastic deformation gradient' , '1' )
2020-10-29 18:20:12 +05:30
case ( 'F_p' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Fp , p )
2020-08-15 19:32:10 +05:30
call results_writeDataset ( group , selected_tensors , output_constituent ( p ) % label ( o ) , &
2019-04-06 15:36:34 +05:30
'plastic deformation gradient' , '1' )
2020-10-29 18:20:12 +05:30
case ( 'F_i' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Fi , p )
2020-08-15 19:32:10 +05:30
call results_writeDataset ( group , selected_tensors , output_constituent ( p ) % label ( o ) , &
2019-04-06 15:36:34 +05:30
'inelastic deformation gradient' , '1' )
2020-10-29 18:20:12 +05:30
case ( 'L_p' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Lp , p )
2020-08-15 19:32:10 +05:30
call results_writeDataset ( group , selected_tensors , output_constituent ( p ) % label ( o ) , &
2019-04-06 15:36:34 +05:30
'plastic velocity gradient' , '1/s' )
2020-10-29 18:20:12 +05:30
case ( 'L_i' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Li , p )
2020-08-15 19:32:10 +05:30
call results_writeDataset ( group , selected_tensors , output_constituent ( p ) % label ( o ) , &
2019-04-06 15:36:34 +05:30
'inelastic velocity gradient' , '1/s' )
2020-08-15 19:32:10 +05:30
case ( 'P' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_P , p )
2020-08-15 19:32:10 +05:30
call results_writeDataset ( group , selected_tensors , output_constituent ( p ) % label ( o ) , &
2020-11-18 19:22:16 +05:30
'First Piola-Kirchhoff stress' , 'Pa' )
2020-08-15 19:32:10 +05:30
case ( 'S' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_S , p )
2020-08-15 19:32:10 +05:30
call results_writeDataset ( group , selected_tensors , output_constituent ( p ) % label ( o ) , &
2020-11-18 19:22:16 +05:30
'Second Piola-Kirchhoff stress' , 'Pa' )
2020-08-15 19:32:10 +05:30
case ( 'O' )
2019-04-07 18:17:21 +05:30
select case ( lattice_structure ( p ) )
2020-03-14 23:41:26 +05:30
case ( lattice_ISO_ID )
structureLabel = 'iso'
case ( lattice_FCC_ID )
structureLabel = 'fcc'
case ( lattice_BCC_ID )
structureLabel = 'bcc'
case ( lattice_BCT_ID )
structureLabel = 'bct'
case ( lattice_HEX_ID )
structureLabel = 'hex'
case ( lattice_ORT_ID )
structureLabel = 'ort'
2019-04-07 18:17:21 +05:30
end select
2019-04-07 17:58:08 +05:30
selected_rotations = select_rotations ( crystallite_orientation , p )
2020-08-15 19:32:10 +05:30
call results_writeDataset ( group , selected_rotations , output_constituent ( p ) % label ( o ) , &
2020-03-14 23:41:26 +05:30
'crystal orientation as quaternion' , structureLabel )
2019-04-06 15:36:34 +05:30
end select
enddo
2019-06-15 19:14:15 +05:30
enddo
2020-02-13 21:40:27 +05:30
2019-06-15 19:14:15 +05:30
contains
2019-04-07 17:58:08 +05:30
2019-06-15 19:14:15 +05:30
!------------------------------------------------------------------------------------------------
!> @brief select tensors for output
!------------------------------------------------------------------------------------------------
2019-04-13 19:04:51 +05:30
function select_tensors ( dataset , instance )
2020-02-13 21:40:27 +05:30
2019-04-13 19:04:51 +05:30
integer , intent ( in ) :: instance
real ( pReal ) , dimension ( : , : , : , : , : ) , intent ( in ) :: dataset
real ( pReal ) , allocatable , dimension ( : , : , : ) :: select_tensors
integer :: e , i , c , j
2020-02-13 21:40:27 +05:30
2020-10-28 01:57:26 +05:30
allocate ( select_tensors ( 3 , 3 , count ( material_phaseAt == instance ) * discretization_nIPs ) )
2019-04-13 19:04:51 +05:30
j = 0
do e = 1 , size ( material_phaseAt , 2 )
2020-10-28 01:57:26 +05:30
do i = 1 , discretization_nIPs
2019-10-19 23:25:00 +05:30
do c = 1 , size ( material_phaseAt , 1 ) !ToDo: this needs to be changed for varying Ngrains
if ( material_phaseAt ( c , e ) == instance ) then
j = j + 1
select_tensors ( 1 : 3 , 1 : 3 , j ) = dataset ( 1 : 3 , 1 : 3 , c , i , e )
endif
2019-04-13 19:04:51 +05:30
enddo
enddo
2019-04-07 17:58:08 +05:30
enddo
2020-02-13 21:40:27 +05:30
2019-04-13 19:04:51 +05:30
end function select_tensors
2020-02-13 21:40:27 +05:30
2019-04-07 17:58:08 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief select rotations for output
2020-02-13 21:40:27 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-13 19:04:51 +05:30
function select_rotations ( dataset , instance )
2020-02-13 21:40:27 +05:30
2019-04-13 19:04:51 +05:30
integer , intent ( in ) :: instance
type ( rotation ) , dimension ( : , : , : ) , intent ( in ) :: dataset
type ( rotation ) , allocatable , dimension ( : ) :: select_rotations
integer :: e , i , c , j
2020-02-13 21:40:27 +05:30
2020-10-28 01:57:26 +05:30
allocate ( select_rotations ( count ( material_phaseAt == instance ) * homogenization_maxNconstituents * discretization_nIPs ) )
2019-04-13 19:04:51 +05:30
j = 0
do e = 1 , size ( material_phaseAt , 2 )
2020-10-28 01:57:26 +05:30
do i = 1 , discretization_nIPs
2019-10-19 23:25:00 +05:30
do c = 1 , size ( material_phaseAt , 1 ) !ToDo: this needs to be changed for varying Ngrains
2019-04-13 19:04:51 +05:30
if ( material_phaseAt ( c , e ) == instance ) then
2019-04-06 15:36:34 +05:30
j = j + 1
2019-04-13 19:04:51 +05:30
select_rotations ( j ) = dataset ( c , i , e )
endif
enddo
enddo
enddo
2020-02-13 21:40:27 +05:30
2019-04-07 17:58:08 +05:30
end function select_rotations
2019-12-19 00:35:51 +05:30
2019-04-06 10:01:02 +05:30
end subroutine crystallite_results
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction
!--------------------------------------------------------------------------------------------------
2020-04-01 12:53:43 +05:30
function integrateStress ( ipc , ip , el , timeFraction ) result ( broken )
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
integer , intent ( in ) :: el , & ! element index
ip , & ! integration point index
ipc ! grain index
real ( pReal ) , optional , intent ( in ) :: timeFraction ! fraction of timestep
2020-02-13 21:40:27 +05:30
2020-02-11 22:20:07 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: F , & ! deformation gradient at end of timestep
2019-04-03 16:24:07 +05:30
Fp_new , & ! plastic deformation gradient at end of timestep
invFp_new , & ! inverse of Fp_new
invFp_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
2020-02-21 13:11:08 +05:30
Fi_new , & ! gradient of intermediate deformation stages
invFi_new , &
invFi_current , & ! inverse of Fi_current
2019-04-03 16:24:07 +05:30
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
2020-02-21 13:11:08 +05:30
Fe , & ! elastic deformation gradient
2019-04-03 16:24:07 +05:30
S , & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
A , &
B , &
temp_33
2020-02-12 10:50:27 +05:30
real ( pReal ) , dimension ( 9 ) :: temp_9 ! needed for matrix inversion by LAPACK
integer , dimension ( 9 ) :: devNull_9 ! needed for matrix inversion by LAPACK
2019-04-03 16:24:07 +05:30
real ( pReal ) , dimension ( 9 , 9 ) :: dRLp_dLp , & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme)
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
2020-02-13 23:13:20 +05:30
real ( pReal ) steplengthLp , &
2019-04-03 16:24:07 +05:30
steplengthLi , &
dt , & ! time increment
2020-03-14 23:41:26 +05:30
atol_Lp , &
atol_Li , &
2020-02-11 22:17:48 +05:30
devNull
2019-04-03 16:24:07 +05:30
integer NiterationStressLp , & ! number of stress integrations
NiterationStressLi , & ! number of inner stress integrations
ierr , & ! error indicator for LAPACK
o , &
p , &
jacoCounterLp , &
jacoCounterLi ! counters to check for Jacobian update
2020-04-01 12:53:43 +05:30
logical :: error , broken
2020-02-13 21:40:27 +05:30
2020-04-01 12:53:43 +05:30
broken = . true .
2014-07-02 17:57:39 +05:30
2019-04-03 16:24:07 +05:30
if ( present ( timeFraction ) ) then
dt = crystallite_subdt ( ipc , ip , el ) * timeFraction
2020-02-11 22:20:07 +05:30
F = 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
2019-04-03 16:24:07 +05:30
else
dt = crystallite_subdt ( ipc , ip , el )
2020-02-11 22:20:07 +05:30
F = crystallite_subF ( 1 : 3 , 1 : 3 , ipc , ip , el )
2019-04-03 16:24:07 +05:30
endif
2014-08-26 20:14:32 +05:30
2020-10-08 01:45:13 +05:30
call constitutive_dependentState ( crystallite_partitionedF ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
2020-04-01 10:48:37 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , ipc , ip , el ) , ipc , ip , el )
2020-02-12 10:50:27 +05:30
Lpguess = crystallite_Lp ( 1 : 3 , 1 : 3 , ipc , ip , el ) ! take as first guess
Liguess = crystallite_Li ( 1 : 3 , 1 : 3 , ipc , ip , el ) ! take as first guess
2020-02-13 21:40:27 +05:30
2020-02-11 22:17:48 +05:30
call math_invert33 ( invFp_current , devNull , error , crystallite_subFp0 ( 1 : 3 , 1 : 3 , ipc , ip , el ) )
2020-02-12 10:50:27 +05:30
if ( error ) return ! error
2020-02-11 22:17:48 +05:30
call math_invert33 ( invFi_current , devNull , error , crystallite_subFi0 ( 1 : 3 , 1 : 3 , ipc , ip , el ) )
2020-02-12 10:50:27 +05:30
if ( error ) return ! error
2020-02-11 22:17:48 +05:30
2020-02-11 22:20:07 +05:30
A = matmul ( F , invFp_current ) ! intermediate tensor needed later to calculate dFe_dLp
2020-02-12 10:32:37 +05:30
jacoCounterLi = 0
steplengthLi = 1.0_pReal
residuumLi_old = 0.0_pReal
Liguess_old = Liguess
2020-02-11 22:17:48 +05:30
2020-02-12 10:32:37 +05:30
NiterationStressLi = 0
LiLoop : do
NiterationStressLi = NiterationStressLi + 1
2020-02-12 10:50:27 +05:30
if ( NiterationStressLi > num % nStress ) return ! error
2020-02-11 22:17:48 +05:30
2019-04-03 16:24:07 +05:30
invFi_new = matmul ( invFi_current , math_I3 - dt * Liguess )
Fi_new = math_inv33 ( invFi_new )
2020-02-13 21:40:27 +05:30
2020-02-12 10:32:37 +05:30
jacoCounterLp = 0
steplengthLp = 1.0_pReal
residuumLp_old = 0.0_pReal
Lpguess_old = Lpguess
2020-02-13 21:40:27 +05:30
2020-02-12 10:32:37 +05:30
NiterationStressLp = 0
LpLoop : do
NiterationStressLp = NiterationStressLp + 1
2020-02-12 10:50:27 +05:30
if ( NiterationStressLp > num % nStress ) return ! error
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
B = math_I3 - dt * Lpguess
Fe = matmul ( matmul ( A , B ) , invFi_new )
call constitutive_SandItsTangents ( S , dS_dFe , dS_dFi , &
2020-02-11 22:17:48 +05:30
Fe , Fi_new , ipc , ip , el )
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
call constitutive_LpAndItsTangents ( Lp_constitutive , dLp_dS , dLp_dFi , &
S , Fi_new , ipc , ip , el )
2014-08-26 20:14:32 +05:30
2019-04-03 16:24:07 +05:30
!* update current residuum and check for convergence of loop
2020-03-14 23:41:26 +05:30
atol_Lp = max ( num % rtol_crystalliteStress * max ( norm2 ( Lpguess ) , norm2 ( Lp_constitutive ) ) , & ! absolute tolerance from largest acceptable relative error
num % atol_crystalliteStress ) ! minimum lower cutoff
2019-04-03 16:24:07 +05:30
residuumLp = Lpguess - Lp_constitutive
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
if ( any ( IEEE_is_NaN ( residuumLp ) ) ) then
2020-02-12 10:50:27 +05:30
return ! error
2020-03-14 23:41:26 +05:30
elseif ( norm2 ( residuumLp ) < atol_Lp ) then ! converged if below absolute tolerance
2020-02-12 10:50:27 +05:30
exit LpLoop
2020-02-11 22:17:48 +05:30
elseif ( NiterationStressLp == 1 . or . norm2 ( residuumLp ) < norm2 ( residuumLp_old ) ) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
2019-04-03 16:24:07 +05:30
residuumLp_old = residuumLp ! ...remember old values and...
Lpguess_old = Lpguess
steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
2019-04-11 10:54:04 +05:30
steplengthLp = num % subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction
2020-02-12 10:50:27 +05:30
Lpguess = Lpguess_old &
+ deltaLp * stepLengthLp
2019-04-03 16:24:07 +05:30
cycle LpLoop
endif
2014-08-26 20:14:32 +05:30
2020-03-31 13:46:14 +05:30
calculateJacobiLi : if ( mod ( jacoCounterLp , num % iJacoLpresiduum ) == 0 ) then
2020-02-12 10:50:27 +05:30
jacoCounterLp = jacoCounterLp + 1
2019-04-03 16:47:21 +05:30
do o = 1 , 3 ; do p = 1 , 3
2020-03-31 13:46:14 +05:30
dFe_dLp ( o , 1 : 3 , p , 1 : 3 ) = - dt * A ( o , p ) * transpose ( invFi_new ) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
2019-04-03 16:47:21 +05:30
enddo ; enddo
2020-09-13 14:48:09 +05:30
dRLp_dLp = math_eye ( 9 ) &
2019-04-03 16:24:07 +05:30
- math_3333to99 ( math_mul3333xx3333 ( math_mul3333xx3333 ( dLp_dS , dS_dFe ) , dFe_dLp ) )
2020-02-12 10:50:27 +05:30
temp_9 = math_33to9 ( residuumLp )
2020-02-13 23:13:20 +05:30
call dgesv ( 9 , 1 , dRLp_dLp , 9 , devNull_9 , temp_9 , 9 , ierr ) ! solve dRLp/dLp * delta Lp = -res for delta Lp
2020-02-12 10:50:27 +05:30
if ( ierr / = 0 ) return ! error
deltaLp = - math_9to33 ( temp_9 )
2020-03-31 13:46:14 +05:30
endif calculateJacobiLi
2014-08-26 20:14:32 +05:30
2020-02-12 10:50:27 +05:30
Lpguess = Lpguess &
+ deltaLp * steplengthLp
2019-04-03 16:24:07 +05:30
enddo LpLoop
2014-09-03 01:41:57 +05:30
2019-04-03 16:24:07 +05:30
call constitutive_LiAndItsTangents ( Li_constitutive , dLi_dS , dLi_dFi , &
S , Fi_new , ipc , ip , el )
2014-09-03 01:41:57 +05:30
2019-04-03 16:24:07 +05:30
!* update current residuum and check for convergence of loop
2020-03-14 23:41:26 +05:30
atol_Li = max ( num % rtol_crystalliteStress * max ( norm2 ( Liguess ) , norm2 ( Li_constitutive ) ) , & ! absolute tolerance from largest acceptable relative error
num % atol_crystalliteStress ) ! minimum lower cutoff
2019-04-03 16:24:07 +05:30
residuumLi = Liguess - Li_constitutive
2020-02-12 10:50:27 +05:30
if ( any ( IEEE_is_NaN ( residuumLi ) ) ) then
return ! error
2020-03-14 23:41:26 +05:30
elseif ( norm2 ( residuumLi ) < atol_Li ) then ! converged if below absolute tolerance
2020-02-12 10:50:27 +05:30
exit LiLoop
2020-02-11 22:17:48 +05:30
elseif ( NiterationStressLi == 1 . or . norm2 ( residuumLi ) < norm2 ( residuumLi_old ) ) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
2019-04-03 16:24:07 +05:30
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...
2020-02-12 10:50:27 +05:30
steplengthLi = num % subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction
Liguess = Liguess_old &
+ deltaLi * steplengthLi
2019-04-03 16:24:07 +05:30
cycle LiLoop
endif
2020-02-13 21:40:27 +05:30
2020-03-31 13:46:14 +05:30
calculateJacobiLp : if ( mod ( jacoCounterLi , num % iJacoLpresiduum ) == 0 ) then
2020-02-12 10:50:27 +05:30
jacoCounterLi = jacoCounterLi + 1
2020-02-13 21:40:27 +05:30
2020-02-12 10:50:27 +05:30
temp_33 = matmul ( matmul ( A , B ) , invFi_current )
2019-04-03 16:47:21 +05:30
do o = 1 , 3 ; do p = 1 , 3
2019-04-03 16:24:07 +05:30
dFe_dLi ( 1 : 3 , o , 1 : 3 , p ) = - dt * math_I3 ( o , p ) * temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
dFi_dLi ( 1 : 3 , o , 1 : 3 , p ) = - dt * math_I3 ( o , p ) * invFi_current
2019-04-03 16:47:21 +05:30
enddo ; enddo
do o = 1 , 3 ; do p = 1 , 3
2019-04-03 16:24:07 +05:30
dFi_dLi ( 1 : 3 , 1 : 3 , o , p ) = matmul ( matmul ( Fi_new , dFi_dLi ( 1 : 3 , 1 : 3 , o , p ) ) , Fi_new )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2020-09-13 14:48:09 +05:30
dRLi_dLi = math_eye ( 9 ) &
2020-02-11 22:17:48 +05:30
- math_3333to99 ( math_mul3333xx3333 ( dLi_dS , math_mul3333xx3333 ( dS_dFe , dFe_dLi ) &
+ math_mul3333xx3333 ( dS_dFi , dFi_dLi ) ) ) &
2019-04-03 16:24:07 +05:30
- math_3333to99 ( math_mul3333xx3333 ( dLi_dFi , dFi_dLi ) )
2020-02-12 10:50:27 +05:30
temp_9 = math_33to9 ( residuumLi )
call dgesv ( 9 , 1 , dRLi_dLi , 9 , devNull_9 , temp_9 , 9 , ierr ) ! solve dRLi/dLp * delta Li = -res for delta Li
if ( ierr / = 0 ) return ! error
deltaLi = - math_9to33 ( temp_9 )
2020-03-31 13:46:14 +05:30
endif calculateJacobiLp
2020-02-13 21:40:27 +05:30
2020-02-12 10:50:27 +05:30
Liguess = Liguess &
+ deltaLi * steplengthLi
2019-04-03 16:24:07 +05:30
enddo LiLoop
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
invFp_new = matmul ( invFp_current , B )
2020-02-11 22:17:48 +05:30
call math_invert33 ( Fp_new , devNull , error , invFp_new )
2020-02-12 10:50:27 +05:30
if ( error ) return ! error
2014-08-26 20:14:32 +05:30
2020-02-21 13:11:08 +05:30
crystallite_P ( 1 : 3 , 1 : 3 , ipc , ip , el ) = matmul ( matmul ( F , invFp_new ) , matmul ( S , transpose ( invFp_new ) ) )
2019-04-03 16:24:07 +05:30
crystallite_S ( 1 : 3 , 1 : 3 , ipc , ip , el ) = S
crystallite_Lp ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Lpguess
crystallite_Li ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Liguess
2020-03-31 13:42:25 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Fp_new / math_det33 ( Fp_new ) ** ( 1.0_pReal / 3.0_pReal ) ! regularize
2019-04-03 16:24:07 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Fi_new
2020-03-31 13:42:25 +05:30
crystallite_Fe ( 1 : 3 , 1 : 3 , ipc , ip , el ) = matmul ( matmul ( F , invFp_new ) , invFi_new )
2020-04-01 12:53:43 +05:30
broken = . false .
2020-02-12 10:50:27 +05:30
2019-01-19 14:05:45 +05:30
end function integrateStress
2010-10-01 17:48:49 +05:30
2012-05-17 20:55:21 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-29 16:43:53 +05:30
subroutine integrateStateFPI ( g , i , e )
2019-01-15 08:57:57 +05:30
2020-09-29 17:06:15 +05:30
integer , intent ( in ) :: &
2020-03-23 11:09:17 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
2020-09-29 17:06:15 +05:30
g !< grain index in grain loop
integer :: &
NiterationState , & !< number of iterations in state loop
2020-03-23 11:09:17 +05:30
p , &
c , &
s , &
2020-04-10 22:23:59 +05:30
size_pl
integer , dimension ( maxval ( phase_Nsources ) ) :: &
size_so
2020-03-23 11:09:17 +05:30
real ( pReal ) :: &
zeta
2020-03-31 12:34:08 +05:30
real ( pReal ) , dimension ( max ( constitutive_plasticity_maxSizeDotState , constitutive_source_maxSizeDotState ) ) :: &
r ! state residuum
2020-04-10 23:06:29 +05:30
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState , 2 ) :: &
plastic_dotState
2020-04-01 00:27:09 +05:30
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , 2 , maxval ( phase_Nsources ) ) :: source_dotState
2020-03-23 11:09:17 +05:30
logical :: &
2020-09-29 16:26:12 +05:30
broken
2020-09-29 17:06:15 +05:30
p = material_phaseAt ( g , e )
c = material_phaseMemberAt ( g , i , e )
broken = constitutive_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedF0 , &
2020-09-29 17:06:15 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedFp0 , &
2020-09-29 17:06:15 +05:30
crystallite_subdt ( g , i , e ) , g , i , e , p , c )
if ( broken ) return
size_pl = plasticState ( p ) % sizeDotState
plasticState ( p ) % state ( 1 : size_pl , c ) = plasticState ( p ) % subState0 ( 1 : size_pl , c ) &
+ plasticState ( p ) % dotState ( 1 : size_pl , c ) &
* crystallite_subdt ( g , i , e )
plastic_dotState ( 1 : size_pl , 2 ) = 0.0_pReal
do s = 1 , phase_Nsources ( p )
size_so ( s ) = sourceState ( p ) % p ( s ) % sizeDotState
sourceState ( p ) % p ( s ) % state ( 1 : size_so ( s ) , c ) = sourceState ( p ) % p ( s ) % subState0 ( 1 : size_so ( s ) , c ) &
+ sourceState ( p ) % p ( s ) % dotState ( 1 : size_so ( s ) , c ) &
* crystallite_subdt ( g , i , e )
source_dotState ( 1 : size_so ( s ) , 2 , s ) = 0.0_pReal
enddo
2014-08-26 20:14:32 +05:30
2020-09-29 17:06:15 +05:30
iteration : do NiterationState = 1 , num % nState
2020-03-23 12:45:33 +05:30
2020-09-29 17:06:15 +05:30
if ( nIterationState > 1 ) plastic_dotState ( 1 : size_pl , 2 ) = plastic_dotState ( 1 : size_pl , 1 )
plastic_dotState ( 1 : size_pl , 1 ) = plasticState ( p ) % dotState ( : , c )
do s = 1 , phase_Nsources ( p )
if ( nIterationState > 1 ) source_dotState ( 1 : size_so ( s ) , 2 , s ) = source_dotState ( 1 : size_so ( s ) , 1 , s )
source_dotState ( 1 : size_so ( s ) , 1 , s ) = sourceState ( p ) % p ( s ) % dotState ( : , c )
enddo
2020-03-24 11:04:42 +05:30
2020-09-29 17:06:15 +05:30
broken = integrateStress ( g , i , e )
if ( broken ) exit iteration
broken = constitutive_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedF0 , &
2020-09-29 17:06:15 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedFp0 , &
2020-09-29 17:06:15 +05:30
crystallite_subdt ( g , i , e ) , g , i , e , p , c )
if ( broken ) exit iteration
zeta = damper ( plasticState ( p ) % dotState ( : , c ) , plastic_dotState ( 1 : size_pl , 1 ) , &
plastic_dotState ( 1 : size_pl , 2 ) )
plasticState ( p ) % dotState ( : , c ) = plasticState ( p ) % dotState ( : , c ) * zeta &
+ plastic_dotState ( 1 : size_pl , 1 ) * ( 1.0_pReal - zeta )
r ( 1 : size_pl ) = plasticState ( p ) % state ( 1 : size_pl , c ) &
- plasticState ( p ) % subState0 ( 1 : size_pl , c ) &
- plasticState ( p ) % dotState ( 1 : size_pl , c ) * crystallite_subdt ( g , i , e )
plasticState ( p ) % state ( 1 : size_pl , c ) = plasticState ( p ) % state ( 1 : size_pl , c ) &
- r ( 1 : size_pl )
crystallite_converged ( g , i , e ) = converged ( r ( 1 : size_pl ) , &
plasticState ( p ) % state ( 1 : size_pl , c ) , &
plasticState ( p ) % atol ( 1 : size_pl ) )
do s = 1 , phase_Nsources ( p )
zeta = damper ( sourceState ( p ) % p ( s ) % dotState ( : , c ) , &
source_dotState ( 1 : size_so ( s ) , 1 , s ) , &
source_dotState ( 1 : size_so ( s ) , 2 , s ) )
sourceState ( p ) % p ( s ) % dotState ( : , c ) = sourceState ( p ) % p ( s ) % dotState ( : , c ) * zeta &
+ source_dotState ( 1 : size_so ( s ) , 1 , s ) * ( 1.0_pReal - zeta )
r ( 1 : size_so ( s ) ) = sourceState ( p ) % p ( s ) % state ( 1 : size_so ( s ) , c ) &
- sourceState ( p ) % p ( s ) % subState0 ( 1 : size_so ( s ) , c ) &
- sourceState ( p ) % p ( s ) % dotState ( 1 : size_so ( s ) , c ) * crystallite_subdt ( g , i , e )
sourceState ( p ) % p ( s ) % state ( 1 : size_so ( s ) , c ) = sourceState ( p ) % p ( s ) % state ( 1 : size_so ( s ) , c ) &
- r ( 1 : size_so ( s ) )
crystallite_converged ( g , i , e ) = &
crystallite_converged ( g , i , e ) . and . converged ( r ( 1 : size_so ( s ) ) , &
sourceState ( p ) % p ( s ) % state ( 1 : size_so ( s ) , c ) , &
sourceState ( p ) % p ( s ) % atol ( 1 : size_so ( s ) ) )
enddo
2020-03-24 14:47:18 +05:30
2020-09-29 17:06:15 +05:30
if ( crystallite_converged ( g , i , e ) ) then
broken = constitutive_deltaState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , g , i , e , p , c )
exit iteration
endif
enddo iteration
2019-01-28 16:19:24 +05:30
2020-03-23 11:09:17 +05:30
contains
2019-01-28 15:58:46 +05:30
2020-03-23 11:09:17 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate the damping for correction of state and dot state
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function damper ( current , previous , previous2 )
2020-02-13 21:40:27 +05:30
2020-03-23 11:09:17 +05:30
real ( pReal ) , dimension ( : ) , intent ( in ) :: &
current , previous , previous2
2020-02-13 21:40:27 +05:30
2020-03-23 11:09:17 +05:30
real ( pReal ) :: dot_prod12 , dot_prod22
2020-02-13 21:40:27 +05:30
2020-03-23 11:09:17 +05:30
dot_prod12 = dot_product ( current - previous , previous - previous2 )
dot_prod22 = dot_product ( previous - previous2 , previous - previous2 )
if ( ( dot_product ( current , previous ) < 0.0_pReal . or . dot_prod12 < 0.0_pReal ) . and . dot_prod22 > 0.0_pReal ) then
damper = 0.75_pReal + 0.25_pReal * tanh ( 2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22 )
else
damper = 1.0_pReal
endif
2020-02-13 21:40:27 +05:30
2020-03-23 11:09:17 +05:30
end function damper
2019-01-28 15:58:46 +05:30
2019-01-15 08:57:57 +05:30
end subroutine integrateStateFPI
2010-10-01 17:48:49 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-29 09:41:29 +05:30
!> @brief integrate state with 1st order explicit Euler method
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-29 16:43:53 +05:30
subroutine integrateStateEuler ( g , i , e )
2014-08-26 20:14:32 +05:30
2020-09-29 17:06:15 +05:30
integer , intent ( in ) :: &
2020-03-24 15:57:53 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
2020-09-29 17:06:15 +05:30
g !< grain index in grain loop
integer :: &
2020-03-24 15:57:53 +05:30
p , &
c , &
s , &
sizeDotState
logical :: &
2020-09-29 16:26:12 +05:30
broken
2020-03-24 15:57:53 +05:30
2020-09-29 17:06:15 +05:30
p = material_phaseAt ( g , e )
c = material_phaseMemberAt ( g , i , e )
broken = constitutive_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedF0 , &
2020-09-29 17:06:15 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedFp0 , &
2020-09-29 17:06:15 +05:30
crystallite_subdt ( g , i , e ) , g , i , e , p , c )
if ( broken ) return
sizeDotState = plasticState ( p ) % sizeDotState
plasticState ( p ) % state ( 1 : sizeDotState , c ) = plasticState ( p ) % subState0 ( 1 : sizeDotState , c ) &
+ plasticState ( p ) % dotState ( 1 : sizeDotState , c ) &
* crystallite_subdt ( g , i , e )
do s = 1 , phase_Nsources ( p )
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) = sourceState ( p ) % p ( s ) % subState0 ( 1 : sizeDotState , c ) &
+ sourceState ( p ) % p ( s ) % dotState ( 1 : sizeDotState , c ) &
* crystallite_subdt ( g , i , e )
enddo
2020-03-24 15:57:53 +05:30
2020-09-29 17:06:15 +05:30
broken = constitutive_deltaState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , g , i , e , p , c )
if ( broken ) return
2020-03-24 16:07:00 +05:30
2020-09-29 17:06:15 +05:30
broken = integrateStress ( g , i , e )
crystallite_converged ( g , i , e ) = . not . broken
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
!--------------------------------------------------------------------------------------------------
2020-09-29 16:43:53 +05:30
subroutine integrateStateAdaptiveEuler ( g , i , e )
2014-07-02 17:57:39 +05:30
2020-09-29 17:06:15 +05:30
integer , intent ( in ) :: &
e , & !< element index in element loop
i , & !< integration point index in ip loop
g !< grain index in grain loop
2020-03-24 17:00:43 +05:30
integer :: &
p , &
c , &
s , &
sizeDotState
logical :: &
2020-09-29 16:26:12 +05:30
broken
2020-03-24 17:00:43 +05:30
2020-04-01 10:32:23 +05:30
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState ) :: residuum_plastic
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , maxval ( phase_Nsources ) ) :: residuum_source
2020-03-24 17:00:43 +05:30
2020-09-29 17:06:15 +05:30
p = material_phaseAt ( g , e )
c = material_phaseMemberAt ( g , i , e )
broken = constitutive_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedF0 , &
2020-09-29 17:06:15 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedFp0 , &
2020-09-29 17:06:15 +05:30
crystallite_subdt ( g , i , e ) , g , i , e , p , c )
if ( broken ) return
sizeDotState = plasticState ( p ) % sizeDotState
residuum_plastic ( 1 : sizeDotState ) = - plasticState ( p ) % dotstate ( 1 : sizeDotState , c ) * 0.5_pReal * crystallite_subdt ( g , i , e )
plasticState ( p ) % state ( 1 : sizeDotState , c ) = plasticState ( p ) % subState0 ( 1 : sizeDotState , c ) &
+ plasticState ( p ) % dotstate ( 1 : sizeDotState , c ) * crystallite_subdt ( g , i , e )
do s = 1 , phase_Nsources ( p )
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
residuum_source ( 1 : sizeDotState , s ) = - sourceState ( p ) % p ( s ) % dotstate ( 1 : sizeDotState , c ) &
* 0.5_pReal * crystallite_subdt ( g , i , e )
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) = sourceState ( p ) % p ( s ) % subState0 ( 1 : sizeDotState , c ) &
+ sourceState ( p ) % p ( s ) % dotstate ( 1 : sizeDotState , c ) * crystallite_subdt ( g , i , e )
enddo
broken = constitutive_deltaState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , g , i , e , p , c )
if ( broken ) return
broken = integrateStress ( g , i , e )
if ( broken ) return
broken = constitutive_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedF0 , &
2020-09-29 17:06:15 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedFp0 , &
2020-09-29 17:06:15 +05:30
crystallite_subdt ( g , i , e ) , g , i , e , p , c )
if ( broken ) return
sizeDotState = plasticState ( p ) % sizeDotState
crystallite_converged ( g , i , e ) = converged ( residuum_plastic ( 1 : sizeDotState ) &
+ 0.5_pReal * plasticState ( p ) % dotState ( : , c ) * crystallite_subdt ( g , i , e ) , &
plasticState ( p ) % state ( 1 : sizeDotState , c ) , &
plasticState ( p ) % atol ( 1 : sizeDotState ) )
do s = 1 , phase_Nsources ( p )
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
crystallite_converged ( g , i , e ) = &
crystallite_converged ( g , i , e ) . and . converged ( residuum_source ( 1 : sizeDotState , s ) &
+ 0.5_pReal * sourceState ( p ) % p ( s ) % dotState ( : , c ) * crystallite_subdt ( g , i , e ) , &
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) , &
sourceState ( p ) % p ( s ) % atol ( 1 : sizeDotState ) )
enddo
2020-02-13 21:40:27 +05:30
2019-01-15 08:57:57 +05:30
end subroutine integrateStateAdaptiveEuler
2010-10-01 17:48:49 +05:30
2020-04-17 13:47:00 +05:30
!---------------------------------------------------------------------------------------------------
2020-04-17 15:21:39 +05:30
!> @brief Integrate state (including stress integration) with the classic Runge Kutta method
2020-04-17 13:47:00 +05:30
!---------------------------------------------------------------------------------------------------
2020-09-29 16:43:53 +05:30
subroutine integrateStateRK4 ( g , i , e )
2020-04-01 12:24:20 +05:30
2020-09-29 17:06:15 +05:30
integer , intent ( in ) :: g , i , e
2014-09-10 14:07:12 +05:30
2020-04-01 15:00:01 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , parameter :: &
A = reshape ( [ &
2020-03-27 02:16:28 +05:30
0.5_pReal , 0.0_pReal , 0.0_pReal , &
0.0_pReal , 0.5_pReal , 0.0_pReal , &
2020-04-17 15:28:03 +05:30
0.0_pReal , 0.0_pReal , 1.0_pReal ] , &
2020-04-17 15:21:39 +05:30
shape ( A ) )
2020-04-01 15:00:01 +05:30
real ( pReal ) , dimension ( 3 ) , parameter :: &
2020-04-17 13:47:00 +05:30
C = [ 0.5_pReal , 0.5_pReal , 1.0_pReal ]
2020-04-01 15:00:01 +05:30
real ( pReal ) , dimension ( 4 ) , parameter :: &
2020-04-17 13:47:00 +05:30
B = [ 1.0_pReal / 6.0_pReal , 1.0_pReal / 3.0_pReal , 1.0_pReal / 3.0_pReal , 1.0_pReal / 6.0_pReal ]
2020-03-27 02:16:28 +05:30
2020-09-29 16:43:53 +05:30
call integrateStateRK ( g , i , e , A , B , C )
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
2020-04-17 13:47:00 +05:30
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the Cash-Carp method
!---------------------------------------------------------------------------------------------------
2020-09-29 16:43:53 +05:30
subroutine integrateStateRKCK45 ( g , i , e )
2020-04-01 12:24:20 +05:30
2020-09-29 17:06:15 +05:30
integer , intent ( in ) :: g , i , e
2014-11-12 22:10:50 +05:30
2020-03-25 15:03:41 +05:30
real ( pReal ) , dimension ( 5 , 5 ) , parameter :: &
A = reshape ( [ &
2020-04-17 15:21:39 +05:30
1._pReal / 5._pReal , . 0_pReal , . 0_pReal , . 0_pReal , . 0_pReal , &
3._pReal / 4 0._pReal , 9._pReal / 4 0._pReal , . 0_pReal , . 0_pReal , . 0_pReal , &
3_pReal / 1 0._pReal , - 9._pReal / 1 0._pReal , 6._pReal / 5._pReal , . 0_pReal , . 0_pReal , &
- 1 1._pReal / 5 4._pReal , 5._pReal / 2._pReal , - 7 0.0_pReal / 2 7.0_pReal , 3 5.0_pReal / 2 7.0_pReal , . 0_pReal , &
163 1._pReal / 5529 6._pReal , 17 5._pReal / 51 2._pReal , 57 5._pReal / 1382 4._pReal , 4427 5._pReal / 11059 2._pReal , 25 3._pReal / 409 6._pReal ] , &
shape ( A ) )
real ( pReal ) , dimension ( 5 ) , parameter :: &
C = [ 0.2_pReal , 0.3_pReal , 0.6_pReal , 1.0_pReal , 0.875_pReal ]
2020-03-25 15:03:41 +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 , &
2020-04-17 13:47:00 +05:30
12 5.0_pReal / 59 4.0_pReal , . 0_pReal , 51 2.0_pReal / 177 1.0_pReal ] , &
2020-03-25 15:03:41 +05:30
DB = B - &
[ 282 5.0_pReal / 2764 8.0_pReal , . 0_pReal , 1857 5.0_pReal / 4838 4.0_pReal , &
2020-04-17 15:21:39 +05:30
1352 5.0_pReal / 5529 6.0_pReal , 27 7.0_pReal / 1433 6.0_pReal , 1._pReal / 4._pReal ]
2020-04-17 11:54:35 +05:30
2020-09-29 16:43:53 +05:30
call integrateStateRK ( g , i , e , A , B , C , DB )
2020-04-17 11:54:35 +05:30
end subroutine integrateStateRKCK45
!--------------------------------------------------------------------------------------------------
2020-04-17 15:21:39 +05:30
!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
!! embedded explicit Runge-Kutta method
2020-04-17 11:54:35 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-29 16:43:53 +05:30
subroutine integrateStateRK ( g , i , e , A , B , CC , DB )
2020-04-17 11:54:35 +05:30
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: A
real ( pReal ) , dimension ( : ) , intent ( in ) :: B , CC
real ( pReal ) , dimension ( : ) , intent ( in ) , optional :: DB
2020-09-29 17:06:15 +05:30
integer , intent ( in ) :: &
e , & !< element index in element loop
i , & !< integration point index in ip loop
g !< grain index in grain loop
2020-03-25 15:03:41 +05:30
integer :: &
stage , & ! stage index in integration stage loop
n , &
p , &
c , &
s , &
sizeDotState
logical :: &
2020-09-29 16:26:12 +05:30
broken
2020-04-17 11:47:32 +05:30
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , size ( B ) , maxval ( phase_Nsources ) ) :: source_RKdotState
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState , size ( B ) ) :: plastic_RKdotState
2019-01-30 13:26:16 +05:30
2020-09-29 17:06:15 +05:30
p = material_phaseAt ( g , e )
c = material_phaseMemberAt ( g , i , e )
broken = constitutive_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedF0 , &
2020-09-29 17:06:15 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedFp0 , &
2020-09-29 17:06:15 +05:30
crystallite_subdt ( g , i , e ) , g , i , e , p , c )
if ( broken ) return
do stage = 1 , size ( A , 1 )
sizeDotState = plasticState ( p ) % sizeDotState
plastic_RKdotState ( 1 : sizeDotState , stage ) = plasticState ( p ) % dotState ( : , c )
plasticState ( p ) % dotState ( : , c ) = A ( 1 , stage ) * plastic_RKdotState ( 1 : sizeDotState , 1 )
do s = 1 , phase_Nsources ( p )
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
source_RKdotState ( 1 : sizeDotState , stage , s ) = sourceState ( p ) % p ( s ) % dotState ( : , c )
sourceState ( p ) % p ( s ) % dotState ( : , c ) = A ( 1 , stage ) * source_RKdotState ( 1 : sizeDotState , 1 , s )
enddo
2020-02-13 21:40:27 +05:30
2020-09-29 17:06:15 +05:30
do n = 2 , stage
sizeDotState = plasticState ( p ) % sizeDotState
plasticState ( p ) % dotState ( : , c ) = plasticState ( p ) % dotState ( : , c ) &
+ A ( n , stage ) * plastic_RKdotState ( 1 : sizeDotState , n )
do s = 1 , phase_Nsources ( p )
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
sourceState ( p ) % p ( s ) % dotState ( : , c ) = sourceState ( p ) % p ( s ) % dotState ( : , c ) &
+ A ( n , stage ) * source_RKdotState ( 1 : sizeDotState , n , s )
enddo
enddo
2020-02-13 21:40:27 +05:30
2020-09-29 17:06:15 +05:30
sizeDotState = plasticState ( p ) % sizeDotState
plasticState ( p ) % state ( 1 : sizeDotState , c ) = plasticState ( p ) % subState0 ( 1 : sizeDotState , c ) &
+ plasticState ( p ) % dotState ( 1 : sizeDotState , c ) &
* crystallite_subdt ( g , i , e )
do s = 1 , phase_Nsources ( p )
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) = sourceState ( p ) % p ( s ) % subState0 ( 1 : sizeDotState , c ) &
+ sourceState ( p ) % p ( s ) % dotState ( 1 : sizeDotState , c ) &
* crystallite_subdt ( g , i , e )
enddo
2020-02-13 21:40:27 +05:30
2020-09-29 17:06:15 +05:30
broken = integrateStress ( g , i , e , CC ( stage ) )
if ( broken ) exit
2020-02-13 21:40:27 +05:30
2020-09-29 17:06:15 +05:30
broken = constitutive_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedF0 , &
2020-09-29 17:06:15 +05:30
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
2020-10-08 01:45:13 +05:30
crystallite_partitionedFp0 , &
2020-09-29 17:06:15 +05:30
crystallite_subdt ( g , i , e ) * CC ( stage ) , g , i , e , p , c )
if ( broken ) exit
2020-02-13 21:40:27 +05:30
2020-09-29 17:06:15 +05:30
enddo
if ( broken ) return
sizeDotState = plasticState ( p ) % sizeDotState
plastic_RKdotState ( 1 : sizeDotState , size ( B ) ) = plasticState ( p ) % dotState ( : , c )
plasticState ( p ) % dotState ( : , c ) = matmul ( plastic_RKdotState ( 1 : sizeDotState , 1 : size ( B ) ) , B )
plasticState ( p ) % state ( 1 : sizeDotState , c ) = plasticState ( p ) % subState0 ( 1 : sizeDotState , c ) &
+ plasticState ( p ) % dotState ( 1 : sizeDotState , c ) &
* crystallite_subdt ( g , i , e )
if ( present ( DB ) ) &
broken = . not . converged ( matmul ( plastic_RKdotState ( 1 : sizeDotState , 1 : size ( DB ) ) , DB ) &
* crystallite_subdt ( g , i , e ) , &
plasticState ( p ) % state ( 1 : sizeDotState , c ) , &
plasticState ( p ) % atol ( 1 : sizeDotState ) )
do s = 1 , phase_Nsources ( p )
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
source_RKdotState ( 1 : sizeDotState , size ( B ) , s ) = sourceState ( p ) % p ( s ) % dotState ( : , c )
sourceState ( p ) % p ( s ) % dotState ( : , c ) = matmul ( source_RKdotState ( 1 : sizeDotState , 1 : size ( B ) , s ) , B )
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) = sourceState ( p ) % p ( s ) % subState0 ( 1 : sizeDotState , c ) &
+ sourceState ( p ) % p ( s ) % dotState ( 1 : sizeDotState , c ) &
* crystallite_subdt ( g , i , e )
if ( present ( DB ) ) &
broken = broken . or . . not . converged ( matmul ( source_RKdotState ( 1 : sizeDotState , 1 : size ( DB ) , s ) , DB ) &
* crystallite_subdt ( g , i , e ) , &
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) , &
sourceState ( p ) % p ( s ) % atol ( 1 : sizeDotState ) )
enddo
if ( broken ) return
2014-08-26 20:14:32 +05:30
2020-09-29 17:06:15 +05:30
broken = constitutive_deltaState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , g , i , e , p , c )
if ( broken ) return
2020-02-13 21:40:27 +05:30
2020-09-29 17:06:15 +05:30
broken = integrateStress ( g , i , e )
crystallite_converged ( g , i , e ) = . not . broken
2020-02-13 21:40:27 +05:30
2020-04-17 11:54:35 +05:30
end subroutine integrateStateRK
2011-05-11 22:08:45 +05:30
2014-08-26 20:14:32 +05:30
2020-03-27 02:32:28 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief determines whether a point is converged
!--------------------------------------------------------------------------------------------------
logical pure function converged ( residuum , state , atol )
2020-02-13 21:40:27 +05:30
2019-01-31 13:44:02 +05:30
real ( pReal ) , intent ( in ) , dimension ( : ) :: &
2020-03-14 23:41:26 +05:30
residuum , state , atol
2019-04-11 10:54:04 +05:30
real ( pReal ) :: &
rTol
rTol = num % rTol_crystalliteState
2019-01-31 13:44:02 +05:30
2020-03-14 23:41:26 +05:30
converged = all ( abs ( residuum ) < = max ( atol , rtol * abs ( state ) ) )
2019-01-31 13:44:02 +05:30
2020-03-27 02:32:28 +05:30
end function converged
2019-01-23 22:34:19 +05:30
2019-01-24 22:29:38 +05:30
2020-02-25 14:10:38 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Write current restart information (Field and constitutive data) to file.
2020-02-25 22:23:15 +05:30
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
2020-02-25 14:10:38 +05:30
!--------------------------------------------------------------------------------------------------
2020-02-25 22:23:15 +05:30
subroutine crystallite_restartWrite
2020-02-25 14:10:38 +05:30
integer :: i
integer ( HID_T ) :: fileHandle , groupHandle
character ( len = pStringLen ) :: fileName , datasetName
2020-09-22 16:39:12 +05:30
print * , ' writing field and constitutive data required for restart to file' ; flush ( IO_STDOUT )
2020-02-25 14:10:38 +05:30
write ( fileName , '(a,i0,a)' ) trim ( getSolverJobName ( ) ) / / '_' , worldrank , '.hdf5'
fileHandle = HDF5_openFile ( fileName , 'a' )
2020-10-08 01:45:13 +05:30
call HDF5_write ( fileHandle , crystallite_partitionedF , 'F' )
2020-10-30 18:06:53 +05:30
call HDF5_write ( fileHandle , crystallite_Fp , 'F_p' )
call HDF5_write ( fileHandle , crystallite_Fi , 'F_i' )
call HDF5_write ( fileHandle , crystallite_Lp , 'L_p' )
call HDF5_write ( fileHandle , crystallite_Li , 'L_i' )
call HDF5_write ( fileHandle , crystallite_S , 'S' )
2020-02-25 14:10:38 +05:30
2020-11-18 19:22:16 +05:30
groupHandle = HDF5_addGroup ( fileHandle , 'phase' )
2020-08-15 19:32:10 +05:30
do i = 1 , size ( material_name_phase )
2020-11-18 19:22:16 +05:30
write ( datasetName , '(i0,a)' ) i , '_omega'
2020-02-25 14:10:38 +05:30
call HDF5_write ( groupHandle , plasticState ( i ) % state , datasetName )
enddo
call HDF5_closeGroup ( groupHandle )
2020-11-18 19:22:16 +05:30
groupHandle = HDF5_addGroup ( fileHandle , 'homogenization' )
2020-10-24 18:26:03 +05:30
do i = 1 , size ( material_name_homogenization )
2020-11-18 19:22:16 +05:30
write ( datasetName , '(i0,a)' ) i , '_omega'
2020-02-25 14:10:38 +05:30
call HDF5_write ( groupHandle , homogState ( i ) % state , datasetName )
enddo
call HDF5_closeGroup ( groupHandle )
call HDF5_closeFile ( fileHandle )
2020-02-25 22:23:15 +05:30
end subroutine crystallite_restartWrite
!--------------------------------------------------------------------------------------------------
!> @brief Read data for restart
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restartRead
integer :: i
integer ( HID_T ) :: fileHandle , groupHandle
character ( len = pStringLen ) :: fileName , datasetName
2020-09-17 22:58:41 +05:30
print '(/,a,i0,a)' , ' reading restart information of increment from file'
2020-02-25 22:23:15 +05:30
write ( fileName , '(a,i0,a)' ) trim ( getSolverJobName ( ) ) / / '_' , worldrank , '.hdf5'
fileHandle = HDF5_openFile ( fileName )
call HDF5_read ( fileHandle , crystallite_F0 , 'F' )
2020-10-30 18:06:53 +05:30
call HDF5_read ( fileHandle , crystallite_Fp0 , 'F_p' )
call HDF5_read ( fileHandle , crystallite_Fi0 , 'F_i' )
call HDF5_read ( fileHandle , crystallite_Lp0 , 'L_p' )
call HDF5_read ( fileHandle , crystallite_Li0 , 'L_i' )
2020-02-25 22:23:15 +05:30
call HDF5_read ( fileHandle , crystallite_S0 , 'S' )
2020-11-18 19:22:16 +05:30
groupHandle = HDF5_openGroup ( fileHandle , 'phase' )
2020-08-15 19:32:10 +05:30
do i = 1 , size ( material_name_phase )
2020-11-18 19:22:16 +05:30
write ( datasetName , '(i0,a)' ) i , '_omega'
2020-02-25 22:23:15 +05:30
call HDF5_read ( groupHandle , plasticState ( i ) % state0 , datasetName )
enddo
call HDF5_closeGroup ( groupHandle )
2020-11-18 19:22:16 +05:30
groupHandle = HDF5_openGroup ( fileHandle , 'homogenization' )
2020-10-24 18:26:03 +05:30
do i = 1 , size ( material_name_homogenization )
2020-11-18 19:22:16 +05:30
write ( datasetName , '(i0,a)' ) i , '_omega'
2020-02-25 22:23:15 +05:30
call HDF5_read ( groupHandle , homogState ( i ) % state0 , datasetName )
enddo
call HDF5_closeGroup ( groupHandle )
call HDF5_closeFile ( fileHandle )
end subroutine crystallite_restartRead
2020-02-25 14:10:38 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
2020-02-25 22:23:15 +05:30
subroutine crystallite_forward
2020-02-25 14:10:38 +05:30
integer :: i , j
2020-10-08 01:45:13 +05:30
crystallite_F0 = crystallite_partitionedF
2020-02-25 14:10:38 +05:30
crystallite_Fp0 = crystallite_Fp
crystallite_Lp0 = crystallite_Lp
crystallite_Fi0 = crystallite_Fi
crystallite_Li0 = crystallite_Li
crystallite_S0 = crystallite_S
do i = 1 , size ( plasticState )
plasticState ( i ) % state0 = plasticState ( i ) % state
enddo
do i = 1 , size ( sourceState )
do j = 1 , phase_Nsources ( i )
sourceState ( i ) % p ( j ) % state0 = sourceState ( i ) % p ( j ) % state
enddo ; enddo
2020-10-24 18:26:03 +05:30
do i = 1 , size ( material_name_homogenization )
2020-02-25 14:10:38 +05:30
homogState ( i ) % state0 = homogState ( i ) % state
thermalState ( i ) % state0 = thermalState ( i ) % state
damageState ( i ) % state0 = damageState ( i ) % state
enddo
2020-02-25 22:23:15 +05:30
end subroutine crystallite_forward
2020-02-25 14:10:38 +05:30
2013-02-22 04:38:36 +05:30
end module crystallite
