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

removed temperature integration and corresponding data structures and debugging options 

temperature is stored in crystallite, but homogeneous on one IP (not an component (grain) quantity and an input value parsed in by the BVP solver.
introduced heat, a component (grain) quantity which is homogenized before returned to the heat transfer solver.
went ahead with removal of dummy functions in homogenization and constitutive, this time mainly reduced function signatures to reflect actually needed quantities.
This commit is contained in:
Martin Diehl 2013-10-16 13:04:59 +00:00
parent 6a1c40d540
commit dc95c82d4a
18 changed files with 757 additions and 1005 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

205
code/CPFEM.f90
View File

@ -1,4 +1,4 @@
! Copyright 2011-13 Max-Planck-Institut für Eisenforschung GmbH !nd Copyright 2011-13 Max-Planck-Institut für Eisenforschung GmbH
! !
! This file is part of DAMASK, ! This file is part of DAMASK,
! the Düsseldorf Advanced Material Simulation Kit. ! the Düsseldorf Advanced Material Simulation Kit.
@ -37,8 +37,8 @@ module CPFEM
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE_knownGood !< known good tangent real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_dcsdE_knownGood !< known good tangent
logical :: CPFEM_init_done = .false., & !< remember whether init has been done already logical :: CPFEM_init_done = .false., & !< remember whether init has been done already
CPFEM_init_inProgress = .false., & !< remember whether first IP is currently performing init CPFEM_init_inProgress = .false., & !< remember whether first ip is currently performing init
CPFEM_calc_done = .false. !< remember whether first IP has already calced the results CPFEM_calc_done = .false. !< remember whether first ip has already calced the results
integer(pInt), parameter, public :: & integer(pInt), parameter, public :: &
CPFEM_CALCRESULTS = 2_pInt**0_pInt, & CPFEM_CALCRESULTS = 2_pInt**0_pInt, &
CPFEM_AGERESULTS = 2_pInt**1_pInt, & CPFEM_AGERESULTS = 2_pInt**1_pInt, &
@ -55,7 +55,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief call (thread safe) all module initializations !> @brief call (thread safe) all module initializations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll(Temperature,element,IP) subroutine CPFEM_initAll(temperature,el,ip)
use prec, only: & use prec, only: &
prec_init prec_init
use numerics, only: & use numerics, only: &
@ -83,9 +83,9 @@ subroutine CPFEM_initAll(Temperature,element,IP)
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: element, & ! FE element number integer(pInt), intent(in) :: el, & ! FE el number
IP ! FE integration point number ip ! FE integration point number
real(pReal), intent(in) :: Temperature ! temperature real(pReal), intent(in) :: temperature ! temperature
real(pReal) rnd real(pReal) rnd
integer(pInt) i,n integer(pInt) i,n
@ -108,12 +108,12 @@ subroutine CPFEM_initAll(Temperature,element,IP)
call debug_init call debug_init
call math_init call math_init
call FE_init call FE_init
call mesh_init(IP, element) ! pass on coordinates to alter calcMode of first ip call mesh_init(ip, el) ! pass on coordinates to alter calcMode of first ip
call lattice_init call lattice_init
call material_init call material_init
call constitutive_init call constitutive_init
call crystallite_init(Temperature) ! (have to) use temperature of first IP for whole model call crystallite_init(temperature) ! (have to) use temperature of first ip for whole model
call homogenization_init(Temperature) call homogenization_init
call CPFEM_init call CPFEM_init
#ifndef Spectral #ifndef Spectral
call DAMASK_interface_init() ! Spectral solver init is already done call DAMASK_interface_init() ! Spectral solver init is already done
@ -172,7 +172,7 @@ subroutine CPFEM_init
homogenization_state0 homogenization_state0
implicit none implicit none
integer(pInt) i,j,k,l,m integer(pInt) :: i,j,k,l,m
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>' write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
@ -245,8 +245,7 @@ subroutine CPFEM_init
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) then if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) then
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs) write(6,'(a32,1x,6(i8,1x))') 'CPFEM_cs: ', shape(CPFEM_cs)
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE) write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE)
write(6,'(a32,1x,6(i8,1x))') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood) write(6,'(a32,1x,6(i8,1x),/)') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
write(6,*)
write(6,*) 'symmetricSolver: ', symmetricSolver write(6,*) 'symmetricSolver: ', symmetricSolver
endif endif
flush(6) flush(6)
@ -257,11 +256,12 @@ end subroutine CPFEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief perform initialization at first call, update variables and call the actual material model !> @brief perform initialization at first call, update variables and call the actual material model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, element, IP, cauchyStress, jacobian) subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, elFE, ip, cauchyStress, jacobian)
! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE use numerics, only: &
use numerics, only: defgradTolerance, & defgradTolerance, &
iJacoStiffness iJacoStiffness
use debug, only: debug_level, & use debug, only: &
debug_level, &
debug_CPFEM, & debug_CPFEM, &
debug_levelBasic, & debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
@ -276,7 +276,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
debug_stressMin, & debug_stressMin, &
debug_jacobianMax, & debug_jacobianMax, &
debug_jacobianMin debug_jacobianMin
use FEsolving, only: outdatedFFN1, & use FEsolving, only: &
outdatedFFN1, &
terminallyIll, & terminallyIll, &
cycleCounter, & cycleCounter, &
theInc, & theInc, &
@ -285,7 +286,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
FEsolving_execElem, & FEsolving_execElem, &
FEsolving_execIP, & FEsolving_execIP, &
restartWrite restartWrite
use math, only: math_identity2nd, & use math, only: &
math_identity2nd, &
math_mul33x33, & math_mul33x33, &
math_det33, & math_det33, &
math_transpose33, & math_transpose33, &
@ -294,15 +296,19 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
math_Mandel66to3333, & math_Mandel66to3333, &
math_Mandel33to6, & math_Mandel33to6, &
math_Mandel6to33 math_Mandel6to33
use mesh, only: mesh_FEasCP, & use mesh, only: &
mesh_FEasCP, &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips, & mesh_maxNips, &
mesh_element mesh_element
use material, only: homogenization_maxNgrains, & use material, only: &
homogenization_maxNgrains, &
microstructure_elemhomo, & microstructure_elemhomo, &
material_phase material_phase
use constitutive, only: constitutive_state0,constitutive_state use constitutive, only: &
use crystallite, only: crystallite_partionedF,& constitutive_state0,constitutive_state
use crystallite, only: &
crystallite_partionedF,&
crystallite_F0, & crystallite_F0, &
crystallite_Fp0, & crystallite_Fp0, &
crystallite_Fp, & crystallite_Fp, &
@ -311,8 +317,10 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
crystallite_dPdF0, & crystallite_dPdF0, &
crystallite_dPdF, & crystallite_dPdF, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_Tstar_v crystallite_Tstar_v, &
use homogenization, only: homogenization_sizeState, & crystallite_temperature
use homogenization, only: &
homogenization_sizeState, &
homogenization_state, & homogenization_state, &
homogenization_state0, & homogenization_state0, &
materialpoint_F, & materialpoint_F, &
@ -321,17 +329,17 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
materialpoint_dPdF, & materialpoint_dPdF, &
materialpoint_results, & materialpoint_results, &
materialpoint_sizeResults, & materialpoint_sizeResults, &
materialpoint_Temperature, &
materialpoint_stressAndItsTangent, & materialpoint_stressAndItsTangent, &
materialpoint_postResults materialpoint_postResults
use IO, only: IO_write_jobRealFile, & use IO, only: &
IO_write_jobRealFile, &
IO_warning IO_warning
use DAMASK_interface use DAMASK_interface
implicit none implicit none
integer(pInt), intent(in) :: element, & !< FE element number integer(pInt), intent(in) :: elFE, & !< FE element number
IP !< FE integration point number ip !< integration point number
real(pReal), intent(inout) :: Temperature !< temperature real(pReal), intent(in) :: temperature !< temperature
real(pReal), intent(in) :: dt !< time increment real(pReal), intent(in) :: dt !< time increment
real(pReal), dimension (3,3), intent(in) :: ffn, & !< deformation gradient for t=t0 real(pReal), dimension (3,3), intent(in) :: ffn, & !< deformation gradient for t=t0
ffn1 !< deformation gradient for t=t1 ffn1 !< deformation gradient for t=t1
@ -339,7 +347,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs
real(pReal), dimension(6), intent(out), optional :: cauchyStress !< stress vector in Mandel notation real(pReal), dimension(6), intent(out), optional :: cauchyStress !< stress vector in Mandel notation
real(pReal), dimension(6,6), intent(out), optional :: jacobian !< jacobian in Mandel notation real(pReal), dimension(6,6), intent(out), optional :: jacobian !< jacobian in Mandel notation (Consistent tangent dcs/dE)
real(pReal) J_inverse, & ! inverse of Jacobian real(pReal) J_inverse, & ! inverse of Jacobian
rnd rnd
@ -348,19 +356,20 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
real(pReal), dimension (3,3,3,3) :: H_sym, & real(pReal), dimension (3,3,3,3) :: H_sym, &
H, & H, &
jacobian3333 ! jacobian in Matrix notation jacobian3333 ! jacobian in Matrix notation
integer(pInt) cp_en, & ! crystal plasticity element number integer(pInt) elCP, & ! crystal plasticity element number
i, j, k, l, m, n i, j, k, l, m, n
logical updateJaco ! flag indicating if JAcobian has to be updated logical updateJaco ! flag indicating if JAcobian has to be updated
cp_en = mesh_FEasCP('elem',element) elCP = mesh_FEasCP('elem',elFE)
!if elCP = 0_pInt return ToDo: needed?
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt & if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt &
.and. cp_en == debug_e .and. IP == debug_i) then .and. elCP == debug_e .and. ip == debug_i) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(/,a)') '#############################################' write(6,'(/,a)') '#############################################'
write(6,'(a1,a22,1x,i8,a13)') '#','element', cp_en, '#' write(6,'(a1,a22,1x,i8,a13)') '#','element', elCP, '#'
write(6,'(a1,a22,1x,i8,a13)') '#','IP', IP, '#' write(6,'(a1,a22,1x,i8,a13)') '#','ip', ip, '#'
write(6,'(a1,a22,1x,f15.7,a6)') '#','theTime', theTime, '#' write(6,'(a1,a22,1x,f15.7,a6)') '#','theTime', theTime, '#'
write(6,'(a1,a22,1x,f15.7,a6)') '#','theDelta', theDelta, '#' write(6,'(a1,a22,1x,f15.7,a6)') '#','theDelta', theDelta, '#'
write(6,'(a1,a22,1x,i8,a13)') '#','theInc', theInc, '#' write(6,'(a1,a22,1x,i8,a13)') '#','theInc', theInc, '#'
@ -398,9 +407,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a)') '<< CPFEM >> aging states' write(6,'(a)') '<< CPFEM >> aging states'
if (debug_e <= mesh_NcpElems .and. debug_i <= mesh_maxNips) then if (debug_e <= mesh_NcpElems .and. debug_i <= mesh_maxNips) then
write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)))') '<< CPFEM >> aged state of el ip grain',& write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)),/)') '<< CPFEM >> aged state of elFE ip grain',&
debug_e, debug_i, 1, constitutive_state(1,debug_i,debug_e)%p debug_e, debug_i, 1, constitutive_state(1,debug_i,debug_e)%p
write(6,*)
endif endif
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -479,18 +487,18 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
!* In case that no parallel execution is required, there is no need to collect FEM input !* In case that no parallel execution is required, there is no need to collect FEM input
if (.not. parallelExecution) then if (.not. parallelExecution) then
materialpoint_Temperature(IP,cp_en) = Temperature crystallite_temperature(ip,elCP) = temperature
materialpoint_F0(1:3,1:3,IP,cp_en) = ffn materialpoint_F0(1:3,1:3,ip,elCP) = ffn
materialpoint_F(1:3,1:3,IP,cp_en) = ffn1 materialpoint_F(1:3,1:3,ip,elCP) = ffn1
elseif (iand(mode, CPFEM_COLLECT) /= 0_pInt) then elseif (iand(mode, CPFEM_COLLECT) /= 0_pInt) then
call random_number(rnd) call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
materialpoint_Temperature(IP,cp_en) = Temperature crystallite_temperature(ip,elCP) = temperature
materialpoint_F0(1:3,1:3,IP,cp_en) = ffn materialpoint_F0(1:3,1:3,ip,elCP) = ffn
materialpoint_F(1:3,1:3,IP,cp_en) = ffn1 materialpoint_F(1:3,1:3,ip,elCP) = ffn1
CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress
CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6)
CPFEM_calc_done = .false. CPFEM_calc_done = .false.
endif endif
@ -502,13 +510,13 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
!*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one !*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one
if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then
! if (.not. terminallyIll .and. .not. outdatedFFN1) then ! if (.not. terminallyIll .and. .not. outdatedFFN1) then
if (any(abs(ffn1 - materialpoint_F(1:3,1:3,IP,cp_en)) > defgradTolerance)) then if (any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at el ip',cp_en,IP write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elCP,ip
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en)) write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',math_transpose33(materialpoint_F(1:3,1:3,ip,elCP))
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1) write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -516,8 +524,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
endif endif
call random_number(rnd) call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,IP,cp_en) = rnd*CPFEM_odd_stress CPFEM_cs(1:6,ip,elCP) = rnd*CPFEM_odd_stress
CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6) CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian*math_identity2nd(6)
!*** deformation gradient is not outdated !*** deformation gradient is not outdated
@ -525,18 +533,18 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
else else
updateJaco = mod(cycleCounter,iJacoStiffness) == 0 updateJaco = mod(cycleCounter,iJacoStiffness) == 0
!* no parallel computation, so we use just one single element and IP for computation !* no parallel computation, so we use just one single elFE and ip for computation
if (.not. parallelExecution) then if (.not. parallelExecution) then
FEsolving_execElem(1) = cp_en FEsolving_execElem(1) = elCP
FEsolving_execElem(2) = cp_en FEsolving_execElem(2) = elCP
if (.not. microstructure_elemhomo(mesh_element(4,cp_en)) .or. & ! calculate unless homogeneous if (.not. microstructure_elemhomo(mesh_element(4,elCP)) .or. & ! calculate unless homogeneous
(microstructure_elemhomo(mesh_element(4,cp_en)) .and. IP == 1_pInt)) then ! and then only first IP (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip == 1_pInt)) then ! and then only first ip
FEsolving_execIP(1,cp_en) = IP FEsolving_execIP(1,elCP) = ip
FEsolving_execIP(2,cp_en) = IP FEsolving_execIP(2,elCP) = ip
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for el ip ',cp_en,IP write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elCP,ip
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent
@ -561,98 +569,81 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, Temperature, dt, el
if ( terminallyIll ) then if ( terminallyIll ) then
call random_number(rnd) call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,IP,cp_en) = rnd * CPFEM_odd_stress CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress
CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6) CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6)
else else
if (microstructure_elemhomo(mesh_element(4,cp_en)) .and. IP > 1_pInt) then ! me homogenous? --> copy from first IP if (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip > 1_pInt) then ! me homogenous? --> copy from first ip
materialpoint_P(1:3,1:3,IP,cp_en) = materialpoint_P(1:3,1:3,1,cp_en) materialpoint_P(1:3,1:3,ip,elCP) = materialpoint_P(1:3,1:3,1,elCP)
materialpoint_F(1:3,1:3,IP,cp_en) = materialpoint_F(1:3,1:3,1,cp_en) materialpoint_F(1:3,1:3,ip,elCP) = materialpoint_F(1:3,1:3,1,elCP)
materialpoint_dPdF(1:3,1:3,1:3,1:3,IP,cp_en) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,cp_en) materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,elCP) = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,elCP)
materialpoint_results(1:materialpoint_sizeResults,IP,cp_en) = materialpoint_results(1:materialpoint_sizeResults,1,cp_en) materialpoint_results(1:materialpoint_sizeResults,ip,elCP) = materialpoint_results(1:materialpoint_sizeResults,1,elCP)
endif endif
! translate from P to CS ! translate from P to CS
Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,IP,cp_en), math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en))) Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,ip,elCP), math_transpose33(materialpoint_F(1:3,1:3,ip,elCP)))
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,IP,cp_en)) J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
CPFEM_cs(1:6,IP,cp_en) = math_Mandel33to6(J_inverse * Kirchhoff) CPFEM_cs(1:6,ip,elCP) = math_Mandel33to6(J_inverse * Kirchhoff)
! translate from dP/dF to dCS/dE ! translate from dP/dF to dCS/dE
H = 0.0_pReal H = 0.0_pReal
do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
H(i,j,k,l) = H(i,j,k,l) + & H(i,j,k,l) = H(i,j,k,l) + &
materialpoint_F(j,m,IP,cp_en) * & materialpoint_F(j,m,ip,elCP) * &
materialpoint_F(l,n,IP,cp_en) * & materialpoint_F(l,n,ip,elCP) * &
materialpoint_dPdF(i,m,k,n,IP,cp_en) - & materialpoint_dPdF(i,m,k,n,ip,elCP) - &
math_I3(j,l) * materialpoint_F(i,m,IP,cp_en) * materialpoint_P(k,m,IP,cp_en) + & math_I3(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) + &
0.5_pReal * (math_I3(i,k) * Kirchhoff(j,l) + math_I3(j,l) * Kirchhoff(i,k) + & 0.5_pReal * (math_I3(i,k) * Kirchhoff(j,l) + math_I3(j,l) * Kirchhoff(i,k) + &
math_I3(i,l) * Kirchhoff(j,k) + math_I3(j,k) * Kirchhoff(i,l)) math_I3(i,l) * Kirchhoff(j,k) + math_I3(j,k) * Kirchhoff(i,l))
enddo; enddo; enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo; enddo; enddo
do i=1,3; do j=1,3; do k=1,3; do l=1,3 !< @ToDo use forall forall(i=1:3, j=1:3,k=1:3,l=1:3) &
H_sym(i,j,k,l) = 0.25_pReal * (H(i,j,k,l) + H(j,i,k,l) + H(i,j,l,k) + H(j,i,l,k)) H_sym(i,j,k,l) = 0.25_pReal * (H(i,j,k,l) + H(j,i,k,l) + H(i,j,l,k) + H(j,i,l,k))
enddo; enddo; enddo; enddo CPFEM_dcsde(1:6,1:6,ip,elCP) = math_Mandel3333to66(J_inverse * H_sym)
CPFEM_dcsde(1:6,1:6,IP,cp_en) = math_Mandel3333to66(J_inverse * H_sym)
endif endif
endif endif
!* remember extreme values of stress and jacobian !* remember extreme values of stress and jacobian
cauchyStress33 = math_Mandel6to33(CPFEM_cs(1:6,ip,elCP))
cauchyStress33 = math_Mandel6to33(CPFEM_cs(1:6,IP,cp_en))
if (maxval(cauchyStress33) > debug_stressMax) then if (maxval(cauchyStress33) > debug_stressMax) then
debug_stressMaxLocation = [cp_en, IP] debug_stressMaxLocation = [elCP, ip]
debug_stressMax = maxval(cauchyStress33) debug_stressMax = maxval(cauchyStress33)
endif endif
if (minval(cauchyStress33) < debug_stressMin) then if (minval(cauchyStress33) < debug_stressMin) then
debug_stressMinLocation = [cp_en, IP] debug_stressMinLocation = [elCP, ip]
debug_stressMin = minval(cauchyStress33) debug_stressMin = minval(cauchyStress33)
endif endif
jacobian3333 = math_Mandel66to3333(CPFEM_dcsdE(1:6,1:6,IP,cp_en)) jacobian3333 = math_Mandel66to3333(CPFEM_dcsdE(1:6,1:6,ip,elCP))
if (maxval(jacobian3333) > debug_jacobianMax) then if (maxval(jacobian3333) > debug_jacobianMax) then
debug_jacobianMaxLocation = [cp_en, IP] debug_jacobianMaxLocation = [elCP, ip]
debug_jacobianMax = maxval(jacobian3333) debug_jacobianMax = maxval(jacobian3333)
endif endif
if (minval(jacobian3333) < debug_jacobianMin) then if (minval(jacobian3333) < debug_jacobianMin) then
debug_jacobianMinLocation = [cp_en, IP] debug_jacobianMinLocation = [elCP, ip]
debug_jacobianMin = minval(jacobian3333) debug_jacobianMin = minval(jacobian3333)
endif endif
!* report stress and stiffness !* report stress and stiffness
if ((iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & if ((iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
.and. ((debug_e == cp_en .and. debug_i == IP) & .and. ((debug_e == elCP .and. debug_i == ip) &
.or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)') '<< CPFEM >> stress/MPa at el ip ', & write(6,'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)') '<< CPFEM >> stress/MPa at elFE ip ', &
cp_en, IP, CPFEM_cs(1:6,IP,cp_en)/1.0e6_pReal elCP, ip, CPFEM_cs(1:6,ip,elCP)/1.0e6_pReal
write(6,'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))') '<< CPFEM >> Jacobian/GPa at el ip ', & write(6,'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))') '<< CPFEM >> Jacobian/GPa at elFE ip ', &
cp_en, IP, transpose(CPFEM_dcsdE(1:6,1:6,IP,cp_en))/1.0e9_pReal elCP, ip, transpose(CPFEM_dcsdE(1:6,1:6,ip,elCP))/1.0e9_pReal
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
endif endif
!*** homogenized result except for potentially non-isothermal starting condition
if (theTime > 0.0_pReal) then
Temperature = materialpoint_Temperature(IP,cp_en)
endif
!*** warn if stiffness close to zero !*** warn if stiffness close to zero
if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip)
if (all(abs(CPFEM_dcsdE(1:6,1:6,IP,cp_en)) < 1e-10_pReal)) then
call IO_warning(601,cp_en,IP)
endif
!*** copy to output if required (FEM solver) !*** copy to output if required (FEM solver)
if(present(cauchyStress)) cauchyStress = CPFEM_cs(1:6,ip,elCP)
if(present(cauchyStress)) cauchyStress = CPFEM_cs(1:6,IP,cp_en) if(present(jacobian)) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP)
if(present(jacobian)) jacobian = CPFEM_dcsdE(1:6,1:6,IP,cp_en)
end subroutine CPFEM_general end subroutine CPFEM_general

2
code/DAMASK_spectral_driver.f90
View File

@ -151,7 +151,7 @@ program DAMASK_spectral_Driver
external :: quit external :: quit
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
call CPFEM_initAll(temperature = 300.0_pReal, element = 1_pInt, IP= 1_pInt) call CPFEM_initAll(temperature = 300.0_pReal, el = 1_pInt, ip = 1_pInt)
write(6,'(/,a)') ' <<<+- DAMASK_spectral_driver init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_driver init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()

5
code/DAMASK_spectral_utilities.f90
View File

@ -819,10 +819,11 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
CPFEM_AGERESULTS, & CPFEM_AGERESULTS, &
CPFEM_BACKUPJACOBIAN, & CPFEM_BACKUPJACOBIAN, &
CPFEM_RESTOREJACOBIAN CPFEM_RESTOREJACOBIAN
use crystallite, only: &
crystallite_temperature
use homogenization, only: & use homogenization, only: &
materialpoint_F0, & materialpoint_F0, &
materialpoint_F, & materialpoint_F, &
materialpoint_Temperature, &
materialpoint_P, & materialpoint_P, &
materialpoint_dPdF materialpoint_dPdF
@ -864,7 +865,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)]) materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)])
materialpoint_F = reshape(F, [3,3,1,product(grid)]) materialpoint_F = reshape(F, [3,3,1,product(grid)])
materialpoint_Temperature = temperature crystallite_temperature = temperature
call debug_reset() call debug_reset()

105
code/constitutive.f90
View File

@ -66,7 +66,7 @@ module constitutive
constitutive_TandItsTangent, & constitutive_TandItsTangent, &
constitutive_collectDotState, & constitutive_collectDotState, &
constitutive_collectDeltaState, & constitutive_collectDeltaState, &
constitutive_dotTemperature, & constitutive_heat, &
constitutive_postResults constitutive_postResults
private :: & private :: &
@ -463,12 +463,24 @@ pure function constitutive_homogenizedC(ipc,ip,el)
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
material_phase material_phase
use constitutive_none use constitutive_none, only: &
use constitutive_j2 CONSTITUTIVE_NONE_label, &
use constitutive_phenopowerlaw constitutive_none_homogenizedC
use constitutive_titanmod use constitutive_j2, only: &
use constitutive_dislotwin CONSTITUTIVE_J2_label, &
use constitutive_nonlocal constitutive_j2_homogenizedC
use constitutive_phenopowerlaw, only: &
CONSTITUTIVE_PHENOPOWERLAW_label, &
constitutive_phenopowerlaw_homogenizedC
use constitutive_titanmod, only: &
CONSTITUTIVE_TITANMOD_label, &
constitutive_titanmod_homogenizedC
use constitutive_dislotwin, only: &
CONSTITUTIVE_DISLOTWIN_label, &
constitutive_dislotwin_homogenizedC
use constitutive_nonlocal, only: &
CONSTITUTIVE_NONLOCAL_label, &
constitutive_nonlocal_homogenizedC
implicit none implicit none
real(pReal), dimension(6,6) :: constitutive_homogenizedC real(pReal), dimension(6,6) :: constitutive_homogenizedC
@ -480,13 +492,16 @@ pure function constitutive_homogenizedC(ipc,ip,el)
select case (phase_plasticity(material_phase(ipc,ip,el))) select case (phase_plasticity(material_phase(ipc,ip,el)))
case (constitutive_none_label) case (constitutive_none_label)
constitutive_homogenizedC = constitutive_none_homogenizedC(constitutive_state,ipc,ip,el) constitutive_homogenizedC = constitutive_none_homogenizedC(ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
constitutive_homogenizedC = constitutive_j2_homogenizedC(constitutive_state,ipc,ip,el) constitutive_homogenizedC = constitutive_j2_homogenizedC(ipc,ip,el)
case (constitutive_phenopowerlaw_label) case (constitutive_phenopowerlaw_label)
constitutive_homogenizedC = constitutive_phenopowerlaw_homogenizedC(constitutive_state,ipc,ip,el) constitutive_homogenizedC = constitutive_phenopowerlaw_homogenizedC(ipc,ip,el)
case (constitutive_nonlocal_label)
constitutive_homogenizedC = constitutive_nonlocal_homogenizedC(ipc,ip,el)
case (constitutive_titanmod_label) case (constitutive_titanmod_label)
constitutive_homogenizedC = constitutive_titanmod_homogenizedC(constitutive_state,ipc,ip,el) constitutive_homogenizedC = constitutive_titanmod_homogenizedC(constitutive_state,ipc,ip,el)
@ -494,9 +509,6 @@ pure function constitutive_homogenizedC(ipc,ip,el)
case (constitutive_dislotwin_label) case (constitutive_dislotwin_label)
constitutive_homogenizedC = constitutive_dislotwin_homogenizedC(constitutive_state,ipc,ip,el) constitutive_homogenizedC = constitutive_dislotwin_homogenizedC(constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
constitutive_homogenizedC = constitutive_nonlocal_homogenizedC(constitutive_state,ipc,ip,el)
end select end select
end function constitutive_homogenizedC end function constitutive_homogenizedC
@ -505,7 +517,7 @@ end function constitutive_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different constitutive models !> @brief calls microstructure function of the different constitutive models
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_microstructure(Temperature, Fe, Fp, ipc, ip, el) subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el)
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
material_phase material_phase
@ -525,7 +537,7 @@ subroutine constitutive_microstructure(Temperature, Fe, Fp, ipc, ip, el)
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature temperature
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fe, & !< elastic deformation gradient Fe, & !< elastic deformation gradient
Fp !< plastic deformation gradient Fp !< plastic deformation gradient
@ -533,13 +545,13 @@ subroutine constitutive_microstructure(Temperature, Fe, Fp, ipc, ip, el)
select case (phase_plasticity(material_phase(ipc,ip,el))) select case (phase_plasticity(material_phase(ipc,ip,el)))
case (constitutive_titanmod_label) case (constitutive_titanmod_label)
call constitutive_titanmod_microstructure(Temperature,constitutive_state,ipc,ip,el) call constitutive_titanmod_microstructure(temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislotwin_label) case (constitutive_dislotwin_label)
call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el) call constitutive_dislotwin_microstructure(temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Fe, Fp, ipc,ip,el) call constitutive_nonlocal_microstructure(constitutive_state,Fe,Fp,ipc,ip,el)
end select end select
@ -549,13 +561,14 @@ end subroutine constitutive_microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ipc, ip, el) subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ipc, ip, el)
use math, only: &
math_identity2nd
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
material_phase material_phase
use constitutive_none, only: & use constitutive_none, only: &
constitutive_none_label, & constitutive_none_label
constitutive_none_LpAndItsTangent
use constitutive_j2, only: & use constitutive_j2, only: &
constitutive_j2_label, & constitutive_j2_label, &
constitutive_j2_LpAndItsTangent constitutive_j2_LpAndItsTangent
@ -589,22 +602,23 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ip
select case (phase_plasticity(material_phase(ipc,ip,el))) select case (phase_plasticity(material_phase(ipc,ip,el)))
case (constitutive_none_label) case (constitutive_none_label)
call constitutive_none_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) Lp = 0.0_pReal
dLp_dTstar = math_identity2nd(9)
case (constitutive_j2_label) case (constitutive_j2_label)
call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label) case (constitutive_phenopowerlaw_label)
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,constitutive_state,ipc,ip,el)
case (constitutive_titanmod_label) case (constitutive_titanmod_label)
call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislotwin_label) case (constitutive_dislotwin_label)
call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state(ipc,ip,el), ipc,ip,el) call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, constitutive_state(ipc,ip,el), ipc,ip,el)
end select end select
@ -744,10 +758,10 @@ subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, sub
constitutive_dotState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway? constitutive_dotState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway?
case (constitutive_j2_label) case (constitutive_j2_label)
constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label) case (constitutive_phenopowerlaw_label)
constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,constitutive_state,ipc,ip,el)
case (constitutive_titanmod_label) case (constitutive_titanmod_label)
constitutive_dotState(ipc,ip,el)%p = constitutive_titanmod_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotState(ipc,ip,el)%p = constitutive_titanmod_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
@ -778,7 +792,7 @@ end subroutine constitutive_collectDotState
!> @brief contains the constitutive equation for calculating the incremental change of !> @brief contains the constitutive equation for calculating the incremental change of
!> microstructure based on the current stress and state !> microstructure based on the current stress and state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDeltaState(Tstar_v, Temperature, ipc, ip, el) subroutine constitutive_collectDeltaState(Tstar_v, ipc, ip, el)
use prec, only: & use prec, only: &
pLongInt pLongInt
use debug, only: & use debug, only: &
@ -799,8 +813,6 @@ subroutine constitutive_collectDeltaState(Tstar_v, Temperature, ipc, ip, el)
ipc, & !< grain number ipc, & !< grain number
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
real(pReal), intent(in) :: &
Temperature
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress Tstar_v !< 2nd Piola-Kirchhoff stress
integer(pLongInt) :: & integer(pLongInt) :: &
@ -814,7 +826,7 @@ subroutine constitutive_collectDeltaState(Tstar_v, Temperature, ipc, ip, el)
select case (phase_plasticity(material_phase(ipc,ip,el))) select case (phase_plasticity(material_phase(ipc,ip,el)))
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),constitutive_state, Tstar_v,Temperature,ipc,ip,el) call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),constitutive_state, Tstar_v,ipc,ip,el)
case default case default
constitutive_deltaState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway? constitutive_deltaState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway?
@ -837,7 +849,7 @@ end subroutine constitutive_collectDeltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el) real(pReal) function constitutive_heat(Tstar_v,Temperature,ipc,ip,el)
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< grain number ipc, & !< grain number
@ -847,14 +859,16 @@ real(pReal) function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
Temperature Temperature
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress Tstar_v !< 2nd Piola-Kirchhoff stress
constitutive_dotTemperature = 0.0_pReal
end function constitutive_dotTemperature constitutive_heat = 0.0_pReal
end function constitutive_heat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results !> @brief returns array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el) function constitutive_postResults(Tstar_v, Fe, temperature, ipc, ip, el)
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
@ -888,8 +902,7 @@ function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el)
real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: & real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: &
constitutive_postResults constitutive_postResults
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature, & temperature
dt !< timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
Fe !< elastic deformation gradient Fe !< elastic deformation gradient
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
@ -902,20 +915,20 @@ function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el)
case (constitutive_none_label) case (constitutive_none_label)
constitutive_postResults = 0.0_pReal constitutive_postResults = 0.0_pReal
case (constitutive_titanmod_label)
constitutive_postResults = constitutive_titanmod_postResults(constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) constitutive_postResults = constitutive_j2_postResults(Tstar_v,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label) case (constitutive_phenopowerlaw_label)
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,constitutive_state,ipc,ip,el)
case (constitutive_titanmod_label)
constitutive_postResults = constitutive_titanmod_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_dislotwin_label) case (constitutive_dislotwin_label)
constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, Fe, Temperature, dt, constitutive_state, & constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, Fe, constitutive_state, &
constitutive_dotstate(ipc,ip,el), ipc, ip, el) constitutive_dotstate(ipc,ip,el), ipc, ip, el)
end select end select

11
code/constitutive_dislotwin.f90
View File

@ -934,10 +934,10 @@ pure function constitutive_dislotwin_homogenizedC(state,ipc,ip,el)
sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0 sumf = sum(state(ipc,ip,el)%p((3_pInt*ns+1_pInt):(3_pInt*ns+nt))) ! safe for nt == 0
!* Homogenized elasticity matrix !* Homogenized elasticity matrix
constitutive_dislotwin_homogenizedC = (1.0_pReal-sumf)*constitutive_dislotwin_Cslip_66(:,:,matID) constitutive_dislotwin_homogenizedC = (1.0_pReal-sumf)*constitutive_dislotwin_Cslip_66(1:6,1:6,matID)
do i=1_pInt,nt do i=1_pInt,nt
constitutive_dislotwin_homogenizedC = & constitutive_dislotwin_homogenizedC = &
constitutive_dislotwin_homogenizedC + state(ipc,ip,el)%p(3_pInt*ns+i)*constitutive_dislotwin_Ctwin_66(:,:,i,matID) constitutive_dislotwin_homogenizedC + state(ipc,ip,el)%p(3_pInt*ns+i)*constitutive_dislotwin_Ctwin_66(1:6,1:6,i,matID)
enddo enddo
end function constitutive_dislotwin_homogenizedC end function constitutive_dislotwin_homogenizedC
@ -945,7 +945,7 @@ pure function constitutive_dislotwin_homogenizedC(state,ipc,ip,el)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_dislotwin_microstructure(Temperature,state,ipc,ip,el) subroutine constitutive_dislotwin_microstructure(temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: & use math, only: &
@ -1529,7 +1529,7 @@ end function constitutive_dislotwin_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results !> @brief return array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) function constitutive_dislotwin_postResults(Tstar_v,Temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: & use math, only: &
@ -1558,8 +1558,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,state,ipc,ip,
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
temperature, & !< temperature at integration point temperature !< temperature at integration point
dt
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point

20
code/constitutive_j2.f90
View File

@ -340,9 +340,7 @@ end function constitutive_j2_aTolState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix !> @brief returns the homogenized elasticity matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_j2_homogenizedC(state,ipc,ip,el) pure function constitutive_j2_homogenizedC(ipc,ip,el)
use prec, only: &
p_vec
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
@ -358,8 +356,6 @@ pure function constitutive_j2_homogenizedC(state,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(1:6,1:6,& constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(1:6,1:6,&
phase_plasticityInstance(material_phase(ipc,ip,el))) phase_plasticityInstance(material_phase(ipc,ip,el)))
@ -370,8 +366,7 @@ end function constitutive_j2_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ipc,ip,el)
temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: & use math, only: &
@ -396,8 +391,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,&
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at IP
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
@ -449,7 +442,7 @@ end subroutine constitutive_j2_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_j2_dotState(Tstar_v,temperature,state,ipc,ip,el) pure function constitutive_j2_dotState(Tstar_v,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: & use math, only: &
@ -467,8 +460,6 @@ pure function constitutive_j2_dotState(Tstar_v,temperature,state,ipc,ip,el)
constitutive_j2_dotState constitutive_j2_dotState
real(pReal), dimension(6), intent(in):: & real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
@ -533,7 +524,7 @@ end function constitutive_j2_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results !> @brief return array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_j2_postResults(Tstar_v,temperature,dt,state,ipc,ip,el) pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: & use math, only: &
@ -550,9 +541,6 @@ pure function constitutive_j2_postResults(Tstar_v,temperature,dt,state,ipc,ip,el
implicit none implicit none
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature, & !< temperature at integration point
dt
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point

46
code/constitutive_none.f90
View File

@ -49,8 +49,7 @@ module constitutive_none
public :: & public :: &
constitutive_none_init, & constitutive_none_init, &
constitutive_none_homogenizedC, & constitutive_none_homogenizedC
constitutive_none_LpAndItsTangent
contains contains
@ -185,7 +184,7 @@ end subroutine constitutive_none_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix !> @brief returns the homogenized elasticity matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_none_homogenizedC(state,ipc,ip,el) pure function constitutive_none_homogenizedC(ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use mesh, only: & use mesh, only: &
@ -203,51 +202,10 @@ pure function constitutive_none_homogenizedC(state,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
constitutive_none_homogenizedC = constitutive_none_Cslip_66(1:6,1:6,& constitutive_none_homogenizedC = constitutive_none_Cslip_66(1:6,1:6,&
phase_plasticityInstance(material_phase(ipc,ip,el))) phase_plasticityInstance(material_phase(ipc,ip,el)))
end function constitutive_none_homogenizedC end function constitutive_none_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @details dummy function, returns 0.0 and Identity
!--------------------------------------------------------------------------------------------------
pure subroutine constitutive_none_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_dev_v, &
temperature, state, ipc, ip, el)
use prec, only: &
p_vec
use math, only: &
math_identity2nd
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: &
homogenization_maxNgrains, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: &
Tstar_dev_v !< deviatoric part of 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at IP
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
Lp = 0.0_pReal ! set Lp to zero
dLp_dTstar99 = math_identity2nd(9) ! set dLp_dTstar to Identity
end subroutine constitutive_none_LpAndItsTangent
end module constitutive_none end module constitutive_none

440
code/constitutive_nonlocal.f90
View File

@ -45,7 +45,7 @@ character (len=*), parameter, public :: &
CONSTITUTIVE_NONLOCAL_LABEL = 'nonlocal' CONSTITUTIVE_NONLOCAL_LABEL = 'nonlocal'
character(len=22), dimension(11), parameter, private :: & character(len=22), dimension(11), parameter, private :: &
BASICSTATES = (/'rhoSglEdgePosMobile ', & BASICSTATES = ['rhoSglEdgePosMobile ', &
'rhoSglEdgeNegMobile ', & 'rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ', & 'rhoSglScrewPosMobile ', &
'rhoSglScrewNegMobile ', & 'rhoSglScrewNegMobile ', &
@ -55,20 +55,20 @@ BASICSTATES = (/'rhoSglEdgePosMobile ', &
'rhoSglScrewNegImmobile', & 'rhoSglScrewNegImmobile', &
'rhoDipEdge ', & 'rhoDipEdge ', &
'rhoDipScrew ', & 'rhoDipScrew ', &
'accumulatedshear ' /) !< list of "basic" microstructural state variables that are independent from other state variables 'accumulatedshear ' ] !< list of "basic" microstructural state variables that are independent from other state variables
character(len=16), dimension(3), parameter, private :: & character(len=16), dimension(3), parameter, private :: &
DEPENDENTSTATES = (/'rhoForest ', & DEPENDENTSTATES = ['rhoForest ', &
'tauThreshold ', & 'tauThreshold ', &
'tauBack ' /) !< list of microstructural state variables that depend on other state variables 'tauBack ' ] !< list of microstructural state variables that depend on other state variables
character(len=20), dimension(6), parameter, private :: & character(len=20), dimension(6), parameter, private :: &
OTHERSTATES = (/'velocityEdgePos ', & OTHERSTATES = ['velocityEdgePos ', &
'velocityEdgeNeg ', & 'velocityEdgeNeg ', &
'velocityScrewPos ', & 'velocityScrewPos ', &
'velocityScrewNeg ', & 'velocityScrewNeg ', &
'maxDipoleHeightEdge ', & 'maxDipoleHeightEdge ', &
'maxDipoleHeightScrew' /) !< list of other dependent state variables that are not updated by microstructure 'maxDipoleHeightScrew' ] !< list of other dependent state variables that are not updated by microstructure
real(pReal), parameter, private :: & real(pReal), parameter, private :: &
KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin
@ -1200,30 +1200,25 @@ end function constitutive_nonlocal_aTolState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix !> @brief returns the homogenized elasticity matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_nonlocal_homogenizedC(state,g,ip,el) pure function constitutive_nonlocal_homogenizedC(ipc,ip,el)
use mesh, only: &
use mesh, only: mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: homogenization_maxNgrains, & use material, only: &
homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_plasticityInstance phase_plasticityInstance
implicit none
!*** input variables implicit none
integer(pInt), intent(in) :: g, & ! current grain ID integer(pInt), intent(in) :: &
ipc, & ! current grain ID
ip, & ! current integration point ip, & ! current integration point
el ! current element el ! current element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! microstructural state real(pReal), dimension(6,6) :: &
constitutive_nonlocal_homogenizedC
!*** output variables constitutive_nonlocal_homogenizedC = &
real(pReal), dimension(6,6) :: constitutive_nonlocal_homogenizedC ! homogenized elasticity matrix Cslip66(1:6,1:6,phase_plasticityInstance(material_phase(ipc,ip,el)))
!*** local variables
integer(pInt) :: matID ! current instance of this plasticity
matID = phase_plasticityInstance(material_phase(g,ip,el))
constitutive_nonlocal_homogenizedC = Cslip66(1:6,1:6,matID)
end function constitutive_nonlocal_homogenizedC end function constitutive_nonlocal_homogenizedC
@ -1231,7 +1226,7 @@ end function constitutive_nonlocal_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates quantities characterizing the microstructure !> @brief calculates quantities characterizing the microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_nonlocal_microstructure(state, Temperature, Fe, Fp, gr, ip, el) subroutine constitutive_nonlocal_microstructure(state, Fe, Fp, gr, ip, el)
use IO, only: & use IO, only: &
IO_error IO_error
@ -1279,7 +1274,6 @@ implicit none
integer(pInt), intent(in) :: gr, & ! current grain ID integer(pInt), intent(in) :: gr, & ! current grain ID
ip, & ! current integration point ip, & ! current integration point
el ! current element el ! current element
real(pReal), intent(in) :: Temperature ! temperature
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Fe, & ! elastic deformation gradient Fe, & ! elastic deformation gradient
Fp ! elastic deformation gradient Fp ! elastic deformation gradient
@ -1547,7 +1541,7 @@ end subroutine constitutive_nonlocal_microstructure
!> @brief calculates kinetics !> @brief calculates kinetics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, & subroutine constitutive_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, &
tauThreshold, c, Temperature, g, ip, el) tauThreshold, c, Temperature, ipc, ip, el)
use debug, only: debug_level, & use debug, only: debug_level, &
debug_constitutive, & debug_constitutive, &
@ -1563,18 +1557,18 @@ use material, only: material_phase, &
implicit none implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: g, & !< current grain number integer(pInt), intent(in) :: ipc, & !< current grain number
ip, & !< current integration point ip, & !< current integration point
el, & !< current element number el, & !< current element number
c !< dislocation character (1:edge, 2:screw) c !< dislocation character (1:edge, 2:screw)
real(pReal), intent(in) :: Temperature !< temperature real(pReal), intent(in) :: Temperature !< temperature
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))), & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))), &
intent(in) :: tau, & !< resolved external shear stress (without non Schmid effects) intent(in) :: tau, & !< resolved external shear stress (without non Schmid effects)
tauNS, & !< resolved external shear stress (including non Schmid effects) tauNS, & !< resolved external shear stress (including non Schmid effects)
tauThreshold !< threshold shear stress tauThreshold !< threshold shear stress
!*** output variables !*** output variables
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))), & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))), &
intent(out) :: v, & !< velocity intent(out) :: v, & !< velocity
dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions) dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions)
dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions) dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions)
@ -1606,7 +1600,7 @@ real(pReal) tauRel_P, &
mobility !< dislocation mobility mobility !< dislocation mobility
instance = phase_plasticityInstance(material_phase(g,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
ns = totalNslip(instance) ns = totalNslip(instance)
v = 0.0_pReal v = 0.0_pReal
@ -1689,10 +1683,10 @@ endif
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,*) write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_kinetics at el ip g',el,ip,g write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_kinetics at el ip ipc',el,ip,ipc
write(6,*) write(6,*)
write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold / 1e6_pReal write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold / 1e6_pReal
write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tau / MPa', tau / 1e6_pReal write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tau / MPa', tau / 1e6_pReal
@ -1709,7 +1703,7 @@ end subroutine constitutive_nonlocal_kinetics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, g, ip, el) subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, ipc, ip, el)
use math, only: math_Plain3333to99, & use math, only: math_Plain3333to99, &
math_mul6x6, & math_mul6x6, &
@ -1733,7 +1727,7 @@ use mesh, only: mesh_ipVolume
implicit none implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: g, & !< current grain number integer(pInt), intent(in) :: ipc, & !< current grain number
ip, & !< current integration point ip, & !< current integration point
el !< current element number el !< current element number
real(pReal), intent(in) :: Temperature !< temperature real(pReal), intent(in) :: Temperature !< temperature
@ -1758,14 +1752,14 @@ integer(pInt) matID, & !< current in
s, & !< index of my current slip system s, & !< index of my current slip system
sLattice !< index of my current slip system according to lattice order sLattice !< index of my current slip system according to lattice order
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 !< derivative of Lp with respect to Tstar (3x3x3x3 matrix) real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 !< derivative of Lp with respect to Tstar (3x3x3x3 matrix)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: &
rhoSgl !< single dislocation densities (including blocked) rhoSgl !< single dislocation densities (including blocked)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: &
v, & !< velocity v, & !< velocity
tauNS, & !< resolved shear stress including non Schmid and backstress terms tauNS, & !< resolved shear stress including non Schmid and backstress terms
dv_dtau, & !< velocity derivative with respect to the shear stress dv_dtau, & !< velocity derivative with respect to the shear stress
dv_dtauNS !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
tau, & !< resolved shear stress including backstress terms tau, & !< resolved shear stress including backstress terms
gdotTotal, & !< shear rate gdotTotal, & !< shear rate
tauBack, & !< back stress from dislocation gradients on same slip system tauBack, & !< back stress from dislocation gradients on same slip system
@ -1777,7 +1771,7 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,e
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal dLp_dTstar3333 = 0.0_pReal
matID = phase_plasticityInstance(material_phase(g,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
structID = constitutive_nonlocal_structure(matID) structID = constitutive_nonlocal_structure(matID)
ns = totalNslip(matID) ns = totalNslip(matID)
@ -1823,7 +1817,7 @@ tau = tau + tauBack
! edges ! edges
call constitutive_nonlocal_kinetics(v(1:ns,1), dv_dtau(1:ns,1), dv_dtauNS(1:ns,1), & call constitutive_nonlocal_kinetics(v(1:ns,1), dv_dtau(1:ns,1), dv_dtauNS(1:ns,1), &
tau(1:ns), tauNS(1:ns,1), tauThreshold(1:ns), & tau(1:ns), tauNS(1:ns,1), tauThreshold(1:ns), &
1_pInt, Temperature, g, ip, el) 1_pInt, Temperature, ipc, ip, el)
v(1:ns,2) = v(1:ns,1) v(1:ns,2) = v(1:ns,1)
dv_dtau(1:ns,2) = dv_dtau(1:ns,1) dv_dtau(1:ns,2) = dv_dtau(1:ns,1)
dv_dtauNS(1:ns,2) = dv_dtauNS(1:ns,1) dv_dtauNS(1:ns,2) = dv_dtauNS(1:ns,1)
@ -1839,7 +1833,7 @@ else
do t = 3_pInt,4_pInt do t = 3_pInt,4_pInt
call constitutive_nonlocal_kinetics(v(1:ns,t), dv_dtau(1:ns,t), dv_dtauNS(1:ns,t), & call constitutive_nonlocal_kinetics(v(1:ns,t), dv_dtau(1:ns,t), dv_dtauNS(1:ns,t), &
tau(1:ns), tauNS(1:ns,t), tauThreshold(1:ns), & tau(1:ns), tauNS(1:ns,t), tauThreshold(1:ns), &
2_pInt , Temperature, g, ip, el) 2_pInt , Temperature, ipc, ip, el)
enddo enddo
endif endif
@ -1892,10 +1886,10 @@ dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
write(6,*) write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip g ',el,ip,g write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip ipc ',el,ip,ipc
write(6,*) write(6,*)
write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal
write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',transpose(Lp) write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',transpose(Lp)
@ -1908,7 +1902,7 @@ end subroutine constitutive_nonlocal_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure !> @brief (instantaneous) incremental change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_nonlocal_deltaState(deltaState, state, Tstar_v, Temperature, g,ip,el) subroutine constitutive_nonlocal_deltaState(deltaState, state, Tstar_v, ipc,ip,el)
use debug, only: debug_level, & use debug, only: debug_level, &
debug_constitutive, & debug_constitutive, &
@ -1931,10 +1925,9 @@ use material, only: homogenization_maxNgrains, &
implicit none implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: g, & ! current grain number integer(pInt), intent(in) :: ipc, & ! current grain number
ip, & ! current integration point ip, & ! current integration point
el ! current element number el ! current element number
real(pReal), intent(in) :: Temperature ! temperature
real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation
!*** input/output variables !*** input/output variables
@ -1952,18 +1945,18 @@ integer(pInt) matID, & ! current insta
t, & ! type of dislocation t, & ! type of dislocation
s, & ! index of my current slip system s, & ! index of my current slip system
sLattice ! index of my current slip system according to lattice order sLattice ! index of my current slip system according to lattice order
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),10) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),10) :: &
deltaRho, & ! density increment deltaRho, & ! density increment
deltaRhoRemobilization, & ! density increment by remobilization deltaRhoRemobilization, & ! density increment by remobilization
deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change) deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: &
rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles) rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: &
v ! dislocation glide velocity v ! dislocation glide velocity
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
tau, & ! current resolved shear stress tau, & ! current resolved shear stress
tauBack ! current back stress from pileups on same slip system tauBack ! current back stress from pileups on same slip system
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),2) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: &
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
dLower, & ! minimum stable dipole distance for edges and screws dLower, & ! minimum stable dipole distance for edges and screws
dUpper, & ! current maximum stable dipole distance for edges and screws dUpper, & ! current maximum stable dipole distance for edges and screws
@ -1973,15 +1966,15 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,e
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,*) write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_deltaState at el ip g ',el,ip,g write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_deltaState at el ip ipc ',el,ip,ipc
write(6,*) write(6,*)
endif endif
#endif #endif
matID = phase_plasticityInstance(material_phase(g,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
structID = constitutive_nonlocal_structure(matID) structID = constitutive_nonlocal_structure(matID)
ns = totalNslip(matID) ns = totalNslip(matID)
@ -1990,15 +1983,15 @@ ns = totalNslip(matID)
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
rhoSgl(s,t) = max(state(g,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities
rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,matID)) rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID))
v(s,t) = state(g,ip,el)%p(iV(s,t,matID)) v(s,t) = state(ipc,ip,el)%p(iV(s,t,matID))
endforall endforall
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
rhoDip(s,c) = max(state(g,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities
dUpperOld(s,c) = state(g,ip,el)%p(iD(s,c,matID)) dUpperOld(s,c) = state(ipc,ip,el)%p(iD(s,c,matID))
endforall endforall
tauBack = state(g,ip,el)%p(iTauB(1:ns,matID)) tauBack = state(ipc,ip,el)%p(iTauB(1:ns,matID))
where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) & where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
.or. abs(rhoSgl) < significantRho(matID)) & .or. abs(rhoSgl) < significantRho(matID)) &
@ -2063,7 +2056,7 @@ forall (t=1_pInt:4_pInt) &
!*** store new maximum dipole height in state !*** store new maximum dipole height in state
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
state(g,ip,el)%p(iD(s,c,matID)) = dUpper(s,c) state(ipc,ip,el)%p(iD(s,c,matID)) = dUpper(s,c)
@ -2084,7 +2077,7 @@ forall (s = 1:ns, c = 1_pInt:2_pInt) &
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation remobilization', deltaRhoRemobilization(1:ns,1:8) write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation remobilization', deltaRhoRemobilization(1:ns,1:8)
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress
@ -2099,7 +2092,7 @@ end subroutine constitutive_nonlocal_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, state, state0, timestep, subfrac, g,ip,el) function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, state, state0, timestep, subfrac, ipc,ip,el)
use prec, only: DAMASK_NaN use prec, only: DAMASK_NaN
use numerics, only: numerics_integrationMode, & use numerics, only: numerics_integrationMode, &
@ -2144,7 +2137,7 @@ use lattice, only: lattice_Sslip_v, &
implicit none implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: g, & !< current grain number integer(pInt), intent(in) :: ipc, & !< current grain number
ip, & !< current integration point ip, & !< current integration point
el !< current element number el !< current element number
real(pReal), intent(in) :: Temperature, & !< temperature real(pReal), intent(in) :: Temperature, & !< temperature
@ -2162,7 +2155,7 @@ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), in
!*** input/output variables !*** input/output variables
!*** output variables !*** output variables
real(pReal), dimension(constitutive_nonlocal_sizeDotState(phase_plasticityInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(constitutive_nonlocal_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_nonlocal_dotState !< evolution of state variables / microstructure constitutive_nonlocal_dotState !< evolution of state variables / microstructure
!*** local variables !*** local variables
@ -2184,38 +2177,38 @@ integer(pInt) matID, & !< current inst
s, & !< index of my current slip system s, & !< index of my current slip system
sLattice, & !< index of my current slip system according to lattice order sLattice, & !< index of my current slip system according to lattice order
deads deads
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),10) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),10) :: &
rhoDot, & !< density evolution rhoDot, & !< density evolution
rhoDotMultiplication, & !< density evolution by multiplication rhoDotMultiplication, & !< density evolution by multiplication
rhoDotFlux, & !< density evolution by flux rhoDotFlux, & !< density evolution by flux
rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
rhoDotThermalAnnihilation !< density evolution by thermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
rhoSglOriginal, & rhoSglOriginal, &
neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
rhoSgl0, & !< single dislocation densities at start of cryst inc (positive/negative screw and edge without dipoles) rhoSgl0, & !< single dislocation densities at start of cryst inc (positive/negative screw and edge without dipoles)
my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: &
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, & !< dislocation glide velocity at start of cryst inc v0, & !< dislocation glide velocity at start of cryst inc
my_v, & !< dislocation glide velocity of central ip my_v, & !< dislocation glide velocity of central ip
neighbor_v, & !< dislocation glide velocity of enighboring ip neighbor_v, & !< dislocation glide velocity of enighboring ip
gdot !< shear rates gdot !< shear rates
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
rhoForest, & !< forest dislocation density rhoForest, & !< forest dislocation density
tauThreshold, & !< threshold shear stress tauThreshold, & !< threshold shear stress
tau, & !< current resolved shear stress tau, & !< current resolved shear stress
tauBack, & !< current back stress from pileups on same slip system tauBack, & !< current back stress from pileups on same slip system
vClimb, & !< climb velocity of edge dipoles vClimb, & !< climb velocity of edge dipoles
nSources nSources
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),2) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: &
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
rhoDipOriginal, & rhoDipOriginal, &
dLower, & !< minimum stable dipole distance for edges and screws dLower, & !< minimum stable dipole distance for edges and screws
dUpper !< current maximum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: &
m !< direction of dislocation motion m !< direction of dislocation motion
real(pReal), dimension(3,3) :: my_F, & !< my total deformation gradient real(pReal), dimension(3,3) :: my_F, & !< my total deformation gradient
neighbor_F, & !< total deformation gradient of my neighbor neighbor_F, & !< total deformation gradient of my neighbor
@ -2237,16 +2230,16 @@ logical considerEnteringFlux, &
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
write(6,*) write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip g ',el,ip,g write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip ipc ',el,ip,ipc
write(6,*) write(6,*)
endif endif
#endif #endif
matID = phase_plasticityInstance(material_phase(g,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
structID = constitutive_nonlocal_structure(matID) structID = constitutive_nonlocal_structure(matID)
ns = totalNslip(matID) ns = totalNslip(matID)
@ -2258,16 +2251,16 @@ gdot = 0.0_pReal
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
rhoSgl(s,t) = max(state(g,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities
rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,matID)) rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID))
v(s,t) = state(g,ip,el)%p(iV(s,t,matID)) v(s,t) = state(ipc,ip,el)%p(iV(s,t,matID))
endforall endforall
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
rhoDip(s,c) = max(state(g,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities
endforall endforall
rhoForest = state(g,ip,el)%p(iRhoF(1:ns,matID)) rhoForest = state(ipc,ip,el)%p(iRhoF(1:ns,matID))
tauThreshold = state(g,ip,el)%p(iTauF(1:ns,matID)) tauThreshold = state(ipc,ip,el)%p(iTauF(1:ns,matID))
tauBack = state(g,ip,el)%p(iTauB(1:ns,matID)) tauBack = state(ipc,ip,el)%p(iTauB(1:ns,matID))
rhoSglOriginal = rhoSgl rhoSglOriginal = rhoSgl
rhoDipOriginal = rhoDip rhoDipOriginal = rhoDip
@ -2280,9 +2273,9 @@ where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
if (numerics_timeSyncing) then if (numerics_timeSyncing) then
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
rhoSgl0(s,t) = max(state0(g,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) rhoSgl0(s,t) = max(state0(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal)
rhoSgl0(s,t+4_pInt) = state0(g,ip,el)%p(iRhoB(s,t,matID)) rhoSgl0(s,t+4_pInt) = state0(ipc,ip,el)%p(iRhoB(s,t,matID))
v0(s,t) = state0(g,ip,el)%p(iV(s,t,matID)) v0(s,t) = state0(ipc,ip,el)%p(iV(s,t,matID))
endforall endforall
where (abs(rhoSgl0) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) & where (abs(rhoSgl0) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
.or. abs(rhoSgl0) < significantRho(matID)) & .or. abs(rhoSgl0) < significantRho(matID)) &
@ -2308,7 +2301,7 @@ forall (t = 1_pInt:4_pInt) &
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot
@ -2364,16 +2357,16 @@ else
endwhere endwhere
do s = 1_pInt,ns do s = 1_pInt,ns
if (nSources(s) < 1.0_pReal) then if (nSources(s) < 1.0_pReal) then
if (sourceProbability(s,g,ip,el) > 1.0_pReal) then if (sourceProbability(s,ipc,ip,el) > 1.0_pReal) then
call random_number(rnd) call random_number(rnd)
sourceProbability(s,g,ip,el) = rnd sourceProbability(s,ipc,ip,el) = rnd
!$OMP FLUSH(sourceProbability) !$OMP FLUSH(sourceProbability)
endif endif
if (sourceProbability(s,g,ip,el) > 1.0_pReal - nSources(s)) then if (sourceProbability(s,ipc,ip,el) > 1.0_pReal - nSources(s)) then
rhoDotMultiplication(s,1:4) = sum(rhoSglOriginal(s,1:4) * abs(v(s,1:4))) / meshlength rhoDotMultiplication(s,1:4) = sum(rhoSglOriginal(s,1:4) * abs(v(s,1:4))) / meshlength
endif endif
else else
sourceProbability(s,g,ip,el) = 2.0_pReal sourceProbability(s,ipc,ip,el) = 2.0_pReal
rhoDotMultiplication(s,1:4) = & rhoDotMultiplication(s,1:4) = &
(sum(abs(gdot(s,1:2))) * fEdgeMultiplication(matID) + sum(abs(gdot(s,3:4)))) & (sum(abs(gdot(s,1:2))) * fEdgeMultiplication(matID) + sum(abs(gdot(s,3:4)))) &
/ burgers(s,matID) * sqrt(rhoForest(s)) / lambda0(s,matID) / burgers(s,matID) * sqrt(rhoForest(s)) / lambda0(s,matID)
@ -2381,7 +2374,7 @@ else
enddo enddo
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
write(6,'(a,/,4(12x,12(f12.5,1x),/))') '<< CONST >> sources', nSources write(6,'(a,/,4(12x,12(f12.5,1x),/))') '<< CONST >> sources', nSources
write(6,*) write(6,*)
@ -2401,7 +2394,7 @@ endif
rhoDotFlux = 0.0_pReal rhoDotFlux = 0.0_pReal
if (.not. phase_localPlasticity(material_phase(g,ip,el))) then ! only for nonlocal plasticity if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then ! only for nonlocal plasticity
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
@ -2433,8 +2426,8 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then
m(1:3,1:ns,3) = -lattice_st(1:3, slipSystemLattice(1:ns,matID), structID) m(1:3,1:ns,3) = -lattice_st(1:3, slipSystemLattice(1:ns,matID), structID)
m(1:3,1:ns,4) = lattice_st(1:3, slipSystemLattice(1:ns,matID), structID) m(1:3,1:ns,4) = lattice_st(1:3, slipSystemLattice(1:ns,matID), structID)
my_Fe = Fe(1:3,1:3,g,ip,el) my_Fe = Fe(1:3,1:3,ipc,ip,el)
my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,g,ip,el)) my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,ipc,ip,el))
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! loop through my neighbors do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! loop through my neighbors
neighbor_el = mesh_ipNeighborhood(1,n,ip,el) neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
@ -2447,9 +2440,9 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then
opposite_n = mesh_ipNeighborhood(3,opposite_neighbor,ip,el) opposite_n = mesh_ipNeighborhood(3,opposite_neighbor,ip,el)
if (neighbor_n > 0_pInt) then ! if neighbor exists, average deformation gradient if (neighbor_n > 0_pInt) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phase(g,neighbor_ip,neighbor_el)) neighbor_instance = phase_plasticityInstance(material_phase(ipc,neighbor_ip,neighbor_el))
neighbor_Fe = Fe(1:3,1:3,g,neighbor_ip,neighbor_el) neighbor_Fe = Fe(1:3,1:3,ipc,neighbor_ip,neighbor_el)
neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,g,neighbor_ip,neighbor_el)) neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,ipc,neighbor_ip,neighbor_el))
Favg = 0.5_pReal * (my_F + neighbor_F) Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average else ! if no neighbor, take my value as average
Favg = my_F Favg = my_F
@ -2471,17 +2464,17 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then
endif endif
if (considerEnteringFlux) then if (considerEnteringFlux) then
if(numerics_timeSyncing .and. (subfrac(g,neighbor_ip,neighbor_el) /= subfrac(g,ip,el))) then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal if(numerics_timeSyncing .and. (subfrac(ipc,neighbor_ip,neighbor_el) /= subfrac(ipc,ip,el))) then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal
forall (s = 1:ns, t = 1_pInt:4_pInt) forall (s = 1:ns, t = 1_pInt:4_pInt)
neighbor_v(s,t) = state0(g,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance)) neighbor_v(s,t) = state0(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance))
neighbor_rhoSgl(s,t) = max(state0(g,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal) neighbor_rhoSgl(s,t) = max(state0(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal)
neighbor_rhoSgl(s,t+4_pInt) = state0(g,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance)) neighbor_rhoSgl(s,t+4_pInt) = state0(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance))
endforall endforall
else else
forall (s = 1:ns, t = 1_pInt:4_pInt) forall (s = 1:ns, t = 1_pInt:4_pInt)
neighbor_v(s,t) = state(g,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance)) neighbor_v(s,t) = state(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance))
neighbor_rhoSgl(s,t) = max(state(g,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal) neighbor_rhoSgl(s,t) = max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal)
neighbor_rhoSgl(s,t+4_pInt) = state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance)) neighbor_rhoSgl(s,t+4_pInt) = state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance))
endforall endforall
endif endif
where (abs(neighbor_rhoSgl) * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal & where (abs(neighbor_rhoSgl) * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal &
@ -2541,11 +2534,11 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then
my_rhoSgl = rhoSgl my_rhoSgl = rhoSgl
my_v = v my_v = v
if(numerics_timeSyncing) then if(numerics_timeSyncing) then
if (subfrac(g,ip,el) == 0.0_pReal) then if (subfrac(ipc,ip,el) == 0.0_pReal) then
my_rhoSgl = rhoSgl0 my_rhoSgl = rhoSgl0
my_v = v0 my_v = v0
elseif (neighbor_n > 0_pInt) then elseif (neighbor_n > 0_pInt) then
if (subfrac(g,neighbor_ip,neighbor_el) == 0.0_pReal) then if (subfrac(ipc,neighbor_ip,neighbor_el) == 0.0_pReal) then
my_rhoSgl = rhoSgl0 my_rhoSgl = rhoSgl0
my_v = v0 my_v = v0
endif endif
@ -2658,18 +2651,18 @@ rhoDot = rhoDotFlux &
+ rhoDotThermalAnnihilation + rhoDotThermalAnnihilation
if (numerics_integrationMode == 1_pInt) then ! save rates for output if in central integration mode if (numerics_integrationMode == 1_pInt) then ! save rates for output if in central integration mode
rhoDotFluxOutput(1:ns,1:8,g,ip,el) = rhoDotFlux(1:ns,1:8) rhoDotFluxOutput(1:ns,1:8,ipc,ip,el) = rhoDotFlux(1:ns,1:8)
rhoDotMultiplicationOutput(1:ns,1:2,g,ip,el) = rhoDotMultiplication(1:ns,[1,3]) rhoDotMultiplicationOutput(1:ns,1:2,ipc,ip,el) = rhoDotMultiplication(1:ns,[1,3])
rhoDotSingle2DipoleGlideOutput(1:ns,1:2,g,ip,el) = rhoDotSingle2DipoleGlide(1:ns,9:10) rhoDotSingle2DipoleGlideOutput(1:ns,1:2,ipc,ip,el) = rhoDotSingle2DipoleGlide(1:ns,9:10)
rhoDotAthermalAnnihilationOutput(1:ns,1:2,g,ip,el) = rhoDotAthermalAnnihilation(1:ns,9:10) rhoDotAthermalAnnihilationOutput(1:ns,1:2,ipc,ip,el) = rhoDotAthermalAnnihilation(1:ns,9:10)
rhoDotThermalAnnihilationOutput(1:ns,1:2,g,ip,el) = rhoDotThermalAnnihilation(1:ns,9:10) rhoDotThermalAnnihilationOutput(1:ns,1:2,ipc,ip,el) = rhoDotThermalAnnihilation(1:ns,9:10)
rhoDotEdgeJogsOutput(1:ns,g,ip,el) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1) rhoDotEdgeJogsOutput(1:ns,ipc,ip,el) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1)
endif endif
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', rhoDotMultiplication(1:ns,1:4) * timestep write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', rhoDotMultiplication(1:ns,1:4) * timestep
write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', rhoDotFlux(1:ns,1:8) * timestep write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', rhoDotFlux(1:ns,1:8) * timestep
@ -2892,7 +2885,7 @@ end subroutine constitutive_nonlocal_updateCompatibility
!********************************************************************* !*********************************************************************
!* calculates quantities characterizing the microstructure * !* calculates quantities characterizing the microstructure *
!********************************************************************* !*********************************************************************
function constitutive_nonlocal_dislocationstress(state, Fe, g, ip, el) pure function constitutive_nonlocal_dislocationstress(state, Fe, ipc, ip, el)
use math, only: math_mul33x33, & use math, only: math_mul33x33, &
math_mul33x3, & math_mul33x3, &
@ -2917,7 +2910,7 @@ implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: g, & ! current grain ID integer(pInt), intent(in) :: ipc, & ! current grain ID
ip, & ! current integration point ip, & ! current integration point
el ! current element el ! current element
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
@ -2975,11 +2968,11 @@ real(pReal), dimension(2,2,maxval(totalNslip)) :: &
neighbor_rhoExcess ! excess density at neighbor material point (edge/screw,mobile/dead,slipsystem) neighbor_rhoExcess ! excess density at neighbor material point (edge/screw,mobile/dead,slipsystem)
real(pReal), dimension(2,maxval(totalNslip)) :: & real(pReal), dimension(2,maxval(totalNslip)) :: &
rhoExcessDead rhoExcessDead
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: &
rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-) rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-)
logical inversionError logical inversionError
phase = material_phase(g,ip,el) phase = material_phase(ipc,ip,el)
instance = phase_plasticityInstance(phase) instance = phase_plasticityInstance(phase)
latticeStruct = constitutive_nonlocal_structure(instance) latticeStruct = constitutive_nonlocal_structure(instance)
ns = totalNslip(instance) ns = totalNslip(instance)
@ -2989,8 +2982,8 @@ ns = totalNslip(instance)
!*** get basic states !*** get basic states
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
rhoSgl(s,t) = max(state(g,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities
rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,instance)) rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance))
endforall endforall
@ -3001,7 +2994,7 @@ endforall
constitutive_nonlocal_dislocationstress = 0.0_pReal constitutive_nonlocal_dislocationstress = 0.0_pReal
if (.not. phase_localPlasticity(phase)) then if (.not. phase_localPlasticity(phase)) then
call math_invert33(Fe(1:3,1:3,g,ip,el), invFe, detFe, inversionError) call math_invert33(Fe(1:3,1:3,ipc,ip,el), invFe, detFe, inversionError)
!* in case of periodic surfaces we have to find out how many periodic images in each direction we need !* in case of periodic surfaces we have to find out how many periodic images in each direction we need
@ -3025,7 +3018,7 @@ if (.not. phase_localPlasticity(phase)) then
do neighbor_el = 1_pInt,mesh_NcpElems do neighbor_el = 1_pInt,mesh_NcpElems
ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))) ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)))
neighbor_phase = material_phase(g,neighbor_ip,neighbor_el) neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el)
if (phase_localPlasticity(neighbor_phase)) then if (phase_localPlasticity(neighbor_phase)) then
cycle cycle
endif endif
@ -3035,10 +3028,10 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
call math_invert33(Fe(1:3,1:3,1,neighbor_ip,neighbor_el), neighbor_invFe, detFe, inversionError) call math_invert33(Fe(1:3,1:3,1,neighbor_ip,neighbor_el), neighbor_invFe, detFe, inversionError)
neighbor_ipVolumeSideLength = mesh_ipVolume(neighbor_ip,neighbor_el) ** (1.0_pReal/3.0_pReal) ! reference volume used here neighbor_ipVolumeSideLength = mesh_ipVolume(neighbor_ip,neighbor_el) ** (1.0_pReal/3.0_pReal) ! reference volume used here
forall (s = 1_pInt:neighbor_ns, c = 1_pInt:2_pInt) forall (s = 1_pInt:neighbor_ns, c = 1_pInt:2_pInt)
neighbor_rhoExcess(c,1,s) = state(g,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)) & ! positive mobiles neighbor_rhoExcess(c,1,s) = state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)) & ! positive mobiles
- state(g,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)) ! negative mobiles - state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)) ! negative mobiles
neighbor_rhoExcess(c,2,s) = abs(state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads neighbor_rhoExcess(c,2,s) = abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads
- abs(state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) ! negative deads - abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) ! negative deads
endforall endforall
Tdislo_neighborLattice = 0.0_pReal Tdislo_neighborLattice = 0.0_pReal
do deltaX = periodicImages(1,1),periodicImages(2,1) do deltaX = periodicImages(1,1),periodicImages(2,1)
@ -3097,8 +3090,8 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
cycle cycle
elseif (j > 1_pInt) then elseif (j > 1_pInt) then
x = connection_neighborSlip(1) + sign(0.5_pReal * segmentLength, & x = connection_neighborSlip(1) + sign(0.5_pReal * segmentLength, &
state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,1,neighbor_instance)) & state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,1,neighbor_instance)) &
- state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,2,neighbor_instance))) - state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2,neighbor_instance)))
xsquare = x * x xsquare = x * x
endif endif
@ -3143,8 +3136,8 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
cycle cycle
elseif (j > 1_pInt) then elseif (j > 1_pInt) then
y = connection_neighborSlip(2) + sign(0.5_pReal * segmentLength, & y = connection_neighborSlip(2) + sign(0.5_pReal * segmentLength, &
state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,3,neighbor_instance)) & state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,3,neighbor_instance)) &
- state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,4,neighbor_instance))) - state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,4,neighbor_instance)))
ysquare = y * y ysquare = y * y
endif endif
@ -3197,8 +3190,8 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
else else
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
rhoExcessDead(c,s) = state(g,ip,el)%p(iRhoB(s,2*c-1,instance)) & ! positive deads (here we use symmetry: if this has negative sign it is treated as negative density at positive position instead of positive density at negative position) rhoExcessDead(c,s) = state(ipc,ip,el)%p(iRhoB(s,2*c-1,instance)) & ! positive deads (here we use symmetry: if this has negative sign it is treated as negative density at positive position instead of positive density at negative position)
+ state(g,ip,el)%p(iRhoB(s,2*c,instance)) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position) + state(ipc,ip,el)%p(iRhoB(s,2*c,instance)) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position)
do s = 1_pInt,ns do s = 1_pInt,ns
if (all(abs(rhoExcessDead(:,s)) < significantRho(instance))) then if (all(abs(rhoExcessDead(:,s)) < significantRho(instance))) then
@ -3240,80 +3233,81 @@ endif
end function constitutive_nonlocal_dislocationstress end function constitutive_nonlocal_dislocationstress
!********************************************************************* !--------------------------------------------------------------------------------------------------
!* return array of constitutive results * !> @brief return array of constitutive results
!********************************************************************* !--------------------------------------------------------------------------------------------------
function constitutive_nonlocal_postResults(Tstar_v, Fe, Temperature, dt, state, dotState, g,ip,el) pure function constitutive_nonlocal_postResults(Tstar_v,Fe,state,dotState,ipc,ip,el)
use math, only: &
use math, only: math_mul6x6, & math_mul6x6, &
math_mul33x3, & math_mul33x3, &
math_mul33x33, & math_mul33x33, &
pi pi
use mesh, only: mesh_NcpElems, & use mesh, only: &
mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: homogenization_maxNgrains, & use material, only: &
homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_plasticityInstance, & phase_plasticityInstance, &
phase_Noutput phase_Noutput
use lattice, only: lattice_Sslip_v, & use lattice, only: &
lattice_Sslip_v, &
lattice_sd, & lattice_sd, &
lattice_st, & lattice_st, &
lattice_sn lattice_sn
implicit none implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe !< elastic deformation gradient
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
type(p_vec), intent(in) :: dotState ! evolution rate of microstructural state
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
!*** input variables real(pReal), dimension(constitutive_nonlocal_sizePostResults(&
integer(pInt), intent(in) :: g, & ! current grain number phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
ip, & ! current integration point
el ! current element number
real(pReal), intent(in) :: Temperature, & ! temperature
dt ! time increment
real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe ! elastic deformation gradient
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state ! current microstructural state
type(p_vec), intent(in) :: dotState ! evolution rate of microstructural state
!*** output variables
real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
constitutive_nonlocal_postResults constitutive_nonlocal_postResults
!*** local variables integer(pInt) :: &
integer(pInt) matID, & ! current instance of this plasticity matID, & !< current instance of this plasticity
structID, & ! current lattice structure structID, & !< current lattice structure
ns, & ! short notation for the total number of active slip systems ns, & !< short notation for the total number of active slip systems
c, & ! character of dislocation c, & !< character of dislocation
cs, & ! constitutive result index cs, & !< constitutive result index
o, & ! index of current output o, & !< index of current output
t, & ! type of dislocation t, & !< type of dislocation
s, & ! index of my current slip system s, & !< index of my current slip system
sLattice ! index of my current slip system according to lattice order sLattice !< index of my current slip system according to lattice order
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: &
rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles) rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
rhoDotSgl ! evolution rate of single dislocation densities (positive/negative screw and edge without dipoles) rhoDotSgl !< evolution rate of single dislocation densities (positive/negative screw and edge without dipoles)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: &
gdot, & ! shear rates gdot, & !< shear rates
v ! velocities v !< velocities
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
rhoForest, & ! forest dislocation density rhoForest, & !< forest dislocation density
tauThreshold, & ! threshold shear stress tauThreshold, & !< threshold shear stress
tau, & ! current resolved shear stress tau, & !< current resolved shear stress
tauBack ! back stress from pileups on same slip system tauBack !< back stress from pileups on same slip system
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),2) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: &
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
rhoDotDip, & ! evolution rate of dipole dislocation densities (screw and edge dipoles) rhoDotDip, & !< evolution rate of dipole dislocation densities (screw and edge dipoles)
dLower, & ! minimum stable dipole distance for edges and screws dLower, & !< minimum stable dipole distance for edges and screws
dUpper ! current maximum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),2) :: & real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: &
m, & ! direction of dislocation motion for edge and screw (unit vector) m, & !< direction of dislocation motion for edge and screw (unit vector)
m_currentconf ! direction of dislocation motion for edge and screw (unit vector) in current configuration m_currentconf !< direction of dislocation motion for edge and screw (unit vector) in current configuration
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
n_currentconf ! slip system normal (unit vector) in current configuration n_currentconf !< slip system normal (unit vector) in current configuration
real(pReal), dimension(3,3) :: sigma real(pReal), dimension(3,3) :: &
sigma
matID = phase_plasticityInstance(material_phase(g,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
structID = constitutive_nonlocal_structure(matID) structID = constitutive_nonlocal_structure(matID)
ns = totalNslip(matID) ns = totalNslip(matID)
@ -3324,19 +3318,19 @@ constitutive_nonlocal_postResults = 0.0_pReal
!* short hand notations for state variables !* short hand notations for state variables
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
rhoSgl(s,t) = state(g,ip,el)%p(iRhoU(s,t,matID)) rhoSgl(s,t) = state(ipc,ip,el)%p(iRhoU(s,t,matID))
rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,matID)) rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID))
v(s,t) = state(g,ip,el)%p(iV(s,t,matID)) v(s,t) = state(ipc,ip,el)%p(iV(s,t,matID))
rhoDotSgl(s,t) = dotState%p(iRhoU(s,t,matID)) rhoDotSgl(s,t) = dotState%p(iRhoU(s,t,matID))
rhoDotSgl(s,t+4_pInt) = dotState%p(iRhoB(s,t,matID)) rhoDotSgl(s,t+4_pInt) = dotState%p(iRhoB(s,t,matID))
endforall endforall
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
rhoDip(s,c) = state(g,ip,el)%p(iRhoD(s,c,matID)) rhoDip(s,c) = state(ipc,ip,el)%p(iRhoD(s,c,matID))
rhoDotDip(s,c) = dotState%p(iRhoD(s,c,matID)) rhoDotDip(s,c) = dotState%p(iRhoD(s,c,matID))
endforall endforall
rhoForest = state(g,ip,el)%p(iRhoF(1:ns,matID)) rhoForest = state(ipc,ip,el)%p(iRhoF(1:ns,matID))
tauThreshold = state(g,ip,el)%p(iTauF(1:ns,matID)) tauThreshold = state(ipc,ip,el)%p(iTauF(1:ns,matID))
tauBack = state(g,ip,el)%p(iTauB(1:ns,matID)) tauBack = state(ipc,ip,el)%p(iTauB(1:ns,matID))
@ -3371,13 +3365,13 @@ dUpper = max(dUpper,dLower)
m(1:3,1:ns,1) = lattice_sd(1:3,slipSystemLattice(1:ns,matID),structID) m(1:3,1:ns,1) = lattice_sd(1:3,slipSystemLattice(1:ns,matID),structID)
m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,matID),structID) m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,matID),structID)
forall (c = 1_pInt:2_pInt, s = 1_pInt:ns) & forall (c = 1_pInt:2_pInt, s = 1_pInt:ns) &
m_currentconf(1:3,s,c) = math_mul33x3(Fe(1:3,1:3,g,ip,el), m(1:3,s,c)) m_currentconf(1:3,s,c) = math_mul33x3(Fe(1:3,1:3,ipc,ip,el), m(1:3,s,c))
forall (s = 1_pInt:ns) & forall (s = 1_pInt:ns) &
n_currentconf(1:3,s) = math_mul33x3(Fe(1:3,1:3,g,ip,el), & n_currentconf(1:3,s) = math_mul33x3(Fe(1:3,1:3,ipc,ip,el), &
lattice_sn(1:3,slipSystemLattice(s,matID),structID)) lattice_sn(1:3,slipSystemLattice(s,matID),structID))
do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_nonlocal_output(o,matID)) select case(constitutive_nonlocal_output(o,matID))
case ('rho') case ('rho')
@ -3557,66 +3551,66 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
cs = cs + ns cs = cs + ns
case ('rho_dot_gen') case ('rho_dot_gen')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,g,ip,el) & constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,ipc,ip,el) &
+ rhoDotMultiplicationOutput(1:ns,2,g,ip,el) + rhoDotMultiplicationOutput(1:ns,2,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_gen_edge') case ('rho_dot_gen_edge')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,g,ip,el) constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_gen_screw') case ('rho_dot_gen_screw')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,2,g,ip,el) constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,2,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_sgl2dip') case ('rho_dot_sgl2dip')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,g,ip,el) & constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,ipc,ip,el) &
+ rhoDotSingle2DipoleGlideOutput(1:ns,2,g,ip,el) + rhoDotSingle2DipoleGlideOutput(1:ns,2,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_sgl2dip_edge') case ('rho_dot_sgl2dip_edge')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,g,ip,el) constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_sgl2dip_screw') case ('rho_dot_sgl2dip_screw')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,2,g,ip,el) constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,2,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_ann_ath') case ('rho_dot_ann_ath')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotAthermalAnnihilationOutput(1:ns,1,g,ip,el) & constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotAthermalAnnihilationOutput(1:ns,1,ipc,ip,el) &
+ rhoDotAthermalAnnihilationOutput(1:ns,2,g,ip,el) + rhoDotAthermalAnnihilationOutput(1:ns,2,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_ann_the') case ('rho_dot_ann_the')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,g,ip,el) & constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,ipc,ip,el) &
+ rhoDotThermalAnnihilationOutput(1:ns,2,g,ip,el) + rhoDotThermalAnnihilationOutput(1:ns,2,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_ann_the_edge') case ('rho_dot_ann_the_edge')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,g,ip,el) constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_ann_the_screw') case ('rho_dot_ann_the_screw')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,2,g,ip,el) constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,2,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_edgejogs') case ('rho_dot_edgejogs')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotEdgeJogsOutput(1:ns,g,ip,el) constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotEdgeJogsOutput(1:ns,ipc,ip,el)
cs = cs + ns cs = cs + ns
case ('rho_dot_flux') case ('rho_dot_flux')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:4,g,ip,el),2) & constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:4,ipc,ip,el),2) &
+ sum(abs(rhoDotFluxOutput(1:ns,5:8,g,ip,el)),2) + sum(abs(rhoDotFluxOutput(1:ns,5:8,ipc,ip,el)),2)
cs = cs + ns cs = cs + ns
case ('rho_dot_flux_edge') case ('rho_dot_flux_edge')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:2,g,ip,el),2) & constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:2,ipc,ip,el),2) &
+ sum(abs(rhoDotFluxOutput(1:ns,5:6,g,ip,el)),2) + sum(abs(rhoDotFluxOutput(1:ns,5:6,ipc,ip,el)),2)
cs = cs + ns cs = cs + ns
case ('rho_dot_flux_screw') case ('rho_dot_flux_screw')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,3:4,g,ip,el),2) & constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,3:4,ipc,ip,el),2) &
+ sum(abs(rhoDotFluxOutput(1:ns,7:8,g,ip,el)),2) + sum(abs(rhoDotFluxOutput(1:ns,7:8,ipc,ip,el)),2)
cs = cs + ns cs = cs + ns
case ('velocity_edge_pos') case ('velocity_edge_pos')
@ -3716,7 +3710,7 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
cs = cs + ns cs = cs + ns
case('dislocationstress') case('dislocationstress')
sigma = constitutive_nonlocal_dislocationstress(state, Fe, g, ip, el) sigma = constitutive_nonlocal_dislocationstress(state, Fe, ipc, ip, el)
constitutive_nonlocal_postResults(cs+1_pInt) = sigma(1,1) constitutive_nonlocal_postResults(cs+1_pInt) = sigma(1,1)
constitutive_nonlocal_postResults(cs+2_pInt) = sigma(2,2) constitutive_nonlocal_postResults(cs+2_pInt) = sigma(2,2)
constitutive_nonlocal_postResults(cs+3_pInt) = sigma(3,3) constitutive_nonlocal_postResults(cs+3_pInt) = sigma(3,3)
@ -3726,7 +3720,7 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
cs = cs + 6_pInt cs = cs + 6_pInt
case('accumulatedshear') case('accumulatedshear')
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(iGamma(1:ns,matID)) constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(ipc,ip,el)%p(iGamma(1:ns,matID))
cs = cs + ns cs = cs + ns
case('boundarylayer') case('boundarylayer')
@ -3742,6 +3736,6 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
end select end select
enddo enddo
endfunction end function constitutive_nonlocal_postResults
END MODULE end module constitutive_nonlocal

18
code/constitutive_phenopowerlaw.f90
View File

@ -646,7 +646,7 @@ end function constitutive_phenopowerlaw_aTolState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix !> @brief returns the homogenized elasticity matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el) pure function constitutive_phenopowerlaw_homogenizedC(ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use mesh, only: & use mesh, only: &
@ -664,8 +664,6 @@ pure function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
constitutive_phenopowerlaw_homogenizedC = constitutive_phenopowerlaw_Cslip_66(1:6,1:6,& constitutive_phenopowerlaw_homogenizedC = constitutive_phenopowerlaw_Cslip_66(1:6,1:6,&
phase_plasticityInstance(material_phase(ipc,ip,el))) phase_plasticityInstance(material_phase(ipc,ip,el)))
@ -676,8 +674,7 @@ end function constitutive_phenopowerlaw_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ipc,ip,el)
temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: & use math, only: &
@ -709,8 +706,6 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at IP
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
@ -831,7 +826,7 @@ end subroutine constitutive_phenopowerlaw_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_dotState(Tstar_v,temperature,state,ipc,ip,el) function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use lattice, only: & use lattice, only: &
@ -854,8 +849,6 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,temperature,state,ipc,ip,el
implicit none implicit none
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
@ -1003,7 +996,7 @@ end function constitutive_phenopowerlaw_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results !> @brief return array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_postResults(Tstar_v,temperature,dt,state,ipc,ip,el) pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use mesh, only: & use mesh, only: &
@ -1029,9 +1022,6 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,temperature,dt,stat
implicit none implicit none
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature, & !< temperature at integration point
dt
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point

21
code/constitutive_titanmod.f90
View File

@ -1144,11 +1144,11 @@ real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! homogenized elasticity matrix ! homogenized elasticity matrix
constitutive_titanmod_homogenizedC = (1.0_pReal-sumf)*constitutive_titanmod_Cslip_66(:,:,matID) constitutive_titanmod_homogenizedC = (1.0_pReal-sumf)*constitutive_titanmod_Cslip_66(1:6,1:6,matID)
do i=1_pInt,nt do i=1_pInt,nt
constitutive_titanmod_homogenizedC = & constitutive_titanmod_homogenizedC = constitutive_titanmod_homogenizedC &
constitutive_titanmod_homogenizedC + volumefraction_PerTwinSys(i)*constitutive_titanmod_Ctwin_66(:,:,i,matID) + volumefraction_PerTwinSys(i)*&
constitutive_titanmod_Ctwin_66(1:6,1:6,i,matID)
enddo enddo
end function constitutive_titanmod_homogenizedC end function constitutive_titanmod_homogenizedC
@ -1717,7 +1717,7 @@ end function constitutive_titanmod_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results !> @brief return array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_titanmod_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) pure function constitutive_titanmod_postResults(state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use mesh, only: & use mesh, only: &
@ -1730,24 +1730,21 @@ pure function constitutive_titanmod_postResults(Tstar_v,Temperature,dt,state,ipc
phase_Noutput phase_Noutput
implicit none implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature, & !< temperature at integration point
dt
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state state !< microstructure state
real(pReal), dimension(constitutive_titanmod_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_titanmod_postResults
integer(pInt) :: & integer(pInt) :: &
matID, structID,& matID, structID,&
ns,nt,& ns,nt,&
o,i,c o,i,c
real(pReal) :: sumf real(pReal) :: sumf
real(pReal), dimension(constitutive_titanmod_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_titanmod_postResults
real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
volumefraction_PerTwinSys volumefraction_PerTwinSys

342
code/crystallite.f90
View File

File diff suppressed because it is too large Load Diff

33
code/debug.f90
View File

@ -61,25 +61,23 @@ module debug
debug_MARC = 12_pInt, & debug_MARC = 12_pInt, &
debug_ABAQUS = 13_pInt debug_ABAQUS = 13_pInt
integer(pInt), parameter, private :: & integer(pInt), parameter, private :: &
debug_MAXNTYPE = debug_ABAQUS ! must be set to the maximum defined debug type debug_MAXNTYPE = debug_ABAQUS !< must be set to the maximum defined debug type
integer(pInt),protected, dimension(debug_maxNtype+2_pInt), public :: & ! specific ones, and 2 for "all" and "other" integer(pInt),protected, dimension(debug_maxNtype+2_pInt), public :: & ! specific ones, and 2 for "all" and "other"
debug_level = 0_pInt debug_level = 0_pInt
integer(pInt), public :: & integer(pInt), public :: &
debug_cumLpCalls = 0_pInt, & ! total number of calls to LpAndItsTangent debug_cumLpCalls = 0_pInt, & !< total number of calls to LpAndItsTangent
debug_cumDeltaStateCalls = 0_pInt, & ! total number of calls to deltaState debug_cumDeltaStateCalls = 0_pInt, & !< total number of calls to deltaState
debug_cumDotStateCalls = 0_pInt, & ! total number of calls to dotState debug_cumDotStateCalls = 0_pInt, & !< total number of calls to dotState
debug_cumDotTemperatureCalls = 0_pInt, & ! total number of calls to dotTemprature
debug_e = 1_pInt, & debug_e = 1_pInt, &
debug_i = 1_pInt, & debug_i = 1_pInt, &
debug_g = 1_pInt debug_g = 1_pInt
integer(pLongInt), public :: & integer(pLongInt), public :: &
debug_cumLpTicks = 0_pLongInt, & ! total cpu ticks spent in LpAndItsTangent debug_cumLpTicks = 0_pLongInt, & !< total cpu ticks spent in LpAndItsTangent
debug_cumDeltaStateTicks = 0_pLongInt, & ! total cpu ticks spent in deltaState debug_cumDeltaStateTicks = 0_pLongInt, & !< total cpu ticks spent in deltaState
debug_cumDotStateTicks = 0_pLongInt, & ! total cpu ticks spent in dotState debug_cumDotStateTicks = 0_pLongInt !< total cpu ticks spent in dotState
debug_cumDotTemperatureTicks = 0_pLongInt ! total cpu ticks spent in dotTemperature
integer(pInt), dimension(2), public :: & integer(pInt), dimension(2), public :: &
debug_stressMaxLocation = 0_pInt, & debug_stressMaxLocation = 0_pInt, &
@ -88,13 +86,13 @@ module debug
debug_jacobianMinLocation = 0_pInt debug_jacobianMinLocation = 0_pInt
integer(pInt), dimension(:), allocatable, public :: & integer(pInt), dimension(:), allocatable, public :: &
debug_CrystalliteLoopDistribution, & ! distribution of crystallite cutbacks debug_CrystalliteLoopDistribution, & !< distribution of crystallite cutbacks
debug_MaterialpointStateLoopDistribution, & debug_MaterialpointStateLoopDistribution, &
debug_MaterialpointLoopDistribution debug_MaterialpointLoopDistribution
integer(pInt), dimension(:,:), allocatable, public :: & integer(pInt), dimension(:,:), allocatable, public :: &
debug_StressLoopDistribution, & ! distribution of stress iterations until convergence debug_StressLoopDistribution, & !< distribution of stress iterations until convergence
debug_StateLoopDistribution ! distribution of state iterations until convergence debug_StateLoopDistribution !< distribution of state iterations until convergence
real(pReal), public :: & real(pReal), public :: &
debug_stressMax = -huge(1.0_pReal), & debug_stressMax = -huge(1.0_pReal), &
@ -103,7 +101,7 @@ module debug
debug_jacobianMin = huge(1.0_pReal) debug_jacobianMin = huge(1.0_pReal)
character(len=64), parameter, private :: & character(len=64), parameter, private :: &
debug_CONFIGFILE = 'debug.config' ! name of configuration file debug_CONFIGFILE = 'debug.config' !< name of configuration file
#ifdef PETSc #ifdef PETSc
character(len=1024), parameter, public :: & character(len=1024), parameter, public :: &
@ -335,11 +333,9 @@ subroutine debug_reset
debug_cumLpTicks = 0_pLongInt debug_cumLpTicks = 0_pLongInt
debug_cumDeltaStateTicks = 0_pLongInt debug_cumDeltaStateTicks = 0_pLongInt
debug_cumDotStateTicks = 0_pLongInt debug_cumDotStateTicks = 0_pLongInt
debug_cumDotTemperatureTicks = 0_pLongInt
debug_cumLpCalls = 0_pInt debug_cumLpCalls = 0_pInt
debug_cumDeltaStateCalls = 0_pInt debug_cumDeltaStateCalls = 0_pInt
debug_cumDotStateCalls = 0_pInt debug_cumDotStateCalls = 0_pInt
debug_cumDotTemperatureCalls = 0_pInt
debug_stressMaxLocation = 0_pInt debug_stressMaxLocation = 0_pInt
debug_stressMinLocation = 0_pInt debug_stressMinLocation = 0_pInt
debug_jacobianMaxLocation = 0_pInt debug_jacobianMaxLocation = 0_pInt
@ -394,13 +390,6 @@ subroutine debug_info
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',& write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
real(debug_cumDeltaStateTicks,pReal)*1.0e6_pReal/real(tickrate*debug_cumDeltaStateCalls,pReal) real(debug_cumDeltaStateTicks,pReal)*1.0e6_pReal/real(tickrate*debug_cumDeltaStateCalls,pReal)
endif endif
write(6,'(/,a33,1x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls
if (debug_cumdotTemperatureCalls > 0_pInt) then
write(6,'(a33,1x,f12.3)') 'total CPU time/s :',&
real(debug_cumDotTemperatureTicks,pReal)/real(tickrate,pReal)
write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',&
real(debug_cumDotTemperatureTicks,pReal)*1.0e6_pReal/real(tickrate*debug_cumDotTemperatureCalls,pReal)
endif
integral = 0_pInt integral = 0_pInt
write(6,'(3/,a)') 'distribution_StressLoop : stress stiffness' write(6,'(3/,a)') 'distribution_StressLoop : stress stiffness'

93
code/homogenization.f90
View File

@ -36,8 +36,6 @@ module homogenization
private private
type(p_vec), dimension(:,:), allocatable, public :: & type(p_vec), dimension(:,:), allocatable, public :: &
homogenization_state0 !< pointer array to homogenization state at start of FE increment homogenization_state0 !< pointer array to homogenization state at start of FE increment
real(pReal), dimension(:,:), allocatable, public :: &
materialpoint_Temperature !< temperature at IP
real(pReal), dimension(:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:), allocatable, public :: &
materialpoint_F0, & !< def grad of IP at start of FE increment materialpoint_F0, & !< def grad of IP at start of FE increment
materialpoint_F, & !< def grad of IP to be reached at end of FE increment materialpoint_F, & !< def grad of IP to be reached at end of FE increment
@ -63,7 +61,8 @@ module homogenization
real(pReal), dimension(:,:), allocatable, private :: & real(pReal), dimension(:,:), allocatable, private :: &
materialpoint_subFrac, & materialpoint_subFrac, &
materialpoint_subStep, & materialpoint_subStep, &
materialpoint_subdt materialpoint_subdt, &
materialpoint_heat
integer(pInt), dimension(:,:), allocatable, private :: & integer(pInt), dimension(:,:), allocatable, private :: &
homogenization_sizePostResults !< size of postResults array per material point homogenization_sizePostResults !< size of postResults array per material point
integer(pInt), private :: & integer(pInt), private :: &
@ -82,7 +81,7 @@ module homogenization
homogenization_partitionDeformation, & homogenization_partitionDeformation, &
homogenization_updateState, & homogenization_updateState, &
homogenization_averageStressAndItsTangent, & homogenization_averageStressAndItsTangent, &
homogenization_averageTemperature, & homogenization_averageHeat, &
homogenization_postResults homogenization_postResults
contains contains
@ -91,14 +90,16 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief module initialization !> @brief module initialization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_init(Temperature) subroutine homogenization_init()
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use math, only: & use math, only: &
math_I3 math_I3
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_homogenization, & debug_homogenization, &
debug_levelBasic debug_levelBasic, &
debug_e, &
debug_g
use IO, only: & use IO, only: &
IO_error, & IO_error, &
IO_open_file, & IO_open_file, &
@ -121,18 +122,16 @@ subroutine homogenization_init(Temperature)
use homogenization_RGC use homogenization_RGC
implicit none implicit none
real(pReal) Temperature
integer(pInt), parameter :: fileunit = 200_pInt integer(pInt), parameter :: fileunit = 200_pInt
integer(pInt) e,i,p,myInstance integer(pInt) :: e,i,p,myInstance
integer(pInt), dimension(:,:), pointer :: thisSize integer(pInt), dimension(:,:), pointer :: thisSize
character(len=64), dimension(:,:), pointer :: thisOutput character(len=64), dimension(:,:), pointer :: thisOutput
logical :: knownHomogenization logical :: knownHomogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! parse homogenization from config file ! parse homogenization from config file
if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) then ! no local material configuration present... if (.not. IO_open_jobFile_stat(fileunit,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(fileunit,material_configFile) ! ... open material.config file call IO_open_file(fileunit,material_configFile) ! ... open material.config file
endif
call homogenization_isostrain_init(fileunit) call homogenization_isostrain_init(fileunit)
call homogenization_RGC_init(fileunit) call homogenization_RGC_init(fileunit)
close(fileunit) close(fileunit)
@ -173,7 +172,8 @@ subroutine homogenization_init(Temperature)
homogenization_sizeState = 0_pInt homogenization_sizeState = 0_pInt
allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems)) allocate(homogenization_sizePostResults(mesh_maxNips,mesh_NcpElems))
homogenization_sizePostResults = 0_pInt homogenization_sizePostResults = 0_pInt
allocate(materialpoint_heat(mesh_maxNips,mesh_NcpElems))
materialpoint_heat = 0.0_pReal
allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems)) allocate(materialpoint_dPdF(3,3,3,3,mesh_maxNips,mesh_NcpElems))
materialpoint_dPdF = 0.0_pReal materialpoint_dPdF = 0.0_pReal
allocate(materialpoint_F0(3,3,mesh_maxNips,mesh_NcpElems)) allocate(materialpoint_F0(3,3,mesh_maxNips,mesh_NcpElems))
@ -185,8 +185,6 @@ subroutine homogenization_init(Temperature)
materialpoint_subF = 0.0_pReal materialpoint_subF = 0.0_pReal
allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems)) allocate(materialpoint_P(3,3,mesh_maxNips,mesh_NcpElems))
materialpoint_P = 0.0_pReal materialpoint_P = 0.0_pReal
allocate(materialpoint_Temperature(mesh_maxNips,mesh_NcpElems))
materialpoint_Temperature = Temperature
allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems)) allocate(materialpoint_subFrac(mesh_maxNips,mesh_NcpElems))
materialpoint_subFrac = 0.0_pReal materialpoint_subFrac = 0.0_pReal
allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems)) allocate(materialpoint_subStep(mesh_maxNips,mesh_NcpElems))
@ -204,7 +202,7 @@ subroutine homogenization_init(Temperature)
materialpoint_F = materialpoint_F0 materialpoint_F = materialpoint_F0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate and initialize global state and postrestuls variables ! allocate and initialize global state and postresutls variables
elementLooping: do e = 1,mesh_NcpElems elementLooping: do e = 1,mesh_NcpElems
myInstance = homogenization_typeInstance(mesh_element(3,e)) myInstance = homogenization_typeInstance(mesh_element(3,e))
IpLooping: do i = 1,FE_Nips(FE_geomtype(mesh_element(2,e))) IpLooping: do i = 1,FE_Nips(FE_geomtype(mesh_element(2,e)))
@ -249,36 +247,36 @@ subroutine homogenization_init(Temperature)
write(6,'(/,a)') ' <<<+- homogenization init -+>>>' write(6,'(/,a)') ' <<<+- homogenization init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
write(6,'(a16,a)') ' Current time : ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0) write(6,'(a32,1x,7(i8,1x))') 'homogenization_state0: ', shape(homogenization_state0)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0) write(6,'(a32,1x,7(i8,1x))') 'homogenization_subState0: ', shape(homogenization_subState0)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state) write(6,'(a32,1x,7(i8,1x))') 'homogenization_state: ', shape(homogenization_state)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_sizeState: ', shape(homogenization_sizeState) write(6,'(a32,1x,7(i8,1x))') 'homogenization_sizeState: ', shape(homogenization_sizeState)
write(6,'(a32,1x,7(i8,1x))') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults) write(6,'(a32,1x,7(i8,1x),/)') 'homogenization_sizePostResults: ', shape(homogenization_sizePostResults)
write(6,*)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F0: ', shape(materialpoint_F0)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_F: ', shape(materialpoint_F)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF0: ', shape(materialpoint_subF0) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF0: ', shape(materialpoint_subF0)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF: ', shape(materialpoint_subF) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subF: ', shape(materialpoint_subF)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_P: ', shape(materialpoint_P) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_P: ', shape(materialpoint_P)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_Temperature: ', shape(materialpoint_Temperature) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_heat: ', shape(materialpoint_heat)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subFrac: ', shape(materialpoint_subFrac)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subStep: ', shape(materialpoint_subStep) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subStep: ', shape(materialpoint_subStep)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subdt: ', shape(materialpoint_subdt) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_subdt: ', shape(materialpoint_subdt)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_requested: ', shape(materialpoint_requested)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged) write(6,'(a32,1x,7(i8,1x))') 'materialpoint_converged: ', shape(materialpoint_converged)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy) write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy: ', shape(materialpoint_doneAndHappy)
write(6,*) write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_results: ', shape(materialpoint_results)
write(6,'(a32,1x,7(i8,1x))') 'materialpoint_results: ', shape(materialpoint_results)
write(6,*)
write(6,'(a32,1x,7(i8,1x))') 'maxSizeState: ', homogenization_maxSizeState write(6,'(a32,1x,7(i8,1x))') 'maxSizeState: ', homogenization_maxSizeState
write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', homogenization_maxSizePostResults
endif endif
flush(6) flush(6)
if (debug_g < 1 .or. debug_g > homogenization_Ngrains(mesh_element(3,debug_e))) &
call IO_error(602_pInt,ext_msg='component (grain)')
end subroutine homogenization_init end subroutine homogenization_init
@ -309,7 +307,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
constitutive_partionedState0, & constitutive_partionedState0, &
constitutive_state constitutive_state
use crystallite, only: & use crystallite, only: &
crystallite_Temperature, & crystallite_heat, &
crystallite_F0, & crystallite_F0, &
crystallite_Fp0, & crystallite_Fp0, &
crystallite_Fp, & crystallite_Fp, &
@ -319,7 +317,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_dPdF0, & crystallite_dPdF0, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_Tstar_v, & crystallite_Tstar_v, &
crystallite_partionedTemperature0, &
crystallite_partionedF0, & crystallite_partionedF0, &
crystallite_partionedF, & crystallite_partionedF, &
crystallite_partionedFp0, & crystallite_partionedFp0, &
@ -357,13 +354,10 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize to starting condition ! initialize to starting condition
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt .and. & if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then
debug_e > 0 .and. debug_e <= mesh_NcpElems .and. debug_i > 0 .and. debug_i <= mesh_maxNips) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i
write(6,'(a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i
write(6,'(a,/,12x,f14.9)') '<< HOMOG >> Temp0', &
materialpoint_Temperature(debug_i,debug_e)
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', & write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', &
math_transpose33(materialpoint_F0(1:3,1:3,debug_i,debug_e)) math_transpose33(materialpoint_F0(1:3,1:3,debug_i,debug_e))
write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', & write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', &
@ -377,7 +371,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains) forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures
crystallite_partionedTemperature0(g,i,e) = materialpoint_Temperature(i,e) ! ...temperatures
crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) ! ...plastic def grads crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) ! ...plastic def grads
crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness
@ -427,7 +420,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then
! wind forward grain starting point of... ! wind forward grain starting point of...
crystallite_partionedTemperature0(1:myNgrains,i,e) = crystallite_Temperature(1:myNgrains,i,e) ! ...temperatures
crystallite_partionedF0(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) ! ...def grads crystallite_partionedF0(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) ! ...def grads
crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fp(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fp(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
@ -476,8 +468,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! restore... ! restore...
crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures
! ...initial def grad unchanged
crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic def grads
crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness
@ -575,7 +565,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping4: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) IpLooping4: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
call homogenization_averageStressAndItsTangent(i,e) call homogenization_averageStressAndItsTangent(i,e)
materialpoint_Temperature(i,e) = homogenization_averageTemperature(i,e) materialpoint_heat(i,e) = homogenization_averageHeat(i,e)
enddo IpLooping4 enddo IpLooping4
enddo elementLooping4 enddo elementLooping4
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -639,7 +629,7 @@ subroutine materialpoint_postResults(dt)
grainLooping :do g = 1,myNgrains grainLooping :do g = 1,myNgrains
theSize = (1 + crystallite_sizePostResults(myCrystallite)) + (1 + constitutive_sizePostResults(g,i,e)) theSize = (1 + crystallite_sizePostResults(myCrystallite)) + (1 + constitutive_sizePostResults(g,i,e))
materialpoint_results(thePos+1:thePos+theSize,i,e) = crystallite_postResults(dt,g,i,e) ! tell crystallite results materialpoint_results(thePos+1:thePos+theSize,i,e) = crystallite_postResults(g,i,e) ! tell crystallite results
thePos = thePos + theSize thePos = thePos + theSize
enddo grainLooping enddo grainLooping
enddo IpLooping enddo IpLooping
@ -677,14 +667,11 @@ subroutine homogenization_partitionDeformation(ip,el)
case (homogenization_isostrain_label) chosenHomogenization case (homogenization_isostrain_label) chosenHomogenization
call homogenization_isostrain_partitionDeformation(& call homogenization_isostrain_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),&
materialpoint_subF(1:3,1:3,ip,el),& materialpoint_subF(1:3,1:3,ip,el),&
homogenization_state(ip,el), &
ip, &
el) el)
case (homogenization_RGC_label) chosenHomogenization case (homogenization_RGC_label) chosenHomogenization
call homogenization_RGC_partitionDeformation(crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), & call homogenization_RGC_partitionDeformation(&
crystallite_partionedF0(1:3,1:3,1:homogenization_maxNgrains,ip,el),& crystallite_partionedF(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
materialpoint_subF(1:3,1:3,ip,el),& materialpoint_subF(1:3,1:3,ip,el),&
homogenization_state(ip,el), & homogenization_state(ip,el), &
ip, & ip, &
@ -765,18 +752,18 @@ subroutine homogenization_averageStressAndItsTangent(ip,el)
chosenHomogenization: select case(homogenization_type(mesh_element(3,el))) chosenHomogenization: select case(homogenization_type(mesh_element(3,el)))
case (homogenization_isostrain_label) chosenHomogenization case (homogenization_isostrain_label) chosenHomogenization
call homogenization_isostrain_averageStressAndItsTangent(materialpoint_P(1:3,1:3,ip,el), & call homogenization_isostrain_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), &
ip, &
el) el)
case (homogenization_RGC_label) chosenHomogenization case (homogenization_RGC_label) chosenHomogenization
call homogenization_RGC_averageStressAndItsTangent( materialpoint_P(1:3,1:3,ip,el), & call homogenization_RGC_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_P(1:3,1:3,1:homogenization_maxNgrains,ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_maxNgrains,ip,el), &
ip, &
el) el)
end select chosenHomogenization end select chosenHomogenization
@ -784,16 +771,16 @@ end subroutine homogenization_averageStressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average temperature from constituent quantities (does not depend on choosen !> @brief derive average heat from constituent quantities (does not depend on choosen
!! homogenization scheme) !! homogenization scheme)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) function homogenization_averageTemperature(ip,el) real(pReal) function homogenization_averageHeat(ip,el)
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use material, only: & use material, only: &
homogenization_Ngrains homogenization_Ngrains
use crystallite, only: & use crystallite, only: &
crystallite_Temperature crystallite_heat
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -803,11 +790,11 @@ real(pReal) function homogenization_averageTemperature(ip,el)
Ngrains Ngrains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! computing the average temperature ! computing the average heat
Ngrains = homogenization_Ngrains(mesh_element(3,el)) Ngrains = homogenization_Ngrains(mesh_element(3,el))
homogenization_averageTemperature= sum(crystallite_Temperature(1:Ngrains,ip,el))/real(Ngrains,pReal) homogenization_averageHeat= sum(crystallite_heat(1:Ngrains,ip,el))/real(Ngrains,pReal)
end function homogenization_averageTemperature end function homogenization_averageHeat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -835,9 +822,9 @@ function homogenization_postResults(ip,el)
homogenization_postResults = 0.0_pReal homogenization_postResults = 0.0_pReal
chosenHomogenization: select case (homogenization_type(mesh_element(3,el))) chosenHomogenization: select case (homogenization_type(mesh_element(3,el)))
case (homogenization_isostrain_label) chosenHomogenization case (homogenization_isostrain_label) chosenHomogenization
homogenization_postResults = homogenization_isostrain_postResults(homogenization_state(ip,el),ip,el) homogenization_postResults = homogenization_isostrain_postResults(el)
case (homogenization_RGC_label) chosenHomogenization case (homogenization_RGC_label) chosenHomogenization
homogenization_postResults = homogenization_RGC_postResults(homogenization_state(ip,el),ip,el) homogenization_postResults = homogenization_RGC_postResults(homogenization_state(ip,el),el)
end select chosenHomogenization end select chosenHomogenization
end function homogenization_postResults end function homogenization_postResults

89
code/homogenization_RGC.f90
View File

@ -77,7 +77,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file !> @brief allocates all neccessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_init(myFile) subroutine homogenization_RGC_init(myUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: & use debug, only: &
debug_level, & debug_level, &
@ -101,7 +101,7 @@ subroutine homogenization_RGC_init(myFile)
use material use material
implicit none implicit none
integer(pInt), intent(in) :: myFile !< file pointer to material configuration integer(pInt), intent(in) :: myUnit !< file pointer to material configuration
integer(pInt), parameter :: MAXNCHUNKS = 4_pInt integer(pInt), parameter :: MAXNCHUNKS = 4_pInt
integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
integer(pInt) ::section=0_pInt, maxNinstance, i,j,e, output=-1_pInt, mySize, myInstance integer(pInt) ::section=0_pInt, maxNinstance, i,j,e, output=-1_pInt, mySize, myInstance
@ -130,14 +130,14 @@ subroutine homogenization_RGC_init(myFile)
allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems)) allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems))
homogenization_RGC_orientation = spread(spread(math_I3,3,mesh_maxNips),4,mesh_NcpElems) ! initialize to identity homogenization_RGC_orientation = spread(spread(math_I3,3,mesh_maxNips),4,mesh_NcpElems) ! initialize to identity
rewind(myFile) rewind(myUnit)
do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to <homogenization> do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to <homogenization>
line = IO_read(myFile) line = IO_read(myUnit)
enddo enddo
do while (trim(line) /= '#EOF#') do while (trim(line) /= '#EOF#')
line = IO_read(myFile) line = IO_read(myUnit)
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
@ -199,8 +199,7 @@ subroutine homogenization_RGC_init(myFile)
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
do i = 1_pInt,maxNinstance do i = 1_pInt,maxNinstance
write(6,'(a15,1x,i4)') 'instance: ', i write(6,'(a15,1x,i4,/)') 'instance: ', i
write(6,*)
write(6,'(a25,3(1x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1_pInt,3_pInt) write(6,'(a25,3(1x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1_pInt,3_pInt)
write(6,'(a25,1x,e10.3)') 'scaling parameter: ', homogenization_RGC_xiAlpha(i) write(6,'(a25,1x,e10.3)') 'scaling parameter: ', homogenization_RGC_xiAlpha(i)
write(6,'(a25,1x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i) write(6,'(a25,1x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i)
@ -228,11 +227,11 @@ subroutine homogenization_RGC_init(myFile)
mySize = 0_pInt mySize = 0_pInt
end select end select
if (mySize > 0_pInt) then ! any meaningful output found outputFound: if (mySize > 0_pInt) then
homogenization_RGC_sizePostResult(j,i) = mySize homogenization_RGC_sizePostResult(j,i) = mySize
homogenization_RGC_sizePostResults(i) = & homogenization_RGC_sizePostResults(i) = &
homogenization_RGC_sizePostResults(i) + mySize homogenization_RGC_sizePostResults(i) + mySize
endif endif outputFound
enddo enddo
homogenization_RGC_sizeState(i) & homogenization_RGC_sizeState(i) &
@ -249,7 +248,7 @@ end subroutine homogenization_RGC_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_partitionDeformation(F,F0,avgF,state,ip,el) subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use debug, only: & use debug, only: &
@ -268,7 +267,6 @@ subroutine homogenization_RGC_partitionDeformation(F,F0,avgF,state,ip,el)
implicit none implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 !< initial partioned F per grain
real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F
type(p_vec), intent(in) :: state type(p_vec), intent(in) :: state
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -427,7 +425,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating volume discrepancy and stress penalty related to overall volume discrepancy ! calculating volume discrepancy and stress penalty related to overall volume discrepancy
call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID) call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! debugging the mismatch, stress and penalties of grains ! debugging the mismatch, stress and penalties of grains
@ -513,6 +511,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
endif endif
homogenization_RGC_updateState = .false. homogenization_RGC_updateState = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! If convergence reached => done and happy ! If convergence reached => done and happy
if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then
@ -521,8 +520,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a55)')'... done and happy' write(6,'(1x,a55,/)')'... done and happy'
write(6,*)' '
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -552,16 +550,14 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a30,1x,e15.8)')'Constitutive work: ',constitutiveWork write(6,'(1x,a30,1x,e15.8)') 'Constitutive work: ',constitutiveWork
write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/real(nGrain,pReal), & write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/real(nGrain,pReal), &
sum(NN(2,:))/real(nGrain,pReal), & sum(NN(2,:))/real(nGrain,pReal), &
sum(NN(3,:))/real(nGrain,pReal) sum(NN(3,:))/real(nGrain,pReal)
write(6,'(1x,a30,1x,e15.8)')'Penalty energy: ',penaltyEnergy write(6,'(1x,a30,1x,e15.8)') 'Penalty energy: ',penaltyEnergy
write(6,'(1x,a30,1x,e15.8)')'Volume discrepancy: ',volDiscrep write(6,'(1x,a30,1x,e15.8,/)') 'Volume discrepancy: ',volDiscrep
write(6,*)'' write(6,'(1x,a30,1x,e15.8)') 'Maximum relaxation rate: ',maxval(abs(drelax))/dt
write(6,'(1x,a30,1x,e15.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt write(6,'(1x,a30,1x,e15.8,/)') 'Average relaxation rate: ',sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal)
write(6,'(1x,a30,1x,e15.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal)
write(6,*)''
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -577,8 +573,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a55)')'... broken' write(6,'(1x,a55,/)')'... broken'
write(6,*)' '
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -589,8 +584,7 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(1x,a55)')'... not yet done' write(6,'(1x,a55,/)')'... not yet done'
write(6,*)' '
flush(6) flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
endif endif
@ -668,9 +662,9 @@ function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el
p_relax = relax p_relax = relax
p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector
state%p(1:3*nIntFaceTot) = p_relax state%p(1:3*nIntFaceTot) = p_relax
call homogenization_RGC_grainDeformation(pF,F0,avgF,state,ip,el) ! compute the grains deformation from perturbed state call homogenization_RGC_grainDeformation(pF,avgF,state,ip,el) ! compute the grains deformation from perturbed state
call homogenization_RGC_stressPenalty(pR,pNN,avgF,pF,ip,el,homID) ! compute stress penalty due to interface mismatch from perturbed state call homogenization_RGC_stressPenalty(pR,pNN,avgF,pF,ip,el,homID) ! compute stress penalty due to interface mismatch from perturbed state
call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el,homID) ! compute stress penalty due to volume discrepancy from perturbed state call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el) ! compute stress penalty due to volume discrepancy from perturbed state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! computing the global stress residual array from the perturbed state ! computing the global stress residual array from the perturbed state
@ -814,9 +808,7 @@ end function homogenization_RGC_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ip,el ) subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,el)
use prec, only: &
p_vec
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_homogenization,& debug_homogenization,&
@ -826,15 +818,14 @@ subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,
use math, only: math_Plain3333to99 use math, only: math_Plain3333to99
implicit none implicit none
real(pReal), dimension (3,3), intent(out) :: avgP ! average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF ! average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P ! array of current grain stresses real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P !< array of current grain stresses
real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF ! array of current grain stiffnesses real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffnesses
integer(pInt), intent(in) :: el !< element number
real(pReal), dimension (9,9) :: dPdF99 real(pReal), dimension (9,9) :: dPdF99
integer(pInt), intent(in) :: &
ip, & ! integration point number integer(pInt) :: homID, i, j, Ngrains, iGrain
el ! element number
integer(pInt) homID, i, j, Ngrains, iGrain
homID = homogenization_typeInstance(mesh_element(3,el)) homID = homogenization_typeInstance(mesh_element(3,el))
Ngrains = homogenization_Ngrains(mesh_element(3,el)) Ngrains = homogenization_Ngrains(mesh_element(3,el))
@ -866,7 +857,7 @@ end subroutine homogenization_RGC_averageStressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion !> @brief return array of homogenization results for post file inclusion
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function homogenization_RGC_postResults(state,ip,el) pure function homogenization_RGC_postResults(state,el)
use prec, only: & use prec, only: &
p_vec p_vec
use mesh, only: & use mesh, only: &
@ -877,9 +868,7 @@ pure function homogenization_RGC_postResults(state,ip,el)
implicit none implicit none
type(p_vec), intent(in) :: state ! my State type(p_vec), intent(in) :: state ! my State
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: el ! element number
ip, & ! integration point number
el ! element number
integer(pInt) homID,o,c,nIntFaceTot integer(pInt) homID,o,c,nIntFaceTot
real(pReal), dimension(homogenization_RGC_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: & real(pReal), dimension(homogenization_RGC_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: &
homogenization_RGC_postResults homogenization_RGC_postResults
@ -943,10 +932,10 @@ subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID)
xSmoo_RGC xSmoo_RGC
implicit none implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen ! stress-like penalty real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen !< stress-like penalty
real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis ! total amount of mismatch real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis !< total amount of mismatch
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef ! deformation gradients real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef !< deformation gradients
real(pReal), dimension (3,3), intent(in) :: avgF ! initial effective stretch tensor real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor
integer(pInt), intent(in) :: ip,el integer(pInt), intent(in) :: ip,el
integer(pInt), dimension (4) :: intFace integer(pInt), dimension (4) :: intFace
integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim
@ -960,7 +949,6 @@ subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID)
real(pReal), parameter :: nDefToler = 1.0e-10_pReal real(pReal), parameter :: nDefToler = 1.0e-10_pReal
nGDim = homogenization_RGC_Ngrains(1:3,homID) nGDim = homogenization_RGC_Ngrains(1:3,homID)
rPen = 0.0_pReal rPen = 0.0_pReal
nMis = 0.0_pReal nMis = 0.0_pReal
@ -1064,7 +1052,7 @@ end subroutine homogenization_RGC_stressPenalty
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to volume discrepancy !> @brief calculate stress-like penalty due to volume discrepancy
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el, homID) subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el)
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_homogenization,& debug_homogenization,&
@ -1092,7 +1080,7 @@ subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el, homID
integer(pInt), intent(in) :: ip,& ! integration point integer(pInt), intent(in) :: ip,& ! integration point
el el
real(pReal), dimension (homogenization_maxNgrains) :: gVol real(pReal), dimension (homogenization_maxNgrains) :: gVol
integer(pInt) :: homID,iGrain,nGrain,i,j integer(pInt) :: iGrain,nGrain,i,j
nGrain = homogenization_Ngrains(mesh_element(3,el)) nGrain = homogenization_Ngrains(mesh_element(3,el))
@ -1434,7 +1422,7 @@ end function homogenization_RGC_interface1to4
!> @brief calculating the grain deformation gradient (the same with !> @brief calculating the grain deformation gradient (the same with
! homogenization_RGC_partionDeformation, but used only for perturbation scheme) ! homogenization_RGC_partionDeformation, but used only for perturbation scheme)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_grainDeformation(F, F0, avgF, state, ip, el) subroutine homogenization_RGC_grainDeformation(F, avgF, state, ip, el)
use prec, only: & use prec, only: &
p_vec p_vec
use mesh, only: & use mesh, only: &
@ -1446,7 +1434,6 @@ subroutine homogenization_RGC_grainDeformation(F, F0, avgF, state, ip, el)
implicit none implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 !< initiatial partioned F per grain
real(pReal), dimension (3,3), intent(in) :: avgF !< real(pReal), dimension (3,3), intent(in) :: avgF !<
type(p_vec), intent(in) :: state type(p_vec), intent(in) :: state
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &

48
code/homogenization_isostrain.f90
View File

@ -55,7 +55,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file !> @brief allocates all neccessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_isostrain_init(myFile) subroutine homogenization_isostrain_init(myUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use math, only: & use math, only: &
math_Mandel3333to66, & math_Mandel3333to66, &
@ -64,7 +64,7 @@ subroutine homogenization_isostrain_init(myFile)
use material use material
implicit none implicit none
integer(pInt), intent(in) :: myFile integer(pInt), intent(in) :: myUnit
integer(pInt), parameter :: MAXNCHUNKS = 2_pInt integer(pInt), parameter :: MAXNCHUNKS = 2_pInt
integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
integer(pInt) :: & integer(pInt) :: &
@ -96,15 +96,15 @@ subroutine homogenization_isostrain_init(myFile)
allocate(homogenization_isostrain_output(maxval(homogenization_Noutput),maxNinstance)) allocate(homogenization_isostrain_output(maxval(homogenization_Noutput),maxNinstance))
homogenization_isostrain_output = '' homogenization_isostrain_output = ''
rewind(myFile) rewind(myUnit)
section = 0_pInt section = 0_pInt
do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to <homogenization> do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to <homogenization>
line = IO_read(myFile) line = IO_read(myUnit)
enddo enddo
do while (trim(line) /= '#EOF#') do while (trim(line) /= '#EOF#')
line = IO_read(myFile) line = IO_read(myUnit)
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
@ -120,7 +120,7 @@ subroutine homogenization_isostrain_init(myFile)
case ('(output)') case ('(output)')
output = output + 1_pInt output = output + 1_pInt
homogenization_isostrain_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt)) homogenization_isostrain_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('ngrains') case ('ngrains','ncomponents')
homogenization_isostrain_Ngrains(i) = IO_intValue(line,positions,2_pInt) homogenization_isostrain_Ngrains(i) = IO_intValue(line,positions,2_pInt)
case ('mapping') case ('mapping')
homogenization_isostrain_mapping(i) = IO_lc(IO_stringValue(line,positions,2_pInt)) homogenization_isostrain_mapping(i) = IO_lc(IO_stringValue(line,positions,2_pInt))
@ -134,17 +134,17 @@ subroutine homogenization_isostrain_init(myFile)
do j = 1_pInt,maxval(homogenization_Noutput) do j = 1_pInt,maxval(homogenization_Noutput)
select case(homogenization_isostrain_output(j,i)) select case(homogenization_isostrain_output(j,i))
case('ngrains') case('ngrains','ncomponents')
mySize = 1_pInt mySize = 1_pInt
case default case default
mySize = 0_pInt mySize = 0_pInt
end select end select
if (mySize > 0_pInt) then ! any meaningful output found outputFound: if (mySize > 0_pInt) then
homogenization_isostrain_sizePostResult(j,i) = mySize homogenization_isostrain_sizePostResult(j,i) = mySize
homogenization_isostrain_sizePostResults(i) = & homogenization_isostrain_sizePostResults(i) = &
homogenization_isostrain_sizePostResults(i) + mySize homogenization_isostrain_sizePostResults(i) + mySize
endif endif outputFound
enddo enddo
enddo enddo
@ -154,10 +154,9 @@ end subroutine homogenization_isostrain_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_isostrain_partitionDeformation(F,F0,avgF,state,ip,el) subroutine homogenization_isostrain_partitionDeformation(F,avgF,el)
use prec, only: & use prec, only: &
pReal, & pReal
p_vec
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use material, only: & use material, only: &
@ -165,12 +164,9 @@ subroutine homogenization_isostrain_partitionDeformation(F,F0,avgF,state,ip,el)
homogenization_Ngrains homogenization_Ngrains
implicit none implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F ! partioned def grad per grain real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned def grad per grain
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: F0 ! initial partioned def grad per grain real(pReal), dimension (3,3), intent(in) :: avgF !< my average def grad
real(pReal), dimension (3,3), intent(in) :: avgF ! my average def grad
type(p_vec), intent(in) :: state ! my state
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ip, & !< integration point number
el !< element number el !< element number
F = spread(avgF,3,homogenization_Ngrains(mesh_element(3,el))) F = spread(avgF,3,homogenization_Ngrains(mesh_element(3,el)))
@ -181,7 +177,7 @@ end subroutine homogenization_isostrain_partitionDeformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ip,el) subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,el)
use prec, only: & use prec, only: &
pReal pReal
use mesh, only: & use mesh, only: &
@ -196,9 +192,7 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P !< array of current grain stresses real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P !< array of current grain stresses
real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffnesses real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffnesses
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: el !< element number
ip, & !< integration point number
el !< element number
integer(pInt) :: & integer(pInt) :: &
homID, & homID, &
Ngrains Ngrains
@ -224,10 +218,9 @@ end subroutine homogenization_isostrain_averageStressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion !> @brief return array of homogenization results for post file inclusion
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function homogenization_isostrain_postResults(state,ip,el) pure function homogenization_isostrain_postResults(el)
use prec, only: & use prec, only: &
pReal,& pReal
p_vec
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use material, only: & use material, only: &
@ -235,10 +228,7 @@ pure function homogenization_isostrain_postResults(state,ip,el)
homogenization_Noutput homogenization_Noutput
implicit none implicit none
type(p_vec), intent(in) :: state integer(pInt), intent(in) :: el !< element number
integer(pInt), intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(homogenization_isostrain_sizePostResults & real(pReal), dimension(homogenization_isostrain_sizePostResults &
(homogenization_typeInstance(mesh_element(3,el)))) :: & (homogenization_typeInstance(mesh_element(3,el)))) :: &
homogenization_isostrain_postResults homogenization_isostrain_postResults
@ -253,7 +243,7 @@ pure function homogenization_isostrain_postResults(state,ip,el)
do o = 1_pInt,homogenization_Noutput(mesh_element(3,el)) do o = 1_pInt,homogenization_Noutput(mesh_element(3,el))
select case(homogenization_isostrain_output(o,homID)) select case(homogenization_isostrain_output(o,homID))
case ('ngrains') case ('ngrains','ncomponents')
homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal) homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal)
c = c + 1_pInt c = c + 1_pInt
end select end select

7
code/material.f90
View File

@ -122,14 +122,17 @@ module material
homogenization_active homogenization_active
public :: material_init public :: &
material_init
private :: material_parseHomogenization, & private :: &
material_parseHomogenization, &
material_parseMicrostructure, & material_parseMicrostructure, &
material_parseCrystallite, & material_parseCrystallite, &
material_parsePhase, & material_parsePhase, &
material_parseTexture, & material_parseTexture, &
material_populateGrains material_populateGrains
contains contains

4
code/mesh.f90
View File

@ -669,7 +669,7 @@ integer(pInt) function mesh_FEasCP(what,myID)
return return
endif endif
do while (upper-lower > 1_pInt) ! binary search in between bounds binarySearch: do while (upper-lower > 1_pInt)
center = (lower+upper)/2_pInt center = (lower+upper)/2_pInt
if (lookupMap(1_pInt,center) < myID) then if (lookupMap(1_pInt,center) < myID) then
lower = center lower = center
@ -679,7 +679,7 @@ integer(pInt) function mesh_FEasCP(what,myID)
mesh_FEasCP = lookupMap(2_pInt,center) mesh_FEasCP = lookupMap(2_pInt,center)
exit exit
endif endif
enddo enddo binarySearch
end function mesh_FEasCP end function mesh_FEasCP

5
code/numerics.f90
View File

@ -60,7 +60,6 @@ module numerics
stepIncreaseCryst = 1.5_pReal, & !< increase of next substep size when previous substep converged in crystallite stepIncreaseCryst = 1.5_pReal, & !< increase of next substep size when previous substep converged in crystallite
stepIncreaseHomog = 1.5_pReal, & !< increase of next substep size when previous substep converged in homogenization stepIncreaseHomog = 1.5_pReal, & !< increase of next substep size when previous substep converged in homogenization
rTol_crystalliteState = 1.0e-6_pReal, & !< relative tolerance in crystallite state loop rTol_crystalliteState = 1.0e-6_pReal, & !< relative tolerance in crystallite state loop
rTol_crystalliteTemperature= 1.0e-6_pReal, & !< relative tolerance in crystallite temperature loop
rTol_crystalliteStress = 1.0e-6_pReal, & !< relative tolerance in crystallite stress loop rTol_crystalliteStress = 1.0e-6_pReal, & !< relative tolerance in crystallite stress loop
aTol_crystalliteStress = 1.0e-8_pReal, & !< absolute tolerance in crystallite stress loop, Default 1.0e-8: residuum is in Lp and hence strain is on this order aTol_crystalliteStress = 1.0e-8_pReal, & !< absolute tolerance in crystallite stress loop, Default 1.0e-8: residuum is in Lp and hence strain is on this order
numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit numerics_unitlength = 1.0_pReal, & !< determines the physical length of one computational length unit
@ -222,8 +221,6 @@ subroutine numerics_init
stepIncreaseHomog = IO_floatValue(line,positions,2_pInt) stepIncreaseHomog = IO_floatValue(line,positions,2_pInt)
case ('rtol_crystallitestate') case ('rtol_crystallitestate')
rTol_crystalliteState = IO_floatValue(line,positions,2_pInt) rTol_crystalliteState = IO_floatValue(line,positions,2_pInt)
case ('rtol_crystallitetemperature')
rTol_crystalliteTemperature = IO_floatValue(line,positions,2_pInt)
case ('rtol_crystallitestress') case ('rtol_crystallitestress')
rTol_crystalliteStress = IO_floatValue(line,positions,2_pInt) rTol_crystalliteStress = IO_floatValue(line,positions,2_pInt)
case ('atol_crystallitestress') case ('atol_crystallitestress')
@ -377,7 +374,6 @@ subroutine numerics_init
write(6,'(a24,1x,i8)') ' nState: ',nState write(6,'(a24,1x,i8)') ' nState: ',nState
write(6,'(a24,1x,i8)') ' nStress: ',nStress write(6,'(a24,1x,i8)') ' nStress: ',nStress
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState write(6,'(a24,1x,es8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteTemp: ',rTol_crystalliteTemperature
write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator
@ -469,7 +465,6 @@ subroutine numerics_init
if (subStepSizeHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeHomog') if (subStepSizeHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='subStepSizeHomog')
if (stepIncreaseHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseHomog') if (stepIncreaseHomog <= 0.0_pReal) call IO_error(301_pInt,ext_msg='stepIncreaseHomog')
if (rTol_crystalliteState <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteState') if (rTol_crystalliteState <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteState')
if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteTemperature')
if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteStress') if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='rTol_crystalliteStress')
if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='aTol_crystalliteStress') if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(301_pInt,ext_msg='aTol_crystalliteStress')
if (any(numerics_integrator <= 0_pInt) .or. any(numerics_integrator >= 6_pInt)) & if (any(numerics_integrator <= 0_pInt) .or. any(numerics_integrator >= 6_pInt)) &