2012-10-02 18:23:25 +05:30
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
2013-03-06 20:11:15 +05:30
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
2020-07-15 18:05:21 +05:30
!> @brief elasticity, plasticity, damage & thermal internal microstructure state
2012-10-02 18:23:25 +05:30
!--------------------------------------------------------------------------------------------------
module constitutive
2019-12-03 01:13:02 +05:30
use prec
2019-06-15 19:10:22 +05:30
use math
2019-12-05 01:20:46 +05:30
use rotations
2019-06-15 19:10:22 +05:30
use IO
use config
use material
use results
use lattice
use discretization
2020-12-20 15:18:13 +05:30
use parallelization
use HDF5_utilities
use DAMASK_interface
use FEsolving
use results
2020-02-07 16:53:22 +05:30
2019-06-15 19:10:22 +05:30
implicit none
private
2020-12-23 22:00:19 +05:30
2020-12-23 12:50:29 +05:30
enum , bind ( c ) ; enumerator :: &
PLASTICITY_UNDEFINED_ID , &
PLASTICITY_NONE_ID , &
PLASTICITY_ISOTROPIC_ID , &
PLASTICITY_PHENOPOWERLAW_ID , &
PLASTICITY_KINEHARDENING_ID , &
PLASTICITY_DISLOTWIN_ID , &
PLASTICITY_DISLOTUNGSTEN_ID , &
PLASTICITY_NONLOCAL_ID , &
SOURCE_UNDEFINED_ID , &
SOURCE_THERMAL_DISSIPATION_ID , &
SOURCE_THERMAL_EXTERNALHEAT_ID , &
SOURCE_DAMAGE_ISOBRITTLE_ID , &
SOURCE_DAMAGE_ISODUCTILE_ID , &
SOURCE_DAMAGE_ANISOBRITTLE_ID , &
SOURCE_DAMAGE_ANISODUCTILE_ID , &
KINEMATICS_UNDEFINED_ID , &
KINEMATICS_CLEAVAGE_OPENING_ID , &
KINEMATICS_SLIPPLANE_OPENING_ID , &
2020-12-23 18:33:15 +05:30
KINEMATICS_THERMAL_EXPANSION_ID
2020-12-23 12:50:29 +05:30
end enum
2020-12-20 15:18:13 +05:30
real ( pReal ) , dimension ( : , : , : ) , allocatable :: &
2020-12-23 18:01:30 +05:30
crystallite_subdt !< substepped time increment of each grain
2020-12-20 15:18:13 +05:30
type ( rotation ) , dimension ( : , : , : ) , allocatable :: &
crystallite_orientation !< current orientation
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: &
crystallite_F0 , & !< def grad at start of FE inc
crystallite_subF , & !< def grad to be reached at end of crystallite inc
crystallite_Fe , & !< current "elastic" def grad (end of converged time step)
crystallite_subFp0 , & !< plastic def grad at start of crystallite inc
crystallite_subFi0 , & !< intermediate def grad at start of crystallite inc
crystallite_Lp0 , & !< plastic velocitiy grad at start of FE inc
crystallite_partitionedLp0 , & !< plastic velocity grad at start of homog inc
crystallite_S0 , & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
crystallite_partitionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc
2020-12-21 19:21:55 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , public :: &
2020-12-20 15:18:13 +05:30
crystallite_P , & !< 1st Piola-Kirchhoff stress per grain
crystallite_Lp , & !< current plastic velocitiy grad (end of converged time step)
crystallite_S , & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
crystallite_partitionedF0 !< def grad at start of homog inc
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable , public :: &
crystallite_partitionedF !< def grad to be reached at end of homog inc
logical , dimension ( : , : , : ) , allocatable :: &
crystallite_converged !< convergence flag
2020-12-22 15:14:43 +05:30
2020-12-21 00:50:39 +05:30
type :: tTensorContainer
real ( pReal ) , dimension ( : , : , : ) , allocatable :: data
end type
type ( tTensorContainer ) , dimension ( : ) , allocatable :: &
constitutive_mech_Fi , &
constitutive_mech_Fi0 , &
2020-12-24 16:24:09 +05:30
constitutive_mech_partitionedFi0 , &
2020-12-21 14:29:13 +05:30
constitutive_mech_Li , &
constitutive_mech_Li0 , &
2020-12-24 16:24:09 +05:30
constitutive_mech_partitionedLi0 , &
2020-12-23 02:51:11 +05:30
constitutive_mech_Fp , &
constitutive_mech_Fp0 , &
2020-12-24 16:24:09 +05:30
constitutive_mech_partitionedFp0
2020-12-23 02:51:11 +05:30
2020-12-21 00:50:39 +05:30
2020-12-20 15:18:13 +05:30
type :: tNumerics
integer :: &
iJacoLpresiduum , & !< frequency of Jacobian update of residuum in Lp
nState , & !< state loop limit
nStress !< stress loop limit
real ( pReal ) :: &
subStepMinCryst , & !< minimum (relative) size of sub-step allowed during cutback
subStepSizeCryst , & !< size of first substep when cutback
subStepSizeLp , & !< size of first substep when cutback in Lp calculation
subStepSizeLi , & !< size of first substep when cutback in Li calculation
stepIncreaseCryst , & !< increase of next substep size when previous substep converged
rtol_crystalliteState , & !< relative tolerance in state loop
rtol_crystalliteStress , & !< relative tolerance in stress loop
atol_crystalliteStress !< absolute tolerance in stress loop
end type tNumerics
type ( tNumerics ) :: num ! numerics parameters. Better name?
type :: tDebugOptions
logical :: &
basic , &
extensive , &
selective
integer :: &
element , &
ip , &
grain
end type tDebugOptions
type ( tDebugOptions ) :: debugCrystallite
2020-12-24 14:52:41 +05:30
2020-02-07 16:53:22 +05:30
2020-12-23 22:00:19 +05:30
integer ( kind ( PLASTICITY_undefined_ID ) ) , dimension ( : ) , allocatable , public :: &
2020-08-15 19:32:10 +05:30
phase_plasticity !< plasticity of each phase
2020-12-15 22:15:11 +05:30
integer ( kind ( SOURCE_undefined_ID ) ) , dimension ( : , : ) , allocatable :: &
2020-08-15 19:32:10 +05:30
phase_source , & !< active sources mechanisms of each phase
2020-12-15 22:15:11 +05:30
phase_kinematics !< active kinematic mechanisms of each phase
2020-08-15 19:32:10 +05:30
2020-12-15 22:15:11 +05:30
integer , dimension ( : ) , allocatable , public :: & !< ToDo: should be protected (bug in Intel compiler)
2020-08-15 19:32:10 +05:30
phase_Nsources , & !< number of source mechanisms active in each phase
phase_Nkinematics , & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations , & !< number of stiffness degradation mechanisms active in each phase
phase_plasticityInstance , & !< instance of particular plasticity of each phase
phase_elasticityInstance !< instance of particular elasticity of each phase
2020-12-15 22:15:11 +05:30
logical , dimension ( : ) , allocatable , public :: & ! ToDo: should be protected (bug in Intel Compiler)
2020-08-15 19:32:10 +05:30
phase_localPlasticity !< flags phases with local constitutive law
type ( tPlasticState ) , allocatable , dimension ( : ) , public :: &
plasticState
type ( tSourceState ) , allocatable , dimension ( : ) , public :: &
sourceState
2019-06-15 19:10:22 +05:30
integer , public , protected :: &
constitutive_plasticity_maxSizeDotState , &
constitutive_source_maxSizeDotState
2020-02-07 16:53:22 +05:30
2019-12-03 01:13:02 +05:30
interface
2020-12-21 16:44:09 +05:30
! == cleaned:begin =================================================================================
2020-11-03 02:09:33 +05:30
module subroutine mech_init
end subroutine mech_init
2020-02-07 16:53:22 +05:30
2020-07-09 04:31:08 +05:30
module subroutine damage_init
end subroutine damage_init
2020-07-03 20:15:11 +05:30
2020-07-09 04:31:08 +05:30
module subroutine thermal_init
end subroutine thermal_init
2020-02-07 16:53:22 +05:30
2020-12-21 16:44:09 +05:30
module subroutine mech_results ( group , ph )
character ( len = * ) , intent ( in ) :: group
integer , intent ( in ) :: ph
end subroutine mech_results
2020-12-22 16:50:00 +05:30
module subroutine damage_results ( group , ph )
character ( len = * ) , intent ( in ) :: group
integer , intent ( in ) :: ph
end subroutine damage_results
2020-12-21 16:44:09 +05:30
module subroutine mech_restart_read ( fileHandle )
integer ( HID_T ) , intent ( in ) :: fileHandle
end subroutine mech_restart_read
2020-12-22 14:33:19 +05:30
module subroutine mech_initializeRestorationPoints ( ph , me )
integer , intent ( in ) :: ph , me
end subroutine mech_initializeRestorationPoints
2020-12-23 03:57:56 +05:30
module subroutine constitutive_mech_windForward ( ph , me )
integer , intent ( in ) :: ph , me
end subroutine constitutive_mech_windForward
2020-12-23 11:28:54 +05:30
module subroutine constitutive_mech_forward
end subroutine constitutive_mech_forward
2020-07-10 20:40:23 +05:30
2020-12-23 11:28:54 +05:30
! == cleaned:end ===================================================================================
2020-12-20 21:41:43 +05:30
2020-12-24 14:52:41 +05:30
module function crystallite_stress ( dt , co , ip , el )
real ( pReal ) , intent ( in ) :: dt
integer , intent ( in ) :: co , ip , el
logical :: crystallite_stress
end function crystallite_stress
2020-12-23 22:00:19 +05:30
module function constitutive_homogenizedC ( co , ip , el ) result ( C )
integer , intent ( in ) :: co , ip , el
real ( pReal ) , dimension ( 6 , 6 ) :: C
end function constitutive_homogenizedC
2020-12-23 11:37:18 +05:30
module subroutine source_damage_anisoBrittle_dotState ( S , co , ip , el )
2020-07-10 20:40:23 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-07-12 20:14:26 +05:30
ip , & !< integration point
el !< element
2020-07-10 20:40:23 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
S
end subroutine source_damage_anisoBrittle_dotState
2020-12-23 11:37:18 +05:30
module subroutine source_damage_anisoDuctile_dotState ( co , ip , el )
2020-07-10 20:40:23 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-07-12 20:14:26 +05:30
ip , & !< integration point
el !< element
2020-07-10 20:40:23 +05:30
end subroutine source_damage_anisoDuctile_dotState
2020-12-23 11:37:18 +05:30
module subroutine source_damage_isoDuctile_dotState ( co , ip , el )
2020-07-10 20:40:23 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-07-12 20:14:26 +05:30
ip , & !< integration point
el !< element
2020-07-10 20:40:23 +05:30
end subroutine source_damage_isoDuctile_dotState
module subroutine source_thermal_externalheat_dotState ( phase , of )
integer , intent ( in ) :: &
phase , &
of
end subroutine source_thermal_externalheat_dotState
2020-07-15 18:05:21 +05:30
module subroutine constitutive_damage_getRateAndItsTangents ( phiDot , dPhiDot_dPhi , phi , ip , el )
integer , intent ( in ) :: &
ip , & !< integration point number
el !< element number
real ( pReal ) , intent ( in ) :: &
phi !< damage parameter
real ( pReal ) , intent ( inout ) :: &
phiDot , &
dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents
module subroutine constitutive_thermal_getRateAndItsTangents ( TDot , dTDot_dT , T , S , Lp , ip , el )
integer , intent ( in ) :: &
ip , & !< integration point number
el !< element number
real ( pReal ) , intent ( in ) :: &
T
real ( pReal ) , intent ( in ) , dimension ( : , : , : , : , : ) :: &
S , & !< current 2nd Piola Kitchoff stress vector
Lp !< plastic velocity gradient
real ( pReal ) , intent ( inout ) :: &
TDot , &
dTDot_dT
end subroutine constitutive_thermal_getRateAndItsTangents
2020-12-23 22:00:19 +05:30
2020-07-11 03:11:56 +05:30
2020-07-10 21:59:36 +05:30
module subroutine plastic_nonlocal_updateCompatibility ( orientation , instance , i , e )
integer , intent ( in ) :: &
instance , &
i , &
e
2020-10-28 01:57:26 +05:30
type ( rotation ) , dimension ( 1 , discretization_nIPs , discretization_Nelems ) , intent ( in ) :: &
2020-07-10 21:59:36 +05:30
orientation !< crystal orientation
end subroutine plastic_nonlocal_updateCompatibility
2020-07-02 00:52:05 +05:30
module subroutine plastic_isotropic_LiAndItsTangent ( Li , dLi_dMi , Mi , instance , of )
2019-12-03 01:13:02 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
Li !< inleastic velocity gradient
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( out ) :: &
dLi_dMi !< derivative of Li with respect to Mandel stress
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
2020-02-07 16:53:22 +05:30
Mi !< Mandel stress
2019-12-03 01:13:02 +05:30
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_isotropic_LiAndItsTangent
2020-02-07 16:53:22 +05:30
2020-12-23 11:37:18 +05:30
module subroutine kinematics_cleavage_opening_LiAndItsTangent ( Ld , dLd_dTstar , S , co , ip , el )
2020-07-09 04:31:08 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< grain number
2020-07-12 20:14:26 +05:30
ip , & !< integration point number
el !< element number
2020-07-09 04:31:08 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
S
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 ) :: &
2020-07-12 20:14:26 +05:30
Ld !< damage velocity gradient
2020-07-09 04:31:08 +05:30
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 , 3 , 3 ) :: &
2020-07-12 20:14:26 +05:30
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
2020-07-09 04:31:08 +05:30
end subroutine kinematics_cleavage_opening_LiAndItsTangent
2020-12-23 11:37:18 +05:30
module subroutine kinematics_slipplane_opening_LiAndItsTangent ( Ld , dLd_dTstar , S , co , ip , el )
2020-07-09 04:31:08 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< grain number
2020-07-12 20:14:26 +05:30
ip , & !< integration point number
el !< element number
2020-07-09 04:31:08 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
S
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 ) :: &
2020-07-12 20:14:26 +05:30
Ld !< damage velocity gradient
2020-07-09 04:31:08 +05:30
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 , 3 , 3 ) :: &
2020-07-12 20:14:26 +05:30
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
2020-07-09 04:31:08 +05:30
end subroutine kinematics_slipplane_opening_LiAndItsTangent
2020-12-23 11:37:18 +05:30
module subroutine kinematics_thermal_expansion_LiAndItsTangent ( Li , dLi_dTstar , co , ip , el )
2020-07-09 04:31:08 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< grain number
2020-07-12 20:14:26 +05:30
ip , & !< integration point number
el !< element number
2020-07-09 04:31:08 +05:30
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 ) :: &
2020-07-12 20:14:26 +05:30
Li !< thermal velocity gradient
2020-07-09 04:31:08 +05:30
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 , 3 , 3 ) :: &
2020-07-12 20:14:26 +05:30
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
2020-09-13 14:09:17 +05:30
end subroutine kinematics_thermal_expansion_LiAndItsTangent
2020-02-07 16:53:22 +05:30
2019-12-03 02:08:41 +05:30
2020-12-23 11:37:18 +05:30
module subroutine source_damage_isoBrittle_deltaState ( C , Fe , co , ip , el )
2019-12-05 01:20:46 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-07-12 20:14:26 +05:30
ip , & !< integration point
el !< element
2020-07-09 04:31:08 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
Fe
real ( pReal ) , intent ( in ) , dimension ( 6 , 6 ) :: &
C
end subroutine source_damage_isoBrittle_deltaState
2020-07-15 18:05:21 +05:30
module subroutine constitutive_plastic_LpAndItsTangents ( Lp , dLp_dS , dLp_dFi , &
2020-12-23 11:37:18 +05:30
S , Fi , co , ip , el )
2020-07-15 18:05:21 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-07-15 18:05:21 +05:30
ip , & !< integration point
el !< element
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
S , & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 ) :: &
Lp !< plastic velocity gradient
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 , 3 , 3 ) :: &
dLp_dS , &
dLp_dFi !< derivative of Lp with respect to Fi
end subroutine constitutive_plastic_LpAndItsTangents
2020-12-23 11:37:18 +05:30
module subroutine constitutive_plastic_dependentState ( F , co , ip , el )
2020-07-15 18:05:21 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-07-15 18:05:21 +05:30
ip , & !< integration point
el !< element
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
2020-12-23 02:51:11 +05:30
F !< elastic deformation gradient
2020-09-13 14:09:17 +05:30
end subroutine constitutive_plastic_dependentState
2020-07-15 18:05:21 +05:30
2020-07-12 16:57:28 +05:30
2020-12-19 22:13:37 +05:30
2020-12-23 11:37:18 +05:30
module subroutine constitutive_hooke_SandItsTangents ( S , dS_dFe , dS_dFi , Fe , Fi , co , ip , el )
2020-11-04 17:13:52 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-11-04 17:13:52 +05:30
ip , & !< integration point
el !< element
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
Fe , & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 ) :: &
S !< 2nd Piola-Kirchhoff stress tensor
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 , 3 , 3 ) :: &
dS_dFe , & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
end subroutine constitutive_hooke_SandItsTangents
2020-12-21 12:35:38 +05:30
end interface
2020-11-04 17:13:52 +05:30
2020-07-02 02:21:21 +05:30
2020-12-23 12:50:29 +05:30
2020-07-02 04:55:24 +05:30
type ( tDebugOptions ) :: debugConstitutive
2020-09-13 14:09:17 +05:30
2019-06-15 19:10:22 +05:30
public :: &
constitutive_init , &
constitutive_homogenizedC , &
constitutive_LiAndItsTangents , &
2020-07-15 18:05:21 +05:30
constitutive_damage_getRateAndItsTangents , &
constitutive_thermal_getRateAndItsTangents , &
2020-08-15 19:32:10 +05:30
constitutive_results , &
constitutive_allocateState , &
2020-12-20 13:57:37 +05:30
constitutive_forward , &
constitutive_restore , &
2020-08-15 19:32:10 +05:30
plastic_nonlocal_updateCompatibility , &
source_active , &
2020-12-21 15:27:18 +05:30
kinematics_active , &
converged , &
2020-12-20 15:18:13 +05:30
crystallite_init , &
crystallite_stress , &
crystallite_stressTangent , &
crystallite_orientations , &
crystallite_push33ToRef , &
crystallite_restartWrite , &
2020-12-24 14:52:41 +05:30
integrateSourceState , &
2020-12-20 15:18:13 +05:30
crystallite_restartRead , &
2020-12-22 04:03:32 +05:30
constitutive_initializeRestorationPoints , &
2020-12-22 14:53:46 +05:30
constitutive_windForward , &
2020-12-23 22:00:19 +05:30
crystallite_restore , &
PLASTICITY_UNDEFINED_ID , &
PLASTICITY_NONE_ID , &
PLASTICITY_ISOTROPIC_ID , &
PLASTICITY_PHENOPOWERLAW_ID , &
PLASTICITY_KINEHARDENING_ID , &
PLASTICITY_DISLOTWIN_ID , &
PLASTICITY_DISLOTUNGSTEN_ID , &
PLASTICITY_NONLOCAL_ID , &
SOURCE_UNDEFINED_ID , &
SOURCE_THERMAL_DISSIPATION_ID , &
SOURCE_THERMAL_EXTERNALHEAT_ID , &
SOURCE_DAMAGE_ISOBRITTLE_ID , &
SOURCE_DAMAGE_ISODUCTILE_ID , &
SOURCE_DAMAGE_ANISOBRITTLE_ID , &
SOURCE_DAMAGE_ANISODUCTILE_ID , &
KINEMATICS_UNDEFINED_ID , &
KINEMATICS_CLEAVAGE_OPENING_ID , &
KINEMATICS_SLIPPLANE_OPENING_ID , &
KINEMATICS_THERMAL_EXPANSION_ID
2020-12-21 15:27:18 +05:30
2012-03-09 01:55:28 +05:30
contains
2012-10-02 18:23:25 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-13 14:09:17 +05:30
!> @brief Initialze constitutive models for individual physics
2012-10-02 18:23:25 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-06 01:56:23 +05:30
subroutine constitutive_init
2014-10-11 02:25:09 +05:30
2019-06-15 19:10:22 +05:30
integer :: &
2020-12-23 19:33:03 +05:30
ph , & !< counter in phase loop
so !< counter in source loop
2020-07-02 02:21:21 +05:30
class ( tNode ) , pointer :: &
2020-08-15 19:32:10 +05:30
debug_constitutive , &
2020-11-03 02:09:33 +05:30
phases
2020-07-02 02:21:21 +05:30
2020-12-23 19:33:03 +05:30
2020-09-13 14:09:17 +05:30
debug_constitutive = > config_debug % get ( 'constitutive' , defaultVal = emptyList )
debugConstitutive % basic = debug_constitutive % contains ( 'basic' )
debugConstitutive % extensive = debug_constitutive % contains ( 'extensive' )
2020-07-03 20:15:11 +05:30
debugConstitutive % selective = debug_constitutive % contains ( 'selective' )
2020-09-13 14:09:17 +05:30
debugConstitutive % element = config_debug % get_asInt ( 'element' , defaultVal = 1 )
debugConstitutive % ip = config_debug % get_asInt ( 'integrationpoint' , defaultVal = 1 )
debugConstitutive % grain = config_debug % get_asInt ( 'grain' , defaultVal = 1 )
2020-07-02 02:21:21 +05:30
2020-08-15 19:32:10 +05:30
!--------------------------------------------------------------------------------------------------
! initialize constitutive laws
2020-12-23 19:33:03 +05:30
print '(/,a)' , ' <<<+- constitutive init -+>>>' ; flush ( IO_STDOUT )
2020-11-03 02:09:33 +05:30
call mech_init
2020-07-09 04:31:08 +05:30
call damage_init
call thermal_init
2020-02-07 16:53:22 +05:30
2020-12-19 22:13:37 +05:30
phases = > config_material % get ( 'phase' )
2019-06-15 19:10:22 +05:30
constitutive_source_maxSizeDotState = 0
2020-12-23 19:33:03 +05:30
PhaseLoop2 : do ph = 1 , phases % length
2016-01-16 22:57:19 +05:30
!--------------------------------------------------------------------------------------------------
2020-06-26 15:14:17 +05:30
! partition and initialize state
2020-12-23 19:33:03 +05:30
plasticState ( ph ) % partitionedState0 = plasticState ( ph ) % state0
plasticState ( ph ) % state = plasticState ( ph ) % partitionedState0
forall ( so = 1 : phase_Nsources ( ph ) )
sourceState ( ph ) % p ( so ) % partitionedState0 = sourceState ( ph ) % p ( so ) % state0
sourceState ( ph ) % p ( so ) % state = sourceState ( ph ) % p ( so ) % partitionedState0
2019-06-15 19:10:22 +05:30
end forall
2020-07-13 18:18:23 +05:30
2020-03-16 23:58:50 +05:30
constitutive_source_maxSizeDotState = max ( constitutive_source_maxSizeDotState , &
2020-12-23 19:33:03 +05:30
maxval ( sourceState ( ph ) % p % sizeDotState ) )
2019-06-15 19:10:22 +05:30
enddo PhaseLoop2
2020-03-16 23:58:50 +05:30
constitutive_plasticity_maxSizeDotState = maxval ( plasticState % sizeDotState )
2014-06-24 14:54:59 +05:30
2012-10-02 18:23:25 +05:30
end subroutine constitutive_init
2009-03-06 15:32:36 +05:30
2020-08-15 19:32:10 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief checks if a source mechanism is active or not
!--------------------------------------------------------------------------------------------------
2020-12-22 16:50:00 +05:30
function source_active ( source_label , src_length ) result ( active_source )
2020-08-15 19:32:10 +05:30
2020-09-13 14:09:17 +05:30
character ( len = * ) , intent ( in ) :: source_label !< name of source mechanism
2020-08-15 19:32:10 +05:30
integer , intent ( in ) :: src_length !< max. number of sources in system
logical , dimension ( : , : ) , allocatable :: active_source
class ( tNode ) , pointer :: &
phases , &
phase , &
sources , &
2020-09-13 14:09:17 +05:30
src
2020-08-15 19:32:10 +05:30
integer :: p , s
2020-09-13 14:09:17 +05:30
phases = > config_material % get ( 'phase' )
2020-08-15 19:32:10 +05:30
allocate ( active_source ( src_length , phases % length ) , source = . false . )
do p = 1 , phases % length
phase = > phases % get ( p )
sources = > phase % get ( 'source' , defaultVal = emptyList )
do s = 1 , sources % length
src = > sources % get ( s )
if ( src % get_asString ( 'type' ) == source_label ) active_source ( s , p ) = . true .
enddo
enddo
end function source_active
!--------------------------------------------------------------------------------------------------
!> @brief checks if a kinematic mechanism is active or not
!--------------------------------------------------------------------------------------------------
2020-12-22 16:50:00 +05:30
function kinematics_active ( kinematics_label , kinematics_length ) result ( active_kinematics )
2020-08-15 19:32:10 +05:30
character ( len = * ) , intent ( in ) :: kinematics_label !< name of kinematic mechanism
integer , intent ( in ) :: kinematics_length !< max. number of kinematics in system
logical , dimension ( : , : ) , allocatable :: active_kinematics
class ( tNode ) , pointer :: &
phases , &
phase , &
kinematics , &
2020-09-13 14:09:17 +05:30
kinematics_type
2020-08-15 19:32:10 +05:30
integer :: p , k
2020-09-13 14:09:17 +05:30
phases = > config_material % get ( 'phase' )
2020-08-15 19:32:10 +05:30
allocate ( active_kinematics ( kinematics_length , phases % length ) , source = . false . )
do p = 1 , phases % length
phase = > phases % get ( p )
kinematics = > phase % get ( 'kinematics' , defaultVal = emptyList )
do k = 1 , kinematics % length
kinematics_type = > kinematics % get ( k )
if ( kinematics_type % get_asString ( 'type' ) == kinematics_label ) active_kinematics ( k , p ) = . true .
enddo
enddo
end function kinematics_active
2020-09-13 14:09:17 +05:30
2020-08-15 19:32:10 +05:30
2014-11-01 00:33:08 +05:30
!--------------------------------------------------------------------------------------------------
2015-10-14 00:22:01 +05:30
!> @brief contains the constitutive equation for calculating the velocity gradient
2018-12-30 17:05:26 +05:30
! ToDo: MD: S is Mi?
2014-11-01 00:33:08 +05:30
!--------------------------------------------------------------------------------------------------
2019-02-18 11:54:56 +05:30
subroutine constitutive_LiAndItsTangents ( Li , dLi_dS , dLi_dFi , &
2020-12-23 11:37:18 +05:30
S , Fi , co , ip , el )
2015-10-14 00:22:01 +05:30
2019-06-15 19:10:22 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2019-06-15 19:10:22 +05:30
ip , & !< integration point
el !< element
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
S !< 2nd Piola-Kirchhoff stress
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
Fi !< intermediate deformation gradient
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 ) :: &
Li !< intermediate velocity gradient
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 , 3 , 3 ) :: &
dLi_dS , & !< derivative of Li with respect to S
dLi_dFi
2020-02-07 16:53:22 +05:30
2019-06-15 19:10:22 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: &
my_Li , & !< intermediate velocity gradient
FiInv , &
temp_33
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: &
my_dLi_dS
real ( pReal ) :: &
detFi
integer :: &
k , i , j , &
instance , of
2020-07-03 20:15:11 +05:30
2019-06-15 19:10:22 +05:30
Li = 0.0_pReal
dLi_dS = 0.0_pReal
dLi_dFi = 0.0_pReal
2020-02-07 16:53:22 +05:30
2020-12-23 11:37:18 +05:30
plasticityType : select case ( phase_plasticity ( material_phaseAt ( co , el ) ) )
2019-06-15 19:10:22 +05:30
case ( PLASTICITY_isotropic_ID ) plasticityType
2020-12-23 11:37:18 +05:30
of = material_phasememberAt ( co , ip , el )
instance = phase_plasticityInstance ( material_phaseAt ( co , el ) )
2020-07-02 00:52:05 +05:30
call plastic_isotropic_LiAndItsTangent ( my_Li , my_dLi_dS , S , instance , of )
2019-06-15 19:10:22 +05:30
case default plasticityType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select plasticityType
2020-02-07 16:53:22 +05:30
2019-06-15 19:10:22 +05:30
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
2020-02-07 16:53:22 +05:30
2020-12-23 11:37:18 +05:30
KinematicsLoop : do k = 1 , phase_Nkinematics ( material_phaseAt ( co , el ) )
kinematicsType : select case ( phase_kinematics ( k , material_phaseAt ( co , el ) ) )
2019-06-15 19:10:22 +05:30
case ( KINEMATICS_cleavage_opening_ID ) kinematicsType
2020-12-23 11:37:18 +05:30
call kinematics_cleavage_opening_LiAndItsTangent ( my_Li , my_dLi_dS , S , co , ip , el )
2019-06-15 19:10:22 +05:30
case ( KINEMATICS_slipplane_opening_ID ) kinematicsType
2020-12-23 11:37:18 +05:30
call kinematics_slipplane_opening_LiAndItsTangent ( my_Li , my_dLi_dS , S , co , ip , el )
2019-06-15 19:10:22 +05:30
case ( KINEMATICS_thermal_expansion_ID ) kinematicsType
2020-12-23 11:37:18 +05:30
call kinematics_thermal_expansion_LiAndItsTangent ( my_Li , my_dLi_dS , co , ip , el )
2019-06-15 19:10:22 +05:30
case default kinematicsType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select kinematicsType
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
enddo KinematicsLoop
2020-02-07 16:53:22 +05:30
2019-06-15 19:10:22 +05:30
FiInv = math_inv33 ( Fi )
detFi = math_det33 ( Fi )
Li = matmul ( matmul ( Fi , Li ) , FiInv ) * detFi !< push forward to intermediate configuration
temp_33 = matmul ( FiInv , Li )
2020-02-07 16:53:22 +05:30
2019-06-15 19:10:22 +05:30
do i = 1 , 3 ; do j = 1 , 3
dLi_dS ( 1 : 3 , 1 : 3 , i , j ) = matmul ( matmul ( Fi , dLi_dS ( 1 : 3 , 1 : 3 , i , j ) ) , FiInv ) * detFi
dLi_dFi ( 1 : 3 , 1 : 3 , i , j ) = dLi_dFi ( 1 : 3 , 1 : 3 , i , j ) + Li * FiInv ( j , i )
dLi_dFi ( 1 : 3 , i , 1 : 3 , j ) = dLi_dFi ( 1 : 3 , i , 1 : 3 , j ) + math_I3 * temp_33 ( j , i ) + Li * FiInv ( j , i )
2019-07-01 10:39:51 +05:30
enddo ; enddo
2015-10-14 00:22:01 +05:30
2018-08-28 18:37:39 +05:30
end subroutine constitutive_LiAndItsTangents
2014-11-01 00:33:08 +05:30
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
2020-12-19 22:52:40 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
2020-12-23 19:33:03 +05:30
function constitutive_damage_collectDotState ( S , co , ip , el , ph , of ) result ( broken )
2020-12-19 22:52:40 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-12-19 22:52:40 +05:30
ip , & !< integration point
el , & !< element
2020-12-23 19:33:03 +05:30
ph , &
2020-12-19 22:52:40 +05:30
of
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
S !< 2nd Piola Kirchhoff stress (vector notation)
integer :: &
2020-12-23 19:33:03 +05:30
so !< counter in source loop
2020-12-19 22:52:40 +05:30
logical :: broken
broken = . false .
2020-12-23 19:33:03 +05:30
SourceLoop : do so = 1 , phase_Nsources ( ph )
2020-07-10 20:40:23 +05:30
2020-12-23 19:33:03 +05:30
sourceType : select case ( phase_source ( so , ph ) )
2020-07-10 20:40:23 +05:30
case ( SOURCE_damage_anisoBrittle_ID ) sourceType
2020-12-23 11:37:18 +05:30
call source_damage_anisoBrittle_dotState ( S , co , ip , el ) ! correct stress?
2020-07-10 20:40:23 +05:30
case ( SOURCE_damage_isoDuctile_ID ) sourceType
2020-12-23 11:37:18 +05:30
call source_damage_isoDuctile_dotState ( co , ip , el )
2020-07-10 20:40:23 +05:30
case ( SOURCE_damage_anisoDuctile_ID ) sourceType
2020-12-23 11:37:18 +05:30
call source_damage_anisoDuctile_dotState ( co , ip , el )
2020-07-10 20:40:23 +05:30
2020-12-22 16:50:00 +05:30
end select sourceType
2020-12-23 19:33:03 +05:30
broken = broken . or . any ( IEEE_is_NaN ( sourceState ( ph ) % p ( so ) % dotState ( : , of ) ) )
2020-12-22 16:50:00 +05:30
enddo SourceLoop
end function constitutive_damage_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
2020-12-23 02:51:11 +05:30
function constitutive_thermal_collectDotState ( ph , me ) result ( broken )
2020-12-22 16:50:00 +05:30
2020-12-23 02:51:11 +05:30
integer , intent ( in ) :: ph , me
2020-12-22 16:50:00 +05:30
logical :: broken
2020-12-23 02:51:11 +05:30
integer :: i
2020-12-22 16:50:00 +05:30
2020-12-23 02:51:11 +05:30
broken = . false .
2020-12-22 16:50:00 +05:30
2020-12-23 02:51:11 +05:30
SourceLoop : do i = 1 , phase_Nsources ( ph )
2020-07-10 20:40:23 +05:30
2020-12-23 02:51:11 +05:30
if ( phase_source ( i , ph ) == SOURCE_thermal_externalheat_ID ) &
call source_thermal_externalheat_dotState ( ph , me )
2020-07-10 20:40:23 +05:30
2020-12-23 02:51:11 +05:30
broken = broken . or . any ( IEEE_is_NaN ( sourceState ( ph ) % p ( i ) % dotState ( : , me ) ) )
2020-07-10 20:40:23 +05:30
enddo SourceLoop
2020-04-01 11:11:55 +05:30
2020-12-22 16:50:00 +05:30
end function constitutive_thermal_collectDotState
2009-03-06 15:32:36 +05:30
2019-06-15 19:10:22 +05:30
2020-12-19 22:52:40 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
2020-12-23 18:33:15 +05:30
function constitutive_damage_deltaState ( Fe , co , ip , el , ph , of ) result ( broken )
2020-12-19 22:52:40 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-12-19 22:52:40 +05:30
ip , & !< integration point
el , & !< element
2020-12-23 18:33:15 +05:30
ph , &
2020-12-19 22:52:40 +05:30
of
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
2020-12-20 14:14:25 +05:30
Fe !< elastic deformation gradient
2020-12-19 22:52:40 +05:30
integer :: &
2020-12-23 18:33:15 +05:30
so , &
2020-12-19 22:52:40 +05:30
myOffset , &
mySize
logical :: &
broken
broken = . false .
2012-05-16 20:13:26 +05:30
2020-12-23 18:33:15 +05:30
sourceLoop : do so = 1 , phase_Nsources ( ph )
2018-08-25 19:29:34 +05:30
2020-12-23 18:33:15 +05:30
sourceType : select case ( phase_source ( so , ph ) )
2018-08-25 19:29:34 +05:30
2019-06-15 19:10:22 +05:30
case ( SOURCE_damage_isoBrittle_ID ) sourceType
2020-12-23 11:37:18 +05:30
call source_damage_isoBrittle_deltaState ( constitutive_homogenizedC ( co , ip , el ) , Fe , &
co , ip , el )
2020-12-23 18:33:15 +05:30
broken = any ( IEEE_is_NaN ( sourceState ( ph ) % p ( so ) % deltaState ( : , of ) ) )
2020-04-01 18:12:38 +05:30
if ( . not . broken ) then
2020-12-23 18:33:15 +05:30
myOffset = sourceState ( ph ) % p ( so ) % offsetDeltaState
mySize = sourceState ( ph ) % p ( so ) % sizeDeltaState
sourceState ( ph ) % p ( so ) % state ( myOffset + 1 : myOffset + mySize , of ) = &
sourceState ( ph ) % p ( so ) % state ( myOffset + 1 : myOffset + mySize , of ) + sourceState ( ph ) % p ( so ) % deltaState ( 1 : mySize , of )
2020-04-01 18:12:38 +05:30
endif
2018-08-25 19:29:34 +05:30
2019-06-15 19:10:22 +05:30
end select sourceType
2018-08-25 19:29:34 +05:30
2019-06-15 19:10:22 +05:30
enddo SourceLoop
2015-05-28 22:32:23 +05:30
2020-12-22 23:54:00 +05:30
end function constitutive_damage_deltaState
2014-09-10 23:56:12 +05:30
2014-09-22 23:45:19 +05:30
2020-08-15 19:32:10 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Allocate the components of the state structure for a given phase
!--------------------------------------------------------------------------------------------------
subroutine constitutive_allocateState ( state , &
2020-10-28 02:03:30 +05:30
Nconstituents , sizeState , sizeDotState , sizeDeltaState )
2020-08-15 19:32:10 +05:30
class ( tState ) , intent ( out ) :: &
state
integer , intent ( in ) :: &
2020-10-28 02:03:30 +05:30
Nconstituents , &
2020-08-15 19:32:10 +05:30
sizeState , &
sizeDotState , &
sizeDeltaState
state % sizeState = sizeState
state % sizeDotState = sizeDotState
state % sizeDeltaState = sizeDeltaState
state % offsetDeltaState = sizeState - sizeDeltaState ! deltaState occupies latter part of state by definition
allocate ( state % atol ( sizeState ) , source = 0.0_pReal )
2020-10-28 02:03:30 +05:30
allocate ( state % state0 ( sizeState , Nconstituents ) , source = 0.0_pReal )
allocate ( state % partitionedState0 ( sizeState , Nconstituents ) , source = 0.0_pReal )
allocate ( state % subState0 ( sizeState , Nconstituents ) , source = 0.0_pReal )
allocate ( state % state ( sizeState , Nconstituents ) , source = 0.0_pReal )
2020-08-15 19:32:10 +05:30
2020-10-28 02:03:30 +05:30
allocate ( state % dotState ( sizeDotState , Nconstituents ) , source = 0.0_pReal )
2020-08-15 19:32:10 +05:30
2020-10-28 02:03:30 +05:30
allocate ( state % deltaState ( sizeDeltaState , Nconstituents ) , source = 0.0_pReal )
2020-08-15 19:32:10 +05:30
end subroutine constitutive_allocateState
2020-12-20 13:57:37 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
subroutine constitutive_restore ( ip , el )
2020-12-20 13:57:37 +05:30
integer , intent ( in ) :: &
2020-12-23 11:44:07 +05:30
ip , & !< integration point number
el !< element number
2020-12-20 13:57:37 +05:30
integer :: &
2020-12-23 11:44:07 +05:30
co , & !< constituent number
2020-12-20 13:57:37 +05:30
s
2020-12-23 11:44:07 +05:30
do co = 1 , homogenization_Nconstituents ( material_homogenizationAt ( el ) )
do s = 1 , phase_Nsources ( material_phaseAt ( co , el ) )
sourceState ( material_phaseAt ( co , el ) ) % p ( s ) % state ( : , material_phasememberAt ( co , ip , el ) ) = &
sourceState ( material_phaseAt ( co , el ) ) % p ( s ) % partitionedState0 ( : , material_phasememberAt ( co , ip , el ) )
2020-12-20 13:57:37 +05:30
enddo
enddo
end subroutine constitutive_restore
!--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
subroutine constitutive_forward
integer :: i , j
2020-12-23 11:28:54 +05:30
crystallite_F0 = crystallite_partitionedF
crystallite_Lp0 = crystallite_Lp
crystallite_S0 = crystallite_S
call constitutive_mech_forward ( )
2020-12-20 13:57:37 +05:30
do i = 1 , size ( sourceState )
do j = 1 , phase_Nsources ( i )
sourceState ( i ) % p ( j ) % state0 = sourceState ( i ) % p ( j ) % state
enddo ; enddo
end subroutine constitutive_forward
2018-12-05 04:25:39 +05:30
!--------------------------------------------------------------------------------------------------
2018-12-12 11:10:57 +05:30
!> @brief writes constitutive results to HDF5 output file
2018-12-05 04:25:39 +05:30
!--------------------------------------------------------------------------------------------------
2019-04-06 10:01:02 +05:30
subroutine constitutive_results
2019-12-19 00:35:51 +05:30
2020-12-22 13:15:01 +05:30
integer :: ph
character ( len = : ) , allocatable :: group
2020-12-22 16:50:00 +05:30
call results_closeGroup ( results_addGroup ( '/current/phase/' ) )
2020-12-22 13:15:01 +05:30
do ph = 1 , size ( material_name_phase )
2020-12-22 16:50:00 +05:30
group = '/current/phase/' / / trim ( material_name_phase ( ph ) ) / / '/'
2020-12-22 13:15:01 +05:30
call results_closeGroup ( results_addGroup ( group ) )
call mech_results ( group , ph )
2020-12-22 16:50:00 +05:30
call damage_results ( group , ph )
2020-12-22 13:15:01 +05:30
enddo
2018-12-05 04:25:39 +05:30
end subroutine constitutive_results
2020-08-15 19:32:10 +05:30
2020-12-20 15:18:13 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief allocates and initialize per grain variables
!--------------------------------------------------------------------------------------------------
subroutine crystallite_init
integer :: &
2020-12-21 00:50:39 +05:30
Nconstituents , &
2020-12-23 17:52:11 +05:30
ph , &
me , &
2020-12-23 13:27:58 +05:30
co , & !< counter in integration point component loop
ip , & !< counter in integration point loop
el , & !< counter in element loop
2020-12-20 15:18:13 +05:30
cMax , & !< maximum number of integration point components
iMax , & !< maximum number of integration points
eMax !< maximum number of elements
class ( tNode ) , pointer :: &
num_crystallite , &
debug_crystallite , & ! pointer to debug options for crystallite
phases , &
phase , &
mech
2020-12-23 11:44:07 +05:30
2020-12-20 15:18:13 +05:30
print '(/,a)' , ' <<<+- crystallite init -+>>>'
debug_crystallite = > config_debug % get ( 'crystallite' , defaultVal = emptyList )
debugCrystallite % extensive = debug_crystallite % contains ( 'extensive' )
cMax = homogenization_maxNconstituents
iMax = discretization_nIPs
eMax = discretization_Nelems
allocate ( crystallite_partitionedF ( 3 , 3 , cMax , iMax , eMax ) , source = 0.0_pReal )
allocate ( crystallite_S0 , &
2020-12-23 02:51:11 +05:30
crystallite_F0 , crystallite_Lp0 , &
2020-12-20 15:18:13 +05:30
crystallite_partitionedS0 , &
2020-12-23 02:51:11 +05:30
crystallite_partitionedF0 , &
2020-12-21 14:29:13 +05:30
crystallite_partitionedLp0 , &
2020-12-20 15:18:13 +05:30
crystallite_S , crystallite_P , &
2020-12-23 02:51:11 +05:30
crystallite_Fe , crystallite_Lp , &
2020-12-24 15:50:34 +05:30
crystallite_subF , &
2020-12-20 15:18:13 +05:30
crystallite_subFp0 , crystallite_subFi0 , &
source = crystallite_partitionedF )
2020-12-23 19:27:53 +05:30
allocate ( crystallite_subdt ( cMax , iMax , eMax ) , source = 0.0_pReal )
2020-12-20 15:18:13 +05:30
allocate ( crystallite_orientation ( cMax , iMax , eMax ) )
allocate ( crystallite_converged ( cMax , iMax , eMax ) , source = . true . )
num_crystallite = > config_numerics % get ( 'crystallite' , defaultVal = emptyDict )
num % subStepMinCryst = num_crystallite % get_asFloat ( 'subStepMin' , defaultVal = 1.0e-3_pReal )
num % subStepSizeCryst = num_crystallite % get_asFloat ( 'subStepSize' , defaultVal = 0.25_pReal )
num % stepIncreaseCryst = num_crystallite % get_asFloat ( 'stepIncrease' , defaultVal = 1.5_pReal )
num % subStepSizeLp = num_crystallite % get_asFloat ( 'subStepSizeLp' , defaultVal = 0.5_pReal )
num % subStepSizeLi = num_crystallite % get_asFloat ( 'subStepSizeLi' , defaultVal = 0.5_pReal )
num % rtol_crystalliteState = num_crystallite % get_asFloat ( 'rtol_State' , defaultVal = 1.0e-6_pReal )
num % rtol_crystalliteStress = num_crystallite % get_asFloat ( 'rtol_Stress' , defaultVal = 1.0e-6_pReal )
num % atol_crystalliteStress = num_crystallite % get_asFloat ( 'atol_Stress' , defaultVal = 1.0e-8_pReal )
num % iJacoLpresiduum = num_crystallite % get_asInt ( 'iJacoLpresiduum' , defaultVal = 1 )
num % nState = num_crystallite % get_asInt ( 'nState' , defaultVal = 20 )
num % nStress = num_crystallite % get_asInt ( 'nStress' , defaultVal = 40 )
if ( num % subStepMinCryst < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepMinCryst' )
if ( num % subStepSizeCryst < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepSizeCryst' )
if ( num % stepIncreaseCryst < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'stepIncreaseCryst' )
if ( num % subStepSizeLp < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepSizeLp' )
if ( num % subStepSizeLi < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'subStepSizeLi' )
if ( num % rtol_crystalliteState < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'rtol_crystalliteState' )
if ( num % rtol_crystalliteStress < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'rtol_crystalliteStress' )
if ( num % atol_crystalliteStress < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'atol_crystalliteStress' )
if ( num % iJacoLpresiduum < 1 ) call IO_error ( 301 , ext_msg = 'iJacoLpresiduum' )
if ( num % nState < 1 ) call IO_error ( 301 , ext_msg = 'nState' )
if ( num % nStress < 1 ) call IO_error ( 301 , ext_msg = 'nStress' )
phases = > config_material % get ( 'phase' )
2020-12-21 00:50:39 +05:30
allocate ( constitutive_mech_Fi ( phases % length ) )
allocate ( constitutive_mech_Fi0 ( phases % length ) )
2020-12-24 16:24:09 +05:30
allocate ( constitutive_mech_partitionedFi0 ( phases % length ) )
2020-12-23 02:51:11 +05:30
allocate ( constitutive_mech_Fp ( phases % length ) )
allocate ( constitutive_mech_Fp0 ( phases % length ) )
2020-12-24 16:24:09 +05:30
allocate ( constitutive_mech_partitionedFp0 ( phases % length ) )
2020-12-21 14:29:13 +05:30
allocate ( constitutive_mech_Li ( phases % length ) )
allocate ( constitutive_mech_Li0 ( phases % length ) )
2020-12-24 16:24:09 +05:30
allocate ( constitutive_mech_partitionedLi0 ( phases % length ) )
2020-12-23 17:52:11 +05:30
do ph = 1 , phases % length
Nconstituents = count ( material_phaseAt == ph ) * discretization_nIPs
allocate ( constitutive_mech_Fi ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( constitutive_mech_Fi0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
2020-12-24 16:24:09 +05:30
allocate ( constitutive_mech_partitionedFi0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
2020-12-23 17:52:11 +05:30
allocate ( constitutive_mech_Fp ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( constitutive_mech_Fp0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
2020-12-24 16:24:09 +05:30
allocate ( constitutive_mech_partitionedFp0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
2020-12-23 17:52:11 +05:30
allocate ( constitutive_mech_Li ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( constitutive_mech_Li0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
2020-12-24 16:24:09 +05:30
allocate ( constitutive_mech_partitionedLi0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
2020-12-20 15:18:13 +05:30
enddo
2020-12-20 22:52:04 +05:30
print '(a42,1x,i10)' , ' # of elements: ' , eMax
print '(a42,1x,i10)' , ' # of integration points/element: ' , iMax
print '(a42,1x,i10)' , 'max # of constituents/integration point: ' , cMax
flush ( IO_STDOUT )
2020-12-20 15:18:13 +05:30
2020-12-23 17:52:11 +05:30
!$OMP PARALLEL DO PRIVATE(ph,me)
2020-12-23 13:27:58 +05:30
do el = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do ip = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 ) ; do co = 1 , homogenization_Nconstituents ( material_homogenizationAt ( el ) )
2020-12-21 00:50:39 +05:30
2020-12-23 17:52:11 +05:30
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
constitutive_mech_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = material_orientation0 ( co , ip , el ) % asMatrix ( ) ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
constitutive_mech_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) &
/ math_det33 ( constitutive_mech_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) ) ** ( 1.0_pReal / 3.0_pReal )
constitutive_mech_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = math_I3
2020-12-21 00:50:39 +05:30
2020-12-23 13:27:58 +05:30
crystallite_F0 ( 1 : 3 , 1 : 3 , co , ip , el ) = math_I3
2020-12-21 00:50:39 +05:30
2020-12-23 17:52:11 +05:30
crystallite_Fe ( 1 : 3 , 1 : 3 , co , ip , el ) = math_inv33 ( matmul ( constitutive_mech_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , &
constitutive_mech_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) ) ) ! assuming that euler angles are given in internal strain free configuration
constitutive_mech_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
constitutive_mech_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-21 00:50:39 +05:30
2020-12-24 16:24:09 +05:30
constitutive_mech_partitionedFi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
constitutive_mech_partitionedFp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-21 00:50:39 +05:30
2020-12-20 15:18:13 +05:30
enddo ; enddo
enddo
!$OMP END PARALLEL DO
crystallite_partitionedF0 = crystallite_F0
crystallite_partitionedF = crystallite_F0
call crystallite_orientations ( )
2020-12-23 17:52:11 +05:30
!$OMP PARALLEL DO PRIVATE(ph,me)
2020-12-23 13:27:58 +05:30
do el = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do ip = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
do co = 1 , homogenization_Nconstituents ( material_homogenizationAt ( el ) )
2020-12-23 17:52:11 +05:30
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
2020-12-24 13:23:02 +05:30
call constitutive_plastic_dependentState ( crystallite_partitionedF0 ( 1 : 3 , 1 : 3 , co , ip , el ) , co , ip , el ) ! update dependent state variables to be consistent with basic states
2020-12-20 15:18:13 +05:30
enddo
enddo
enddo
!$OMP END PARALLEL DO
end subroutine crystallite_init
2020-12-23 14:35:02 +05:30
2020-12-20 15:18:13 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Backup data for homog cutback.
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
subroutine constitutive_initializeRestorationPoints ( ip , el )
2020-12-20 15:18:13 +05:30
integer , intent ( in ) :: &
2020-12-23 11:44:07 +05:30
ip , & !< integration point number
el !< element number
2020-12-20 15:18:13 +05:30
integer :: &
2020-12-23 11:44:07 +05:30
co , & !< constituent number
2020-12-22 14:33:19 +05:30
s , ph , me
2020-12-20 15:18:13 +05:30
2020-12-23 11:44:07 +05:30
do co = 1 , homogenization_Nconstituents ( material_homogenizationAt ( el ) )
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
crystallite_partitionedLp0 ( 1 : 3 , 1 : 3 , co , ip , el ) = crystallite_Lp0 ( 1 : 3 , 1 : 3 , co , ip , el )
crystallite_partitionedF0 ( 1 : 3 , 1 : 3 , co , ip , el ) = crystallite_F0 ( 1 : 3 , 1 : 3 , co , ip , el )
crystallite_partitionedS0 ( 1 : 3 , 1 : 3 , co , ip , el ) = crystallite_S0 ( 1 : 3 , 1 : 3 , co , ip , el )
2020-12-20 15:18:13 +05:30
2020-12-22 14:33:19 +05:30
call mech_initializeRestorationPoints ( ph , me )
2020-12-23 03:57:56 +05:30
2020-12-23 11:44:07 +05:30
do s = 1 , phase_Nsources ( material_phaseAt ( co , el ) )
sourceState ( material_phaseAt ( co , el ) ) % p ( s ) % partitionedState0 ( : , material_phasememberAt ( co , ip , el ) ) = &
sourceState ( material_phaseAt ( co , el ) ) % p ( s ) % state0 ( : , material_phasememberAt ( co , ip , el ) )
2020-12-20 15:18:13 +05:30
enddo
enddo
2020-12-22 04:03:32 +05:30
end subroutine constitutive_initializeRestorationPoints
2020-12-20 15:18:13 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward.
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
subroutine constitutive_windForward ( ip , el )
2020-12-20 15:18:13 +05:30
integer , intent ( in ) :: &
2020-12-23 11:44:07 +05:30
ip , & !< integration point number
el !< element number
2020-12-20 15:18:13 +05:30
integer :: &
2020-12-23 11:44:07 +05:30
co , & !< constituent number
2020-12-23 03:57:56 +05:30
s , ph , me
2020-12-23 11:44:07 +05:30
do co = 1 , homogenization_Nconstituents ( material_homogenizationAt ( el ) )
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
crystallite_partitionedF0 ( 1 : 3 , 1 : 3 , co , ip , el ) = crystallite_partitionedF ( 1 : 3 , 1 : 3 , co , ip , el )
crystallite_partitionedLp0 ( 1 : 3 , 1 : 3 , co , ip , el ) = crystallite_Lp ( 1 : 3 , 1 : 3 , co , ip , el )
crystallite_partitionedS0 ( 1 : 3 , 1 : 3 , co , ip , el ) = crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el )
2020-12-20 15:18:13 +05:30
2020-12-23 03:57:56 +05:30
call constitutive_mech_windForward ( ph , me )
2020-12-23 11:44:07 +05:30
do s = 1 , phase_Nsources ( material_phaseAt ( co , el ) )
2020-12-23 03:57:56 +05:30
sourceState ( ph ) % p ( s ) % partitionedState0 ( : , me ) = sourceState ( ph ) % p ( s ) % state ( : , me )
2020-12-20 15:18:13 +05:30
enddo
enddo
2020-12-22 14:53:46 +05:30
end subroutine constitutive_windForward
2020-12-20 15:18:13 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
subroutine crystallite_restore ( ip , el , includeL )
2020-12-20 15:18:13 +05:30
integer , intent ( in ) :: &
2020-12-23 11:44:07 +05:30
ip , & !< integration point number
el !< element number
2020-12-20 15:18:13 +05:30
logical , intent ( in ) :: &
includeL !< protect agains fake cutback
integer :: &
2020-12-23 11:44:07 +05:30
co , p , m !< constituent number
2020-12-20 15:18:13 +05:30
2020-12-23 11:44:07 +05:30
do co = 1 , homogenization_Nconstituents ( material_homogenizationAt ( el ) )
p = material_phaseAt ( co , el )
m = material_phaseMemberAt ( co , ip , el )
2020-12-20 15:18:13 +05:30
if ( includeL ) then
2020-12-23 11:44:07 +05:30
crystallite_Lp ( 1 : 3 , 1 : 3 , co , ip , el ) = crystallite_partitionedLp0 ( 1 : 3 , 1 : 3 , co , ip , el )
2020-12-24 16:24:09 +05:30
constitutive_mech_Li ( p ) % data ( 1 : 3 , 1 : 3 , m ) = constitutive_mech_partitionedLi0 ( p ) % data ( 1 : 3 , 1 : 3 , m )
2020-12-20 15:18:13 +05:30
endif ! maybe protecting everything from overwriting makes more sense
2020-12-21 14:29:13 +05:30
2020-12-24 16:24:09 +05:30
constitutive_mech_Fp ( p ) % data ( 1 : 3 , 1 : 3 , m ) = constitutive_mech_partitionedFp0 ( p ) % data ( 1 : 3 , 1 : 3 , m )
constitutive_mech_Fi ( p ) % data ( 1 : 3 , 1 : 3 , m ) = constitutive_mech_partitionedFi0 ( p ) % data ( 1 : 3 , 1 : 3 , m )
2020-12-23 11:44:07 +05:30
crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) = crystallite_partitionedS0 ( 1 : 3 , 1 : 3 , co , ip , el )
2020-12-20 15:18:13 +05:30
2020-12-23 11:44:07 +05:30
plasticState ( material_phaseAt ( co , el ) ) % state ( : , material_phasememberAt ( co , ip , el ) ) = &
plasticState ( material_phaseAt ( co , el ) ) % partitionedState0 ( : , material_phasememberAt ( co , ip , el ) )
2020-12-20 15:18:13 +05:30
enddo
end subroutine crystallite_restore
!--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF).
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
function crystallite_stressTangent ( co , ip , el ) result ( dPdF )
2020-12-20 15:18:13 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dPdF
integer , intent ( in ) :: &
2020-12-23 11:44:07 +05:30
co , & !< counter in constituent loop
ip , & !< counter in integration point loop
el !< counter in element loop
2020-12-20 15:18:13 +05:30
integer :: &
o , &
2020-12-21 00:50:39 +05:30
p , pp , m
2020-12-20 15:18:13 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: devNull , &
invSubFp0 , invSubFi0 , invFp , invFi , &
temp_33_1 , temp_33_2 , temp_33_3 , temp_33_4
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dSdFe , &
dSdF , &
dSdFi , &
dLidS , & ! tangent in lattice configuration
dLidFi , &
dLpdS , &
dLpdFi , &
dFidS , &
dFpinvdF , &
rhs_3333 , &
lhs_3333 , &
temp_3333
real ( pReal ) , dimension ( 9 , 9 ) :: temp_99
logical :: error
2020-12-23 11:44:07 +05:30
pp = material_phaseAt ( co , el )
m = material_phaseMemberAt ( co , ip , el )
2020-12-20 15:18:13 +05:30
2020-12-23 12:42:56 +05:30
call constitutive_hooke_SandItsTangents ( devNull , dSdFe , dSdFi , &
crystallite_Fe ( 1 : 3 , 1 : 3 , co , ip , el ) , &
constitutive_mech_Fi ( pp ) % data ( 1 : 3 , 1 : 3 , m ) , co , ip , el )
call constitutive_LiAndItsTangents ( devNull , dLidS , dLidFi , &
crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) , &
constitutive_mech_Fi ( pp ) % data ( 1 : 3 , 1 : 3 , m ) , &
co , ip , el )
invFp = math_inv33 ( constitutive_mech_Fp ( pp ) % data ( 1 : 3 , 1 : 3 , m ) )
invFi = math_inv33 ( constitutive_mech_Fi ( pp ) % data ( 1 : 3 , 1 : 3 , m ) )
invSubFp0 = math_inv33 ( crystallite_subFp0 ( 1 : 3 , 1 : 3 , co , ip , el ) )
invSubFi0 = math_inv33 ( crystallite_subFi0 ( 1 : 3 , 1 : 3 , co , ip , el ) )
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 ( co , ip , el ) * 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 ) &
+ invFi * invFi ( p , o )
rhs_3333 ( 1 : 3 , 1 : 3 , o , p ) = rhs_3333 ( 1 : 3 , 1 : 3 , o , p ) &
- crystallite_subdt ( co , ip , el ) * matmul ( invSubFi0 , dLidS ( 1 : 3 , 1 : 3 , o , p ) )
enddo ; enddo
call math_invert ( temp_99 , error , math_3333to99 ( lhs_3333 ) )
if ( error ) then
call IO_warning ( warning_ID = 600 , el = el , ip = ip , g = co , &
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_plastic_LpAndItsTangents ( devNull , dLpdS , dLpdFi , &
crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) , &
constitutive_mech_Fi ( pp ) % data ( 1 : 3 , 1 : 3 , m ) , co , ip , el )
dLpdS = math_mul3333xx3333 ( dLpdFi , dFidS ) + dLpdS
2020-12-20 15:18:13 +05:30
!--------------------------------------------------------------------------------------------------
! calculate dSdF
2020-12-23 12:42:56 +05:30
temp_33_1 = transpose ( matmul ( invFp , invFi ) )
temp_33_2 = matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , co , ip , el ) , invSubFp0 )
temp_33_3 = matmul ( matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , co , ip , el ) , invFp ) , invSubFi0 )
do o = 1 , 3 ; do p = 1 , 3
rhs_3333 ( p , o , 1 : 3 , 1 : 3 ) = matmul ( dSdFe ( p , o , 1 : 3 , 1 : 3 ) , temp_33_1 )
temp_3333 ( 1 : 3 , 1 : 3 , p , o ) = matmul ( matmul ( temp_33_2 , dLpdS ( 1 : 3 , 1 : 3 , p , o ) ) , invFi ) &
+ matmul ( temp_33_3 , dLidS ( 1 : 3 , 1 : 3 , p , o ) )
enddo ; enddo
lhs_3333 = crystallite_subdt ( co , ip , el ) * math_mul3333xx3333 ( dSdFe , temp_3333 ) &
+ math_mul3333xx3333 ( dSdFi , dFidS )
call math_invert ( temp_99 , error , math_eye ( 9 ) + math_3333to99 ( lhs_3333 ) )
if ( error ) then
call IO_warning ( warning_ID = 600 , el = el , ip = ip , g = co , &
ext_msg = 'inversion error in analytic tangent calculation' )
dSdF = rhs_3333
else
dSdF = math_mul3333xx3333 ( math_99to3333 ( temp_99 ) , rhs_3333 )
endif
2020-12-20 15:18:13 +05:30
!--------------------------------------------------------------------------------------------------
! calculate dFpinvdF
2020-12-23 12:42:56 +05:30
temp_3333 = math_mul3333xx3333 ( dLpdS , dSdF )
do o = 1 , 3 ; do p = 1 , 3
dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) = - crystallite_subdt ( co , ip , el ) &
* matmul ( invSubFp0 , matmul ( temp_3333 ( 1 : 3 , 1 : 3 , p , o ) , invFi ) )
enddo ; enddo
2020-12-20 15:18:13 +05:30
!--------------------------------------------------------------------------------------------------
! assemble dPdF
2020-12-23 12:42:56 +05:30
temp_33_1 = matmul ( crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) , transpose ( invFp ) )
temp_33_2 = matmul ( invFp , temp_33_1 )
temp_33_3 = matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , co , ip , el ) , invFp )
temp_33_4 = matmul ( temp_33_3 , crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) )
dPdF = 0.0_pReal
do p = 1 , 3
dPdF ( p , 1 : 3 , p , 1 : 3 ) = transpose ( temp_33_2 )
enddo
do o = 1 , 3 ; do p = 1 , 3
dPdF ( 1 : 3 , 1 : 3 , p , o ) = dPdF ( 1 : 3 , 1 : 3 , p , o ) &
+ matmul ( matmul ( crystallite_subF ( 1 : 3 , 1 : 3 , co , ip , el ) , &
dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) , temp_33_1 ) &
+ matmul ( matmul ( temp_33_3 , dSdF ( 1 : 3 , 1 : 3 , p , o ) ) , &
transpose ( invFp ) ) &
+ matmul ( temp_33_4 , transpose ( dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) )
enddo ; enddo
2020-12-20 15:18:13 +05:30
end function crystallite_stressTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates orientations
!--------------------------------------------------------------------------------------------------
subroutine crystallite_orientations
integer &
2020-12-23 11:44:07 +05:30
co , & !< counter in integration point component loop
ip , & !< counter in integration point loop
el !< counter in element loop
2020-12-20 15:18:13 +05:30
!$OMP PARALLEL DO
2020-12-23 11:44:07 +05:30
do el = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
do ip = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
do co = 1 , homogenization_Nconstituents ( material_homogenizationAt ( el ) )
call crystallite_orientation ( co , ip , el ) % fromMatrix ( transpose ( math_rotationalPart ( crystallite_Fe ( 1 : 3 , 1 : 3 , co , ip , el ) ) ) )
2020-12-20 15:18:13 +05:30
enddo ; enddo ; enddo
!$OMP END PARALLEL DO
nonlocalPresent : if ( any ( plasticState % nonlocal ) ) then
!$OMP PARALLEL DO
2020-12-23 11:44:07 +05:30
do el = FEsolving_execElem ( 1 ) , FEsolving_execElem ( 2 )
if ( plasticState ( material_phaseAt ( 1 , el ) ) % nonlocal ) then
do ip = FEsolving_execIP ( 1 ) , FEsolving_execIP ( 2 )
2020-12-20 15:18:13 +05:30
call plastic_nonlocal_updateCompatibility ( crystallite_orientation , &
2020-12-23 11:44:07 +05:30
phase_plasticityInstance ( material_phaseAt ( 1 , el ) ) , ip , el )
2020-12-20 15:18:13 +05:30
enddo
endif
enddo
!$OMP END PARALLEL DO
endif nonlocalPresent
end subroutine crystallite_orientations
!--------------------------------------------------------------------------------------------------
!> @brief Map 2nd order tensor to reference config
!--------------------------------------------------------------------------------------------------
2020-12-23 11:37:18 +05:30
function crystallite_push33ToRef ( co , ip , el , tensor33 )
2020-12-20 15:18:13 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: tensor33
real ( pReal ) , dimension ( 3 , 3 ) :: T
integer , intent ( in ) :: &
el , &
ip , &
2020-12-23 11:37:18 +05:30
co
2020-12-23 12:44:42 +05:30
2020-12-23 12:42:56 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: crystallite_push33ToRef
2020-12-20 15:18:13 +05:30
2020-12-23 12:44:42 +05:30
2020-12-23 11:37:18 +05:30
T = matmul ( material_orientation0 ( co , ip , el ) % asMatrix ( ) , & ! ToDo: initial orientation correct?
transpose ( math_inv33 ( crystallite_subF ( 1 : 3 , 1 : 3 , co , ip , el ) ) ) )
2020-12-20 15:18:13 +05:30
crystallite_push33ToRef = matmul ( transpose ( T ) , matmul ( tensor33 , T ) )
end function crystallite_push33ToRef
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
subroutine integrateSourceState ( co , ip , el )
2020-12-20 15:18:13 +05:30
integer , intent ( in ) :: &
2020-12-23 11:44:07 +05:30
el , & !< element index in element loop
ip , & !< integration point index in ip loop
co !< grain index in grain loop
2020-12-23 12:44:42 +05:30
2020-12-20 15:18:13 +05:30
integer :: &
NiterationState , & !< number of iterations in state loop
2020-12-23 11:44:07 +05:30
ph , &
2020-12-23 12:44:42 +05:30
me , &
2020-12-23 18:33:15 +05:30
so
2020-12-20 15:18:13 +05:30
integer , dimension ( maxval ( phase_Nsources ) ) :: &
size_so
real ( pReal ) :: &
zeta
2020-12-23 18:33:15 +05:30
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState ) :: &
2020-12-20 15:18:13 +05:30
r ! state residuum
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , 2 , maxval ( phase_Nsources ) ) :: source_dotState
logical :: &
broken
2020-12-23 12:42:56 +05:30
2020-12-23 11:44:07 +05:30
ph = material_phaseAt ( co , el )
2020-12-23 12:44:42 +05:30
me = material_phaseMemberAt ( co , ip , el )
2020-12-20 15:18:13 +05:30
2020-12-23 12:44:42 +05:30
broken = constitutive_thermal_collectDotState ( ph , me )
broken = broken . or . constitutive_damage_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) , co , ip , el , ph , me )
2020-12-20 15:18:13 +05:30
if ( broken ) return
2020-12-23 12:44:42 +05:30
do so = 1 , phase_Nsources ( ph )
size_so ( so ) = sourceState ( ph ) % p ( so ) % sizeDotState
sourceState ( ph ) % p ( so ) % state ( 1 : size_so ( so ) , me ) = sourceState ( ph ) % p ( so ) % subState0 ( 1 : size_so ( so ) , me ) &
+ sourceState ( ph ) % p ( so ) % dotState ( 1 : size_so ( so ) , me ) &
2020-12-23 11:44:07 +05:30
* crystallite_subdt ( co , ip , el )
2020-12-23 12:44:42 +05:30
source_dotState ( 1 : size_so ( so ) , 2 , so ) = 0.0_pReal
2020-12-20 15:18:13 +05:30
enddo
iteration : do NiterationState = 1 , num % nState
2020-12-23 12:44:42 +05:30
do so = 1 , phase_Nsources ( ph )
if ( nIterationState > 1 ) source_dotState ( 1 : size_so ( so ) , 2 , so ) = source_dotState ( 1 : size_so ( so ) , 1 , so )
source_dotState ( 1 : size_so ( so ) , 1 , so ) = sourceState ( ph ) % p ( so ) % dotState ( : , me )
2020-12-20 15:18:13 +05:30
enddo
2020-12-23 12:44:42 +05:30
broken = constitutive_thermal_collectDotState ( ph , me )
broken = broken . or . constitutive_damage_collectDotState ( crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) , co , ip , el , ph , me )
2020-12-20 15:18:13 +05:30
if ( broken ) exit iteration
2020-12-23 12:44:42 +05:30
do so = 1 , phase_Nsources ( ph )
zeta = damper ( sourceState ( ph ) % p ( so ) % dotState ( : , me ) , &
source_dotState ( 1 : size_so ( so ) , 1 , so ) , &
source_dotState ( 1 : size_so ( so ) , 2 , so ) )
sourceState ( ph ) % p ( so ) % dotState ( : , me ) = sourceState ( ph ) % p ( so ) % dotState ( : , me ) * zeta &
+ source_dotState ( 1 : size_so ( so ) , 1 , so ) * ( 1.0_pReal - zeta )
r ( 1 : size_so ( so ) ) = sourceState ( ph ) % p ( so ) % state ( 1 : size_so ( so ) , me ) &
- sourceState ( ph ) % p ( so ) % subState0 ( 1 : size_so ( so ) , me ) &
- sourceState ( ph ) % p ( so ) % dotState ( 1 : size_so ( so ) , me ) * crystallite_subdt ( co , ip , el )
sourceState ( ph ) % p ( so ) % state ( 1 : size_so ( so ) , me ) = sourceState ( ph ) % p ( so ) % state ( 1 : size_so ( so ) , me ) &
- r ( 1 : size_so ( so ) )
2020-12-23 11:44:07 +05:30
crystallite_converged ( co , ip , el ) = &
2020-12-23 12:44:42 +05:30
crystallite_converged ( co , ip , el ) . and . converged ( r ( 1 : size_so ( so ) ) , &
sourceState ( ph ) % p ( so ) % state ( 1 : size_so ( so ) , me ) , &
sourceState ( ph ) % p ( so ) % atol ( 1 : size_so ( so ) ) )
2020-12-20 15:18:13 +05:30
enddo
2020-12-23 11:44:07 +05:30
if ( crystallite_converged ( co , ip , el ) ) then
2020-12-23 12:44:42 +05:30
broken = constitutive_damage_deltaState ( crystallite_Fe ( 1 : 3 , 1 : 3 , co , ip , el ) , co , ip , el , ph , me )
2020-12-20 15:18:13 +05:30
exit iteration
endif
enddo iteration
contains
!--------------------------------------------------------------------------------------------------
!> @brief calculate the damping for correction of state and dot state
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function damper ( current , previous , previous2 )
real ( pReal ) , dimension ( : ) , intent ( in ) :: &
current , previous , previous2
real ( pReal ) :: dot_prod12 , dot_prod22
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
end function damper
end subroutine integrateSourceState
!--------------------------------------------------------------------------------------------------
!> @brief determines whether a point is converged
!--------------------------------------------------------------------------------------------------
logical pure function converged ( residuum , state , atol )
real ( pReal ) , intent ( in ) , dimension ( : ) :: &
residuum , state , atol
real ( pReal ) :: &
rTol
rTol = num % rTol_crystalliteState
converged = all ( abs ( residuum ) < = max ( atol , rtol * abs ( state ) ) )
end function converged
!--------------------------------------------------------------------------------------------------
!> @brief Write current restart information (Field and constitutive data) to file.
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restartWrite
integer :: i
integer ( HID_T ) :: fileHandle , groupHandle
character ( len = pStringLen ) :: fileName , datasetName
print * , ' writing field and constitutive data required for restart to file' ; flush ( IO_STDOUT )
write ( fileName , '(a,i0,a)' ) trim ( getSolverJobName ( ) ) / / '_' , worldrank , '.hdf5'
fileHandle = HDF5_openFile ( fileName , 'a' )
call HDF5_write ( fileHandle , crystallite_partitionedF , 'F' )
call HDF5_write ( fileHandle , crystallite_Lp , 'L_p' )
call HDF5_write ( fileHandle , crystallite_S , 'S' )
groupHandle = HDF5_addGroup ( fileHandle , 'phase' )
do i = 1 , size ( material_name_phase )
write ( datasetName , '(i0,a)' ) i , '_omega'
call HDF5_write ( groupHandle , plasticState ( i ) % state , datasetName )
2020-12-21 00:50:39 +05:30
write ( datasetName , '(i0,a)' ) i , '_F_i'
call HDF5_write ( groupHandle , constitutive_mech_Fi ( i ) % data , datasetName )
2020-12-21 14:29:13 +05:30
write ( datasetName , '(i0,a)' ) i , '_L_i'
call HDF5_write ( groupHandle , constitutive_mech_Li ( i ) % data , datasetName )
2020-12-23 02:51:11 +05:30
write ( datasetName , '(i0,a)' ) i , '_F_p'
call HDF5_write ( groupHandle , constitutive_mech_Fp ( i ) % data , datasetName )
2020-12-20 15:18:13 +05:30
enddo
call HDF5_closeGroup ( groupHandle )
groupHandle = HDF5_addGroup ( fileHandle , 'homogenization' )
do i = 1 , size ( material_name_homogenization )
write ( datasetName , '(i0,a)' ) i , '_omega'
call HDF5_write ( groupHandle , homogState ( i ) % state , datasetName )
enddo
call HDF5_closeGroup ( groupHandle )
call HDF5_closeFile ( fileHandle )
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
print '(/,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_Lp0 , 'L_p' )
call HDF5_read ( fileHandle , crystallite_S0 , 'S' )
groupHandle = HDF5_openGroup ( fileHandle , 'phase' )
do i = 1 , size ( material_name_phase )
write ( datasetName , '(i0,a)' ) i , '_omega'
call HDF5_read ( groupHandle , plasticState ( i ) % state0 , datasetName )
2020-12-21 00:50:39 +05:30
write ( datasetName , '(i0,a)' ) i , '_F_i'
call HDF5_read ( groupHandle , constitutive_mech_Fi0 ( i ) % data , datasetName )
2020-12-21 14:29:13 +05:30
write ( datasetName , '(i0,a)' ) i , '_L_i'
call HDF5_read ( groupHandle , constitutive_mech_Li0 ( i ) % data , datasetName )
2020-12-23 02:51:11 +05:30
write ( datasetName , '(i0,a)' ) i , '_F_p'
call HDF5_read ( groupHandle , constitutive_mech_Fp0 ( i ) % data , datasetName )
2020-12-20 15:18:13 +05:30
enddo
call HDF5_closeGroup ( groupHandle )
groupHandle = HDF5_openGroup ( fileHandle , 'homogenization' )
do i = 1 , size ( material_name_homogenization )
write ( datasetName , '(i0,a)' ) i , '_omega'
call HDF5_read ( groupHandle , homogState ( i ) % state0 , datasetName )
enddo
call HDF5_closeGroup ( groupHandle )
call HDF5_closeFile ( fileHandle )
end subroutine crystallite_restartRead
2013-02-11 16:13:45 +05:30
end module constitutive
