223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'development' into new-gmsh-version

This commit is contained in:
Sharan Roongta 2021-01-12 14:13:39 +01:00
parent dc55aa002a 1e50fcc770
commit c5ec5dce03
23 changed files with 1036 additions and 728 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
v3.0.0-alpha2-173-g584c7cc3a
v3.0.0-alpha2-230-g0fc670d01

3
src/DAMASK_marc.f90
View File

@ -364,7 +364,8 @@ subroutine flux(f,ts,n,time)
real(pReal), dimension(2), intent(out) :: &
f
call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1)))
f(2) = 0.0_pReal
call thermal_conduction_getSource(f(1), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1)))
end subroutine flux

5
src/commercialFEM_fileList.f90
View File

@ -32,8 +32,8 @@
#include "constitutive_plastic_disloTungsten.f90"
#include "constitutive_plastic_nonlocal.f90"
#include "constitutive_thermal.f90"
#include "source_thermal_dissipation.f90"
#include "source_thermal_externalheat.f90"
#include "constitutive_thermal_dissipation.f90"
#include "constitutive_thermal_externalheat.f90"
#include "kinematics_thermal_expansion.f90"
#include "constitutive_damage.f90"
#include "source_damage_isoBrittle.f90"
@ -51,4 +51,5 @@
#include "homogenization_mech_none.f90"
#include "homogenization_mech_isostrain.f90"
#include "homogenization_mech_RGC.f90"
#include "homogenization_thermal.f90"
#include "CPFEM.f90"

568
src/constitutive.f90
View File

@ -21,21 +21,6 @@ module constitutive
private
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, &
@ -81,15 +66,11 @@ module constitutive
type(tDebugOptions) :: debugCrystallite
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public :: &
phase_plasticity !< plasticity of each phase
integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: &
phase_source, & !< active sources mechanisms of each phase
integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: &
phase_kinematics !< active kinematic mechanisms of each phase
integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler)
thermal_Nsources, &
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
@ -102,7 +83,7 @@ module constitutive
type(tPlasticState), allocatable, dimension(:), public :: &
plasticState
type(tSourceState), allocatable, dimension(:), public :: &
sourceState
damageState, thermalState
integer, public, protected :: &
@ -112,13 +93,15 @@ module constitutive
interface
! == cleaned:begin =================================================================================
module subroutine mech_init
module subroutine mech_init(phases)
class(tNode), pointer :: phases
end subroutine mech_init
module subroutine damage_init
end subroutine damage_init
module subroutine thermal_init
module subroutine thermal_init(phases)
class(tNode), pointer :: phases
end subroutine thermal_init
@ -137,21 +120,37 @@ module constitutive
integer, intent(in) :: ph, me
end subroutine mech_initializeRestorationPoints
module subroutine constitutive_mech_windForward(ph,me)
module subroutine thermal_initializeRestorationPoints(ph,me)
integer, intent(in) :: ph, me
end subroutine constitutive_mech_windForward
end subroutine thermal_initializeRestorationPoints
module subroutine mech_windForward(ph,me)
integer, intent(in) :: ph, me
end subroutine mech_windForward
module subroutine thermal_windForward(ph,me)
integer, intent(in) :: ph, me
end subroutine thermal_windForward
module subroutine mech_forward()
end subroutine mech_forward
module subroutine thermal_forward()
end subroutine thermal_forward
module subroutine constitutive_mech_forward
end subroutine constitutive_mech_forward
module subroutine mech_restore(ip,el,includeL)
integer, intent(in) :: &
ip, &
el
logical, intent(in) :: &
includeL
integer, intent(in) :: ip, el
logical, intent(in) :: includeL
end subroutine mech_restore
module subroutine thermal_restore(ip,el)
integer, intent(in) :: ip, el
end subroutine thermal_restore
module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
@ -172,43 +171,67 @@ module constitutive
end subroutine mech_restartRead
module function constitutive_mech_getS(co,ip,el) result(S)
integer, intent(in) :: co, ip, el
module function mech_S(ph,me) result(S)
integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: S
end function constitutive_mech_getS
end function mech_S
module function constitutive_mech_getLp(co,ip,el) result(Lp)
integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: Lp
end function constitutive_mech_getLp
module function mech_L_p(ph,me) result(L_p)
integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: L_p
end function mech_L_p
module function constitutive_mech_getF(co,ip,el) result(F)
integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: F
end function constitutive_mech_getF
module function constitutive_mech_getF_e(co,ip,el) result(F_e)
integer, intent(in) :: co, ip, el
module function mech_F_e(ph,me) result(F_e)
integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: F_e
end function constitutive_mech_getF_e
end function mech_F_e
module function constitutive_mech_getP(co,ip,el) result(P)
integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: P
end function constitutive_mech_getP
module function constitutive_thermal_T(co,ip,el) result(T)
integer, intent(in) :: co, ip, el
module function thermal_T(ph,me) result(T)
integer, intent(in) :: ph,me
real(pReal) :: T
end function constitutive_thermal_T
end function thermal_T
module subroutine constitutive_mech_setF(F,co,ip,el)
real(pReal), dimension(3,3), intent(in) :: F
integer, intent(in) :: co, ip, el
end subroutine constitutive_mech_setF
module subroutine constitutive_thermal_setT(T,co,ip,el)
real(pReal), intent(in) :: T
integer, intent(in) :: co, ip, el
end subroutine constitutive_thermal_setT
! == cleaned:end ===================================================================================
module function integrateThermalState(Delta_t,co,ip,el) result(broken)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: broken
end function integrateThermalState
module function integrateDamageState(dt,co,ip,el) result(broken)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: broken
end function integrateDamageState
module function crystallite_stress(dt,co,ip,el) result(converged_)
real(pReal), intent(in) :: dt
integer, intent(in) :: co, ip, el
@ -260,16 +283,15 @@ module constitutive
dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el)
module subroutine constitutive_thermal_getRate(TDot, T,ip,el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
T
real(pReal), intent(inout) :: &
TDot, &
dTDot_dT
end subroutine constitutive_thermal_getRateAndItsTangents
real(pReal), intent(out) :: &
TDot
end subroutine constitutive_thermal_getRate
@ -282,6 +304,7 @@ module constitutive
orientation !< crystal orientation
end subroutine plastic_nonlocal_updateCompatibility
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
@ -344,23 +367,6 @@ module constitutive
end subroutine source_damage_isoBrittle_deltaState
module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
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
module subroutine constitutive_plastic_dependentState(co,ip,el)
integer, intent(in) :: &
co, & !< component-ID of integration point
@ -368,23 +374,6 @@ module constitutive
el !< element
end subroutine constitutive_plastic_dependentState
module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
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
end interface
@ -394,15 +383,13 @@ module constitutive
public :: &
constitutive_init, &
constitutive_homogenizedC, &
constitutive_LiAndItsTangents, &
constitutive_damage_getRateAndItsTangents, &
constitutive_thermal_getRateAndItsTangents, &
constitutive_thermal_getRate, &
constitutive_results, &
constitutive_allocateState, &
constitutive_forward, &
constitutive_restore, &
plastic_nonlocal_updateCompatibility, &
source_active, &
kinematics_active, &
converged, &
crystallite_init, &
@ -412,29 +399,14 @@ module constitutive
crystallite_push33ToRef, &
constitutive_restartWrite, &
constitutive_restartRead, &
integrateSourceState, &
constitutive_mech_setF, &
integrateThermalState, &
integrateDamageState, &
constitutive_thermal_setT, &
constitutive_mech_getP, &
constitutive_mech_getLp, &
constitutive_mech_setF, &
constitutive_mech_getF, &
constitutive_mech_getS, &
constitutive_initializeRestorationPoints, &
constitutive_windForward, &
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, &
@ -456,6 +428,8 @@ subroutine constitutive_init
phases
print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT)
debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList)
debugConstitutive%basic = debug_constitutive%contains('basic')
debugConstitutive%extensive = debug_constitutive%contains('extensive')
@ -464,15 +438,14 @@ subroutine constitutive_init
debugConstitutive%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1)
debugConstitutive%grain = config_debug%get_asInt('grain',defaultVal = 1)
!--------------------------------------------------------------------------------------------------
! initialize constitutive laws
print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT)
call mech_init
call damage_init
call thermal_init
phases => config_material%get('phase')
call mech_init(phases)
call damage_init
call thermal_init(phases)
constitutive_source_maxSizeDotState = 0
PhaseLoop2:do ph = 1,phases%length
!--------------------------------------------------------------------------------------------------
@ -480,48 +453,18 @@ subroutine constitutive_init
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
damageState(ph)%p(so)%partitionedState0 = damageState(ph)%p(so)%state0
damageState(ph)%p(so)%state = damageState(ph)%p(so)%partitionedState0
end forall
constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, &
maxval(sourceState(ph)%p%sizeDotState))
maxval(damageState(ph)%p%sizeDotState))
enddo PhaseLoop2
constitutive_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
end subroutine constitutive_init
!--------------------------------------------------------------------------------------------------
!> @brief checks if a source mechanism is active or not
!--------------------------------------------------------------------------------------------------
function source_active(source_label,src_length) result(active_source)
character(len=*), intent(in) :: source_label !< name of source mechanism
integer, intent(in) :: src_length !< max. number of sources in system
logical, dimension(:,:), allocatable :: active_source
class(tNode), pointer :: &
phases, &
phase, &
sources, &
src
integer :: p,s
phases => config_material%get('phase')
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
@ -554,197 +497,6 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki
end function kinematics_active
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi?
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
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
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
Li = 0.0_pReal
dLi_dS = 0.0_pReal
dLi_dFi = 0.0_pReal
plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_isotropic_ID) plasticityType
of = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(material_phaseAt(co,el))
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of)
case default plasticityType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select plasticityType
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el))
kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
case (KINEMATICS_slipplane_opening_ID) kinematicsType
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
case (KINEMATICS_thermal_expansion_ID) kinematicsType
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el)
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
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)
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)
enddo; enddo
end subroutine constitutive_LiAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
of
integer :: &
so !< counter in source loop
logical :: broken
broken = .false.
SourceLoop: do so = 1, phase_Nsources(ph)
sourceType: select case (phase_source(so,ph))
case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState(constitutive_mech_getS(co,ip,el), co, ip, el) ! correct stress?
case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_dotState(co, ip, el)
case (SOURCE_damage_anisoDuctile_ID) sourceType
call source_damage_anisoDuctile_dotState(co, ip, el)
end select sourceType
broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of)))
enddo SourceLoop
end function constitutive_damage_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function constitutive_thermal_collectDotState(ph,me) result(broken)
integer, intent(in) :: ph, me
logical :: broken
integer :: i
broken = .false.
SourceLoop: do i = 1, phase_Nsources(ph)
if (phase_source(i,ph) == SOURCE_thermal_externalheat_ID) &
call source_thermal_externalheat_dotState(ph,me)
broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(i)%dotState(:,me)))
enddo SourceLoop
end function constitutive_thermal_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
of
real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient
integer :: &
so, &
myOffset, &
mySize
logical :: &
broken
broken = .false.
sourceLoop: do so = 1, phase_Nsources(ph)
sourceType: select case (phase_source(so,ph))
case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, &
co, ip, el)
broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of)))
if(.not. broken) then
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)
endif
end select sourceType
enddo SourceLoop
end function constitutive_damage_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief Allocate the components of the state structure for a given phase
!--------------------------------------------------------------------------------------------------
@ -795,12 +547,13 @@ subroutine constitutive_restore(ip,el,includeL)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
do so = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = &
sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el))
damageState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = &
damageState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el))
enddo
enddo
call mech_restore(ip,el,includeL)
call thermal_restore(ip,el)
end subroutine constitutive_restore
@ -809,16 +562,17 @@ end subroutine constitutive_restore
!> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
subroutine constitutive_forward
subroutine constitutive_forward()
integer :: ph, so
call constitutive_mech_forward()
call mech_forward()
call thermal_forward()
do ph = 1, size(sourceState)
do ph = 1, size(damageState)
do so = 1,phase_Nsources(ph)
sourceState(ph)%p(so)%state0 = sourceState(ph)%p(so)%state
damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state
enddo; enddo
end subroutine constitutive_forward
@ -827,7 +581,7 @@ end subroutine constitutive_forward
!--------------------------------------------------------------------------------------------------
!> @brief writes constitutive results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine constitutive_results
subroutine constitutive_results()
integer :: ph
character(len=:), allocatable :: group
@ -851,7 +605,7 @@ end subroutine constitutive_results
!--------------------------------------------------------------------------------------------------
!> @brief allocates and initialize per grain variables
!--------------------------------------------------------------------------------------------------
subroutine crystallite_init
subroutine crystallite_init()
integer :: &
ph, &
@ -919,7 +673,10 @@ subroutine crystallite_init
do ph = 1, phases%length
do so = 1, phase_Nsources(ph)
allocate(sourceState(ph)%p(so)%subState0,source=sourceState(ph)%p(so)%state0) ! ToDo: hack
allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0) ! ToDo: hack
enddo
do so = 1, thermal_Nsources(ph)
allocate(thermalState(ph)%p(so)%subState0,source=thermalState(ph)%p(so)%state0) ! ToDo: hack
enddo
enddo
@ -963,10 +720,12 @@ subroutine constitutive_initializeRestorationPoints(ip,el)
me = material_phaseMemberAt(co,ip,el)
call mech_initializeRestorationPoints(ph,me)
call thermal_initializeRestorationPoints(ph,me)
do so = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state0(:,me)
do so = 1, size(damageState(ph)%p)
damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state0(:,me)
enddo
enddo
end subroutine constitutive_initializeRestorationPoints
@ -990,10 +749,13 @@ subroutine constitutive_windForward(ip,el)
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call constitutive_mech_windForward(ph,me)
call mech_windForward(ph,me)
call thermal_windForward(ph,me)
do so = 1, phase_Nsources(material_phaseAt(co,el))
sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state(:,me)
damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state(:,me)
enddo
enddo
end subroutine constitutive_windForward
@ -1011,7 +773,7 @@ subroutine crystallite_orientations(co,ip,el)
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(&
constitutive_mech_getF_e(co,ip,el))))
mech_F_e(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)))))
if (plasticState(material_phaseAt(1,el))%nonlocal) &
call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
@ -1043,109 +805,7 @@ function crystallite_push33ToRef(co,ip,el, tensor33)
end function crystallite_push33ToRef
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
function integrateSourceState(dt,co,ip,el) result(broken)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
integer :: &
NiterationState, & !< number of iterations in state loop
ph, &
me, &
so
integer, dimension(maxval(phase_Nsources)) :: &
size_so
real(pReal) :: &
zeta
real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
r ! state residuum
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: &
broken, converged_
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
converged_ = .true.
broken = constitutive_thermal_collectDotState(ph,me)
broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) return
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) * dt
source_dotState(1:size_so(so),2,so) = 0.0_pReal
enddo
iteration: do NiterationState = 1, num%nState
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)
enddo
broken = constitutive_thermal_collectDotState(ph,me)
broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) exit iteration
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) * dt
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))
converged_ = converged_ .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)))
enddo
if(converged_) then
broken = constitutive_damage_deltaState(constitutive_mech_getF_e(co,ip,el),co,ip,el,ph,me)
exit iteration
endif
enddo iteration
broken = broken .or. .not. converged_
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 function integrateSourceState
!--------------------------------------------------------------------------------------------------

274
src/constitutive_damage.f90
View File

@ -2,6 +2,16 @@
!> @brief internal microstructure state for all damage sources and kinematics constitutive models
!----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_damage
enum, bind(c); enumerator :: &
DAMAGE_UNDEFINED_ID, &
DAMAGE_ISOBRITTLE_ID, &
DAMAGE_ISODUCTILE_ID, &
DAMAGE_ANISOBRITTLE_ID, &
DAMAGE_ANISODUCTILE_ID
end enum
integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: &
phase_source !< active sources mechanisms of each phase
interface
@ -119,24 +129,24 @@ module subroutine damage_init
phases => config_material%get('phase')
allocate(sourceState (phases%length))
allocate(damageState (phases%length))
allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics
do ph = 1,phases%length
phase => phases%get(ph)
sources => phase%get('source',defaultVal=emptyList)
phase_Nsources(ph) = sources%length
allocate(sourceState(ph)%p(phase_Nsources(ph)))
allocate(damageState(ph)%p(phase_Nsources(ph)))
enddo
allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID)
allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID)
! initialize source mechanisms
if(maxval(phase_Nsources) /= 0) then
where(source_damage_isoBrittle_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_isoBrittle_ID
where(source_damage_isoDuctile_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_isoDuctile_ID
where(source_damage_anisoBrittle_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_anisoBrittle_ID
where(source_damage_anisoDuctile_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_anisoDuctile_ID
where(source_damage_isoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID
where(source_damage_isoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID
where(source_damage_anisoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID
where(source_damage_anisoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID
endif
!--------------------------------------------------------------------------------------------------
@ -189,16 +199,16 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)
case (DAMAGE_ISOBRITTLE_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_isoDuctile_ID)
case (DAMAGE_ISODUCTILE_ID)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoBrittle_ID)
case (DAMAGE_ANISOBRITTLE_ID)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoDuctile_ID)
case (DAMAGE_ANISODUCTILE_ID)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case default
@ -214,6 +224,111 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
module function integrateDamageState(dt,co,ip,el) result(broken)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: broken
integer :: &
NiterationState, & !< number of iterations in state loop
ph, &
me, &
so
integer, dimension(maxval(phase_Nsources)) :: &
size_so
real(pReal) :: &
zeta
real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
r ! state residuum
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: &
converged_
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
converged_ = .true.
broken = constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) return
do so = 1, phase_Nsources(ph)
size_so(so) = damageState(ph)%p(so)%sizeDotState
damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%subState0(1:size_so(so),me) &
+ damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt
source_dotState(1:size_so(so),2,so) = 0.0_pReal
enddo
iteration: do NiterationState = 1, num%nState
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) = damageState(ph)%p(so)%dotState(:,me)
enddo
broken = constitutive_damage_collectDotState(co,ip,el,ph,me)
if(broken) exit iteration
do so = 1, phase_Nsources(ph)
zeta = damper(damageState(ph)%p(so)%dotState(:,me), &
source_dotState(1:size_so(so),1,so),&
source_dotState(1:size_so(so),2,so))
damageState(ph)%p(so)%dotState(:,me) = damageState(ph)%p(so)%dotState(:,me) * zeta &
+ source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta)
r(1:size_so(so)) = damageState(ph)%p(so)%state (1:size_so(so),me) &
- damageState(ph)%p(so)%subState0(1:size_so(so),me) &
- damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt
damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%state(1:size_so(so),me) &
- r(1:size_so(so))
converged_ = converged_ .and. converged(r(1:size_so(so)), &
damageState(ph)%p(so)%state(1:size_so(so),me), &
damageState(ph)%p(so)%atol(1:size_so(so)))
enddo
if(converged_) then
broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me)
exit iteration
endif
enddo iteration
broken = broken .or. .not. converged_
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 function integrateDamageState
!----------------------------------------------------------------------------------------------
!< @brief writes damage sources results to HDF5 output file
!----------------------------------------------------------------------------------------------
@ -226,23 +341,23 @@ module subroutine damage_results(group,ph)
sourceLoop: do so = 1, phase_Nsources(ph)
if (phase_source(so,ph) /= SOURCE_UNDEFINED_ID) &
if (phase_source(so,ph) /= DAMAGE_UNDEFINED_ID) &
call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage'
sourceType: select case (phase_source(so,ph))
case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_results(ph,group//'sources/')
case (SOURCE_damage_anisoDuctile_ID) sourceType
call source_damage_anisoDuctile_results(ph,group//'sources/')
case (SOURCE_damage_isoBrittle_ID) sourceType
case (DAMAGE_ISOBRITTLE_ID) sourceType
call source_damage_isoBrittle_results(ph,group//'sources/')
case (SOURCE_damage_isoDuctile_ID) sourceType
case (DAMAGE_ISODUCTILE_ID) sourceType
call source_damage_isoDuctile_results(ph,group//'sources/')
case (DAMAGE_ANISOBRITTLE_ID) sourceType
call source_damage_anisoBrittle_results(ph,group//'sources/')
case (DAMAGE_ANISODUCTILE_ID) sourceType
call source_damage_anisoDuctile_results(ph,group//'sources/')
end select sourceType
enddo SourceLoop
@ -250,4 +365,123 @@ module subroutine damage_results(group,ph)
end subroutine damage_results
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
of
integer :: &
so !< counter in source loop
logical :: broken
broken = .false.
SourceLoop: do so = 1, phase_Nsources(ph)
sourceType: select case (phase_source(so,ph))
case (DAMAGE_ISODUCTILE_ID) sourceType
call source_damage_isoDuctile_dotState(co, ip, el)
case (DAMAGE_ANISODUCTILE_ID) sourceType
call source_damage_anisoDuctile_dotState(co, ip, el)
case (DAMAGE_ANISOBRITTLE_ID) sourceType
call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),&
co, ip, el) ! correct stress?
end select sourceType
broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,of)))
enddo SourceLoop
end function constitutive_damage_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken)
integer, intent(in) :: &
co, & !< component-ID of integration point
ip, & !< integration point
el, & !< element
ph, &
of
real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient
integer :: &
so, &
myOffset, &
mySize
logical :: &
broken
broken = .false.
sourceLoop: do so = 1, phase_Nsources(ph)
sourceType: select case (phase_source(so,ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, &
co, ip, el)
broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,of)))
if(.not. broken) then
myOffset = damageState(ph)%p(so)%offsetDeltaState
mySize = damageState(ph)%p(so)%sizeDeltaState
damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = &
damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + damageState(ph)%p(so)%deltaState(1:mySize,of)
endif
end select sourceType
enddo SourceLoop
end function constitutive_damage_deltaState
!--------------------------------------------------------------------------------------------------
!> @brief checks if a source mechanism is active or not
!--------------------------------------------------------------------------------------------------
function source_active(source_label,src_length) result(active_source)
character(len=*), intent(in) :: source_label !< name of source mechanism
integer, intent(in) :: src_length !< max. number of sources in system
logical, dimension(:,:), allocatable :: active_source
class(tNode), pointer :: &
phases, &
phase, &
sources, &
src
integer :: p,s
phases => config_material%get('phase')
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
end submodule constitutive_damage

186
src/constitutive_mech.f90
View File

@ -7,12 +7,20 @@ submodule(constitutive) constitutive_mech
ELASTICITY_UNDEFINED_ID, &
ELASTICITY_HOOKE_ID, &
STIFFNESS_DEGRADATION_UNDEFINED_ID, &
STIFFNESS_DEGRADATION_DAMAGE_ID
STIFFNESS_DEGRADATION_DAMAGE_ID, &
PLASTICITY_UNDEFINED_ID, &
PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_PHENOPOWERLAW_ID, &
PLASTICITY_KINEHARDENING_ID, &
PLASTICITY_DISLOTWIN_ID, &
PLASTICITY_DISLOTUNGSTEN_ID, &
PLASTICITY_NONLOCAL_ID
end enum
integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable :: &
integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: &
phase_elasticity !< elasticity of each phase
integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: &
integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: &
phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase
type(tTensorContainer), dimension(:), allocatable :: &
@ -41,6 +49,12 @@ submodule(constitutive) constitutive_mech
constitutive_mech_partitionedS0
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: &
phase_plasticity !< plasticity of each phase
interface
module function plastic_none_init() result(myPlasticity)
@ -319,7 +333,10 @@ contains
!> @brief Initialize mechanical field related constitutive models
!> @details Initialize elasticity, plasticity and stiffness degradation models.
!--------------------------------------------------------------------------------------------------
module subroutine mech_init
module subroutine mech_init(phases)
class(tNode), pointer :: &
phases
integer :: &
el, &
@ -331,7 +348,6 @@ module subroutine mech_init
Nconstituents
class(tNode), pointer :: &
num_crystallite, &
phases, &
phase, &
mech, &
elastic, &
@ -341,7 +357,6 @@ module subroutine mech_init
!-------------------------------------------------------------------------------------------------
! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC?
phases => config_material%get('phase')
allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
allocate(phase_elasticityInstance(phases%length), source = 0)
allocate(phase_NstiffnessDegradations(phases%length),source=0)
@ -529,7 +544,7 @@ end function plastic_active
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermediate deformation gradients using Hooke's law
!--------------------------------------------------------------------------------------------------
module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, co, ip, el)
integer, intent(in) :: &
@ -615,7 +630,7 @@ end subroutine constitutive_plastic_dependentState
! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e.
! Mp in, dLp_dMp out
!--------------------------------------------------------------------------------------------------
module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
@ -1485,7 +1500,7 @@ end subroutine mech_initializeRestorationPoints
!--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward.
!--------------------------------------------------------------------------------------------------
module subroutine constitutive_mech_windForward(ph,me)
module subroutine mech_windForward(ph,me)
integer, intent(in) :: ph, me
@ -1499,14 +1514,14 @@ module subroutine constitutive_mech_windForward(ph,me)
plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state(:,me)
end subroutine constitutive_mech_windForward
end subroutine mech_windForward
!--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
module subroutine constitutive_mech_forward()
module subroutine mech_forward()
integer :: ph
@ -1521,7 +1536,7 @@ module subroutine constitutive_mech_forward()
plasticState(ph)%state0 = plasticState(ph)%state
enddo
end subroutine constitutive_mech_forward
end subroutine mech_forward
@ -1585,7 +1600,10 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
do so = 1, phase_Nsources(ph)
sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%partitionedState0(:,me)
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%partitionedState0(:,me)
enddo
do so = 1, thermal_Nsources(ph)
thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me)
enddo
subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)
subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)
@ -1613,7 +1631,10 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
subFi0 = constitutive_mech_Fi(ph)%data(1:3,1:3,me)
subState0 = plasticState(ph)%state(:,me)
do so = 1, phase_Nsources(ph)
sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%state(:,me)
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me)
enddo
do so = 1, thermal_Nsources(ph)
thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%state(:,me)
enddo
endif
!--------------------------------------------------------------------------------------------------
@ -1629,7 +1650,10 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
endif
plasticState(ph)%state(:,me) = subState0
do so = 1, phase_Nsources(ph)
sourceState(ph)%p(so)%state(:,me) = sourceState(ph)%p(so)%subState0(:,me)
damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me)
enddo
do so = 1, thermal_Nsources(ph)
thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me)
enddo
todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair)
@ -1643,7 +1667,8 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el)
endif
enddo cutbackLooping
@ -1678,8 +1703,7 @@ module subroutine mech_restore(ip,el,includeL)
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)
constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me)
plasticState (material_phaseAt(co,el))%state( :,material_phasememberAt(co,ip,el)) = &
plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phasememberAt(co,ip,el))
plasticState(ph)%state(:,me) = plasticState(ph)%partitionedState0(:,me)
enddo
end subroutine mech_restore
@ -1846,31 +1870,37 @@ module subroutine mech_restartRead(groupHandle,ph)
end subroutine mech_restartRead
! getter for non-mech (e.g. thermal)
module function constitutive_mech_getS(co,ip,el) result(S)
!----------------------------------------------------------------------------------------------
!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics)
!----------------------------------------------------------------------------------------------
module function mech_S(ph,me) result(S)
integer, intent(in) :: co, ip, el
integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: S
S = constitutive_mech_S(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
S = constitutive_mech_S(ph)%data(1:3,1:3,me)
end function constitutive_mech_getS
end function mech_S
! getter for non-mech (e.g. thermal)
module function constitutive_mech_getLp(co,ip,el) result(Lp)
!----------------------------------------------------------------------------------------------
!< @brief Get plastic velocity gradient (for use by non-mech physics)
!----------------------------------------------------------------------------------------------
module function mech_L_p(ph,me) result(L_p)
integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: Lp
integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: L_p
Lp = constitutive_mech_Lp(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
L_p = constitutive_mech_Lp(ph)%data(1:3,1:3,me)
end function constitutive_mech_getLp
end function mech_L_p
! getter for non-mech (e.g. thermal)
!----------------------------------------------------------------------------------------------
!< @brief Get deformation gradient (for use by homogenization)
!----------------------------------------------------------------------------------------------
module function constitutive_mech_getF(co,ip,el) result(F)
integer, intent(in) :: co, ip, el
@ -1882,20 +1912,24 @@ module function constitutive_mech_getF(co,ip,el) result(F)
end function constitutive_mech_getF
! getter for non-mech (e.g. thermal)
module function constitutive_mech_getF_e(co,ip,el) result(F_e)
!----------------------------------------------------------------------------------------------
!< @brief Get elastic deformation gradient (for use by non-mech physics)
!----------------------------------------------------------------------------------------------
module function mech_F_e(ph,me) result(F_e)
integer, intent(in) :: co, ip, el
integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: F_e
F_e = constitutive_mech_Fe(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
F_e = constitutive_mech_Fe(ph)%data(1:3,1:3,me)
end function constitutive_mech_getF_e
end function mech_F_e
! getter for non-mech (e.g. thermal)
!----------------------------------------------------------------------------------------------
!< @brief Get second Piola-Kichhoff stress (for use by homogenization)
!----------------------------------------------------------------------------------------------
module function constitutive_mech_getP(co,ip,el) result(P)
integer, intent(in) :: co, ip, el
@ -1918,5 +1952,85 @@ module subroutine constitutive_mech_setF(F,co,ip,el)
end subroutine constitutive_mech_setF
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi?
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, co, ip, el)
integer, intent(in) :: &
co, & !< component-ID of integration point
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
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
Li = 0.0_pReal
dLi_dS = 0.0_pReal
dLi_dFi = 0.0_pReal
plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_isotropic_ID) plasticityType
of = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(material_phaseAt(co,el))
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of)
case default plasticityType
my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select plasticityType
Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS
KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el))
kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
case (KINEMATICS_slipplane_opening_ID) kinematicsType
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
case (KINEMATICS_thermal_expansion_ID) kinematicsType
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el)
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
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)
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)
enddo; enddo
end subroutine constitutive_LiAndItsTangents
end submodule constitutive_mech

316
src/constitutive_thermal.f90
View File

@ -3,6 +3,21 @@
!----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_thermal
enum, bind(c); enumerator :: &
THERMAL_UNDEFINED_ID ,&
THERMAL_DISSIPATION_ID, &
THERMAL_EXTERNALHEAT_ID
end enum
type :: tDataContainer
real(pReal), dimension(:), allocatable :: T
end type tDataContainer
integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: &
thermal_source
type(tDataContainer), dimension(:), allocatable :: current
integer :: thermal_source_maxSizeDotState
interface
module function source_thermal_dissipation_init(source_length) result(mySources)
@ -21,7 +36,7 @@ submodule(constitutive) constitutive_thermal
end function kinematics_thermal_expansion_init
module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase)
module subroutine thermal_dissipation_getRate(TDot, Tstar,Lp,phase)
integer, intent(in) :: &
phase !< phase ID of element
real(pReal), intent(in), dimension(3,3) :: &
@ -29,18 +44,16 @@ submodule(constitutive) constitutive_thermal
real(pReal), intent(in), dimension(3,3) :: &
Lp !< plastic velocuty gradient for a given element
real(pReal), intent(out) :: &
TDot, &
dTDot_dT
end subroutine source_thermal_dissipation_getRateAndItsTangent
TDot
end subroutine thermal_dissipation_getRate
module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of)
module subroutine thermal_externalheat_getRate(TDot, phase,of)
integer, intent(in) :: &
phase, &
of
real(pReal), intent(out) :: &
TDot, &
dTDot_dT
end subroutine source_thermal_externalheat_getRateAndItsTangent
TDot
end subroutine thermal_externalheat_getRate
end interface
@ -49,13 +62,59 @@ contains
!----------------------------------------------------------------------------------------------
!< @brief initializes thermal sources and kinematics mechanism
!----------------------------------------------------------------------------------------------
module subroutine thermal_init
module subroutine thermal_init(phases)
! initialize source mechanisms
if(maxval(phase_Nsources) /= 0) then
where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID
where(source_thermal_externalheat_init(maxval(phase_Nsources))) phase_source = SOURCE_thermal_externalheat_ID
class(tNode), pointer :: &
phases
class(tNode), pointer :: &
phase, thermal, sources
integer :: &
ph, so, &
Nconstituents
print'(/,a)', ' <<<+- constitutive_thermal init -+>>>'
allocate(current(phases%length))
allocate(thermalState (phases%length))
allocate(thermal_Nsources(phases%length),source = 0)
do ph = 1, phases%length
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
allocate(current(ph)%T(Nconstituents))
phase => phases%get(ph)
if(phase%contains('thermal')) then
thermal => phase%get('thermal')
sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length
endif
allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
enddo
allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
if(maxval(thermal_Nsources) /= 0) then
where(source_thermal_dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID
where(source_thermal_externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID
endif
thermal_source_maxSizeDotState = 0
PhaseLoop2:do ph = 1,phases%length
do so = 1,thermal_Nsources(ph)
thermalState(ph)%p(so)%partitionedState0 = thermalState(ph)%p(so)%state0
thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%partitionedState0
enddo
thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, &
maxval(thermalState(ph)%p%sizeDotState))
enddo PhaseLoop2
!--------------------------------------------------------------------------------------------------
!initialize kinematic mechanisms
@ -68,75 +127,236 @@ end subroutine thermal_init
!----------------------------------------------------------------------------------------------
!< @brief calculates thermal dissipation rate
!----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el)
module subroutine constitutive_thermal_getRate(TDot, T, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
T !< plastic velocity gradient
real(pReal), intent(inout) :: &
TDot, &
dTDot_dT
real(pReal), intent(out) :: &
TDot
real(pReal) :: &
my_Tdot, &
my_dTdot_dT
real(pReal), dimension(3,3) :: Lp, S
my_Tdot
integer :: &
phase, &
ph, &
homog, &
instance, &
grain, &
source, &
constituent
me, &
so, &
co
homog = material_homogenizationAt(el)
instance = thermal_typeInstance(homog)
do grain = 1, homogenization_Nconstituents(homog)
phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_thermal_dissipation_ID)
Lp = constitutive_mech_getLp(grain,ip,el)
S = constitutive_mech_getS(grain,ip,el)
call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
S, Lp, phase)
TDot = 0.0_pReal
do co = 1, homogenization_Nconstituents(homog)
ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el)
do so = 1, thermal_Nsources(ph)
select case(thermal_source(so,ph))
case (THERMAL_DISSIPATION_ID)
call thermal_dissipation_getRate(my_Tdot, mech_S(ph,me),mech_L_p(ph,me),ph)
case (SOURCE_thermal_externalheat_ID)
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
phase, constituent)
case (THERMAL_EXTERNALHEAT_ID)
call thermal_externalheat_getRate(my_Tdot, ph,me)
case default
my_Tdot = 0.0_pReal
my_dTdot_dT = 0.0_pReal
end select
Tdot = Tdot + my_Tdot
dTdot_dT = dTdot_dT + my_dTdot_dT
enddo
enddo
end subroutine constitutive_thermal_getRateAndItsTangents
end subroutine constitutive_thermal_getRate
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function constitutive_thermal_collectDotState(ph,me) result(broken)
integer, intent(in) :: ph, me
logical :: broken
integer :: i
broken = .false.
SourceLoop: do i = 1, thermal_Nsources(ph)
if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) &
call source_thermal_externalheat_dotState(ph,me)
broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me)))
enddo SourceLoop
end function constitutive_thermal_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
module function integrateThermalState(Delta_t,co,ip,el) result(broken)
! getter for non-thermal (e.g. mech)
module function constitutive_thermal_T(co,ip,el) result(T)
real(pReal), intent(in) :: Delta_t
integer, intent(in) :: &
el, & !< element index in element loop
ip, & !< integration point index in ip loop
co !< grain index in grain loop
logical :: &
broken
integer, intent(in) :: co, ip, el
integer :: &
ph, &
me, &
so, &
sizeDotState
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
broken = constitutive_thermal_collectDotState(ph,me)
if(broken) return
do so = 1, thermal_Nsources(ph)
sizeDotState = thermalState(ph)%p(so)%sizeDotState
thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%subState0(1:sizeDotState,me) &
+ thermalState(ph)%p(so)%dotState(1:sizeDotState,me) * Delta_t
enddo
end function integrateThermalState
module subroutine thermal_initializeRestorationPoints(ph,me)
integer, intent(in) :: ph, me
integer :: so
do so = 1, size(thermalState(ph)%p)
thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me)
enddo
end subroutine thermal_initializeRestorationPoints
module subroutine thermal_windForward(ph,me)
integer, intent(in) :: ph, me
integer :: so
do so = 1, size(thermalState(ph)%p)
thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state(:,me)
enddo
end subroutine thermal_windForward
module subroutine thermal_forward()
integer :: ph, so
do ph = 1, size(thermalState)
do so = 1, size(thermalState(ph)%p)
thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state
enddo
enddo
end subroutine thermal_forward
module subroutine thermal_restore(ip,el)
integer, intent(in) :: ip, el
integer :: co, ph, me, so
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
do so = 1, size(thermalState(ph)%p)
thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me)
enddo
enddo
end subroutine thermal_restore
!----------------------------------------------------------------------------------------------
!< @brief Get temperature (for use by non-thermal physics)
!----------------------------------------------------------------------------------------------
module function thermal_T(ph,me) result(T)
integer, intent(in) :: ph, me
real(pReal) :: T
integer :: ho, tme
ho = material_homogenizationAt(el)
tme = material_homogenizationMemberAt(ip,el)
T = current(ph)%T(me)
T = temperature(ho)%p(tme)
end function thermal_T
end function constitutive_thermal_T
!----------------------------------------------------------------------------------------------
!< @brief Set temperature
!----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_setT(T,co,ip,el)
real(pReal), intent(in) :: T
integer, intent(in) :: co, ip, el
current(material_phaseAt(co,el))%T(material_phaseMemberAt(co,ip,el)) = T
end subroutine constitutive_thermal_setT
!--------------------------------------------------------------------------------------------------
!> @brief checks if a source mechanism is active or not
!--------------------------------------------------------------------------------------------------
function thermal_active(source_label,src_length) result(active_source)
character(len=*), intent(in) :: source_label !< name of source mechanism
integer, intent(in) :: src_length !< max. number of sources in system
logical, dimension(:,:), allocatable :: active_source
class(tNode), pointer :: &
phases, &
phase, &
sources, thermal, &
src
integer :: p,s
phases => config_material%get('phase')
allocate(active_source(src_length,phases%length), source = .false. )
do p = 1, phases%length
phase => phases%get(p)
if (phase%contains('thermal')) then
thermal => phase%get('thermal',defaultVal=emptyList)
sources => thermal%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
endif
enddo
end function thermal_active
end submodule constitutive_thermal

28
src/source_thermal_dissipation.f90 → src/constitutive_thermal_dissipation.f90
View File

@ -4,7 +4,7 @@
!> @brief material subroutine for thermal source due to plastic dissipation
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_thermal) source_thermal_dissipation
submodule(constitutive:constitutive_thermal) source_dissipation
integer, dimension(:), allocatable :: &
source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism?
@ -33,13 +33,14 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
class(tNode), pointer :: &
phases, &
phase, &
sources, &
sources, thermal, &
src
integer :: Ninstances,sourceOffset,Nconstituents,p
print'(/,a)', ' <<<+- source_thermal_dissipation init -+>>>'
print'(/,a)', ' <<<+- thermal_dissipation init -+>>>'
mySources = thermal_active('dissipation',source_length)
mySources = source_active('thermal_dissipation',source_length)
Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstances == 0) return
@ -51,18 +52,19 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
do p = 1, phases%length
phase => phases%get(p)
if(count(mySources(:,p)) == 0) cycle
if(any(mySources(:,p))) source_thermal_dissipation_instance(p) = count(mySources(:,1:p))
sources => phase%get('source')
if(count(mySources(:,p)) == 0) cycle
thermal => phase%get('thermal')
sources => thermal%get('source')
do sourceOffset = 1, sources%length
if(mySources(sourceOffset,p)) then
source_thermal_dissipation_offset(p) = sourceOffset
associate(prm => param(source_thermal_dissipation_instance(p)))
src => sources%get(sourceOffset)
prm%kappa = src%get_asFloat('kappa')
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,0,0,0)
call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,0,0,0)
end associate
endif
@ -76,7 +78,7 @@ end function source_thermal_dissipation_init
!--------------------------------------------------------------------------------------------------
!> @brief Ninstancess dissipation rate
!--------------------------------------------------------------------------------------------------
module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase)
module subroutine thermal_dissipation_getRate(TDot, Tstar, Lp, phase)
integer, intent(in) :: &
phase
@ -86,14 +88,12 @@ module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT
Lp
real(pReal), intent(out) :: &
TDot, &
dTDot_dT
TDot
associate(prm => param(source_thermal_dissipation_instance(phase)))
TDot = prm%kappa*sum(abs(Tstar*Lp))
dTDot_dT = 0.0_pReal
end associate
end subroutine source_thermal_dissipation_getRateAndItsTangent
end subroutine thermal_dissipation_getRate
end submodule source_thermal_dissipation
end submodule source_dissipation

29
src/source_thermal_externalheat.f90 → src/constitutive_thermal_externalheat.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Michigan State University
!> @brief material subroutine for variable heat source
!--------------------------------------------------------------------------------------------------
submodule(constitutive:constitutive_thermal) source_thermal_externalheat
submodule(constitutive:constitutive_thermal) source_externalheat
integer, dimension(:), allocatable :: &
@ -37,13 +37,14 @@ module function source_thermal_externalheat_init(source_length) result(mySources
class(tNode), pointer :: &
phases, &
phase, &
sources, &
sources, thermal, &
src
integer :: Ninstances,sourceOffset,Nconstituents,p
print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>'
print'(/,a)', ' <<<+- thermal_externalheat init -+>>>'
mySources = thermal_active('externalheat',source_length)
mySources = source_active('thermal_externalheat',source_length)
Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstances == 0) return
@ -57,7 +58,8 @@ module function source_thermal_externalheat_init(source_length) result(mySources
phase => phases%get(p)
if(any(mySources(:,p))) source_thermal_externalheat_instance(p) = count(mySources(:,1:p))
if(count(mySources(:,p)) == 0) cycle
sources => phase%get('source')
thermal => phase%get('thermal')
sources => thermal%get('source')
do sourceOffset = 1, sources%length
if(mySources(sourceOffset,p)) then
source_thermal_externalheat_offset(p) = sourceOffset
@ -70,9 +72,8 @@ module function source_thermal_externalheat_init(source_length) result(mySources
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,1,1,0)
end associate
endif
enddo
enddo
@ -95,7 +96,7 @@ module subroutine source_thermal_externalheat_dotState(phase, of)
sourceOffset = source_thermal_externalheat_offset(phase)
sourceState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time
thermalState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time
end subroutine source_thermal_externalheat_dotState
@ -103,14 +104,13 @@ end subroutine source_thermal_externalheat_dotState
!--------------------------------------------------------------------------------------------------
!> @brief returns local heat generation rate
!--------------------------------------------------------------------------------------------------
module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of)
module subroutine thermal_externalheat_getRate(TDot, phase, of)
integer, intent(in) :: &
phase, &
of
real(pReal), intent(out) :: &
TDot, &
dTDot_dT
TDot
integer :: &
sourceOffset, interval
@ -121,7 +121,7 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d
associate(prm => param(source_thermal_externalheat_instance(phase)))
do interval = 1, prm%nIntervals ! scan through all rate segments
frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) &
frac_time = (thermalState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) &
/ (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment
if ( (frac_time < 0.0_pReal .and. interval == 1) &
.or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
@ -130,9 +130,8 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d
prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries...
! ...or extrapolate if outside of bounds
enddo
dTDot_dT = 0.0
end associate
end subroutine source_thermal_externalheat_getRateAndItsTangent
end subroutine thermal_externalheat_getRate
end submodule source_thermal_externalheat
end submodule source_externalheat

8
src/damage_none.f90
View File

@ -25,10 +25,10 @@ subroutine damage_none_init
if (damage_type(h) /= DAMAGE_NONE_ID) cycle
Nmaterialpoints = count(material_homogenizationAt == h)
damageState(h)%sizeState = 0
allocate(damageState(h)%state0 (0,Nmaterialpoints))
allocate(damageState(h)%subState0(0,Nmaterialpoints))
allocate(damageState(h)%state (0,Nmaterialpoints))
damageState_h(h)%sizeState = 0
allocate(damageState_h(h)%state0 (0,Nmaterialpoints))
allocate(damageState_h(h)%subState0(0,Nmaterialpoints))
allocate(damageState_h(h)%state (0,Nmaterialpoints))
allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal)

10
src/damage_nonlocal.f90
View File

@ -76,12 +76,12 @@ subroutine damage_nonlocal_init
#endif
Nmaterialpoints = count(material_homogenizationAt == h)
damageState(h)%sizeState = 1
allocate(damageState(h)%state0 (1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal)
damageState_h(h)%sizeState = 1
allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState_h(h)%subState0(1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal)
damage(h)%p => damageState(h)%state(1,:)
damage(h)%p => damageState_h(h)%state(1,:)
end associate
enddo

4
src/grid/grid_thermal_spectral.f90
View File

@ -256,7 +256,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
PetscObject :: dummy
PetscErrorCode :: ierr
integer :: i, j, k, cell
real(pReal) :: Tdot, dTdot_dT
real(pReal) :: Tdot
T_current = x_scal
!--------------------------------------------------------------------------------------------------
@ -278,7 +278,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T_current(i,j,k), 1, cell)
call thermal_conduction_getSource(Tdot, T_current(i,j,k), 1, cell)
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
+ thermal_conduction_getMassDensity (1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)*(T_lastInc(i,j,k) - &

25
src/homogenization.f90
View File

@ -27,6 +27,8 @@ module homogenization
!--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point
real(pReal), dimension(:), allocatable, public :: &
homogenization_T
real(pReal), dimension(:,:,:), allocatable, public :: &
homogenization_F0, & !< def grad of IP at start of FE increment
homogenization_F !< def grad of IP to be reached at end of FE increment
@ -56,6 +58,9 @@ module homogenization
num_homog !< pointer to mechanical homogenization numerics data
end subroutine mech_init
module subroutine thermal_init
end subroutine thermal_init
module subroutine mech_partition(subF,ip,el)
real(pReal), intent(in), dimension(3,3) :: &
subF
@ -64,6 +69,13 @@ module homogenization
el !< element number
end subroutine mech_partition
module subroutine thermal_partition(T,ip,el)
real(pReal), intent(in) :: T
integer, intent(in) :: &
ip, & !< integration point
el !< element number
end subroutine thermal_partition
module subroutine mech_homogenize(dt,ip,el)
real(pReal), intent(in) :: dt
integer, intent(in) :: &
@ -126,9 +138,10 @@ subroutine homogenization_init
call mech_init(num_homog)
call thermal_init()
if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init
if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init
if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init(homogenization_T)
if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init(homogenization_T)
if (any(damage_type == DAMAGE_none_ID)) call damage_none_init
if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init
@ -173,7 +186,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
if (homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me)
if (damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me)
if (damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State0(:,me)
cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog)
@ -187,7 +200,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
call constitutive_windForward(ip,el)
if(homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me)
if(damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State(:,me)
endif steppingNeeded
elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
@ -202,7 +215,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
call constitutive_restore(ip,el,subStep < 1.0_pReal)
if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
if(damageState(ho)%sizeState > 0) damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%subState0(:,me)
endif
if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
@ -313,7 +326,7 @@ subroutine homogenization_forward
do ho = 1, size(material_name_homogenization)
homogState (ho)%state0 = homogState (ho)%state
damageState(ho)%state0 = damageState(ho)%state
damageState_h(ho)%state0 = damageState_h(ho)%state
enddo
end subroutine homogenization_forward

10
src/homogenization_mech.f90
View File

@ -113,24 +113,24 @@ module subroutine mech_partition(subF,ip,el)
el !< element number
integer :: co
real(pReal) :: F(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt(el))) :: Fs
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
F(1:3,1:3,1) = subF
Fs(1:3,1:3,1) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(F,subF)
call mech_isostrain_partitionDeformation(Fs,subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(F,subF,ip,el)
call mech_RGC_partitionDeformation(Fs,subF,ip,el)
end select chosenHomogenization
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
call constitutive_mech_setF(F(1:3,1:3,co),co,ip,el)
call constitutive_mech_setF(Fs(1:3,1:3,co),co,ip,el)
enddo

39
src/homogenization_thermal.f90 Normal file
View File

@ -0,0 +1,39 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_thermal
contains
!--------------------------------------------------------------------------------------------------
!> @brief Allocate variables and set parameters.
!--------------------------------------------------------------------------------------------------
module subroutine thermal_init()
print'(/,a)', ' <<<+- homogenization_thermal init -+>>>'
allocate(homogenization_T(discretization_nIPs*discretization_Nelems))
end subroutine thermal_init
!--------------------------------------------------------------------------------------------------
!> @brief Partition T onto the individual constituents.
!--------------------------------------------------------------------------------------------------
module subroutine thermal_partition(T,ip,el)
real(pReal), intent(in) :: T
integer, intent(in) :: &
ip, & !< integration point
el !< element number
call constitutive_thermal_setT(T,1,ip,el)
end subroutine thermal_partition
end submodule homogenization_thermal

83
src/lattice.f90
View File

@ -453,12 +453,13 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine lattice_init
integer :: Nphases, p,i
integer :: Nphases, ph,i
class(tNode), pointer :: &
phases, &
phase, &
mech, &
elasticity
elasticity, &
thermal
print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT)
@ -476,67 +477,71 @@ subroutine lattice_init
lattice_mu, lattice_nu,&
source=[(0.0_pReal,i=1,Nphases)])
do p = 1, phases%length
phase => phases%get(p)
do ph = 1, phases%length
phase => phases%get(ph)
mech => phase%get('mechanics')
elasticity => mech%get('elasticity')
lattice_C66(1,1,p) = elasticity%get_asFloat('C_11')
lattice_C66(1,2,p) = elasticity%get_asFloat('C_12')
lattice_C66(1,1,ph) = elasticity%get_asFloat('C_11')
lattice_C66(1,2,ph) = elasticity%get_asFloat('C_12')
lattice_C66(1,3,p) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal)
lattice_C66(2,2,p) = elasticity%get_asFloat('C_22',defaultVal=0.0_pReal)
lattice_C66(2,3,p) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal)
lattice_C66(3,3,p) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal)
lattice_C66(4,4,p) = elasticity%get_asFloat('C_44',defaultVal=0.0_pReal)
lattice_C66(5,5,p) = elasticity%get_asFloat('C_55',defaultVal=0.0_pReal)
lattice_C66(6,6,p) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal)
lattice_C66(1,3,ph) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal)
lattice_C66(2,2,ph) = elasticity%get_asFloat('C_22',defaultVal=0.0_pReal)
lattice_C66(2,3,ph) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal)
lattice_C66(3,3,ph) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal)
lattice_C66(4,4,ph) = elasticity%get_asFloat('C_44',defaultVal=0.0_pReal)
lattice_C66(5,5,ph) = elasticity%get_asFloat('C_55',defaultVal=0.0_pReal)
lattice_C66(6,6,ph) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal)
select case(phase%get_asString('lattice'))
case('cF')
lattice_structure(p) = lattice_FCC_ID
lattice_structure(ph) = lattice_FCC_ID
case('cI')
lattice_structure(p) = lattice_BCC_ID
lattice_structure(ph) = lattice_BCC_ID
case('hP')
lattice_structure(p) = lattice_HEX_ID
lattice_structure(ph) = lattice_HEX_ID
case('tI')
lattice_structure(p) = lattice_BCT_ID
lattice_structure(ph) = lattice_BCT_ID
case('oP')
lattice_structure(p) = lattice_ORT_ID
lattice_structure(ph) = lattice_ORT_ID
case('aP')
lattice_structure(p) = lattice_ISO_ID
lattice_structure(ph) = lattice_ISO_ID
case default
call IO_error(130,ext_msg='lattice_init: '//phase%get_asString('lattice'))
end select
lattice_C66(1:6,1:6,p) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,p),phase%get_asString('lattice'))
lattice_C66(1:6,1:6,ph) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,ph),phase%get_asString('lattice'))
lattice_nu(p) = lattice_equivalent_nu(lattice_C66(1:6,1:6,p),'voigt')
lattice_mu(p) = lattice_equivalent_mu(lattice_C66(1:6,1:6,p),'voigt')
lattice_nu(ph) = lattice_equivalent_nu(lattice_C66(1:6,1:6,ph),'voigt')
lattice_mu(ph) = lattice_equivalent_mu(lattice_C66(1:6,1:6,ph),'voigt')
lattice_C66(1:6,1:6,p) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,p))) ! Literature data is in Voigt notation
lattice_C66(1:6,1:6,ph) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,ph))) ! Literature data is in Voigt notation
do i = 1, 6
if (abs(lattice_C66(i,i,p))<tol_math_check) &
call IO_error(135,el=i,ip=p,ext_msg='matrix diagonal "el"ement of phase "ip"')
if (abs(lattice_C66(i,i,ph))<tol_math_check) &
call IO_error(135,el=i,ip=ph,ext_msg='matrix diagonal "el"ement of phase "ip"')
enddo
lattice_rho(ph) = phase%get_asFloat('rho', defaultVal=0.0_pReal)
! SHOULD NOT BE PART OF LATTICE BEGIN
lattice_K(1,1,p) = phase%get_asFloat('K_11',defaultVal=0.0_pReal)
lattice_K(2,2,p) = phase%get_asFloat('K_22',defaultVal=0.0_pReal)
lattice_K(3,3,p) = phase%get_asFloat('K_33',defaultVal=0.0_pReal)
lattice_K(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,p), &
if (phase%contains('thermal')) then
thermal => phase%get('thermal')
lattice_K(1,1,ph) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal)
lattice_K(2,2,ph) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal)
lattice_K(3,3,ph) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal)
lattice_K(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,ph), &
phase%get_asString('lattice'))
lattice_c_p(ph) = thermal%get_asFloat('c_p', defaultVal=0.0_pReal)
endif
lattice_D(1,1,ph) = phase%get_asFloat('D_11',defaultVal=0.0_pReal)
lattice_D(2,2,ph) = phase%get_asFloat('D_22',defaultVal=0.0_pReal)
lattice_D(3,3,ph) = phase%get_asFloat('D_33',defaultVal=0.0_pReal)
lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), &
phase%get_asString('lattice'))
lattice_c_p(p) = phase%get_asFloat('c_p', defaultVal=0.0_pReal)
lattice_rho(p) = phase%get_asFloat('rho', defaultVal=0.0_pReal)
lattice_D(1,1,p) = phase%get_asFloat('D_11',defaultVal=0.0_pReal)
lattice_D(2,2,p) = phase%get_asFloat('D_22',defaultVal=0.0_pReal)
lattice_D(3,3,p) = phase%get_asFloat('D_33',defaultVal=0.0_pReal)
lattice_D(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,p), &
phase%get_asString('lattice'))
lattice_M(p) = phase%get_asFloat('M',defaultVal=0.0_pReal)
lattice_M(ph) = phase%get_asFloat('M',defaultVal=0.0_pReal)
! SHOULD NOT BE PART OF LATTICE END
call selfTest

4
src/material.f90
View File

@ -61,7 +61,7 @@ module material
type(tState), allocatable, dimension(:), public :: &
homogState, &
damageState
damageState_h
type(Rotation), dimension(:,:,:), allocatable, public, protected :: &
material_orientation0 !< initial orientation of each grain,IP,element
@ -101,7 +101,7 @@ subroutine material_init(restart)
allocate(homogState (size(material_name_homogenization)))
allocate(damageState (size(material_name_homogenization)))
allocate(damageState_h (size(material_name_homogenization)))
allocate(temperature (size(material_name_homogenization)))
allocate(damage (size(material_name_homogenization)))

16
src/source_damage_anisoBrittle.f90
View File

@ -101,9 +101,9 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
end associate
@ -146,7 +146,7 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_anisoBrittle_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
damageState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
do i = 1, prm%sum_N_cl
traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
@ -154,8 +154,8 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
damageState(phase)%p(sourceOffset)%dotState(1,constituent) &
= damageState(phase)%p(sourceOffset)%dotState(1,constituent) &
+ prm%dot_o / prm%s_crit(i) &
* ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
(max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
@ -185,7 +185,7 @@ module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, d
sourceOffset = source_damage_anisoBrittle_offset(phase)
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
@ -204,7 +204,7 @@ module subroutine source_damage_anisoBrittle_results(phase,group)
integer :: o
associate(prm => param(source_damage_anisoBrittle_instance(phase)), &
stt => sourceState(phase)%p(source_damage_anisoBrittle_offset(phase))%state)
stt => damageState(phase)%p(source_damage_anisoBrittle_offset(phase))%state)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('f_phi')

12
src/source_damage_anisoDuctile.f90
View File

@ -87,9 +87,9 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
end associate
@ -128,7 +128,7 @@ module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_anisoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &
damageState(phase)%p(sourceOffset)%dotState(1,constituent) &
= sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit)
end associate
@ -154,7 +154,7 @@ module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, d
sourceOffset = source_damage_anisoDuctile_offset(phase)
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
@ -173,7 +173,7 @@ module subroutine source_damage_anisoDuctile_results(phase,group)
integer :: o
associate(prm => param(source_damage_anisoDuctile_instance(phase)), &
stt => sourceState(phase)%p(source_damage_anisoDuctile_offset(phase))%state)
stt => damageState(phase)%p(source_damage_anisoDuctile_offset(phase))%state)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('f_phi')

24
src/source_damage_isoBrittle.f90
View File

@ -74,9 +74,9 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,1)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
end associate
@ -124,13 +124,13 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el)
strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit
! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit
if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then
sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
strainenergy - sourceState(phase)%p(sourceOffset)%state(1,constituent)
if (strainenergy > damageState(phase)%p(sourceOffset)%subState0(1,constituent)) then
damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
strainenergy - damageState(phase)%p(sourceOffset)%state(1,constituent)
else
sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - &
sourceState(phase)%p(sourceOffset)%state(1,constituent)
damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
damageState(phase)%p(sourceOffset)%subState0(1,constituent) - &
damageState(phase)%p(sourceOffset)%state(1,constituent)
endif
end associate
@ -158,8 +158,8 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo
associate(prm => param(source_damage_isoBrittle_instance(phase)))
localphiDot = 1.0_pReal &
- phi*sourceState(phase)%p(sourceOffset)%state(1,constituent)
dLocalphiDot_dPhi = - sourceState(phase)%p(sourceOffset)%state(1,constituent)
- phi*damageState(phase)%p(sourceOffset)%state(1,constituent)
dLocalphiDot_dPhi = - damageState(phase)%p(sourceOffset)%state(1,constituent)
end associate
end subroutine source_damage_isoBrittle_getRateAndItsTangent
@ -176,7 +176,7 @@ module subroutine source_damage_isoBrittle_results(phase,group)
integer :: o
associate(prm => param(source_damage_isoBrittle_instance(phase)), &
stt => sourceState(phase)%p(source_damage_isoBrittle_offset(phase))%state)
stt => damageState(phase)%p(source_damage_isoBrittle_offset(phase))%state)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('f_phi')

12
src/source_damage_isoDuctile.f90
View File

@ -78,9 +78,9 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
end associate
@ -119,7 +119,7 @@ module subroutine source_damage_isoDuctile_dotState(co, ip, el)
damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_isoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
damageState(phase)%p(sourceOffset)%dotState(1,constituent) = &
sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit
end associate
@ -145,7 +145,7 @@ module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLo
sourceOffset = source_damage_isoDuctile_offset(phase)
dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent)
localphiDot = 1.0_pReal &
+ dLocalphiDot_dPhi*phi
@ -164,7 +164,7 @@ module subroutine source_damage_isoDuctile_results(phase,group)
integer :: o
associate(prm => param(source_damage_isoDuctile_instance(phase)), &
stt => sourceState(phase)%p(source_damage_isoDuctile_offset(phase))%state)
stt => damageState(phase)%p(source_damage_isoDuctile_offset(phase))%state)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('f_phi')

46
src/thermal_conduction.f90
View File

@ -10,6 +10,7 @@ module thermal_conduction
use results
use constitutive
use YAML_types
use discretization
implicit none
private
@ -24,7 +25,7 @@ module thermal_conduction
public :: &
thermal_conduction_init, &
thermal_conduction_getSourceAndItsTangent, &
thermal_conduction_getSource, &
thermal_conduction_getConductivity, &
thermal_conduction_getSpecificHeat, &
thermal_conduction_getMassDensity, &
@ -38,25 +39,28 @@ contains
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_init
subroutine thermal_conduction_init(T)
integer :: Ninstances,Nmaterialpoints,h
real(pReal), dimension(:), intent(inout) :: T
integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce
class(tNode), pointer :: &
material_homogenization, &
homog, &
homogThermal
print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6)
Ninstances = count(thermal_type == THERMAL_conduction_ID)
allocate(param(Ninstances))
material_homogenization => config_material%get('homogenization')
do h = 1, size(material_name_homogenization)
if (thermal_type(h) /= THERMAL_conduction_ID) cycle
homog => material_homogenization%get(h)
do ho = 1, size(material_name_homogenization)
if (thermal_type(ho) /= THERMAL_conduction_ID) cycle
homog => material_homogenization%get(ho)
homogThermal => homog%get('thermal')
associate(prm => param(thermal_typeInstance(h)))
associate(prm => param(thermal_typeInstance(ho)))
#if defined (__GFORTRAN__)
prm%output = output_asStrings(homogThermal)
@ -64,21 +68,30 @@ subroutine thermal_conduction_init
prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray)
#endif
Nmaterialpoints=count(material_homogenizationAt==h)
Nmaterialpoints=count(material_homogenizationAt==ho)
allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h))
allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal)
allocate (temperature (ho)%p(Nmaterialpoints), source=thermal_initialT(ho))
allocate (temperatureRate(ho)%p(Nmaterialpoints), source=0.0_pReal)
end associate
enddo
ce = 0
do el = 1, discretization_Nelems
do ip = 1, discretization_nIPs
ce = ce + 1
ho = material_homogenizationAt(el)
if (thermal_type(ho) == THERMAL_conduction_ID) T(ce) = thermal_initialT(ho)
enddo
enddo
end subroutine thermal_conduction_init
!--------------------------------------------------------------------------------------------------
!> @brief return heat generation rate
!--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
subroutine thermal_conduction_getSource(Tdot, T,ip,el)
integer, intent(in) :: &
ip, & !< integration point number
@ -86,20 +99,17 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
real(pReal), intent(in) :: &
T
real(pReal), intent(out) :: &
Tdot, dTdot_dT
Tdot
integer :: &
homog
Tdot = 0.0_pReal
dTdot_dT = 0.0_pReal
homog = material_homogenizationAt(el)
call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el)
call constitutive_thermal_getRate(TDot, T,ip,el)
Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal)
end subroutine thermal_conduction_getSourceAndItsTangent
end subroutine thermal_conduction_getSource
!--------------------------------------------------------------------------------------------------

26
src/thermal_isothermal.f90
View File

@ -6,6 +6,7 @@ module thermal_isothermal
use prec
use config
use material
use discretization
implicit none
public
@ -15,22 +16,33 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
subroutine thermal_isothermal_init
subroutine thermal_isothermal_init(T)
integer :: h,Nmaterialpoints
real(pReal), dimension(:), intent(inout) :: T
integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce
print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6)
do h = 1, size(material_name_homogenization)
if (thermal_type(h) /= THERMAL_isothermal_ID) cycle
do ho = 1, size(thermal_type)
if (thermal_type(ho) /= THERMAL_isothermal_ID) cycle
Nmaterialpoints = count(material_homogenizationAt == h)
Nmaterialpoints = count(material_homogenizationAt == ho)
allocate(temperature (h)%p(Nmaterialpoints),source=thermal_initialT(h))
allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal)
allocate(temperature (ho)%p(Nmaterialpoints),source=thermal_initialT(ho))
allocate(temperatureRate(ho)%p(Nmaterialpoints),source = 0.0_pReal)
enddo
ce = 0
do el = 1, discretization_Nelems
do ip = 1, discretization_nIPs
ce = ce + 1
ho = material_homogenizationAt(el)
if (thermal_type(ho) == THERMAL_isothermal_ID) T(ce) = thermal_initialT(ho)
enddo
enddo
end subroutine thermal_isothermal_init
end module thermal_isothermal