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
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 debug
use numerics
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
2019-04-03 16:24:07 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , public , protected :: &
crystallite_Fe , & !< current "elastic" def grad (end of converged time step)
2020-02-25 14:20:21 +05:30
crystallite_P , & !< 1st Piola-Kirchhoff stress per grain
crystallite_S0 , & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_Fp0 , & !< plastic def grad at start of FE inc
crystallite_Fi0 , & !< intermediate def grad at start of FE inc
crystallite_F0 , & !< def grad at start of FE inc
crystallite_Lp0 , & !< plastic velocitiy grad at start of FE inc
crystallite_Li0 !< intermediate velocitiy grad at start of FE inc
2019-04-03 16:24:07 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , public :: &
crystallite_S , & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
crystallite_partionedS0 , & !< 2nd Piola-Kirchhoff stress vector at start of homog inc
crystallite_Fp , & !< current plastic def grad (end of converged time step)
crystallite_partionedFp0 , & !< plastic def grad at start of homog inc
crystallite_Fi , & !< current intermediate def grad (end of converged time step)
crystallite_partionedFi0 , & !< intermediate def grad at start of homog inc
crystallite_partionedF , & !< def grad to be reached at end of homog inc
crystallite_partionedF0 , & !< def grad at start of homog inc
crystallite_Lp , & !< current plastic velocitiy grad (end of converged time step)
crystallite_partionedLp0 , & !< plastic velocity grad at start of homog inc
crystallite_Li , & !< current intermediate velocitiy grad (end of converged time step)
crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc
2019-05-17 02:26:48 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: &
2019-04-03 16:24:07 +05:30
crystallite_subFp0 , & !< plastic def grad at start of crystallite inc
crystallite_subFi0 , & !< intermediate def grad at start of crystallite inc
crystallite_subF , & !< def grad to be reached at end of crystallite inc
crystallite_subF0 , & !< def grad at start of crystallite inc
crystallite_subLp0 , & !< plastic velocity grad at start of crystallite inc
crystallite_subLi0 !< intermediate velocity grad at start of crystallite inc
2020-02-15 02:01:03 +05:30
real ( pReal ) , dimension ( : , : , : , : , : , : , : ) , allocatable , public , protected :: &
2019-04-03 16:24:07 +05:30
crystallite_dPdF !< current individual dPdF per grain (end of converged time step)
logical , dimension ( : , : , : ) , allocatable , public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
2019-05-17 02:26:48 +05:30
logical , dimension ( : , : , : ) , allocatable :: &
2019-04-03 16:24:07 +05:30
crystallite_converged , & !< convergence flag
crystallite_todo , & !< flag to indicate need for further computation
crystallite_localPlasticity !< indicates this grain to have purely local constitutive law
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-02-13 21:40:27 +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
2019-04-03 16:24:07 +05:30
procedure ( ) , 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 , &
crystallite_forward
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-01-03 18:03:32 +05:30
logical , dimension ( discretization_nIP , discretization_nElem ) :: 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
eMax , & !< maximum number of elements
2019-11-30 20:33:18 +05:30
myNcomponents !< number of components at current IP
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
write ( 6 , '(/,a)' ) ' <<<+- crystallite init -+>>>'
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
cMax = homogenization_maxNgrains
2019-06-07 02:19:17 +05:30
iMax = discretization_nIP
eMax = discretization_nElem
2020-02-13 21:40:27 +05:30
2020-02-29 18:23:04 +05:30
allocate ( crystallite_partionedF ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_S0 , &
crystallite_F0 , crystallite_Fi0 , crystallite_Fp0 , &
crystallite_Li0 , crystallite_Lp0 , &
crystallite_partionedS0 , &
crystallite_partionedF0 , crystallite_partionedFp0 , crystallite_partionedFi0 , &
crystallite_partionedLp0 , crystallite_partionedLi0 , &
crystallite_S , crystallite_P , &
crystallite_Fe , crystallite_Fi , crystallite_Fp , &
crystallite_Li , crystallite_Lp , &
crystallite_subF , crystallite_subF0 , &
crystallite_subFp0 , crystallite_subFi0 , &
crystallite_subLi0 , crystallite_subLp0 , &
source = crystallite_partionedF )
allocate ( crystallite_dPdF ( 3 , 3 , 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_dt ( cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_subdt , 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_localPlasticity ( cMax , iMax , eMax ) , source = . true . )
allocate ( crystallite_requested ( cMax , iMax , eMax ) , source = . false . )
allocate ( crystallite_todo ( cMax , iMax , eMax ) , source = . false . )
allocate ( crystallite_converged ( cMax , iMax , eMax ) , source = . true . )
2020-02-13 21:40:27 +05:30
2019-04-11 14:57:03 +05:30
num % subStepMinCryst = config_numerics % getFloat ( 'substepmincryst' , defaultVal = 1.0e-3_pReal )
num % subStepSizeCryst = config_numerics % getFloat ( 'substepsizecryst' , defaultVal = 0.25_pReal )
num % stepIncreaseCryst = config_numerics % getFloat ( 'stepincreasecryst' , defaultVal = 1.5_pReal )
2020-02-13 21:40:27 +05:30
2019-04-11 14:57:03 +05:30
num % subStepSizeLp = config_numerics % getFloat ( 'substepsizelp' , defaultVal = 0.5_pReal )
num % subStepSizeLi = config_numerics % getFloat ( 'substepsizeli' , defaultVal = 0.5_pReal )
2020-03-14 23:41:26 +05:30
num % rtol_crystalliteState = config_numerics % getFloat ( 'rtol_crystallitestate' , defaultVal = 1.0e-6_pReal )
num % rtol_crystalliteStress = config_numerics % getFloat ( 'rtol_crystallitestress' , defaultVal = 1.0e-6_pReal )
num % atol_crystalliteStress = config_numerics % getFloat ( 'atol_crystallitestress' , defaultVal = 1.0e-8_pReal )
2020-02-13 21:40:27 +05:30
2019-04-11 14:57:03 +05:30
num % iJacoLpresiduum = config_numerics % getInt ( 'ijacolpresiduum' , defaultVal = 1 )
2020-02-13 21:40:27 +05:30
2019-04-11 14:57:03 +05:30
num % nState = config_numerics % getInt ( 'nstate' , defaultVal = 20 )
num % nStress = config_numerics % getInt ( '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
2019-04-03 16:24:07 +05:30
select case ( numerics_integrator )
case ( 1 )
integrateState = > integrateStateFPI
case ( 2 )
integrateState = > integrateStateEuler
case ( 3 )
integrateState = > integrateStateAdaptiveEuler
case ( 4 )
integrateState = > integrateStateRK4
case ( 5 )
integrateState = > integrateStateRKCK45
end select
2020-02-13 21:40:27 +05:30
2019-04-06 15:36:34 +05:30
allocate ( output_constituent ( size ( config_phase ) ) )
2019-04-06 10:01:02 +05:30
do c = 1 , size ( config_phase )
#if defined(__GFORTRAN__)
2020-02-13 21:40:27 +05:30
allocate ( output_constituent ( c ) % label ( 1 ) )
2019-04-06 15:36:34 +05:30
output_constituent ( c ) % label ( 1 ) = 'GfortranBug86277'
output_constituent ( c ) % label = config_phase ( c ) % getStrings ( '(output)' , defaultVal = output_constituent ( c ) % label )
if ( output_constituent ( c ) % label ( 1 ) == 'GfortranBug86277' ) output_constituent ( c ) % label = [ character ( len = pStringLen ) :: ]
2019-04-06 10:01:02 +05:30
#else
2019-04-06 15:36:34 +05:30
output_constituent ( c ) % label = config_phase ( c ) % getStrings ( '(output)' , defaultVal = [ character ( len = pStringLen ) :: ] )
2019-04-06 10:01:02 +05:30
#endif
enddo
2019-04-03 16:24:07 +05:30
2019-04-06 10:01:02 +05:30
call config_deallocate ( 'material.config/phase' )
2018-06-27 00:03:41 +05:30
2013-04-29 16:47:30 +05:30
!--------------------------------------------------------------------------------------------------
! initialize
2019-01-18 16:46:26 +05:30
!$OMP PARALLEL DO PRIVATE(myNcomponents,i,c)
2019-04-03 16:24:07 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2019-06-05 13:35:59 +05:30
myNcomponents = homogenization_Ngrains ( material_homogenizationAt ( e ) )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 ) ; do c = 1 , myNcomponents
2020-02-13 13:19:49 +05:30
crystallite_Fp0 ( 1 : 3 , 1 : 3 , c , i , e ) = material_orientation0 ( c , i , e ) % asMatrix ( ) ! plastic def gradient reflects init orientation
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
2019-06-14 12:19:05 +05:30
crystallite_localPlasticity ( c , i , e ) = phase_localPlasticity ( material_phaseAt ( c , e ) )
2019-04-03 16:24:07 +05:30
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
2019-04-03 16:24:07 +05:30
if ( any ( . not . crystallite_localPlasticity ) . and . . not . usePingPong ) call IO_error ( 601 ) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal?
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
crystallite_partionedFp0 = crystallite_Fp0
crystallite_partionedFi0 = crystallite_Fi0
crystallite_partionedF0 = crystallite_F0
crystallite_partionedF = 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 )
2019-06-05 13:35:59 +05:30
do c = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2020-02-11 22:06:43 +05:30
call constitutive_dependentState ( crystallite_partionedF0 ( 1 : 3 , 1 : 3 , c , i , e ) , &
2020-02-07 17:11:01 +05:30
crystallite_partionedFp0 ( 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 ( )
call crystallite_stressTangent
2019-01-18 19:12:44 +05:30
2018-09-20 09:39:02 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelBasic ) / = 0 ) then
write ( 6 , '(a42,1x,i10)' ) ' # of elements: ' , eMax
write ( 6 , '(a42,1x,i10)' ) 'max # of integration points/element: ' , iMax
write ( 6 , '(a42,1x,i10)' ) 'max # of constituents/integration point: ' , cMax
write ( 6 , '(a42,1x,i10)' ) ' # of nonlocal constituents: ' , count ( . not . crystallite_localPlasticity )
flush ( 6 )
endif
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
call debug_info
call debug_reset
2018-09-20 09:39:02 +05:30
#endif
2010-11-03 22:52:48 +05:30
2012-03-09 01:55:28 +05:30
end subroutine crystallite_init
2009-05-07 21:57:36 +05:30
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
2019-02-17 16:45:46 +05:30
function crystallite_stress ( dummyArgumentToPreventInternalCompilerErrorWithGCC )
2020-02-13 21:40:27 +05:30
2019-06-07 02:19:17 +05:30
logical , dimension ( discretization_nIP , discretization_nElem ) :: crystallite_stress
2019-04-03 16:24:07 +05:30
real ( pReal ) , intent ( in ) , optional :: &
dummyArgumentToPreventInternalCompilerErrorWithGCC
real ( pReal ) :: &
formerSubStep
integer :: &
NiterationCrystallite , & ! number of iterations in crystallite loop
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e , & !< counter in element loop
startIP , endIP , &
s
2020-02-13 21:40:27 +05:30
2019-01-18 19:12:44 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 &
. and . FEsolving_execElem ( 1 ) < = debug_e &
. and . debug_e < = FEsolving_execElem ( 2 ) ) then
write ( 6 , '(/,a,i8,1x,i2,1x,i3)' ) '<< CRYST stress >> boundary and initial values at el ip ipc ' , &
debug_e , debug_i , debug_g
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> F ' , &
transpose ( crystallite_partionedF ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> F0 ' , &
transpose ( crystallite_partionedF0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> Fp0' , &
transpose ( crystallite_partionedFp0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> Fi0' , &
transpose ( crystallite_partionedFi0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> Lp0' , &
transpose ( crystallite_partionedLp0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
write ( 6 , '(a,/,3(12x,3(f14.9,1x)/))' ) '<< CRYST stress >> Li0' , &
transpose ( crystallite_partionedLi0 ( 1 : 3 , 1 : 3 , debug_g , debug_i , debug_e ) )
endif
2019-01-18 19:12:44 +05:30
#endif
!--------------------------------------------------------------------------------------------------
! initialize to starting condition
2019-04-03 16:24:07 +05:30
crystallite_subStep = 0.0_pReal
!$OMP PARALLEL DO
elementLooping1 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 ) ; do c = 1 , homogenization_Ngrains ( 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 ) ) = &
plasticState ( material_phaseAt ( c , e ) ) % partionedState0 ( : , 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 ) ) = &
sourceState ( material_phaseAt ( c , e ) ) % p ( s ) % partionedState0 ( : , material_phaseMemberAt ( c , i , e ) )
2019-04-03 16:24:07 +05:30
enddo
crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedFp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subLp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedLp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedFi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subLi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedLi0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_partionedF0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subFrac ( c , i , e ) = 0.0_pReal
2019-04-11 10:54:04 +05:30
crystallite_subStep ( c , i , e ) = 1.0_pReal / num % subStepSizeCryst
2019-04-03 16:24:07 +05:30
crystallite_todo ( c , i , e ) = . true .
crystallite_converged ( c , i , e ) = . false . ! pretend failed step of 1/subStepSizeCryst
endif homogenizationRequestsCalculation
enddo ; enddo
enddo elementLooping1
!$OMP END PARALLEL DO
singleRun : if ( FEsolving_execELem ( 1 ) == FEsolving_execElem ( 2 ) . and . &
2020-01-25 13:54:42 +05:30
FEsolving_execIP ( 1 ) == FEsolving_execIP ( 2 ) ) then
startIP = FEsolving_execIP ( 1 )
2019-04-03 16:24:07 +05:30
endIP = startIP
else singleRun
startIP = 1
2020-01-25 13:54:42 +05:30
endIP = discretization_nIP
2019-04-03 16:24:07 +05:30
endif singleRun
NiterationCrystallite = 0
cutbackLooping : do while ( any ( crystallite_todo ( : , startIP : endIP , FEsolving_execELem ( 1 ) : FEsolving_execElem ( 2 ) ) ) )
NiterationCrystallite = NiterationCrystallite + 1
2019-02-27 00:17:46 +05:30
2019-01-18 19:12:44 +05:30
#ifdef DEBUG
2019-04-03 16:24:07 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 ) &
write ( 6 , '(a,i6)' ) '<< CRYST stress >> crystallite iteration ' , NiterationCrystallite
2019-01-18 19:12:44 +05:30
#endif
2019-04-03 16:24:07 +05:30
!$OMP PARALLEL DO PRIVATE(formerSubStep)
elementLooping3 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do c = 1 , homogenization_Ngrains ( 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
crystallite_todo ( c , i , e ) = crystallite_subStep ( c , i , e ) > 0.0_pReal ! still time left to integrate on?
if ( crystallite_todo ( c , i , e ) ) then
crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subLp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subLi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subFp0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_subFi0 ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e )
!if abbrevation, make c and p private in omp
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
crystallite_Lp ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subLp0 ( 1 : 3 , 1 : 3 , c , i , e )
crystallite_Li ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subLi0 ( 1 : 3 , 1 : 3 , c , i , e )
endif
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
2019-04-11 10:54:04 +05:30
crystallite_todo ( c , i , e ) = crystallite_subStep ( c , i , e ) > num % subStepMinCryst ! still on track or already done (beyond repair)
2019-04-03 16:24:07 +05:30
endif
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
! prepare for integration
2019-04-03 16:24:07 +05:30
if ( crystallite_todo ( c , i , e ) ) then
crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) = crystallite_subF0 ( 1 : 3 , 1 : 3 , c , i , e ) &
2020-02-21 13:15:11 +05:30
+ crystallite_subStep ( c , i , e ) * ( crystallite_partionedF ( 1 : 3 , 1 : 3 , c , i , e ) &
- crystallite_partionedF0 ( 1 : 3 , 1 : 3 , c , i , e ) )
crystallite_Fe ( 1 : 3 , 1 : 3 , c , i , e ) = matmul ( matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
2020-02-21 13:11:08 +05:30
math_inv33 ( crystallite_Fp ( 1 : 3 , 1 : 3 , c , i , e ) ) ) , &
math_inv33 ( crystallite_Fi ( 1 : 3 , 1 : 3 , c , i , e ) ) )
2019-04-03 16:24:07 +05:30
crystallite_subdt ( c , i , e ) = crystallite_subStep ( c , i , e ) * crystallite_dt ( c , i , e )
crystallite_converged ( c , i , e ) = . false .
endif
enddo
enddo
enddo elementLooping3
!$OMP END PARALLEL DO
2019-01-18 19:12:44 +05:30
!--------------------------------------------------------------------------------------------------
! integrate --- requires fully defined state array (basic + dependent state)
2019-04-03 16:24:07 +05:30
if ( any ( crystallite_todo ) ) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation
2019-04-11 10:54:04 +05:30
where ( . not . crystallite_converged . and . crystallite_subStep > num % subStepMinCryst ) & ! do not try non-converged but fully cutbacked any further
2019-04-03 16:24:07 +05:30
crystallite_todo = . true . ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation
2019-01-18 19:12:44 +05:30
2019-04-03 16:24:07 +05:30
enddo cutbackLooping
2019-01-18 19:12:44 +05:30
! return whether converged or not
2019-04-03 16:24:07 +05:30
crystallite_stress = . false .
elementLooping5 : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
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
2019-01-18 16:46:26 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate tangent (dPdF)
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine crystallite_stressTangent
2019-04-03 16:24:07 +05:30
integer :: &
c , & !< counter in integration point component loop
i , & !< counter in integration point loop
e , & !< counter in element loop
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
2020-02-21 13:11:08 +05:30
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,o,p, &
!$OMP invSubFp0,invSubFi0,invFp,invFi, &
2019-04-03 16:24:07 +05:30
!$OMP rhs_3333,lhs_3333,temp_99,temp_33_1,temp_33_2,temp_33_3,temp_33_4,temp_3333,error)
elementLooping : do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do c = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-04-03 16:24:07 +05:30
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 )
2019-09-21 06:50:33 +05:30
call math_invert ( temp_99 , error , math_identity2nd ( 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
crystallite_dPdF ( 1 : 3 , 1 : 3 , 1 : 3 , 1 : 3 , c , i , e ) = 0.0_pReal
2019-04-03 16:47:21 +05:30
do p = 1 , 3
2020-02-11 22:11:30 +05:30
crystallite_dPdF ( p , 1 : 3 , p , 1 : 3 , c , i , e ) = 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-02-14 10:52:38 +05:30
crystallite_dPdF ( 1 : 3 , 1 : 3 , p , o , c , i , e ) = crystallite_dPdF ( 1 : 3 , 1 : 3 , p , o , c , i , e ) &
+ matmul ( matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , c , i , e ) , &
dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) , temp_33_1 ) &
+ matmul ( matmul ( temp_33_3 , dSdF ( 1 : 3 , 1 : 3 , p , o ) ) , &
2020-02-21 13:11:08 +05:30
transpose ( invFp ) ) &
2020-02-14 10:52:38 +05:30
+ matmul ( temp_33_4 , transpose ( dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) )
2019-04-03 16:47:21 +05:30
enddo ; enddo
2019-04-03 16:24:07 +05:30
enddo ; enddo
enddo elementLooping
!$OMP END PARALLEL DO
2019-01-18 16:46:26 +05:30
end subroutine crystallite_stressTangent
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-19 14:05:45 +05:30
!> @brief calculates orientations
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-19 14:05:45 +05:30
subroutine crystallite_orientations
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 )
2019-06-05 13:35:59 +05:30
do c = 1 , homogenization_Ngrains ( 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
2019-04-03 16:24:07 +05:30
nonlocalPresent : if ( any ( plasticState % nonLocal ) ) then
!$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-03-16 14:19:59 +05:30
if ( plasticState ( material_phaseAt ( 1 , e ) ) % nonLocal ) &
call plastic_nonlocal_updateCompatibility ( crystallite_orientation , &
phase_plasticityInstance ( material_phaseAt ( i , e ) ) , i , e )
2019-04-03 16:24:07 +05:30
enddo ; enddo
!$OMP END PARALLEL DO
endif nonlocalPresent
2014-08-26 20:14:32 +05:30
2019-01-19 14:05:45 +05:30
end subroutine crystallite_orientations
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Map 2nd order tensor to reference config
!--------------------------------------------------------------------------------------------------
function crystallite_push33ToRef ( ipc , ip , el , tensor33 )
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-03-14 23:41:26 +05:30
character ( len = pStringLen ) :: group , structureLabel
2020-02-13 21:40:27 +05:30
2019-04-06 10:01:02 +05:30
do p = 1 , size ( config_name_phase )
2019-04-18 15:25:50 +05:30
group = trim ( 'current/constituent' ) / / '/' / / trim ( config_name_phase ( p ) ) / / '/generic'
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 ) )
case ( 'f' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_partionedF , p )
call results_writeDataset ( group , selected_tensors , 'F' , &
2019-04-06 15:36:34 +05:30
'deformation gradient' , '1' )
case ( 'fe' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Fe , p )
call results_writeDataset ( group , selected_tensors , 'Fe' , &
2019-04-06 15:36:34 +05:30
'elastic deformation gradient' , '1' )
case ( 'fp' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Fp , p )
call results_writeDataset ( group , selected_tensors , 'Fp' , &
2019-04-06 15:36:34 +05:30
'plastic deformation gradient' , '1' )
case ( 'fi' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Fi , p )
call results_writeDataset ( group , selected_tensors , 'Fi' , &
2019-04-06 15:36:34 +05:30
'inelastic deformation gradient' , '1' )
case ( 'lp' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Lp , p )
call results_writeDataset ( group , selected_tensors , 'Lp' , &
2019-04-06 15:36:34 +05:30
'plastic velocity gradient' , '1/s' )
case ( 'li' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_Li , p )
call results_writeDataset ( group , selected_tensors , 'Li' , &
2019-04-06 15:36:34 +05:30
'inelastic velocity gradient' , '1/s' )
case ( 'p' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_P , p )
call results_writeDataset ( group , selected_tensors , 'P' , &
2020-03-19 16:00:36 +05:30
'First Piola-Kirchoff stress' , 'Pa' )
2019-04-06 15:36:34 +05:30
case ( 's' )
2019-04-07 17:58:08 +05:30
selected_tensors = select_tensors ( crystallite_S , p )
call results_writeDataset ( group , selected_tensors , 'S' , &
2020-03-19 16:00:36 +05:30
'Second Piola-Kirchoff stress' , 'Pa' )
2019-04-07 17:58:08 +05:30
case ( 'orientation' )
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 )
call results_writeDataset ( group , selected_rotations , 'orientation' , &
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-02-22 15:56:12 +05:30
allocate ( select_tensors ( 3 , 3 , count ( material_phaseAt == instance ) * discretization_nIP ) )
2019-04-13 19:04:51 +05:30
j = 0
do e = 1 , size ( material_phaseAt , 2 )
2019-10-19 23:25:00 +05:30
do i = 1 , discretization_nIP
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
2019-10-19 23:25:00 +05:30
allocate ( select_rotations ( count ( material_phaseAt == instance ) * homogenization_maxNgrains * discretization_nIP ) )
2019-04-13 19:04:51 +05:30
j = 0
do e = 1 , size ( material_phaseAt , 2 )
2019-10-19 23:25:00 +05:30
do i = 1 , discretization_nIP
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
!--------------------------------------------------------------------------------------------------
2019-04-03 16:24:07 +05:30
logical function integrateStress ( ipc , ip , el , timeFraction )
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
Fe_new , &
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-02-11 22:17:48 +05:30
logical :: error
2019-04-03 16:24:07 +05:30
external :: &
dgesv
2020-02-13 21:40:27 +05:30
2019-04-03 16:24:07 +05:30
integrateStress = . false .
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-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
2019-01-15 08:57:57 +05:30
!* calculate Jacobian for correction term
2019-04-11 14:57:03 +05:30
if ( mod ( jacoCounterLp , num % iJacoLpresiduum ) == 0 ) then
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
dFe_dLp ( o , 1 : 3 , p , 1 : 3 ) = A ( o , p ) * transpose ( invFi_new ) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
enddo ; enddo
2019-04-03 16:24:07 +05:30
dFe_dLp = - dt * dFe_dLp
dRLp_dLp = math_identity2nd ( 9 ) &
- math_3333to99 ( math_mul3333xx3333 ( math_mul3333xx3333 ( dLp_dS , dS_dFe ) , dFe_dLp ) )
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 )
2019-04-03 16:24:07 +05:30
endif
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
2019-04-03 16:24:07 +05:30
!* calculate Jacobian for correction term
2019-04-11 14:57:03 +05:30
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
2019-04-03 16:24:07 +05:30
dRLi_dLi = math_identity2nd ( 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 )
2019-04-03 16:24:07 +05:30
endif
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
2020-02-13 15:47:31 +05:30
Fp_new = Fp_new / math_det33 ( Fp_new ) ** ( 1.0_pReal / 3.0_pReal ) ! regularize
2020-02-11 22:20:07 +05:30
Fe_new = matmul ( matmul ( F , invFp_new ) , invFi_new )
2014-08-26 20:14:32 +05:30
2019-04-03 16:24:07 +05:30
integrateStress = . true .
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
crystallite_Fp ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Fp_new
crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Fi_new
crystallite_Fe ( 1 : 3 , 1 : 3 , ipc , ip , el ) = Fe_new
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
!--------------------------------------------------------------------------------------------------
2019-04-11 11:05:58 +05:30
subroutine integrateStateFPI
2019-01-15 08:57:57 +05:30
2020-03-23 11:09:17 +05:30
integer :: &
NiterationState , & !< number of iterations in state loop
e , & !< element index in element loop
i , & !< integration point index in ip loop
g , & !< grain index in grain loop
p , &
c , &
s , &
sizeDotState
real ( pReal ) :: &
zeta
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState ) :: &
residuum_plastic ! residuum for plastic state
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState ) :: &
residuum_source ! residuum for source state
logical :: &
doneWithIntegration , &
nonlocalBroken
2014-08-26 20:14:32 +05:30
2020-03-23 11:09:17 +05:30
! --+>> PREGUESS FOR STATE <<+--
call update_dotState ( 1.0_pReal )
call update_state ( 1.0_pReal )
2015-10-14 00:22:01 +05:30
2020-03-23 11:09:17 +05:30
NiterationState = 0
doneWithIntegration = . false .
nonlocalBroken = . false .
iteration : do while ( . not . doneWithIntegration . and . NiterationState < num % nState )
NiterationState = NiterationState + 1
2019-01-15 08:57:57 +05:30
2019-02-27 02:01:47 +05:30
#ifdef DEBUG
2019-04-03 16:02:30 +05:30
if ( iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 ) &
2019-02-27 00:17:46 +05:30
write ( 6 , '(a,i6)' ) '<< CRYST stateFPI >> state iteration ' , NiterationState
2019-02-27 02:01:47 +05:30
#endif
2019-02-27 00:17:46 +05:30
2020-03-23 11:09:17 +05:30
!$OMP PARALLEL DO PRIVATE(p,c)
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) . and . &
( . not . nonlocalBroken . or . crystallite_localPlasticity ( g , i , e ) ) ) then
p = material_phaseAt ( g , e ) ; c = material_phaseMemberAt ( g , i , e )
plasticState ( p ) % previousDotState2 ( : , c ) = merge ( plasticState ( p ) % previousDotState ( : , c ) , &
0.0_pReal , &
NiterationState > 1 )
plasticState ( p ) % previousDotState ( : , c ) = plasticState ( p ) % dotState ( : , c )
do s = 1 , phase_Nsources ( p )
sourceState ( p ) % p ( s ) % previousDotState2 ( : , c ) = merge ( sourceState ( p ) % p ( s ) % previousDotState ( : , c ) , &
0.0_pReal , &
NiterationState > 1 )
sourceState ( p ) % p ( s ) % previousDotState ( : , c ) = sourceState ( p ) % p ( s ) % dotState ( : , c )
enddo
call constitutive_dependentState ( crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fp ( 1 : 3 , 1 : 3 , g , i , e ) , &
g , i , e )
crystallite_todo ( g , i , e ) = integrateStress ( g , i , e , 1.0_pReal )
if ( . not . ( crystallite_todo ( g , i , e ) . or . crystallite_localPlasticity ( g , i , e ) ) ) &
nonlocalBroken = . true .
endif
enddo
2020-03-23 05:02:41 +05:30
enddo
enddo
2020-03-23 11:09:17 +05:30
!$OMP END PARALLEL DO
2020-02-13 21:40:27 +05:30
2020-03-23 11:09:17 +05:30
if ( nonlocalBroken ) where ( . not . crystallite_localPlasticity ) crystallite_todo = . false .
2014-05-27 20:16:03 +05:30
2020-03-23 11:09:17 +05:30
call update_dotState ( 1.0_pReal )
2020-02-13 21:40:27 +05:30
2020-03-23 11:09:17 +05:30
!$OMP PARALLEL
!$OMP DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c)
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
p = material_phaseAt ( g , e ) ; c = material_phaseMemberAt ( g , i , e )
sizeDotState = plasticState ( p ) % sizeDotState
zeta = damper ( plasticState ( p ) % dotState ( : , c ) , &
plasticState ( p ) % previousDotState ( : , c ) , &
plasticState ( p ) % previousDotState2 ( : , c ) )
residuum_plastic ( 1 : SizeDotState ) = plasticState ( p ) % state ( 1 : sizeDotState , c ) &
- plasticState ( p ) % subState0 ( 1 : sizeDotState , c ) &
- ( plasticState ( p ) % dotState ( : , c ) * zeta &
+ plasticState ( p ) % previousDotState ( : , c ) * ( 1.0_pReal - zeta ) &
) * crystallite_subdt ( g , i , e )
2019-01-15 08:57:57 +05:30
2020-03-23 11:09:17 +05:30
plasticState ( p ) % state ( 1 : sizeDotState , c ) = plasticState ( p ) % state ( 1 : sizeDotState , c ) &
- residuum_plastic ( 1 : sizeDotState )
plasticState ( p ) % dotState ( : , c ) = plasticState ( p ) % dotState ( : , c ) * zeta &
+ plasticState ( p ) % previousDotState ( : , c ) * ( 1.0_pReal - zeta )
2020-02-13 21:40:27 +05:30
2020-03-23 11:09:17 +05:30
crystallite_converged ( g , i , e ) = converged ( residuum_plastic ( 1 : sizeDotState ) , &
plasticState ( p ) % state ( 1 : sizeDotState , c ) , &
plasticState ( p ) % atol ( 1 : sizeDotState ) )
2020-02-13 21:40:27 +05:30
2014-05-27 20:16:03 +05:30
2020-03-23 11:09:17 +05:30
do s = 1 , phase_Nsources ( p )
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
2020-02-13 21:40:27 +05:30
2020-03-23 11:09:17 +05:30
zeta = damper ( sourceState ( p ) % p ( s ) % dotState ( : , c ) , &
sourceState ( p ) % p ( s ) % previousDotState ( : , c ) , &
sourceState ( p ) % p ( s ) % previousDotState2 ( : , c ) )
2019-01-30 04:01:26 +05:30
2020-03-23 11:09:17 +05:30
residuum_source ( 1 : sizeDotState ) = sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) &
- sourceState ( p ) % p ( s ) % subState0 ( 1 : sizeDotState , c ) &
- ( sourceState ( p ) % p ( s ) % dotState ( : , c ) * zeta &
+ sourceState ( p ) % p ( s ) % previousDotState ( : , c ) * ( 1.0_pReal - zeta ) &
) * crystallite_subdt ( g , i , e )
2019-01-30 02:59:36 +05:30
2020-03-23 11:09:17 +05:30
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) = sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) &
- residuum_source ( 1 : sizeDotState )
sourceState ( p ) % p ( s ) % dotState ( : , c ) = sourceState ( p ) % p ( s ) % dotState ( : , c ) * zeta &
+ sourceState ( p ) % p ( s ) % previousDotState ( : , c ) * ( 1.0_pReal - zeta )
2019-01-30 02:59:36 +05:30
2020-03-23 11:09:17 +05:30
crystallite_converged ( g , i , e ) = &
crystallite_converged ( g , i , e ) . and . converged ( residuum_source ( 1 : sizeDotState ) , &
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) , &
sourceState ( p ) % p ( s ) % atol ( 1 : sizeDotState ) )
enddo
endif
enddo ; enddo ; enddo
!$OMP ENDDO
2019-01-15 08:57:57 +05:30
2020-03-23 11:09:17 +05:30
!$OMP DO
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
!$OMP FLUSH(crystallite_todo)
if ( crystallite_todo ( g , i , e ) . and . crystallite_converged ( g , i , e ) ) then ! converged and still alive...
crystallite_todo ( g , i , e ) = stateJump ( g , i , e )
!$OMP FLUSH(crystallite_todo)
if ( . not . crystallite_todo ( g , i , e ) ) then ! if state jump fails, then convergence is broken
crystallite_converged ( g , i , e ) = . false .
if ( . not . crystallite_localPlasticity ( g , i , e ) ) then ! if broken non-local...
!$OMP CRITICAL (checkTodo)
crystallite_todo = crystallite_todo . and . crystallite_localPlasticity ! ...all non-locals skipped
!$OMP END CRITICAL (checkTodo)
endif
endif
endif
enddo ; enddo ; enddo
!$OMP ENDDO
!$OMP END PARALLEL
2019-01-15 08:57:57 +05:30
2014-08-26 20:14:32 +05:30
2020-03-23 11:09:17 +05:30
if ( any ( plasticState ( : ) % nonlocal ) ) call nonlocalConvergenceCheck
2014-08-26 20:14:32 +05:30
2014-05-27 20:16:03 +05:30
2020-03-23 11:09:17 +05:30
! --- CHECK IF DONE WITH INTEGRATION ---
doneWithIntegration = . true .
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
doneWithIntegration = . false .
exit
endif
enddo ; enddo
enddo
2019-01-15 08:57:57 +05:30
2020-03-23 11:09:17 +05:30
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
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine integrateStateEuler
2014-08-26 20:14:32 +05:30
2019-01-23 10:44:19 +05:30
call update_dotState ( 1.0_pReal )
2019-01-29 09:41:29 +05:30
call update_state ( 1.0_pReal )
2019-01-24 16:03:04 +05:30
call update_deltaState
call update_dependentState
2019-01-24 22:29:38 +05:30
call update_stress ( 1.0_pReal )
2019-01-25 11:50:05 +05:30
call setConvergenceFlag
2019-01-29 09:41:29 +05:30
if ( any ( plasticState ( : ) % nonlocal ) ) call nonlocalConvergenceCheck
2014-08-26 20:14:32 +05:30
2018-09-20 09:57:53 +05:30
end subroutine integrateStateEuler
2010-10-01 17:48:49 +05:30
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine integrateStateAdaptiveEuler
2014-07-02 17:57:39 +05:30
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-15 08:57:57 +05:30
e , & ! element index in element loop
i , & ! integration point index in ip loop
g , & ! grain index in grain loop
2014-06-25 04:51:25 +05:30
p , &
c , &
2019-01-30 03:34:50 +05:30
s , &
2019-01-30 03:16:21 +05:30
sizeDotState
2020-02-13 21:40:27 +05:30
2019-01-30 14:58:47 +05:30
! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState , &
2019-06-07 02:19:17 +05:30
homogenization_maxNgrains , discretization_nIP , discretization_nElem ) :: &
2019-01-30 14:58:47 +05:30
residuum_plastic
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , &
maxval ( phase_Nsources ) , &
2019-06-07 02:19:17 +05:30
homogenization_maxNgrains , discretization_nIP , discretization_nElem ) :: &
2019-01-30 04:15:41 +05:30
residuum_source
2014-05-27 20:16:03 +05:30
2019-01-23 10:44:19 +05:30
!--------------------------------------------------------------------------------------------------
! contribution to state and relative residui and from Euler integration
call update_dotState ( 1.0_pReal )
2014-08-26 20:14:32 +05:30
2019-01-30 03:16:21 +05:30
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
2019-01-30 04:01:26 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-06-14 12:19:05 +05:30
p = material_phaseAt ( g , e ) ; c = material_phaseMemberAt ( g , i , e )
2019-01-30 03:16:21 +05:30
sizeDotState = plasticState ( p ) % sizeDotState
2020-02-13 21:40:27 +05:30
2019-01-30 04:15:41 +05:30
residuum_plastic ( 1 : sizeDotState , g , i , e ) = plasticState ( p ) % dotstate ( 1 : sizeDotState , c ) &
* ( - 0.5_pReal * crystallite_subdt ( g , i , e ) )
2019-02-27 02:01:47 +05:30
plasticState ( p ) % state ( 1 : sizeDotState , c ) = &
plasticState ( p ) % state ( 1 : sizeDotState , c ) + plasticState ( p ) % dotstate ( 1 : sizeDotState , c ) * crystallite_subdt ( g , i , e ) !ToDo: state, partitioned state?
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 03:34:50 +05:30
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
2020-02-13 21:40:27 +05:30
2019-01-30 04:15:41 +05:30
residuum_source ( 1 : sizeDotState , s , g , i , e ) = sourceState ( p ) % p ( s ) % dotstate ( 1 : sizeDotState , c ) &
* ( - 0.5_pReal * crystallite_subdt ( g , i , e ) )
2019-02-27 02:01:47 +05:30
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) = &
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) + sourceState ( p ) % p ( s ) % dotstate ( 1 : sizeDotState , c ) * crystallite_subdt ( g , i , e ) !ToDo: state, partitioned state?
2019-01-15 08:57:57 +05:30
enddo
endif
enddo ; enddo ; enddo
2019-01-29 05:24:02 +05:30
!$OMP END PARALLEL DO
2019-01-15 08:57:57 +05:30
2019-01-30 14:58:47 +05:30
call update_deltaState
call update_dependentState
call update_stress ( 1.0_pReal )
call update_dotState ( 1.0_pReal )
2019-01-15 08:57:57 +05:30
2019-01-30 04:15:41 +05:30
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,c)
2019-01-30 04:01:26 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-06-14 12:19:05 +05:30
p = material_phaseAt ( g , e ) ; c = material_phaseMemberAt ( g , i , e )
2019-01-30 04:01:26 +05:30
sizeDotState = plasticState ( p ) % sizeDotState
2020-02-13 21:40:27 +05:30
2019-01-30 04:15:41 +05:30
residuum_plastic ( 1 : sizeDotState , g , i , e ) = residuum_plastic ( 1 : sizeDotState , g , i , e ) &
+ 0.5_pReal * plasticState ( p ) % dotState ( : , c ) * crystallite_subdt ( g , i , e )
2020-02-13 21:40:27 +05:30
2019-01-31 13:42:44 +05:30
crystallite_converged ( g , i , e ) = converged ( residuum_plastic ( 1 : sizeDotState , g , i , e ) , &
2019-02-27 02:01:47 +05:30
plasticState ( p ) % state ( 1 : sizeDotState , c ) , &
2020-03-15 14:21:40 +05:30
plasticState ( p ) % atol ( 1 : sizeDotState ) )
2015-05-28 22:32:23 +05:30
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 03:34:50 +05:30
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
2020-02-13 21:40:27 +05:30
2019-02-27 02:01:47 +05:30
residuum_source ( 1 : sizeDotState , s , g , i , e ) = &
residuum_source ( 1 : sizeDotState , s , g , i , e ) + 0.5_pReal * sourceState ( p ) % p ( s ) % dotState ( : , c ) * crystallite_subdt ( g , i , e )
2019-01-30 17:06:02 +05:30
2019-02-27 02:01:47 +05:30
crystallite_converged ( g , i , e ) = &
crystallite_converged ( g , i , e ) . and . converged ( residuum_source ( 1 : sizeDotState , s , g , i , e ) , &
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , c ) , &
2020-03-15 14:21:40 +05:30
sourceState ( p ) % p ( s ) % atol ( 1 : sizeDotState ) )
2019-01-30 17:06:02 +05:30
enddo
2020-02-13 21:40:27 +05:30
2013-02-22 04:38:36 +05:30
endif
2019-01-30 04:01:26 +05:30
enddo ; enddo ; enddo
2019-01-29 05:24:02 +05:30
!$OMP END PARALLEL DO
2014-01-22 00:15:41 +05:30
2019-01-29 09:41:29 +05:30
if ( any ( plasticState ( : ) % nonlocal ) ) call nonlocalConvergenceCheck
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
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-15 08:57:57 +05:30
!> @brief integrate stress, state with 4th order explicit Runge Kutta method
2019-01-30 04:47:04 +05:30
! ToDo: This is totally BROKEN: RK4dotState is never used!!!
2013-02-22 04:38:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine integrateStateRK4
2014-09-10 14:07:12 +05:30
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( 4 ) , parameter :: &
TIMESTEPFRACTION = [ 0.5_pReal , 0.5_pReal , 1.0_pReal , 1.0_pReal ] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
real ( pReal ) , dimension ( 4 ) , parameter :: &
WEIGHT = [ 1.0_pReal , 2.0_pReal , 2.0_pReal , 1.0_pReal / 6.0_pReal ] ! weight of slope used for Runge Kutta integration (final weight divided by 6)
2014-09-10 14:07:12 +05:30
2019-04-03 16:02:30 +05:30
integer :: e , & ! element index in element loop
2019-01-15 08:57:57 +05:30
i , & ! integration point index in ip loop
g , & ! grain index in grain loop
p , & ! phase loop
c , &
n , &
2019-01-30 04:31:40 +05:30
s
2019-01-15 08:57:57 +05:30
2019-01-23 10:44:19 +05:30
call update_dotState ( 1.0_pReal )
2018-05-09 20:24:06 +05:30
2014-08-26 20:14:32 +05:30
2019-04-03 16:02:30 +05:30
do n = 1 , 4
2013-04-29 16:47:30 +05:30
2019-01-30 04:41:10 +05:30
!$OMP PARALLEL DO PRIVATE(p,c)
2019-01-29 09:41:29 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-06-14 12:19:05 +05:30
p = material_phaseAt ( g , e ) ; c = material_phaseMemberAt ( g , i , e )
2019-01-30 04:31:40 +05:30
2019-01-30 04:41:10 +05:30
plasticState ( p ) % RK4dotState ( : , c ) = WEIGHT ( n ) * plasticState ( p ) % dotState ( : , c ) &
2019-04-03 16:02:30 +05:30
+ merge ( plasticState ( p ) % RK4dotState ( : , c ) , 0.0_pReal , n > 1 )
do s = 1 , phase_Nsources ( p )
2019-01-30 04:41:10 +05:30
sourceState ( p ) % p ( s ) % RK4dotState ( : , c ) = WEIGHT ( n ) * sourceState ( p ) % p ( s ) % dotState ( : , c ) &
2019-04-03 16:02:30 +05:30
+ merge ( sourceState ( p ) % p ( s ) % RK4dotState ( : , c ) , 0.0_pReal , n > 1 )
2019-01-15 08:57:57 +05:30
enddo
endif
enddo ; enddo ; enddo
2019-01-30 04:41:10 +05:30
!$OMP END PARALLEL DO
2019-01-15 08:57:57 +05:30
2019-01-23 16:21:43 +05:30
call update_state ( TIMESTEPFRACTION ( n ) )
2019-01-24 16:03:04 +05:30
call update_deltaState
2019-01-23 22:34:19 +05:30
call update_dependentState
2019-01-24 22:29:38 +05:30
call update_stress ( TIMESTEPFRACTION ( n ) )
2019-01-15 08:57:57 +05:30
! --- dot state and RK dot state---
first3steps : if ( n < 4 ) then
2019-01-30 04:41:10 +05:30
call update_dotState ( TIMESTEPFRACTION ( n ) )
2019-01-15 08:57:57 +05:30
endif first3steps
enddo
2019-01-29 09:41:29 +05:30
call setConvergenceFlag
if ( any ( plasticState ( : ) % nonlocal ) ) call nonlocalConvergenceCheck
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
end subroutine integrateStateRK4
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with
!> adaptive step size (use 5th order solution to advance = "local extrapolation")
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine integrateStateRKCK45
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( 5 , 5 ) , parameter :: &
A = reshape ( [ &
2019-02-27 02:01:47 +05:30
. 2_pReal , . 075_pReal , . 3_pReal , - 1 1.0_pReal / 5 4.0_pReal , 163 1.0_pReal / 5529 6.0_pReal , &
. 0_pReal , . 225_pReal , - . 9_pReal , 2.5_pReal , 17 5.0_pReal / 51 2.0_pReal , &
. 0_pReal , . 0_pReal , 1.2_pReal , - 7 0.0_pReal / 2 7.0_pReal , 57 5.0_pReal / 1382 4.0_pReal , &
. 0_pReal , . 0_pReal , . 0_pReal , 3 5.0_pReal / 2 7.0_pReal , 4427 5.0_pReal / 11059 2.0_pReal , &
. 0_pReal , . 0_pReal , . 0_pReal , . 0_pReal , 25 3.0_pReal / 409 6.0_pReal ] , &
2019-01-15 08:57:57 +05:30
[ 5 , 5 ] , order = [ 2 , 1 ] ) !< coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( 6 ) , parameter :: &
B = &
[ 3 7.0_pReal / 37 8.0_pReal , . 0_pReal , 25 0.0_pReal / 62 1.0_pReal , &
12 5.0_pReal / 59 4.0_pReal , . 0_pReal , 51 2.0_pReal / 177 1.0_pReal ] , & !< coefficients in Butcher tableau (used for final integration and error estimate)
DB = B - &
2019-02-27 02:01:47 +05:30
[ 282 5.0_pReal / 2764 8.0_pReal , . 0_pReal , 1857 5.0_pReal / 4838 4.0_pReal , &
1352 5.0_pReal / 5529 6.0_pReal , 27 7.0_pReal / 1433 6.0_pReal , 0.25_pReal ] !< coefficients in Butcher tableau (used for final integration and error estimate)
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( 5 ) , parameter :: &
C = [ 0.2_pReal , 0.3_pReal , 0.6_pReal , 1.0_pReal , 0.875_pReal ] !< coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-15 08:57:57 +05:30
e , & ! element index in element loop
i , & ! integration point index in ip loop
g , & ! grain index in grain loop
stage , & ! stage index in integration stage loop
n , &
p , &
cc , &
2019-01-30 15:07:18 +05:30
s , &
2019-01-30 15:18:59 +05:30
sizeDotState
2019-01-29 09:41:29 +05:30
2019-01-30 15:34:49 +05:30
! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of RKCK45
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState , &
2019-06-07 02:19:17 +05:30
homogenization_maxNgrains , discretization_nIP , discretization_nElem ) :: &
2019-01-30 15:34:49 +05:30
residuum_plastic ! relative residuum from evolution in microstructure
2019-01-15 08:57:57 +05:30
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , &
maxval ( phase_Nsources ) , &
2019-06-07 02:19:17 +05:30
homogenization_maxNgrains , discretization_nIP , discretization_nElem ) :: &
2019-01-30 15:34:49 +05:30
residuum_source ! relative residuum from evolution in microstructure
2014-11-01 00:33:08 +05:30
2014-11-12 22:10:50 +05:30
2019-01-23 10:44:19 +05:30
call update_dotState ( 1.0_pReal )
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
! --- SECOND TO SIXTH RUNGE KUTTA STEP ---
2014-11-01 00:33:08 +05:30
2019-04-03 16:02:30 +05:30
do stage = 1 , 5
2014-11-12 22:10:50 +05:30
2019-01-15 08:57:57 +05:30
! --- state update ---
2019-01-30 13:26:16 +05:30
!$OMP PARALLEL DO PRIVATE(p,cc)
2019-01-29 09:41:29 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-06-14 12:19:05 +05:30
p = material_phaseAt ( g , e ) ; cc = material_phaseMemberAt ( g , i , e )
2014-11-12 22:10:50 +05:30
2019-01-30 13:26:16 +05:30
plasticState ( p ) % RKCK45dotState ( stage , : , cc ) = plasticState ( p ) % dotState ( : , cc )
2019-01-15 08:57:57 +05:30
plasticState ( p ) % dotState ( : , cc ) = A ( 1 , stage ) * plasticState ( p ) % RKCK45dotState ( 1 , : , cc )
2019-01-30 13:26:16 +05:30
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 15:07:18 +05:30
sourceState ( p ) % p ( s ) % RKCK45dotState ( stage , : , cc ) = sourceState ( p ) % p ( s ) % dotState ( : , cc )
sourceState ( p ) % p ( s ) % dotState ( : , cc ) = A ( 1 , stage ) * sourceState ( p ) % p ( s ) % RKCK45dotState ( 1 , : , cc )
2019-01-15 08:57:57 +05:30
enddo
2019-01-30 13:26:16 +05:30
2019-04-03 16:02:30 +05:30
do n = 2 , stage
2019-01-30 13:26:16 +05:30
plasticState ( p ) % dotState ( : , cc ) = plasticState ( p ) % dotState ( : , cc ) &
+ A ( n , stage ) * plasticState ( p ) % RKCK45dotState ( n , : , cc )
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 15:07:18 +05:30
sourceState ( p ) % p ( s ) % dotState ( : , cc ) = sourceState ( p ) % p ( s ) % dotState ( : , cc ) &
2019-01-30 15:16:53 +05:30
+ A ( n , stage ) * sourceState ( p ) % p ( s ) % RKCK45dotState ( n , : , cc )
2019-01-15 08:57:57 +05:30
enddo
enddo
2019-01-30 13:26:16 +05:30
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
2019-01-30 13:26:16 +05:30
!$OMP END PARALLEL DO
2014-11-01 00:33:08 +05:30
2019-01-23 16:21:43 +05:30
call update_state ( 1.0_pReal ) !MD: 1.0 correct?
2019-01-24 16:03:04 +05:30
call update_deltaState
call update_dependentState
2019-01-24 22:29:38 +05:30
call update_stress ( C ( stage ) )
call update_dotState ( C ( stage ) )
2014-08-26 20:14:32 +05:30
2019-01-15 08:57:57 +05:30
enddo
2014-08-26 20:14:32 +05:30
2015-03-06 18:39:00 +05:30
2019-01-15 08:57:57 +05:30
!--------------------------------------------------------------------------------------------------
! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE ---
2019-01-30 15:18:59 +05:30
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc)
2019-01-29 09:41:29 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-01-15 08:57:57 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-06-14 12:19:05 +05:30
p = material_phaseAt ( g , e ) ; cc = material_phaseMemberAt ( g , i , e )
2020-02-13 21:40:27 +05:30
2019-01-30 15:18:59 +05:30
sizeDotState = plasticState ( p ) % sizeDotState
2020-02-13 21:40:27 +05:30
2019-01-30 13:51:33 +05:30
plasticState ( p ) % RKCK45dotState ( 6 , : , cc ) = plasticState ( p ) % dotState ( : , cc )
2020-02-13 21:40:27 +05:30
2019-01-30 15:21:24 +05:30
residuum_plastic ( 1 : sizeDotState , g , i , e ) = &
2019-01-31 09:47:46 +05:30
matmul ( transpose ( plasticState ( p ) % RKCK45dotState ( 1 : 6 , 1 : sizeDotState , cc ) ) , DB ) & ! why transpose? Better to transpose constant DB
2019-01-15 08:57:57 +05:30
* crystallite_subdt ( g , i , e )
2020-02-13 21:40:27 +05:30
2019-01-30 13:51:33 +05:30
plasticState ( p ) % dotState ( : , cc ) = &
2019-01-31 09:47:46 +05:30
matmul ( transpose ( plasticState ( p ) % RKCK45dotState ( 1 : 6 , 1 : sizeDotState , cc ) ) , B ) ! why transpose? Better to transpose constant B
2020-02-13 21:40:27 +05:30
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 15:18:59 +05:30
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
2020-02-13 21:40:27 +05:30
2019-01-30 19:22:12 +05:30
sourceState ( p ) % p ( s ) % RKCK45dotState ( 6 , : , cc ) = sourceState ( p ) % p ( s ) % dotState ( : , cc )
2020-02-13 21:40:27 +05:30
2019-01-30 15:21:24 +05:30
residuum_source ( 1 : sizeDotState , s , g , i , e ) = &
2019-01-30 15:18:59 +05:30
matmul ( transpose ( sourceState ( p ) % p ( s ) % RKCK45dotState ( 1 : 6 , 1 : sizeDotState , cc ) ) , DB ) &
2019-01-15 08:57:57 +05:30
* crystallite_subdt ( g , i , e )
2014-08-26 20:14:32 +05:30
2019-01-30 15:07:18 +05:30
sourceState ( p ) % p ( s ) % dotState ( : , cc ) = &
2019-01-30 15:18:59 +05:30
matmul ( transpose ( sourceState ( p ) % p ( s ) % RKCK45dotState ( 1 : 6 , 1 : sizeDotState , cc ) ) , B )
2019-01-15 08:57:57 +05:30
enddo
2020-02-13 21:40:27 +05:30
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
2019-01-30 13:51:33 +05:30
!$OMP END PARALLEL DO
2014-08-26 20:14:32 +05:30
2019-01-23 16:21:43 +05:30
call update_state ( 1.0_pReal )
2020-02-13 21:40:27 +05:30
2019-01-15 08:57:57 +05:30
! --- relative residui and state convergence ---
2014-08-26 20:14:32 +05:30
2019-01-30 15:34:49 +05:30
!$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc)
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-01-30 15:34:49 +05:30
if ( crystallite_todo ( g , i , e ) ) then
2019-06-14 12:19:05 +05:30
p = material_phaseAt ( g , e ) ; cc = material_phaseMemberAt ( g , i , e )
2020-02-13 21:40:27 +05:30
2019-01-30 15:34:49 +05:30
sizeDotState = plasticState ( p ) % sizeDotState
2020-02-13 21:40:27 +05:30
2019-01-31 13:42:44 +05:30
crystallite_todo ( g , i , e ) = converged ( residuum_plastic ( 1 : sizeDotState , g , i , e ) , &
2019-01-30 21:34:58 +05:30
plasticState ( p ) % state ( 1 : sizeDotState , cc ) , &
2020-03-15 14:21:40 +05:30
plasticState ( p ) % atol ( 1 : sizeDotState ) )
2014-08-26 20:14:32 +05:30
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-30 15:34:49 +05:30
sizeDotState = sourceState ( p ) % p ( s ) % sizeDotState
2020-02-13 21:40:27 +05:30
2019-02-27 02:01:47 +05:30
crystallite_todo ( g , i , e ) = &
crystallite_todo ( g , i , e ) . and . converged ( residuum_source ( 1 : sizeDotState , s , g , i , e ) , &
sourceState ( p ) % p ( s ) % state ( 1 : sizeDotState , cc ) , &
2020-03-15 14:21:40 +05:30
sourceState ( p ) % p ( s ) % atol ( 1 : sizeDotState ) )
2019-01-30 19:22:12 +05:30
enddo
2019-01-15 08:57:57 +05:30
endif
enddo ; enddo ; enddo
2019-01-30 15:34:49 +05:30
!$OMP END PARALLEL DO
2014-08-26 20:14:32 +05:30
2019-01-24 16:03:04 +05:30
call update_deltaState
2019-01-23 22:34:19 +05:30
call update_dependentState
2019-01-24 22:29:38 +05:30
call update_stress ( 1.0_pReal )
2019-01-25 11:50:05 +05:30
call setConvergenceFlag
2019-01-29 09:41:29 +05:30
if ( any ( plasticState ( : ) % nonlocal ) ) call nonlocalConvergenceCheck
2020-02-13 21:40:27 +05:30
2019-01-29 09:41:29 +05:30
end subroutine integrateStateRKCK45
2011-05-11 22:08:45 +05:30
2014-08-26 20:14:32 +05:30
2019-01-29 09:41:29 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief sets convergence flag for nonlocal calculations
2020-03-15 17:16:10 +05:30
!> @details one non-converged nonlocal sets all other nonlocals to non-converged to trigger cut back
2019-01-29 09:41:29 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-11 11:05:58 +05:30
subroutine nonlocalConvergenceCheck
2020-02-13 21:40:27 +05:30
2019-01-29 09:41:29 +05:30
if ( any ( . not . crystallite_converged . and . . not . crystallite_localPlasticity ) ) & ! any non-local not yet converged (or broken)...
where ( . not . crystallite_localPlasticity ) crystallite_converged = . false .
2014-06-25 04:51:25 +05:30
2019-01-29 09:41:29 +05:30
end subroutine nonlocalConvergenceCheck
2009-05-07 21:57:36 +05:30
2019-01-19 14:05:45 +05:30
2019-01-25 11:50:05 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Sets convergence flag based on "todo": every point that survived the integration (todo is
! still .true. is considered as converged
!> @details: For explicitEuler, RK4 and RKCK45, adaptive Euler and FPI have their on criteria
!--------------------------------------------------------------------------------------------------
2019-04-11 11:05:58 +05:30
subroutine setConvergenceFlag
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-25 11:50:05 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
g !< grain index in grain loop
2020-02-13 21:40:27 +05:30
2019-04-03 16:47:21 +05:30
!OMP DO PARALLEL PRIVATE
2019-01-25 11:50:05 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-04-03 16:47:21 +05:30
crystallite_converged ( g , i , e ) = crystallite_todo ( g , i , e ) . or . crystallite_converged ( g , i , e ) ! if still "to do" then converged per definition
enddo ; enddo ; enddo
2019-01-25 11:50:05 +05:30
!OMP END DO PARALLEL
end subroutine setConvergenceFlag
2019-01-31 13:44:02 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief determines whether a point is converged
!--------------------------------------------------------------------------------------------------
2020-03-14 23:41:26 +05:30
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
end function converged
2019-01-24 22:29:38 +05:30
!--------------------------------------------------------------------------------------------------
2019-10-11 18:51:29 +05:30
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t) comment seems wrong!
2019-01-24 22:29:38 +05:30
!--------------------------------------------------------------------------------------------------
subroutine update_stress ( timeFraction )
2019-04-11 11:05:58 +05:30
2020-03-23 04:45:00 +05:30
real ( pReal ) , intent ( in ) :: &
timeFraction
integer :: &
2020-03-23 11:09:17 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
2020-03-23 04:45:00 +05:30
g
logical :: &
2020-03-23 11:09:17 +05:30
nonlocalBroken
2019-01-24 22:29:38 +05:30
2020-03-23 11:09:17 +05:30
nonlocalBroken = . false .
2020-03-23 04:45:00 +05:30
!$OMP PARALLEL DO
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2020-03-23 11:09:17 +05:30
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) . and . &
( . not . nonlocalBroken . or . crystallite_localPlasticity ( g , i , e ) ) ) then
2020-03-23 04:45:00 +05:30
crystallite_todo ( g , i , e ) = integrateStress ( g , i , e , timeFraction )
if ( . not . ( crystallite_todo ( g , i , e ) . or . crystallite_localPlasticity ( g , i , e ) ) ) &
2020-03-23 11:09:17 +05:30
nonlocalBroken = . true .
2020-03-23 04:45:00 +05:30
endif
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
2020-03-23 11:09:17 +05:30
if ( nonlocalBroken ) where ( . not . crystallite_localPlasticity ) crystallite_todo = . false .
2019-01-24 22:29:38 +05:30
end subroutine update_stress
2020-03-23 11:09:17 +05:30
2019-01-23 22:34:19 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief tbd
!--------------------------------------------------------------------------------------------------
2019-04-11 14:57:03 +05:30
subroutine update_dependentState
2019-04-03 16:02:30 +05:30
integer :: e , & ! element index in element loop
2019-01-23 22:34:19 +05:30
i , & ! integration point index in ip loop
g ! grain index in grain loop
2019-01-24 03:48:14 +05:30
!$OMP PARALLEL DO
2019-01-23 22:34:19 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-01-23 22:34:19 +05:30
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) &
2019-02-02 00:58:07 +05:30
call constitutive_dependentState ( crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
2019-01-23 22:34:19 +05:30
crystallite_Fp ( 1 : 3 , 1 : 3 , g , i , e ) , &
g , i , e )
enddo ; enddo ; enddo
2019-01-24 03:48:14 +05:30
!$OMP END PARALLEL DO
2019-01-23 22:34:19 +05:30
end subroutine update_dependentState
2019-01-24 22:29:38 +05:30
2019-01-23 16:21:43 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Standard forwarding of state as state = state0 + dotState * (delta t)
!--------------------------------------------------------------------------------------------------
subroutine update_state ( timeFraction )
real ( pReal ) , intent ( in ) :: &
timeFraction
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-23 16:21:43 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
g , & !< grain index in grain loop
p , &
c , &
s , &
mySize
!$OMP PARALLEL DO PRIVATE(mySize,p,c)
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-01-23 16:21:43 +05:30
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) then
2019-06-14 12:19:05 +05:30
p = material_phaseAt ( g , e ) ; c = material_phaseMemberAt ( g , i , e )
2019-01-23 16:21:43 +05:30
mySize = plasticState ( p ) % sizeDotState
plasticState ( p ) % state ( 1 : mySize , c ) = plasticState ( p ) % subState0 ( 1 : mySize , c ) &
+ plasticState ( p ) % dotState ( 1 : mySize , c ) &
* crystallite_subdt ( g , i , e ) * timeFraction
2019-04-03 16:02:30 +05:30
do s = 1 , phase_Nsources ( p )
2019-01-23 16:21:43 +05:30
mySize = sourceState ( p ) % p ( s ) % sizeDotState
sourceState ( p ) % p ( s ) % state ( 1 : mySize , c ) = sourceState ( p ) % p ( s ) % subState0 ( 1 : mySize , c ) &
+ sourceState ( p ) % p ( s ) % dotState ( 1 : mySize , c ) &
* crystallite_subdt ( g , i , e ) * timeFraction
enddo
endif
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
end subroutine update_state
2020-03-23 11:09:17 +05:30
!---------------------------------------------------------------------------------------------------
!> @brief Trigger calculation of all new rates
2019-01-23 10:44:19 +05:30
!> if NaN occurs, crystallite_todo is set to FALSE. Any NaN in a nonlocal propagates to all others
2020-03-23 11:09:17 +05:30
!---------------------------------------------------------------------------------------------------
2019-01-23 10:44:19 +05:30
subroutine update_dotState ( timeFraction )
2020-03-23 11:09:17 +05:30
real ( pReal ) , intent ( in ) :: &
timeFraction
integer :: &
e , & !< element index in element loop
i , & !< integration point index in ip loop
g , & !< grain index in grain loop
p , &
c , &
s
logical :: &
nonlocalBroken
2019-01-23 10:44:19 +05:30
2020-03-23 11:09:17 +05:30
nonlocalBroken = . false .
!$OMP PARALLEL DO PRIVATE (p,c)
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
if ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) . and . &
( . not . nonlocalBroken . or . crystallite_localPlasticity ( g , i , e ) ) ) then
call constitutive_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_partionedF0 , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_partionedFp0 , &
crystallite_subdt ( g , i , e ) * timeFraction , g , i , e )
p = material_phaseAt ( g , e ) ; c = material_phaseMemberAt ( g , i , e )
crystallite_todo ( g , i , e ) = all ( . not . IEEE_is_NaN ( plasticState ( p ) % dotState ( : , c ) ) )
do s = 1 , phase_Nsources ( p )
crystallite_todo ( g , i , e ) = crystallite_todo ( g , i , e ) . and . all ( . not . IEEE_is_NaN ( sourceState ( p ) % p ( s ) % dotState ( : , c ) ) )
enddo
if ( . not . ( crystallite_todo ( g , i , e ) . or . crystallite_localPlasticity ( g , i , e ) ) ) &
nonlocalBroken = . true .
endif
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
2019-01-24 18:45:26 +05:30
2020-03-23 11:09:17 +05:30
if ( nonlocalBroken ) where ( . not . crystallite_localPlasticity ) crystallite_todo = . false .
2019-01-23 10:44:19 +05:30
end subroutine update_DotState
2020-03-23 11:09:17 +05:30
!---------------------------------------------------------------------------------------------------
!> @brief Trigger calculation of all new sudden state change
!> if NaN occurs, crystallite_todo is set to FALSE. Any NaN in a nonlocal propagates to all others
!---------------------------------------------------------------------------------------------------
2019-01-24 11:42:20 +05:30
subroutine update_deltaState
2019-05-17 02:26:48 +05:30
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-24 11:42:20 +05:30
e , & !< element index in element loop
i , & !< integration point index in ip loop
g , & !< grain index in grain loop
p , &
2019-01-30 04:41:10 +05:30
mySize , &
2019-01-24 12:04:30 +05:30
myOffset , &
2019-01-24 11:42:20 +05:30
c , &
2020-02-13 21:40:27 +05:30
s
2019-01-24 18:45:26 +05:30
logical :: &
NaN , &
nonlocalStop
2020-02-13 21:40:27 +05:30
2019-01-24 18:45:26 +05:30
nonlocalStop = . false .
2019-01-24 11:42:20 +05:30
2019-01-30 04:41:10 +05:30
!$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize,NaN)
2019-01-24 11:42:20 +05:30
do e = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
2020-01-25 13:54:42 +05:30
do i = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2019-06-05 13:35:59 +05:30
do g = 1 , homogenization_Ngrains ( material_homogenizationAt ( e ) )
2019-01-24 18:45:26 +05:30
!$OMP FLUSH(nonlocalStop)
2019-01-30 20:36:14 +05:30
if ( ( crystallite_todo ( g , i , e ) . and . . not . crystallite_converged ( g , i , e ) ) . and . . not . nonlocalStop ) then
2019-03-09 05:03:11 +05:30
call constitutive_collectDeltaState ( crystallite_S ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fe ( 1 : 3 , 1 : 3 , g , i , e ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , g , i , e ) , &
g , i , e )
2019-06-14 12:19:05 +05:30
p = material_phaseAt ( g , e ) ; c = material_phaseMemberAt ( g , i , e )
2019-01-24 12:04:30 +05:30
myOffset = plasticState ( p ) % offsetDeltaState
mySize = plasticState ( p ) % sizeDeltaState
NaN = any ( IEEE_is_NaN ( plasticState ( p ) % deltaState ( 1 : mySize , c ) ) )
2020-02-13 21:40:27 +05:30
2019-01-24 12:04:30 +05:30
if ( . not . NaN ) then
2020-02-13 21:40:27 +05:30
2019-04-03 16:02:30 +05:30
plasticState ( p ) % state ( myOffset + 1 : myOffset + mySize , c ) = &
plasticState ( p ) % state ( myOffset + 1 : myOffset + mySize , c ) + plasticState ( p ) % deltaState ( 1 : mySize , c )
do s = 1 , phase_Nsources ( p )
2019-01-30 04:41:10 +05:30
myOffset = sourceState ( p ) % p ( s ) % offsetDeltaState
mySize = sourceState ( p ) % p ( s ) % sizeDeltaState
NaN = NaN . or . any ( IEEE_is_NaN ( sourceState ( p ) % p ( s ) % deltaState ( 1 : mySize , c ) ) )
2020-02-13 21:40:27 +05:30
2019-01-24 16:03:04 +05:30
if ( . not . NaN ) then
2019-04-03 16:02:30 +05:30
sourceState ( p ) % p ( s ) % state ( myOffset + 1 : myOffset + mySize , c ) = &
sourceState ( p ) % p ( s ) % state ( myOffset + 1 : myOffset + mySize , c ) + sourceState ( p ) % p ( s ) % deltaState ( 1 : mySize , c )
2019-01-24 16:03:04 +05:30
endif
enddo
2019-01-24 12:04:30 +05:30
endif
2020-02-13 21:40:27 +05:30
2019-01-24 16:03:04 +05:30
crystallite_todo ( g , i , e ) = . not . NaN
2019-01-24 11:42:20 +05:30
if ( . not . crystallite_todo ( g , i , e ) ) then ! if state jump fails, then convergence is broken
crystallite_converged ( g , i , e ) = . false .
2019-01-24 18:45:26 +05:30
if ( . not . crystallite_localPlasticity ( g , i , e ) ) nonlocalStop = . true .
2019-01-24 11:42:20 +05:30
endif
endif
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
2019-01-24 18:45:26 +05:30
if ( nonlocalStop ) crystallite_todo = crystallite_todo . and . crystallite_localPlasticity
2020-02-13 21:40:27 +05:30
2019-01-24 11:42:20 +05:30
end subroutine update_deltaState
2019-01-19 14:05:45 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates a jump in the state according to the current state and the current stress
!> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
!--------------------------------------------------------------------------------------------------
logical function stateJump ( ipc , ip , el )
2019-04-03 16:02:30 +05:30
integer , intent ( in ) :: &
2019-01-19 14:05:45 +05:30
el , & ! element index
ip , & ! integration point index
ipc ! grain index
2019-04-03 16:02:30 +05:30
integer :: &
2019-01-19 14:05:45 +05:30
c , &
p , &
mySource , &
2019-01-24 11:42:20 +05:30
myOffset , &
mySize
2019-01-19 14:05:45 +05:30
2019-06-14 12:12:12 +05:30
c = material_phaseMemberAt ( ipc , ip , el )
2019-06-14 12:19:05 +05:30
p = material_phaseAt ( ipc , el )
2019-01-19 14:05:45 +05:30
2019-03-09 05:03:11 +05:30
call constitutive_collectDeltaState ( crystallite_S ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
crystallite_Fe ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
crystallite_Fi ( 1 : 3 , 1 : 3 , ipc , ip , el ) , &
ipc , ip , el )
2019-01-19 14:05:45 +05:30
2019-01-24 11:42:20 +05:30
myOffset = plasticState ( p ) % offsetDeltaState
mySize = plasticState ( p ) % sizeDeltaState
2019-01-19 14:05:45 +05:30
2019-01-24 11:42:20 +05:30
if ( any ( IEEE_is_NaN ( plasticState ( p ) % deltaState ( 1 : mySize , c ) ) ) ) then ! NaN occured in deltaState
2019-01-19 14:05:45 +05:30
stateJump = . false .
return
endif
2019-04-03 16:02:30 +05:30
plasticState ( p ) % state ( myOffset + 1 : myOffset + mySize , c ) = &
plasticState ( p ) % state ( myOffset + 1 : myOffset + mySize , c ) + plasticState ( p ) % deltaState ( 1 : mySize , c )
2019-01-19 14:05:45 +05:30
2019-04-03 16:02:30 +05:30
do mySource = 1 , phase_Nsources ( p )
2019-01-24 11:42:20 +05:30
myOffset = sourceState ( p ) % p ( mySource ) % offsetDeltaState
2019-01-25 11:50:05 +05:30
mySize = sourceState ( p ) % p ( mySource ) % sizeDeltaState
2019-01-24 11:42:20 +05:30
if ( any ( IEEE_is_NaN ( sourceState ( p ) % p ( mySource ) % deltaState ( 1 : mySize , c ) ) ) ) then ! NaN occured in deltaState
2019-01-19 14:05:45 +05:30
stateJump = . false .
return
endif
2019-04-03 16:02:30 +05:30
sourceState ( p ) % p ( mySource ) % state ( myOffset + 1 : myOffset + mySize , c ) = &
sourceState ( p ) % p ( mySource ) % state ( myOffset + 1 : myOffset + mySize , c ) + sourceState ( p ) % p ( mySource ) % deltaState ( 1 : mySize , c )
2019-01-19 14:05:45 +05:30
enddo
#ifdef DEBUG
2019-01-24 11:42:20 +05:30
if ( any ( dNeq0 ( plasticState ( p ) % deltaState ( 1 : mySize , c ) ) ) &
2019-04-03 16:02:30 +05:30
. and . iand ( debug_level ( debug_crystallite ) , debug_levelExtensive ) / = 0 &
2019-01-19 14:05:45 +05:30
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
2019-04-03 16:02:30 +05:30
. or . . not . iand ( debug_level ( debug_crystallite ) , debug_levelSelective ) / = 0 ) ) then
2019-01-19 14:05:45 +05:30
write ( 6 , '(a,i8,1x,i2,1x,i3, /)' ) '<< CRYST >> update state at el ip ipc ' , el , ip , ipc
2019-01-24 11:42:20 +05:30
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> deltaState' , plasticState ( p ) % deltaState ( 1 : mySize , c )
2019-01-19 14:05:45 +05:30
write ( 6 , '(a,/,(12x,12(e12.5,1x)),/)' ) '<< CRYST >> new state' , &
2019-04-03 16:02:30 +05:30
plasticState ( p ) % state ( myOffset + 1 : &
2019-01-24 11:42:20 +05:30
myOffset + mySize , c )
2019-01-19 14:05:45 +05:30
endif
#endif
stateJump = . true .
end function stateJump
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
write ( 6 , '(a)' ) ' writing field and constitutive data required for restart to file' ; flush ( 6 )
write ( fileName , '(a,i0,a)' ) trim ( getSolverJobName ( ) ) / / '_' , worldrank , '.hdf5'
fileHandle = HDF5_openFile ( fileName , 'a' )
call HDF5_write ( fileHandle , crystallite_partionedF , 'F' )
call HDF5_write ( fileHandle , crystallite_Fp , 'Fp' )
call HDF5_write ( fileHandle , crystallite_Fi , 'Fi' )
call HDF5_write ( fileHandle , crystallite_Lp , 'Lp' )
call HDF5_write ( fileHandle , crystallite_Li , 'Li' )
call HDF5_write ( fileHandle , crystallite_S , 'S' )
groupHandle = HDF5_addGroup ( fileHandle , 'constituent' )
do i = 1 , size ( phase_plasticity )
write ( datasetName , '(i0,a)' ) i , '_omega_plastic'
call HDF5_write ( groupHandle , plasticState ( i ) % state , datasetName )
enddo
call HDF5_closeGroup ( groupHandle )
groupHandle = HDF5_addGroup ( fileHandle , 'materialpoint' )
do i = 1 , material_Nhomogenization
write ( datasetName , '(i0,a)' ) i , '_omega_homogenization'
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
write ( 6 , '(/,a,i0,a)' ) ' reading restart information of increment from file'
write ( fileName , '(a,i0,a)' ) trim ( getSolverJobName ( ) ) / / '_' , worldrank , '.hdf5'
fileHandle = HDF5_openFile ( fileName )
call HDF5_read ( fileHandle , crystallite_F0 , 'F' )
call HDF5_read ( fileHandle , crystallite_Fp0 , 'Fp' )
call HDF5_read ( fileHandle , crystallite_Fi0 , 'Fi' )
call HDF5_read ( fileHandle , crystallite_Lp0 , 'Lp' )
call HDF5_read ( fileHandle , crystallite_Li0 , 'Li' )
call HDF5_read ( fileHandle , crystallite_S0 , 'S' )
groupHandle = HDF5_openGroup ( fileHandle , 'constituent' )
do i = 1 , size ( phase_plasticity )
write ( datasetName , '(i0,a)' ) i , '_omega_plastic'
call HDF5_read ( groupHandle , plasticState ( i ) % state0 , datasetName )
enddo
call HDF5_closeGroup ( groupHandle )
groupHandle = HDF5_openGroup ( fileHandle , 'materialpoint' )
do i = 1 , material_Nhomogenization
write ( datasetName , '(i0,a)' ) i , '_omega_homogenization'
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
crystallite_F0 = crystallite_partionedF
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
do i = 1 , material_Nhomogenization
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
